11 #include <pxd/calibration/PXDAnalyticGainCalibrationAlgorithm.h>
12 #include <pxd/calibration/PXDCalibrationUtilities.h>
13 #include <pxd/dbobjects/PXDGainMapPar.h>
15 #include <boost/format.hpp>
40 PXDAnalyticGainCalibrationAlgorithm::PXDAnalyticGainCalibrationAlgorithm():
42 minClusters(1000), safetyFactor(2.0), forceContinue(false), strategy(0), useChargeRatioHistogram(true), correctForward(false)
45 " -------------------------- PXDAnalyticGainCalibrationAlgorithm ---------------------------------\n"
47 " Algorithm for estimating pxd gains (conversion factor from charge to ADU) \n"
48 " ----------------------------------------------------------------------------------------\n"
57 auto cluster_counter = getObjectPtr<TH1I>(
"PXDClusterCounter");
63 unsigned short nBinsU = 0;
64 unsigned short nBinsV = 0;
70 B2WARNING(
"Not enough Data: Only " << cluster_counter->GetEntries() <<
" hits were collected but " <<
int(
72 nSensors * nBinsU * nBinsV) <<
" needed!");
75 B2WARNING(
"Continue despite low statistics: Only " << cluster_counter->GetEntries() <<
" hits were collected but" <<
int(
77 nSensors * nBinsU * nBinsV) <<
" would be desirable!");
81 B2INFO(
"Start calibration using a " << nBinsU <<
"x" << nBinsV <<
" grid per sensor.");
82 B2INFO(
"Number of collected clusters is " << cluster_counter->GetEntries());
86 set<VxdID> pxdSensors;
89 for (
auto histoBin = 1; histoBin <= cluster_counter->GetXaxis()->GetNbins(); histoBin++) {
91 string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
94 istringstream stream(label);
96 getline(stream, token,
'_');
97 VxdID sensorID(token);
99 getline(stream, token,
'_');
100 unsigned short uBin = std::stoi(token);
102 getline(stream, token,
'_');
103 unsigned short vBin = std::stoi(token);
106 int numberOfHits = cluster_counter->GetBinContent(histoBin);
114 auto hClusterChargeRatio = getObjectPtr<TH2F>(
"PXDClusterChargeRatio");
115 TH1D* hRatios = hClusterChargeRatio->ProjectionY(
"proj", histoBin, histoBin);
123 B2WARNING(label <<
": Number of hits is too small (" << numberOfHits <<
" < " <<
minClusters <<
124 "). Use default gain.");
127 pxdSensors.insert(sensorID);
134 for (
const auto& sensorID : pxdSensors) {
139 for (
unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
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);
146 gainForwardRegion = temp;
147 vBinToCheck = vBinGood + 1;
152 if (gainForwardRegion != 1.0)
153 for (
unsigned short vBin = nBinsV - 1; vBin >= vBinToCheck; --vBin) {
154 auto gain = gainMapPar->
getContent(sensorID.getID(), uBin, vBin);
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);
165 for (
unsigned short vBin = 0; vBin < nBinsV; ++vBin) {
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);
180 B2RESULT(
"Gain calibration on sensor=" << sensorID <<
" and vBin=" << vBin <<
" was successful on " << nGood <<
"/" << nBinsU <<
184 if (nGood > 0 && nBad > 0) {
186 for (
unsigned short uBin = 0; uBin < nBinsU; ++uBin) {
187 auto gain = gainMapPar->
getContent(sensorID.getID(), uBin, vBin);
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.");
204 B2INFO(
"PXD Gain Calibration Successful");
214 double (*estimateGainFromHist)(TH1*);
215 double (*estimateGainFromVec)(std::vector<double>&);
223 B2FATAL(
"strategy unavailable, use 0 for medians or 1 for landau fit!");
228 gain = estimateGainFromHist(hist);
234 const string treename = str(format(
"tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
236 vector<double> ratios;
238 auto tree = getObjectPtr<TTree>(treename);
239 tree->SetBranchAddress(
"signal", &m_signal);
240 tree->SetBranchAddress(
"estimated", &m_estimated);
243 const auto nEntries_MC = tree->GetEntries();
244 for (
int i = 0; i < nEntries_MC; ++i) {
246 ratios.push_back(
double(m_signal) / m_estimated);
248 gain = estimateGainFromVec(ratios);
253 B2WARNING(
"Retrieved negative median/MPV for sensor=" << sensorID <<
" uBin=" << uBin <<
" vBin=" << vBin <<
254 ". Set gain to default value (=1.0) as well.");
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);
264 return gain * gainFromDB;
272 auto dbtree = getObjectPtr<TTree>(
"dbtree");
273 dbtree->SetBranchAddress(
"run", &m_run);
274 dbtree->SetBranchAddress(
"exp", &m_exp);
275 dbtree->SetBranchAddress(
"gainMap", &gainMapPtr);
282 const auto nEntries = dbtree->GetEntries();
283 for (
int i = 0; i < nEntries; ++i) {
289 gainMapPtr =
nullptr;
291 return sum / counter;
298 if (m_boundaries.empty()) {
299 B2INFO(
"This is the first run encountered, let's say it is a boundary.");