Belle II Software  release-05-02-19
PXDCalibrationUtilities.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2021 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Qingyuan Liu *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <framework/logging/Logger.h>
12 #include <pxd/calibration/PXDCalibrationUtilities.h>
13 #include <vxd/dataobjects/VxdID.h>
14 
15 #include <string>
16 #include <algorithm>
17 #include <set>
18 
19 #include <sstream>
20 #include <iostream>
21 
22 #include <boost/format.hpp>
23 
24 #include <TF1.h>
25 
26 using namespace std;
27 using boost::format;
28 
29 namespace Belle2 {
34  namespace PXD {
35 
37  void getNumberOfBins(const std::shared_ptr<TH1I>& histo_ptr, unsigned short& nBinsU, unsigned short& nBinsV)
38  {
39  set<unsigned short> uBinSet;
40  set<unsigned short> vBinSet;
41 
42  // Loop over all bins of input histo
43  for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
44  // The bin label contains the vxdid, uBin and vBin
45  string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
46 
47  // Parse label string format to read sensorID, uBin and vBin
48  istringstream stream(label);
49  string token;
50  getline(stream, token, '_');
51  getline(stream, token, '_');
52  unsigned short uBin = std::stoi(token);
53 
54  getline(stream, token, '_');
55  unsigned short vBin = std::stoi(token);
56 
57  uBinSet.insert(uBin);
58  vBinSet.insert(vBin);
59  }
60 
61  if (uBinSet.empty() || vBinSet.empty()) {
62  B2FATAL("Not able to determine the grid size. Something is wrong with collected data.");
63  } else {
64  nBinsU = *uBinSet.rbegin() + 1;
65  nBinsV = *vBinSet.rbegin() + 1;
66  }
67  }
68 
70  unsigned short getNumberOfSensors(const std::shared_ptr<TH1I>& histo_ptr)
71  {
72  set<unsigned short> sensorSet;
73 
74  // Loop over all bins of input histo
75  for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
76  // The bin label contains the vxdid, uBin and vBin
77  string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
78 
79  // Parse label string format to read sensorID, uBin and vBin
80  istringstream stream(label);
81  string token;
82  getline(stream, token, '_');
83  VxdID sensorID(token);
84  sensorSet.insert(sensorID.getID());
85  }
86  return sensorSet.size();
87  }
88 
90  double CalculateMedian(std::vector<double>& signals)
91  {
92  auto size = signals.size();
93 
94  if (size == 0) {
95  return 0.0; // Undefined, really.
96  } else if (size <= 100) {
97  // sort() or partial_sort is in O(NlogN)
98  sort(signals.begin(), signals.end());
99  if (size % 2 == 0) {
100  return (signals[size / 2 - 1] + signals[size / 2]) / 2;
101  } else {
102  return signals[size / 2];
103  }
104  } else {
105  // nth_element or max_element in O(N) only
106  // All elements before the nth are guanranteed smaller
107  auto n = size / 2;
108  nth_element(signals.begin(), signals.begin() + n, signals.end());
109  auto med = signals[n];
110  if (!(size & 1)) { // if size is even
111  auto max_it = max_element(signals.begin(), signals.begin() + n);
112  med = (*max_it + med) / 2.0;
113  }
114  return med;
115  }
116  }
118  double CalculateMedian(TH1* hist)
119  {
120  double quantiles[1]; // One element just for median
121  double probSums[1] = {0.5}; // median definiton
122  hist->GetQuantiles(1, quantiles, probSums);
123  return quantiles[0];
124  }
125 
127  double FitLandau(TH1* hist)
128  {
129  auto size = hist->GetEntries();
130  if (size == 0) return 0.0; // Undefined.
131 
132  int max = hist->GetBinLowEdge(hist->GetNbinsX() + 1);
133  int min = hist->GetBinLowEdge(1);
134 
135  // create fit function
136  TF1* landau = new TF1("landau", "TMath::Landau(x,[0],[1])*[2]", min, max);
137  landau->SetParNames("MPV", "sigma", "scale");
138  landau->SetParameters(1., 0.1, 1000);
139  landau->SetParLimits(0, 0., 3.);
140 
141  Int_t status = hist->Fit("landau", "Lq", "", min, max);
142  double MPV = landau->GetParameter("MPV");
143 
144  B2INFO("Fit result: " << status << " MPV " << MPV << " sigma " << landau->GetParameter("sigma")
145  << " scale " << landau->GetParameter("scale") << " chi2 " << landau->GetChisquare());
146 
147  // clean up
148  delete landau;
149 
150  // check fit status
151  if (status == 0) return MPV;
152  else {
153  B2WARNING("Fit failed!. using default value.");
154  return 0.0;
155  }
156  }
157 
159  double FitLandau(std::vector<double>& signals)
160  {
161  auto size = signals.size();
162  if (size == 0) return 0.0; // Undefined, really.
163 
164  // get max and min values of signal vector
165  double max = *max_element(signals.begin(), signals.end());
166  double min = *min_element(signals.begin(), signals.end());
167 
168  // create histogram to hold signals and fill it
169  TH1D* hist_signals = new TH1D("", "", max - min, min, max);
170  for (auto it = signals.begin(); it != signals.end(); ++it) {
171  hist_signals->Fill(*it);
172  }
173 
174  double MPV = FitLandau(hist_signals);
175  delete hist_signals;
176 
177  return MPV;
178  }
179  }// namespace PXD
181 } // namespace Belle2
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::PXD::FitLandau
double FitLandau(TH1 *hist)
Helper function to estimate MPV from 1D histogram.
Definition: PXDCalibrationUtilities.cc:127
Belle2::VxdID::getID
baseType getID() const
Get the unique id.
Definition: VxdID.h:104
Belle2::PXD::getNumberOfBins
void getNumberOfBins(const std::shared_ptr< TH1I > &histo_ptr, unsigned short &nBinsU, unsigned short &nBinsV)
Helper function to extract number of bins along u side and v side from counter histogram labels.
Definition: PXDCalibrationUtilities.cc:37
Belle2::PXD::CalculateMedian
double CalculateMedian(std::vector< double > &signals)
Helper function to calculate a median from unsorted signal vector.
Definition: PXDCalibrationUtilities.cc:90
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::PXD::getNumberOfSensors
unsigned short getNumberOfSensors(const std::shared_ptr< TH1I > &histo_ptr)
Helper function to extract number of sensors from counter histogram labels.
Definition: PXDCalibrationUtilities.cc:70