11 #include <pxd/calibration/PXDDataMCGainCalibrationAlgorithm.h>
12 #include <pxd/dbobjects/PXDClusterChargeMapPar.h>
13 #include <pxd/dbobjects/PXDGainMapPar.h>
23 #include <boost/format.hpp>
49 void getNumberOfBins(
const std::shared_ptr<TH1I>& histo_ptr,
unsigned short& nBinsU,
unsigned short& nBinsV)
51 set<unsigned short> uBinSet;
52 set<unsigned short> vBinSet;
55 for (
auto histoBin = 1; histoBin <= histo_ptr->GetNbinsX(); histoBin++) {
57 string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
60 istringstream stream(label);
62 getline(stream, token,
'_');
63 getline(stream, token,
'_');
64 unsigned short uBin = std::stoi(token);
66 getline(stream, token,
'_');
67 unsigned short vBin = std::stoi(token);
73 if (uBinSet.empty() || vBinSet.empty()) {
74 B2FATAL(
"Not able to determine the grid size. Something is wrong with collected data.");
76 nBinsU = *uBinSet.rbegin() + 1;
77 nBinsV = *vBinSet.rbegin() + 1;
84 set<unsigned short> sensorSet;
87 for (
auto histoBin = 1; histoBin <= histo_ptr->GetNbinsX(); histoBin++) {
89 string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
92 istringstream stream(label);
94 getline(stream, token,
'_');
95 VxdID sensorID(token);
96 sensorSet.insert(sensorID.getID());
99 return sensorSet.size();
105 PXDDataMCGainCalibrationAlgorithm::PXDDataMCGainCalibrationAlgorithm():
107 minClusters(5000), noiseSigma(0.6), safetyFactor(2.0), forceContinue(false), strategy(0),
108 doCalibration(false), useChargeHistogram(false), chargePayloadName(
"PXDClusterChargeMapPar")
111 " ------------------------------ PXDDataMCGainCalibrationAlgorithm -----------------------------------\n"
113 " Algorithm for estimating pxd median/MPV cluster charges for different position on sensor in ADU, and optionally calculate gain \n"
114 " ------------------------------------------------------------------------------------------------\n"
124 auto cluster_counter = getObjectPtr<TH1I>(
"PXDClusterCounter");
125 if (!cluster_counter) {
126 B2INFO(
"Not enough Data: cluster counter does not exist ");
134 unsigned short nBinsU = 0;
135 unsigned short nBinsV = 0;
141 B2INFO(
"Not enough Data: Only " << cluster_counter->GetEntries() <<
" hits were collected but " <<
int(
safetyFactor *
minClusters *
142 nSensors * nBinsU * nBinsV) <<
" needed!");
145 B2INFO(
"Continue despite low statistics: Only " << cluster_counter->GetEntries() <<
" hits were collected but" <<
int(
147 nSensors * nBinsU * nBinsV) <<
" would be desirable!");
151 B2INFO(
"Start calibration using a " << nBinsU <<
"x" << nBinsV <<
" grid per sensor.");
157 set<VxdID> pxdSensors;
162 auto dbtree = getObjectPtr<TTree>(
"dbtree");
163 dbtree->SetBranchAddress(
"run", &m_run);
164 dbtree->SetBranchAddress(
"exp", &m_exp);
165 dbtree->SetBranchAddress(
"chargeMap", &chargeMapPtr);
166 dbtree->SetBranchAddress(
"gainMap", &gainMapPtr);
169 for (
auto histoBin = 1; histoBin <= cluster_counter->GetNbinsX(); histoBin++) {
172 string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
175 istringstream stream(label);
177 getline(stream, token,
'_');
178 VxdID sensorID(token);
180 getline(stream, token,
'_');
181 unsigned short uBin = std::stoi(token);
183 getline(stream, token,
'_');
184 unsigned short vBin = std::stoi(token);
187 int numberOfDataHits = cluster_counter->GetBinContent(histoBin);
192 B2INFO(
"start EstimateCharge");
199 B2INFO(
"EstimateCharge: sensor " << sensorID.
getID() <<
" U " << uBin <<
" V " << vBin <<
" Charge " << Charge);
207 B2WARNING(
"Retrieved negative charge for data for sensor=" << sensorID <<
" uBin=" << uBin <<
" vBin=" << vBin <<
208 ". Set gain to default value (=1.0).");
212 double mcCharge = chargeMapPtr->
getContent(sensorID.
getID(), uBin, vBin);
214 if (mcCharge <= 0.0) {
215 B2WARNING(
"Retrieved negative charge for MC from DB for sensor=" << sensorID <<
" uBin=" << uBin <<
" vBin=" << vBin <<
216 ". Set gain to default value (=1.0).");
219 Gain = Charge / mcCharge;
220 B2INFO(
"Estimated Gain: sensor " << sensorID.
getID() <<
" U " << uBin <<
" V " << vBin <<
" Gain " << Gain <<
" = " << Charge <<
232 B2WARNING(label <<
": Number of data hits too small for fitting (" << numberOfDataHits <<
" < " <<
minClusters <<
233 "). Use default value of 0 for charge, 1 for gain.");
239 pxdSensors.insert(sensorID);
243 chargeMapPtr =
nullptr;
244 gainMapPtr =
nullptr;
251 B2INFO(
"PXD Cluster Charge Calibration Successful");
260 for (
const auto& sensorID : pxdSensors) {
261 for (
unsigned short vBin = 0; vBin < nBinsV; ++vBin) {
263 unsigned short nGood = 0;
264 unsigned short nBad = 0;
265 for (
unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
266 auto gain = gainMapPar->
getContent(sensorID.getID(), uBin, vBin);
275 B2RESULT(
"Gain calibration on sensor=" << sensorID <<
" and vBin=" << vBin <<
" was successful on " << nGood <<
"/" << nBinsU <<
279 if (nGood > 0 && nBad > 0) {
281 for (
unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
282 auto gain = gainMapPar->
getContent(sensorID.getID(), uBin, vBin);
284 gainMapPar->
setContent(sensorID.getID(), uBin, vBin, meanGain);
285 B2RESULT(
"Gain calibration on sensor=" << sensorID <<
", vBin=" << vBin <<
" uBin " << uBin <<
": Replace default gain wih average "
297 B2INFO(
"PXD Gain Calibration Successful");
306 unsigned short histoBin)
315 const string treename = str(format(
"tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
317 vector<double> signals;
319 auto tree = getObjectPtr<TTree>(treename);
320 tree->SetBranchAddress(
"signal", &m_signal);
323 const auto nEntries = tree->GetEntries();
326 for (
int i = 0; i < nEntries; i += incr) {
328 double noise = gRandom->Gaus(0.0,
noiseSigma);
329 signals.push_back(m_signal + noise);
335 B2INFO(
"EstimateCharge: sensor " << sensorID.
getID() <<
"(" << layerNumber <<
"," << ladderNumber <<
"," << sensorNumber
336 <<
") U " << uBin <<
" V " << vBin
337 <<
" Charge " << median);
342 double diff = (landaumpv - median);
344 if (landaumpv > 0.) difff = (landaumpv - median) / landaumpv;
345 B2INFO(
"EstimateCharge: sensor " << sensorID.
getID() <<
"(" << layerNumber <<
"," << ladderNumber <<
"," << sensorNumber
346 <<
") U " << uBin <<
" V " << vBin
347 <<
" Charge " << landaumpv <<
" Median " << median <<
" diff = " << diff <<
"/" << difff);
351 for (
auto& each : signals)
353 if (signals.size() > 0) mean = mean / signals.size();
356 B2FATAL(
"strategy unavailable, use 0 for medians, 1 for landau fit and 2 for mean!");
361 auto cluster_counter = getObjectPtr<TH2I>(
"PXDClusterCharge");
362 TH1D* hist_signals = cluster_counter->ProjectionY(
"proj", histoBin, histoBin);
366 B2INFO(
"EstimateCharge: sensor " << sensorID.
getID() <<
"(" << layerNumber <<
"," << ladderNumber <<
"," << sensorNumber
367 <<
") U " << uBin <<
" V " << vBin
368 <<
" Charge " << median);
373 double landaumpv =
FitLandau(hist_signals);
374 double diff = (landaumpv - median);
376 if (landaumpv > 0.) difff = (landaumpv - median) / landaumpv;
377 B2INFO(
"EstimateCharge: sensor " << sensorID.
getID() <<
"(" << layerNumber <<
"," << ladderNumber <<
"," << sensorNumber
378 <<
") U " << uBin <<
" V " << vBin
379 <<
" Charge " << landaumpv <<
" Median " << median <<
" diff = " << diff <<
"/" << difff);
382 return hist_signals->GetMean();
384 B2FATAL(
"strategy unavailable, use 0 for medians, 1 for landau fit and 2 for mean!");
395 auto size = signals.size();
400 sort(signals.begin(), signals.end());
402 return (signals[size / 2 - 1] + signals[size / 2]) / 2;
404 return signals[size / 2];
411 auto size = hist_signals->GetEntries();
413 if (size == 0)
return 0.0;
416 for (
int ibin = 0; ibin < hist_signals->GetNbinsX(); ++ibin) {
417 sum += hist_signals->GetBinContent(ibin + 1);
418 if (sum > size / 2) {
419 return hist_signals->GetBinLowEdge(ibin + 1);
423 B2WARNING(
"Could not find median! using default value 0.0!");
429 auto size = signals.size();
430 if (size == 0)
return 0.0;
433 int max = *max_element(signals.begin(), signals.end());
434 int min = *min_element(signals.begin(), signals.end());
436 if ((max - min) % 2) max--;
437 int nbin = max - min;
440 TH1D* hist_signals =
new TH1D(
"",
"", nbin, min, max);
441 for (
auto it = signals.begin(); it != signals.end(); ++it) {
442 hist_signals->Fill(*it);
446 TF1* landau =
new TF1(
"landau",
"TMath::Landau(x,[0],[1])*[2]", min, max);
447 landau->SetParNames(
"MPV",
"sigma",
"scale");
448 landau->SetParameters(35, 8, 1000);
449 landau->SetParLimits(0, 0., 80.);
454 Int_t status = hist_signals->Fit(
"landau",
"Lq",
"", fitmin, fitmax);
455 double MPV = landau->GetParameter(
"MPV");
457 B2INFO(
"Fit result: " << status <<
" MPV " << MPV <<
" sigma " << landau->GetParameter(
"sigma")
458 <<
" scale " << landau->GetParameter(
"scale") <<
" chi2 " << landau->GetChisquare());
465 if (status == 0)
return MPV;
467 B2WARNING(
"Fit failed! using default value 0.0!");
474 auto size = hist_signals->GetEntries();
475 if (size == 0)
return 0.0;
477 int max = hist_signals->GetBinLowEdge(hist_signals->GetNbinsX() + 1);
478 int min = hist_signals->GetBinLowEdge(1);
481 TF1* landau =
new TF1(
"landau",
"TMath::Landau(x,[0],[1])*[2]", min, max);
482 landau->SetParNames(
"MPV",
"sigma",
"scale");
483 landau->SetParameters(35, 8, 1000);
484 landau->SetParLimits(0, 0., 80.);
489 Int_t status = hist_signals->Fit(
"landau",
"Lq",
"", fitmin, fitmax);
490 double MPV = landau->GetParameter(
"MPV");
492 B2INFO(
"Fit result: " << status <<
" MPV " << MPV <<
" sigma " << landau->GetParameter(
"sigma")
493 <<
" scale " << landau->GetParameter(
"scale") <<
" chi2 " << landau->GetChisquare());
499 if (status == 0)
return MPV;
501 B2WARNING(
"Fit failed! using default value 0.0!");