Belle II Software  release-05-02-19
PXDGainCalibrationAlgorithm.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, Jonas Roetter *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <pxd/calibration/PXDGainCalibrationAlgorithm.h>
12 #include <pxd/dbobjects/PXDGainMapPar.h>
13 #include <pxd/dbobjects/PXDClusterChargeMapPar.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 //ROOT
25 #include <TRandom.h>
26 #include <TH1.h>
27 #include <TF1.h>
28 
29 using namespace std;
30 using boost::format;
31 using namespace Belle2;
32 
33 
34 // Anonymous namespace for data objects used by PXDGainCalibrationAlgorithm class
35 namespace {
36 
38  int m_signal;
40  int m_run;
42  int m_exp;
43 
45  void getNumberOfBins(const std::shared_ptr<TH1I>& histo_ptr, unsigned short& nBinsU, unsigned short& nBinsV)
46  {
47  set<unsigned short> uBinSet;
48  set<unsigned short> vBinSet;
49 
50  // Loop over all bins of input histo
51  for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
52  // The bin label contains the vxdid, uBin and vBin
53  string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
54 
55  // Parse label string format to read sensorID, uBin and vBin
56  istringstream stream(label);
57  string token;
58  getline(stream, token, '_');
59  getline(stream, token, '_');
60  unsigned short uBin = std::stoi(token);
61 
62  getline(stream, token, '_');
63  unsigned short vBin = std::stoi(token);
64 
65  uBinSet.insert(uBin);
66  vBinSet.insert(vBin);
67  }
68 
69  if (uBinSet.empty() || vBinSet.empty()) {
70  B2FATAL("Not able to determine the grid size. Something is wrong with collected data.");
71  } else {
72  nBinsU = *uBinSet.rbegin() + 1;
73  nBinsV = *vBinSet.rbegin() + 1;
74  }
75  }
76 
78  unsigned short getNumberOfSensors(const std::shared_ptr<TH1I>& histo_ptr)
79  {
80  set<unsigned short> sensorSet;
81 
82  // Loop over all bins of input histo
83  for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
84  // The bin label contains the vxdid, uBin and vBin
85  string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
86 
87  // Parse label string format to read sensorID, uBin and vBin
88  istringstream stream(label);
89  string token;
90  getline(stream, token, '_');
91  VxdID sensorID(token);
92  sensorSet.insert(sensorID.getID());
93  }
94  return sensorSet.size();
95  }
96 
97 }
98 
99 
100 PXDGainCalibrationAlgorithm::PXDGainCalibrationAlgorithm():
101  CalibrationAlgorithm("PXDClusterChargeCollector"),
102 
103  minClusters(1000), noiseSigma(0.6), safetyFactor(2.0), forceContinue(false), strategy(0)
104 {
106  " -------------------------- PXDGainCalibrationAlgorithm ---------------------------------\n"
107  " \n"
108  " Algorithm for estimating pxd gains (conversion factor from charge to ADU) \n"
109  " ----------------------------------------------------------------------------------------\n"
110  );
111 }
112 
113 
114 
115 
117 {
118 
119  // Get counter histograms
120  auto cluster_counter = getObjectPtr<TH1I>("PXDClusterCounter");
121 
122  // Extract number of sensors from counter histograms
123  auto nSensors = getNumberOfSensors(cluster_counter);
124 
125  // Extract the number of grid bins from counter histograms
126  unsigned short nBinsU = 0;
127  unsigned short nBinsV = 0;
128  getNumberOfBins(cluster_counter, nBinsU, nBinsV);
129 
130  // Check that we have collected enough Data
131  if (cluster_counter->GetEntries() < int(safetyFactor * minClusters * nSensors * nBinsU * nBinsV)) {
132  if (not forceContinue) {
133  B2WARNING("Not enough Data: Only " << cluster_counter->GetEntries() << " hits were collected but " << int(
135  nSensors * nBinsU * nBinsV) << " needed!");
136  return c_NotEnoughData;
137  } else {
138  B2WARNING("Continue despite low statistics: Only " << cluster_counter->GetEntries() << " hits were collected but" << int(
140  nSensors * nBinsU * nBinsV) << " would be desirable!");
141  }
142  }
143 
144  B2INFO("Start calibration using a " << nBinsU << "x" << nBinsV << " grid per sensor.");
145  B2INFO("Number of collected clusters is " << cluster_counter->GetEntries());
146 
147  // This is the PXD gain correction payload for conditions DB
148  PXDGainMapPar* gainMapPar = new PXDGainMapPar(nBinsU, nBinsV);
149  set<VxdID> pxdSensors;
150 
151  // Loop over all bins of input histo
152  for (auto histoBin = 1; histoBin <= cluster_counter->GetXaxis()->GetNbins(); histoBin++) {
153  // The bin label contains the vxdid, uBin and vBin
154  string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
155 
156  // Parse label string format to read sensorID, uBin and vBin
157  istringstream stream(label);
158  string token;
159  getline(stream, token, '_');
160  VxdID sensorID(token);
161 
162  getline(stream, token, '_');
163  unsigned short uBin = std::stoi(token);
164 
165  getline(stream, token, '_');
166  unsigned short vBin = std::stoi(token);
167 
168  // Read back the counters for number of collected clusters
169  int numberOfHits = cluster_counter->GetBinContent(histoBin);
170 
171  // Only perform fitting, when enough data is available
172  if (numberOfHits >= minClusters) {
173 
174  // Estimate the gain on a certain part of PXD
175  auto gain = EstimateGain(sensorID, uBin, vBin);
176 
177  // Store the gain
178  gainMapPar->setContent(sensorID.getID(), uBin, vBin, gain);
179  } else {
180  B2WARNING(label << ": Number of mc hits too small for fitting (" << numberOfHits << " < " << minClusters <<
181  "). Use default gain.");
182  gainMapPar->setContent(sensorID.getID(), uBin, vBin, 1.0);
183  }
184  pxdSensors.insert(sensorID);
185  }
186 
187  // Post processing of gain map. It is possible that the gain
188  // computation failed on some parts. Here, we replace default
189  // values (1.0) by local averages of neighboring sensor parts.
190 
191  for (const auto& sensorID : pxdSensors) {
192  for (unsigned short vBin = 0; vBin < nBinsV; ++vBin) {
193  float meanGain = 0;
194  unsigned short nGood = 0;
195  unsigned short nBad = 0;
196  for (unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
197  auto gain = gainMapPar->getContent(sensorID.getID(), uBin, vBin);
198  // Filter default gains
199  if (gain != 1.0) {
200  nGood += 1;
201  meanGain += gain;
202  } else {
203  nBad += 1;
204  }
205  }
206  B2RESULT("Gain calibration on sensor=" << sensorID << " and vBin=" << vBin << " was successful on " << nGood << "/" << nBinsU <<
207  " uBins.");
208 
209  // Check if we can repair bad calibrations with a local avarage
210  if (nGood > 0 && nBad > 0) {
211  meanGain /= nGood;
212  for (unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
213  auto gain = gainMapPar->getContent(sensorID.getID(), uBin, vBin);
214  if (gain == 1.0) {
215  gainMapPar->setContent(sensorID.getID(), uBin, vBin, meanGain);
216  B2RESULT("Gain calibration on sensor=" << sensorID << ", vBin=" << vBin << " uBin " << uBin << ": Replace default gain wih average "
217  << meanGain);
218  }
219  }
220  }
221  }
222  }
223 
224 
225  // Save the gain map to database. Note that this will set the database object name to the same as the collector but you
226  // are free to change it.
227  saveCalibration(gainMapPar, "PXDGainMapPar");
228 
229  B2INFO("PXD Gain Calibration Successful");
230  return c_OK;
231 }
232 
233 
234 double PXDGainCalibrationAlgorithm::EstimateGain(VxdID sensorID, unsigned short uBin, unsigned short vBin)
235 {
236  // Construct a tree name for requested part of PXD
237  auto layerNumber = sensorID.getLayerNumber();
238  auto ladderNumber = sensorID.getLadderNumber();
239  auto sensorNumber = sensorID.getSensorNumber();
240  const string treename = str(format("tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
241  // Vector with cluster signals from collected mc
242  vector<double> mc_signals;
243 
244  auto tree_MC = getObjectPtr<TTree>(treename);
245  tree_MC->SetBranchAddress("signal", &m_signal);
246 
247  // Loop over tree_MC
248  const auto nEntries_MC = tree_MC->GetEntries();
249  for (int i = 0; i < nEntries_MC; ++i) {
250  tree_MC->GetEntry(i);
251 
252  double noise = gRandom->Gaus(0.0, noiseSigma);
253  mc_signals.push_back(m_signal + noise);
254  }
255 
256  double dataMedian = GetChargeMedianFromDB(sensorID, uBin, vBin);
257  double mcMedian = -1;
258  // check if dataMedian makes sense
259  if (dataMedian <= 0.0) {
260  B2WARNING("Retrieved negative charge median from DB for sensor=" << sensorID << " uBin=" << uBin << " vBin=" << vBin <<
261  ". Set gain to default value (=1.0) as well.");
262  return 1.0;
263  }
264  if (strategy == 0) mcMedian = CalculateMedian(mc_signals);
265  else if (strategy == 1) mcMedian = FitLandau(mc_signals);
266  else {
267  B2FATAL("strategy unavailable, use 0 for medians or 1 for landau fit!");
268  }
269 
270  // check if mcMedian makes sense
271  if (mcMedian <= 0.0) {
272  B2WARNING("Retrieved negative charge median from DB for sensor=" << sensorID << " uBin=" << uBin << " vBin=" << vBin <<
273  ". Set gain to default value (=1.0) as well.");
274  return 1.0;
275  }
276 
277 
278  double gain = dataMedian / mcMedian;
279  double gainFromDB = GetCurrentGainFromDB(sensorID, uBin, vBin);
280  B2DEBUG(10, "Gain from db used in PXDDigitizer is " << gainFromDB);
281  B2DEBUG(10, "New gain correction derived is " << gain);
282  B2DEBUG(10, "The total gain we should return is " << gain * gainFromDB);
283 
284  return gain * gainFromDB;
285 }
286 
287 double PXDGainCalibrationAlgorithm::GetChargeMedianFromDB(VxdID sensorID, unsigned short uBin, unsigned short vBin)
288 {
289  // Read back db payloads
290  PXDClusterChargeMapPar* chargeMapPtr = 0;
291  PXDGainMapPar* gainMapPtr = 0;
292 
293  auto dbtree = getObjectPtr<TTree>("dbtree");
294  dbtree->SetBranchAddress("run", &m_run);
295  dbtree->SetBranchAddress("exp", &m_exp);
296  dbtree->SetBranchAddress("chargeMap", &chargeMapPtr);
297  dbtree->SetBranchAddress("gainMap", &gainMapPtr);
298 
299  // Compute running average of charge medians from db
300  double sum = 0;
301  int counter = 0;
302 
303  // Loop over dbtree
304  const auto nEntries = dbtree->GetEntries();
305  for (int i = 0; i < nEntries; ++i) {
306  dbtree->GetEntry(i);
307  sum += chargeMapPtr->getContent(sensorID.getID(), uBin, vBin);
308  counter += 1;
309  }
310  delete chargeMapPtr;
311  chargeMapPtr = 0;
312  delete gainMapPtr;
313  gainMapPtr = 0;
314 
315  return sum / counter;
316 }
317 
318 double PXDGainCalibrationAlgorithm::GetCurrentGainFromDB(VxdID sensorID, unsigned short uBin, unsigned short vBin)
319 {
320  // Read back db payloads
321  PXDClusterChargeMapPar* chargeMapPtr = 0;
322  PXDGainMapPar* gainMapPtr = 0;
323 
324  auto dbtree = getObjectPtr<TTree>("dbtree");
325  dbtree->SetBranchAddress("run", &m_run);
326  dbtree->SetBranchAddress("exp", &m_exp);
327  dbtree->SetBranchAddress("chargeMap", &chargeMapPtr);
328  dbtree->SetBranchAddress("gainMap", &gainMapPtr);
329 
330  // Compute running average of gains from db
331  double sum = 0;
332  int counter = 0;
333 
334  // Loop over dbtree
335  const auto nEntries = dbtree->GetEntries();
336  for (int i = 0; i < nEntries; ++i) {
337  dbtree->GetEntry(i);
338  sum += gainMapPtr->getContent(sensorID.getID(), uBin, vBin);
339  counter += 1;
340  }
341  delete chargeMapPtr;
342  chargeMapPtr = 0;
343  delete gainMapPtr;
344  gainMapPtr = 0;
345 
346  return sum / counter;
347 }
348 
349 
350 double PXDGainCalibrationAlgorithm::CalculateMedian(vector<double>& signals)
351 {
352  auto size = signals.size();
353 
354  if (size == 0) {
355  return 0.0; // Undefined, really.
356  } else {
357  sort(signals.begin(), signals.end());
358  if (size % 2 == 0) {
359  return (signals[size / 2 - 1] + signals[size / 2]) / 2;
360  } else {
361  return signals[size / 2];
362  }
363  }
364 }
365 
366 double PXDGainCalibrationAlgorithm::FitLandau(vector<double>& signals)
367 {
368  auto size = signals.size();
369  if (size == 0) return 0.0; // Undefined, really.
370 
371  // get max and min values of signal vector
372  int max = *max_element(signals.begin(), signals.end());
373  int min = *min_element(signals.begin(), signals.end());
374 
375  // create histogram to hold signals and fill it
376  TH1D* hist_signals = new TH1D("", "", max - min, min, max);
377  for (auto it = signals.begin(); it != signals.end(); ++it) {
378  hist_signals->Fill(*it);
379  }
380 
381  // create fit function
382  TF1* landau = new TF1("landau", "TMath::Landau(x,[0],[1])*[2]", min, max);
383  landau->SetParNames("MPV", "sigma", "scale");
384  landau->SetParameters(100, 1, 1000);
385 
386  // do fit and get results
387  Int_t status = hist_signals->Fit("landau", "Lq", "", 0, 350);
388  double MPV = landau->GetParameter("MPV");
389 
390  // clean up
391  delete hist_signals;
392  delete landau;
393 
394  // check fit status
395  if (status == 0) return MPV;
396  else {
397  B2WARNING("Fit failed!. using default value.");
398  return 0.0;
399  }
400 }
Belle2::PXDGainCalibrationAlgorithm::EstimateGain
double EstimateGain(VxdID sensorID, unsigned short uBin, unsigned short vBin)
Estimate gain as ratio of medians from MC and data for a part of PXD.
Definition: PXDGainCalibrationAlgorithm.cc:234
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::PXDGainCalibrationAlgorithm::minClusters
int minClusters
Minimum number of collected clusters for estimating gains.
Definition: PXDGainCalibrationAlgorithm.h:40
Belle2::PXDGainCalibrationAlgorithm::FitLandau
double FitLandau(std::vector< double > &signals)
Calculate MPV from signal vector using a landau fit.
Definition: PXDGainCalibrationAlgorithm.cc:366
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::VxdID::getID
baseType getID() const
Get the unique id.
Definition: VxdID.h:104
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::PXDGainMapPar
The payload class for PXD gain corrections.
Definition: PXDGainMapPar.h:53
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
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::PXDGainCalibrationAlgorithm::noiseSigma
float noiseSigma
Artificial noise sigma for smearing cluster charge.
Definition: PXDGainCalibrationAlgorithm.h:43
Belle2::PXDClusterChargeMapPar::getContent
float getContent(unsigned short sensorID, unsigned short globalID) const
Get content.
Definition: PXDClusterChargeMapPar.h:100
Belle2::PXDGainMapPar::getContent
float getContent(unsigned short sensorID, unsigned short globalID) const
Get content.
Definition: PXDGainMapPar.h:113
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::PXDGainCalibrationAlgorithm::calibrate
virtual EResult calibrate() override
Run algo on data.
Definition: PXDGainCalibrationAlgorithm.cc:116
Belle2::PXDGainCalibrationAlgorithm::safetyFactor
float safetyFactor
Safety factor for determining whether the collected number of clusters is enough.
Definition: PXDGainCalibrationAlgorithm.h:46
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::PXDGainMapPar::setContent
void setContent(unsigned short sensorID, unsigned short globalID, float value)
Set map content.
Definition: PXDGainMapPar.h:78
Belle2::PXDGainCalibrationAlgorithm::strategy
int strategy
strategy to used for gain calibration, 0 for medians, 1 for landau fit
Definition: PXDGainCalibrationAlgorithm.h:52
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::PXDGainCalibrationAlgorithm::forceContinue
bool forceContinue
Force continue in low statistics runs instead of returning c_NotEnoughData.
Definition: PXDGainCalibrationAlgorithm.h:49
Belle2::PXDGainCalibrationAlgorithm::GetChargeMedianFromDB
double GetChargeMedianFromDB(VxdID sensorID, unsigned short uBin, unsigned short vBin)
Retrive charge median value from pulled in data base payload.
Definition: PXDGainCalibrationAlgorithm.cc:287
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::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
Belle2::PXDGainCalibrationAlgorithm::GetCurrentGainFromDB
double GetCurrentGainFromDB(VxdID sensorID, unsigned short uBin, unsigned short vBin)
Retrive current gain value from pulled in data base payload.
Definition: PXDGainCalibrationAlgorithm.cc:318
Belle2::PXDGainCalibrationAlgorithm::CalculateMedian
double CalculateMedian(std::vector< double > &signals)
Calculate a median from unsorted signal vector. The input vector gets sorted.
Definition: PXDGainCalibrationAlgorithm.cc:350
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