Belle II Software  release-05-01-25
PXDAnalyticGainCalibrationAlgorithm.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 <pxd/calibration/PXDAnalyticGainCalibrationAlgorithm.h>
12 #include <pxd/calibration/PXDCalibrationUtilities.h>
13 #include <pxd/dbobjects/PXDGainMapPar.h>
14 
15 #include <boost/format.hpp>
16 
17 //ROOT
18 #include <TH2F.h>
19 
20 using namespace std;
21 using boost::format;
22 using namespace Belle2;
23 using namespace Belle2::PXD;
24 
25 
26 // Anonymous namespace for data objects used by PXDAnalyticGainCalibrationAlgorithm class
27 namespace {
28 
30  int m_signal;
32  float m_estimated;
34  int m_run;
36  int m_exp;
37 }
38 
39 
40 PXDAnalyticGainCalibrationAlgorithm::PXDAnalyticGainCalibrationAlgorithm():
41  CalibrationAlgorithm("PXDPerformanceCollector"),
42  minClusters(1000), safetyFactor(2.0), forceContinue(false), strategy(0), useChargeRatioHistogram(true), correctForward(false)
43 {
45  " -------------------------- PXDAnalyticGainCalibrationAlgorithm ---------------------------------\n"
46  " \n"
47  " Algorithm for estimating pxd gains (conversion factor from charge to ADU) \n"
48  " ----------------------------------------------------------------------------------------\n"
49  );
50 }
51 
52 
54 {
55 
56  // Get counter histograms
57  auto cluster_counter = getObjectPtr<TH1I>("PXDClusterCounter");
58 
59  // Extract number of sensors from counter histograms
60  auto nSensors = getNumberOfSensors(cluster_counter);
61 
62  // Extract the number of grid bins from counter histograms
63  unsigned short nBinsU = 0;
64  unsigned short nBinsV = 0;
65  getNumberOfBins(cluster_counter, nBinsU, nBinsV);
66 
67  // Check that we have collected enough Data
68  if (cluster_counter->GetEntries() < int(safetyFactor * minClusters * nSensors * nBinsU * nBinsV)) {
69  if (not forceContinue) {
70  B2WARNING("Not enough Data: Only " << cluster_counter->GetEntries() << " hits were collected but " << int(
72  nSensors * nBinsU * nBinsV) << " needed!");
73  return c_NotEnoughData;
74  } else {
75  B2WARNING("Continue despite low statistics: Only " << cluster_counter->GetEntries() << " hits were collected but" << int(
77  nSensors * nBinsU * nBinsV) << " would be desirable!");
78  }
79  }
80 
81  B2INFO("Start calibration using a " << nBinsU << "x" << nBinsV << " grid per sensor.");
82  B2INFO("Number of collected clusters is " << cluster_counter->GetEntries());
83 
84  // This is the PXD gain correction payload for conditions DB
85  PXDGainMapPar* gainMapPar = new PXDGainMapPar(nBinsU, nBinsV);
86  set<VxdID> pxdSensors;
87 
88  // Loop over all bins of input histo
89  for (auto histoBin = 1; histoBin <= cluster_counter->GetXaxis()->GetNbins(); histoBin++) {
90  // The bin label contains the vxdid, uBin and vBin
91  string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
92 
93  // Parse label string format to read sensorID, uBin and vBin
94  istringstream stream(label);
95  string token;
96  getline(stream, token, '_');
97  VxdID sensorID(token);
98 
99  getline(stream, token, '_');
100  unsigned short uBin = std::stoi(token);
101 
102  getline(stream, token, '_');
103  unsigned short vBin = std::stoi(token);
104 
105  // Read back the counters for number of collected clusters
106  int numberOfHits = cluster_counter->GetBinContent(histoBin);
107 
108  // Only perform estimation, when enough data is available
109  if (numberOfHits >= minClusters) {
110 
111  double gain = 1.0;
112  // Estimate the gain on a certain part of PXD
114  auto hClusterChargeRatio = getObjectPtr<TH2F>("PXDClusterChargeRatio");
115  TH1D* hRatios = hClusterChargeRatio->ProjectionY("proj", histoBin, histoBin);
116  gain = EstimateGain(sensorID, uBin, vBin, hRatios);
117  } else {
118  gain = EstimateGain(sensorID, uBin, vBin);
119  }
120  // Store the gain
121  gainMapPar->setContent(sensorID.getID(), uBin, vBin, gain);
122  } else {
123  B2WARNING(label << ": Number of hits is too small (" << numberOfHits << " < " << minClusters <<
124  "). Use default gain.");
125  gainMapPar->setContent(sensorID.getID(), uBin, vBin, 1.0);
126  }
127  pxdSensors.insert(sensorID);
128  }
129 
130  // Post processing of gain map. It is possible that the gain
131  // computation failed on some parts. Here, we replace default
132  // values (1.0) by local averages of neighboring sensor parts.
133 
134  for (const auto& sensorID : pxdSensors) {
135 
136  // Special treatement for the last 2 vBin as Bhabha 2-track events
137  // have no enough statistics there if nBinsV = 6
138  if (correctForward && sensorID.getSensorNumber() == 1)
139  for (unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
140  // Search for a vaid gain along v
141  double gainForwardRegion = 1.0;
142  unsigned short vBinToCheck = nBinsV - 1;
143  for (unsigned short vBinGood = nBinsV - 2; vBinGood >= 1; --vBinGood) {
144  auto temp = gainMapPar->getContent(sensorID.getID(), uBin, vBinGood);
145  if (temp != 1.0) {
146  gainForwardRegion = temp;
147  vBinToCheck = vBinGood + 1;
148  break;
149  }
150  }
151  // loop part of the forward regions and check values
152  if (gainForwardRegion != 1.0)
153  for (unsigned short vBin = nBinsV - 1; vBin >= vBinToCheck; --vBin) {
154  auto gain = gainMapPar->getContent(sensorID.getID(), uBin, vBin);
155  if (gain == 1.0) {
156  gainMapPar->setContent(sensorID.getID(), uBin, vBin, gainForwardRegion);
157  B2RESULT("Gain calibration on sensor=" << sensorID << ", vBin=" << vBin << " uBin " << uBin <<
158  ": Replace default gain with that from the closest vBin with non default value "
159  << gainForwardRegion);
160  }
161  }
162  }
163 
164  // general value check
165  for (unsigned short vBin = 0; vBin < nBinsV; ++vBin) {
166 
167  float meanGain = 0;
168  unsigned short nGood = 0;
169  unsigned short nBad = 0;
170  for (unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
171  auto gain = gainMapPar->getContent(sensorID.getID(), uBin, vBin);
172  // Filter default gains
173  if (gain != 1.0) {
174  nGood += 1;
175  meanGain += gain;
176  } else {
177  nBad += 1;
178  }
179  }
180  B2RESULT("Gain calibration on sensor=" << sensorID << " and vBin=" << vBin << " was successful on " << nGood << "/" << nBinsU <<
181  " uBins.");
182 
183  // Check if we can repair bad calibrations with a local avarage
184  if (nGood > 0 && nBad > 0) {
185  meanGain /= nGood;
186  for (unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
187  auto gain = gainMapPar->getContent(sensorID.getID(), uBin, vBin);
188  if (gain == 1.0) {
189  gainMapPar->setContent(sensorID.getID(), uBin, vBin, meanGain);
190  B2RESULT("Gain calibration on sensor=" << sensorID << ", vBin=" << vBin << " uBin " << uBin <<
191  ": Replace default gain with average "
192  << meanGain << " of uBins with non-default gains.");
193  }
194  }
195  }
196  }
197  }
198 
199 
200  // Save the gain map to database. Note that this will set the database object name to the same as the collector but you
201  // are free to change it.
202  saveCalibration(gainMapPar, "PXDGainMapPar");
203 
204  B2INFO("PXD Gain Calibration Successful");
205  return c_OK;
206 }
207 
208 
209 double PXDAnalyticGainCalibrationAlgorithm::EstimateGain(VxdID sensorID, unsigned short uBin, unsigned short vBin, TH1* hist)
210 {
211  double gain = 1.0; // default gain
212 
213  // function pointers for different strategy
214  double (*estimateGainFromHist)(TH1*);
215  double (*estimateGainFromVec)(std::vector<double>&);
216  if (strategy == 0) {
217  estimateGainFromVec = &CalculateMedian;
218  estimateGainFromHist = &CalculateMedian;
219  } else if (strategy == 1) {
220  estimateGainFromVec = &FitLandau;
221  estimateGainFromHist = &FitLandau;
222  } else {
223  B2FATAL("strategy unavailable, use 0 for medians or 1 for landau fit!");
224  }
225 
226  // Do estimation
227  if (hist) { // estimate gain from existing histogram
228  gain = estimateGainFromHist(hist);
229  } else { // estimate from TTree
230  // Construct a tree name for requested part of PXD
231  auto layerNumber = sensorID.getLayerNumber();
232  auto ladderNumber = sensorID.getLadderNumber();
233  auto sensorNumber = sensorID.getSensorNumber();
234  const string treename = str(format("tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
235  // Vector with ratios (cluster charge to its estimation)
236  vector<double> ratios;
237 
238  auto tree = getObjectPtr<TTree>(treename);
239  tree->SetBranchAddress("signal", &m_signal);
240  tree->SetBranchAddress("estimated", &m_estimated);
241 
242  // Loop over tree
243  const auto nEntries_MC = tree->GetEntries();
244  for (int i = 0; i < nEntries_MC; ++i) {
245  tree->GetEntry(i);
246  ratios.push_back(double(m_signal) / m_estimated);
247  }
248  gain = estimateGainFromVec(ratios);
249  }
250 
251  // check if gain makes sense
252  if (gain <= 0.0) {
253  B2WARNING("Retrieved negative median/MPV for sensor=" << sensorID << " uBin=" << uBin << " vBin=" << vBin <<
254  ". Set gain to default value (=1.0) as well.");
255  return 1.0;
256  }
257 
258  // calculate and return the absolute gain
259  double gainFromDB = GetCurrentGainFromDB(sensorID, uBin, vBin);
260  B2DEBUG(10, "Gain from db used in PXDDigitizer is " << gainFromDB);
261  B2DEBUG(10, "New gain correction derived is " << gain);
262  B2DEBUG(10, "The total gain we should return is " << gain * gainFromDB);
263 
264  return gain * gainFromDB;
265 }
266 
267 double PXDAnalyticGainCalibrationAlgorithm::GetCurrentGainFromDB(VxdID sensorID, unsigned short uBin, unsigned short vBin)
268 {
269  // Read back db payloads
270  PXDGainMapPar* gainMapPtr = nullptr;
271 
272  auto dbtree = getObjectPtr<TTree>("dbtree");
273  dbtree->SetBranchAddress("run", &m_run);
274  dbtree->SetBranchAddress("exp", &m_exp);
275  dbtree->SetBranchAddress("gainMap", &gainMapPtr);
276 
277  // Compute running average of gains from db
278  double sum = 0;
279  int counter = 0;
280 
281  // Loop over dbtree
282  const auto nEntries = dbtree->GetEntries();
283  for (int i = 0; i < nEntries; ++i) {
284  dbtree->GetEntry(i);
285  sum += gainMapPtr->getContent(sensorID.getID(), uBin, vBin);
286  counter += 1;
287  }
288  delete gainMapPtr;
289  gainMapPtr = nullptr;
290 
291  return sum / counter;
292 }
293 
294 bool PXDAnalyticGainCalibrationAlgorithm::isBoundaryRequired(const Calibration::ExpRun& /*currentRun*/)
295 {
296  // First run in data as we iterate, but our boundaries weren't set manually already?
297  // Just set the first run to be a boundary and we are done.
298  if (m_boundaries.empty()) {
299  B2INFO("This is the first run encountered, let's say it is a boundary.");
300  return true;
301  } else {
302  return false;
303  }
304 }
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
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::PXD::FitLandau
double FitLandau(TH1 *hist)
Helper function to estimate MPV from 1D histogram.
Definition: PXDCalibrationUtilities.cc:127
Belle2::PXDAnalyticGainCalibrationAlgorithm::strategy
int strategy
strategy to used for gain calibration, 0 for medians, 1 for landau fit
Definition: PXDAnalyticGainCalibrationAlgorithm.h:55
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:108
Belle2::PXDAnalyticGainCalibrationAlgorithm::correctForward
bool correctForward
Flag to update default gains in forward region due to low statistics.
Definition: PXDAnalyticGainCalibrationAlgorithm.h:61
Belle2::PXDAnalyticGainCalibrationAlgorithm::isBoundaryRequired
virtual bool isBoundaryRequired(const Calibration::ExpRun &) override
Decide if a run should be a payload boundary. Only used in certain Python Algorithm Starategies.
Definition: PXDAnalyticGainCalibrationAlgorithm.cc:294
Belle2::VxdID::getID
baseType getID() const
Get the unique id.
Definition: VxdID.h:104
Belle2::PXDAnalyticGainCalibrationAlgorithm::GetCurrentGainFromDB
double GetCurrentGainFromDB(VxdID sensorID, unsigned short uBin, unsigned short vBin)
Retrive current gain value from pulled in data base payload.
Definition: PXDAnalyticGainCalibrationAlgorithm.cc:267
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::PXDAnalyticGainCalibrationAlgorithm::minClusters
int minClusters
Minimum number of collected clusters for estimating gains.
Definition: PXDAnalyticGainCalibrationAlgorithm.h:46
Belle2::PXDAnalyticGainCalibrationAlgorithm::calibrate
virtual EResult calibrate() override
Run algo on data.
Definition: PXDAnalyticGainCalibrationAlgorithm.cc:53
Belle2::PXDAnalyticGainCalibrationAlgorithm::safetyFactor
float safetyFactor
Safety factor for determining whether the collected number of clusters is enough.
Definition: PXDAnalyticGainCalibrationAlgorithm.h:49
Belle2::PXD::CalculateMedian
double CalculateMedian(std::vector< double > &signals)
Helper function to calculate a median from unsorted signal vector.
Definition: PXDCalibrationUtilities.cc:90
Belle2::PXD
Namespace to encapsulate code needed for simulation and reconstrucion of the PXD.
Definition: PXDCalibrationUtilities.h:28
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::PXDAnalyticGainCalibrationAlgorithm::forceContinue
bool forceContinue
Force continue in low statistics runs instead of returning c_NotEnoughData.
Definition: PXDAnalyticGainCalibrationAlgorithm.h:52
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::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
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::PXDAnalyticGainCalibrationAlgorithm::EstimateGain
double EstimateGain(VxdID sensorID, unsigned short uBin, unsigned short vBin, TH1 *hist=nullptr)
Estimate gain as ratio of medians from MC and data for a part of PXD.
Definition: PXDAnalyticGainCalibrationAlgorithm.cc:209
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
Belle2::PXDAnalyticGainCalibrationAlgorithm::useChargeRatioHistogram
bool useChargeRatioHistogram
Flag to use histogram of charge ratio (relative to expected MPV)
Definition: PXDAnalyticGainCalibrationAlgorithm.h:58
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