Belle II Software  release-05-02-19
PXDChargeCalibrationAlgorithm.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Benjamin Schwenker *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <pxd/calibration/PXDChargeCalibrationAlgorithm.h>
12 #include <pxd/dbobjects/PXDClusterChargeMapPar.h>
13 
14 #include <string>
15 #include <algorithm>
16 #include <set>
17 
18 #include <sstream>
19 #include <iostream>
20 
21 #include <boost/format.hpp>
22 
23 //ROOT
24 #include <TRandom.h>
25 #include <TH1I.h>
26 #include <TF1.h>
27 
28 
29 using namespace std;
30 using boost::format;
31 using namespace Belle2;
32 
33 
34 // Anonymous namespace for data objects used by PXDChargeCalibrationAlgorithm class
35 namespace {
36 
38  int m_signal;
39 
41  void getNumberOfBins(const std::shared_ptr<TH1I>& histo_ptr, unsigned short& nBinsU, unsigned short& nBinsV)
42  {
43  set<unsigned short> uBinSet;
44  set<unsigned short> vBinSet;
45 
46  // Loop over all bins of input histo
47  for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
48  // The bin label contains the vxdid, uBin and vBin
49  string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
50 
51  // Parse label string format to read sensorID, uBin and vBin
52  istringstream stream(label);
53  string token;
54  getline(stream, token, '_');
55  getline(stream, token, '_');
56  unsigned short uBin = std::stoi(token);
57 
58  getline(stream, token, '_');
59  unsigned short vBin = std::stoi(token);
60 
61  uBinSet.insert(uBin);
62  vBinSet.insert(vBin);
63  }
64 
65  if (uBinSet.empty() || vBinSet.empty()) {
66  B2FATAL("Not able to determine the grid size. Something is wrong with collected data.");
67  } else {
68  nBinsU = *uBinSet.rbegin() + 1;
69  nBinsV = *vBinSet.rbegin() + 1;
70  }
71  }
72 
74  unsigned short getNumberOfSensors(const std::shared_ptr<TH1I>& histo_ptr)
75  {
76  set<unsigned short> sensorSet;
77 
78  // Loop over all bins of input histo
79  for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
80  // The bin label contains the vxdid, uBin and vBin
81  string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
82 
83  // Parse label string format to read sensorID, uBin and vBin
84  istringstream stream(label);
85  string token;
86  getline(stream, token, '_');
87  VxdID sensorID(token);
88  sensorSet.insert(sensorID.getID());
89  }
90 
91  return sensorSet.size();
92  }
93 
94 }
95 
96 
97 PXDChargeCalibrationAlgorithm::PXDChargeCalibrationAlgorithm():
98  CalibrationAlgorithm("PXDClusterChargeCollector"),
99  minClusters(5000), noiseSigma(0.6), safetyFactor(2.0), forceContinue(false), strategy(0)
100 {
102  " ------------------------------ PXDChargeCalibrationAlgorithm -----------------------------------\n"
103  " \n"
104  " Algorithm for estimating pxd median/MPV cluster charges for different position on sensor in ADU \n"
105  " ------------------------------------------------------------------------------------------------\n"
106  );
107 }
108 
109 
110 
111 
113 {
114 
115  // Get counter histogram
116  auto cluster_counter = getObjectPtr<TH1I>("PXDClusterCounter");
117 
118  // Extract number of sensors from counter histograms
119  auto nSensors = getNumberOfSensors(cluster_counter);
120 
121  // Extract the number of grid bins from counter histograms
122  unsigned short nBinsU = 0;
123  unsigned short nBinsV = 0;
124  getNumberOfBins(cluster_counter, nBinsU, nBinsV);
125 
126  // Check that we have collected enough Data
127  if (cluster_counter->GetEntries() < int(safetyFactor * minClusters * nSensors * nBinsU * nBinsV)) {
128  if (not forceContinue) {
129  B2INFO("Not enough Data: Only " << cluster_counter->GetEntries() << " hits were collected but " << int(safetyFactor * minClusters *
130  nSensors * nBinsU * nBinsV) << " needed!");
131  return c_NotEnoughData;
132  } else {
133  B2INFO("Continue despite low statistics: Only " << cluster_counter->GetEntries() << " hits were collected but" << int(
135  nSensors * nBinsU * nBinsV) << " would be desirable!");
136  }
137  }
138 
139  B2INFO("Start calibration using a " << nBinsU << "x" << nBinsV << " grid per sensor.");
140 
141  // This is the PXD charge calibration payload for conditions DB
142  PXDClusterChargeMapPar* chargeMapPar = new PXDClusterChargeMapPar(nBinsU, nBinsV);
143 
144  // Loop over all bins of input histo
145  for (auto histoBin = 1; histoBin <= cluster_counter->GetXaxis()->GetNbins(); histoBin++) {
146  // The bin label contains the vxdid, uBin and vBin
147  string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
148 
149  // Parse label string format to read sensorID, uBin and vBin
150  istringstream stream(label);
151  string token;
152  getline(stream, token, '_');
153  VxdID sensorID(token);
154 
155  getline(stream, token, '_');
156  unsigned short uBin = std::stoi(token);
157 
158  getline(stream, token, '_');
159  unsigned short vBin = std::stoi(token);
160 
161  // Read back the counters for number of collected clusters
162  int numberOfDataHits = cluster_counter->GetBinContent(histoBin);
163 
164  // Only perform fitting, when enough data is available
165  if (numberOfDataHits >= minClusters) {
166 
167  // Compute the median cluster charge for the part of PXD
168  auto Charge = EstimateCharge(sensorID, uBin, vBin);
169 
170  // Store the charge
171  chargeMapPar->setContent(sensorID.getID(), uBin, vBin, Charge);
172  } else {
173  B2WARNING(label << ": Number of data hits too small for fitting (" << numberOfDataHits << " < " << minClusters <<
174  "). Use default value of 0.");
175  chargeMapPar->setContent(sensorID.getID(), uBin, vBin, 0.0);
176  }
177  }
178 
179  // Save the charge map to database. Note that this will set the database object name to the same as the collector but you
180  // are free to change it.
181  saveCalibration(chargeMapPar, "PXDClusterChargeMapPar");
182 
183  B2INFO("PXD Cluster Charge Calibration Successful");
184  return c_OK;
185 }
186 
187 
188 double PXDChargeCalibrationAlgorithm::EstimateCharge(VxdID sensorID, unsigned short uBin, unsigned short vBin)
189 {
190 
191  // Construct a tree name for requested part of PXD
192  auto layerNumber = sensorID.getLayerNumber();
193  auto ladderNumber = sensorID.getLadderNumber();
194  auto sensorNumber = sensorID.getSensorNumber();
195  const string treename = str(format("tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
196  // Vector with cluster signals from collected data
197  vector<double> signals;
198 
199  // Fill data_signal vector from input data
200  auto tree = getObjectPtr<TTree>(treename);
201  tree->SetBranchAddress("signal", &m_signal);
202 
203  // Loop over tree
204  const auto nEntries = tree->GetEntries();
205  for (int i = 0; i < nEntries; ++i) {
206  tree->GetEntry(i);
207 
208  double noise = gRandom->Gaus(0.0, noiseSigma);
209  signals.push_back(m_signal + noise);
210  }
211  if (strategy == 0) return CalculateMedian(signals);
212  if (strategy == 1) return FitLandau(signals);
213  else {
214  B2FATAL("strategy unavailable, use 0 for medians or 1 for landau fit!");
215  }
216 }
217 
218 
219 double PXDChargeCalibrationAlgorithm::CalculateMedian(vector<double>& signals)
220 {
221  auto size = signals.size();
222 
223  if (size == 0) {
224  return 0.0; // Undefined, really.
225  } else {
226  sort(signals.begin(), signals.end());
227  if (size % 2 == 0) {
228  return (signals[size / 2 - 1] + signals[size / 2]) / 2;
229  } else {
230  return signals[size / 2];
231  }
232  }
233 }
234 
235 double PXDChargeCalibrationAlgorithm::FitLandau(vector<double>& signals)
236 {
237  auto size = signals.size();
238  if (size == 0) return 0.0; // Undefined, really.
239 
240  // get max and min values of vector
241  int max = *max_element(signals.begin(), signals.end());
242  int min = *min_element(signals.begin(), signals.end());
243 
244 
245  // create histogram to hold signals and fill it
246  TH1D* hist_signals = new TH1D("", "", max - min, min, max);
247  for (auto it = signals.begin(); it != signals.end(); ++it) {
248  hist_signals->Fill(*it);
249  }
250 
251  // create fit function
252  TF1* landau = new TF1("landau", "TMath::Landau(x,[0],[1])*[2]", min, max);
253  landau->SetParNames("MPV", "sigma", "scale");
254  landau->SetParameters(100, 1, 1000);
255 
256  //do fit and get results, fit range restricted to exclude low charge peak
257  Int_t status = hist_signals->Fit("landau", "Lq", "", 30, 350);
258  double MPV = landau->GetParameter("MPV");
259 
260  // clean up
261  delete hist_signals;
262  delete landau;
263 
264  // check fit status
265  if (status == 0) return MPV;
266  else {
267  B2WARNING("Fit failed! using default value 0.0!");
268  return 0.0;
269  }
270 }
Belle2::PXDChargeCalibrationAlgorithm::noiseSigma
float noiseSigma
Artificial noise sigma for smearing cluster charge.
Definition: PXDChargeCalibrationAlgorithm.h:43
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::PXDChargeCalibrationAlgorithm::calibrate
virtual EResult calibrate() override
Run algo on data.
Definition: PXDChargeCalibrationAlgorithm.cc:112
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:290
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:108
Belle2::PXDChargeCalibrationAlgorithm::strategy
int strategy
strategy to used for gain calibration, 0 for medians, 1 for landau fit
Definition: PXDChargeCalibrationAlgorithm.h:52
Belle2::VxdID::getID
baseType getID() const
Get the unique id.
Definition: VxdID.h:104
Belle2::PXDChargeCalibrationAlgorithm::minClusters
int minClusters
Minimum number of collected clusters for estimating median charge.
Definition: PXDChargeCalibrationAlgorithm.h:40
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
Belle2::PXDClusterChargeMapPar::setContent
void setContent(unsigned short sensorID, unsigned short globalID, float value)
Set map content.
Definition: PXDClusterChargeMapPar.h:65
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::PXDChargeCalibrationAlgorithm::forceContinue
bool forceContinue
Force continue in low statistics runs instead of returning c_NotEnoughData.
Definition: PXDChargeCalibrationAlgorithm.h:49
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::PXDChargeCalibrationAlgorithm::safetyFactor
float safetyFactor
Safety factor for determining whether the collected number of clusters is enough.
Definition: PXDChargeCalibrationAlgorithm.h:46
Belle2::PXDChargeCalibrationAlgorithm::CalculateMedian
double CalculateMedian(std::vector< double > &signals)
Calculate a median from unsorted signal vector. The input vector gets sorted.
Definition: PXDChargeCalibrationAlgorithm.cc:219
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::PXDChargeCalibrationAlgorithm::FitLandau
double FitLandau(std::vector< double > &signals)
calculate MPV of unsorted signal vector using a Landau fit
Definition: PXDChargeCalibrationAlgorithm.cc:235
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::PXDClusterChargeMapPar
The payload class for PXD cluster charge calibrations.
Definition: PXDClusterChargeMapPar.h:40
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::PXDChargeCalibrationAlgorithm::EstimateCharge
double EstimateCharge(VxdID sensorID, unsigned short uBin, unsigned short vBin)
Estimate median charge form collected clusters on part of PXD.
Definition: PXDChargeCalibrationAlgorithm.cc:188
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
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