11 #include <pxd/calibration/PXDGainCalibrationAlgorithm.h>
12 #include <pxd/dbobjects/PXDGainMapPar.h>
13 #include <pxd/dbobjects/PXDClusterChargeMapPar.h>
22 #include <boost/format.hpp>
45 void getNumberOfBins(
const std::shared_ptr<TH1I>& histo_ptr,
unsigned short& nBinsU,
unsigned short& nBinsV)
47 set<unsigned short> uBinSet;
48 set<unsigned short> vBinSet;
51 for (
auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
53 string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
56 istringstream stream(label);
58 getline(stream, token,
'_');
59 getline(stream, token,
'_');
60 unsigned short uBin = std::stoi(token);
62 getline(stream, token,
'_');
63 unsigned short vBin = std::stoi(token);
69 if (uBinSet.empty() || vBinSet.empty()) {
70 B2FATAL(
"Not able to determine the grid size. Something is wrong with collected data.");
72 nBinsU = *uBinSet.rbegin() + 1;
73 nBinsV = *vBinSet.rbegin() + 1;
80 set<unsigned short> sensorSet;
83 for (
auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
85 string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
88 istringstream stream(label);
90 getline(stream, token,
'_');
91 VxdID sensorID(token);
92 sensorSet.insert(sensorID.getID());
94 return sensorSet.size();
100 PXDGainCalibrationAlgorithm::PXDGainCalibrationAlgorithm():
103 minClusters(1000), noiseSigma(0.6), safetyFactor(2.0), forceContinue(false), strategy(0)
106 " -------------------------- PXDGainCalibrationAlgorithm ---------------------------------\n"
108 " Algorithm for estimating pxd gains (conversion factor from charge to ADU) \n"
109 " ----------------------------------------------------------------------------------------\n"
120 auto cluster_counter = getObjectPtr<TH1I>(
"PXDClusterCounter");
126 unsigned short nBinsU = 0;
127 unsigned short nBinsV = 0;
133 B2WARNING(
"Not enough Data: Only " << cluster_counter->GetEntries() <<
" hits were collected but " <<
int(
135 nSensors * nBinsU * nBinsV) <<
" needed!");
138 B2WARNING(
"Continue despite low statistics: Only " << cluster_counter->GetEntries() <<
" hits were collected but" <<
int(
140 nSensors * nBinsU * nBinsV) <<
" would be desirable!");
144 B2INFO(
"Start calibration using a " << nBinsU <<
"x" << nBinsV <<
" grid per sensor.");
145 B2INFO(
"Number of collected clusters is " << cluster_counter->GetEntries());
149 set<VxdID> pxdSensors;
152 for (
auto histoBin = 1; histoBin <= cluster_counter->GetXaxis()->GetNbins(); histoBin++) {
154 string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
157 istringstream stream(label);
159 getline(stream, token,
'_');
160 VxdID sensorID(token);
162 getline(stream, token,
'_');
163 unsigned short uBin = std::stoi(token);
165 getline(stream, token,
'_');
166 unsigned short vBin = std::stoi(token);
169 int numberOfHits = cluster_counter->GetBinContent(histoBin);
180 B2WARNING(label <<
": Number of mc hits too small for fitting (" << numberOfHits <<
" < " <<
minClusters <<
181 "). Use default gain.");
184 pxdSensors.insert(sensorID);
191 for (
const auto& sensorID : pxdSensors) {
192 for (
unsigned short vBin = 0; vBin < nBinsV; ++vBin) {
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);
206 B2RESULT(
"Gain calibration on sensor=" << sensorID <<
" and vBin=" << vBin <<
" was successful on " << nGood <<
"/" << nBinsU <<
210 if (nGood > 0 && nBad > 0) {
212 for (
unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
213 auto gain = gainMapPar->
getContent(sensorID.getID(), uBin, vBin);
215 gainMapPar->
setContent(sensorID.getID(), uBin, vBin, meanGain);
216 B2RESULT(
"Gain calibration on sensor=" << sensorID <<
", vBin=" << vBin <<
" uBin " << uBin <<
": Replace default gain wih average "
229 B2INFO(
"PXD Gain Calibration Successful");
240 const string treename = str(format(
"tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
242 vector<double> mc_signals;
244 auto tree_MC = getObjectPtr<TTree>(treename);
245 tree_MC->SetBranchAddress(
"signal", &m_signal);
248 const auto nEntries_MC = tree_MC->GetEntries();
249 for (
int i = 0; i < nEntries_MC; ++i) {
250 tree_MC->GetEntry(i);
252 double noise = gRandom->Gaus(0.0,
noiseSigma);
253 mc_signals.push_back(m_signal + noise);
257 double mcMedian = -1;
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.");
267 B2FATAL(
"strategy unavailable, use 0 for medians or 1 for landau fit!");
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.");
278 double gain = dataMedian / mcMedian;
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);
284 return gain * gainFromDB;
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);
304 const auto nEntries = dbtree->GetEntries();
305 for (
int i = 0; i < nEntries; ++i) {
315 return sum / counter;
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);
335 const auto nEntries = dbtree->GetEntries();
336 for (
int i = 0; i < nEntries; ++i) {
346 return sum / counter;
352 auto size = signals.size();
357 sort(signals.begin(), signals.end());
359 return (signals[size / 2 - 1] + signals[size / 2]) / 2;
361 return signals[size / 2];
368 auto size = signals.size();
369 if (size == 0)
return 0.0;
372 int max = *max_element(signals.begin(), signals.end());
373 int min = *min_element(signals.begin(), signals.end());
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);
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);
387 Int_t status = hist_signals->Fit(
"landau",
"Lq",
"", 0, 350);
388 double MPV = landau->GetParameter(
"MPV");
395 if (status == 0)
return MPV;
397 B2WARNING(
"Fit failed!. using default value.");