11 #include <pxd/calibration/PXDChargeCalibrationAlgorithm.h>
12 #include <pxd/dbobjects/PXDClusterChargeMapPar.h>
21 #include <boost/format.hpp>
41 void getNumberOfBins(
const std::shared_ptr<TH1I>& histo_ptr,
unsigned short& nBinsU,
unsigned short& nBinsV)
43 set<unsigned short> uBinSet;
44 set<unsigned short> vBinSet;
47 for (
auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
49 string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
52 istringstream stream(label);
54 getline(stream, token,
'_');
55 getline(stream, token,
'_');
56 unsigned short uBin = std::stoi(token);
58 getline(stream, token,
'_');
59 unsigned short vBin = std::stoi(token);
65 if (uBinSet.empty() || vBinSet.empty()) {
66 B2FATAL(
"Not able to determine the grid size. Something is wrong with collected data.");
68 nBinsU = *uBinSet.rbegin() + 1;
69 nBinsV = *vBinSet.rbegin() + 1;
76 set<unsigned short> sensorSet;
79 for (
auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
81 string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
84 istringstream stream(label);
86 getline(stream, token,
'_');
87 VxdID sensorID(token);
88 sensorSet.insert(sensorID.getID());
91 return sensorSet.size();
97 PXDChargeCalibrationAlgorithm::PXDChargeCalibrationAlgorithm():
99 minClusters(5000), noiseSigma(0.6), safetyFactor(2.0), forceContinue(false), strategy(0)
102 " ------------------------------ PXDChargeCalibrationAlgorithm -----------------------------------\n"
104 " Algorithm for estimating pxd median/MPV cluster charges for different position on sensor in ADU \n"
105 " ------------------------------------------------------------------------------------------------\n"
116 auto cluster_counter = getObjectPtr<TH1I>(
"PXDClusterCounter");
122 unsigned short nBinsU = 0;
123 unsigned short nBinsV = 0;
129 B2INFO(
"Not enough Data: Only " << cluster_counter->GetEntries() <<
" hits were collected but " <<
int(
safetyFactor *
minClusters *
130 nSensors * nBinsU * nBinsV) <<
" needed!");
133 B2INFO(
"Continue despite low statistics: Only " << cluster_counter->GetEntries() <<
" hits were collected but" <<
int(
135 nSensors * nBinsU * nBinsV) <<
" would be desirable!");
139 B2INFO(
"Start calibration using a " << nBinsU <<
"x" << nBinsV <<
" grid per sensor.");
145 for (
auto histoBin = 1; histoBin <= cluster_counter->GetXaxis()->GetNbins(); histoBin++) {
147 string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
150 istringstream stream(label);
152 getline(stream, token,
'_');
153 VxdID sensorID(token);
155 getline(stream, token,
'_');
156 unsigned short uBin = std::stoi(token);
158 getline(stream, token,
'_');
159 unsigned short vBin = std::stoi(token);
162 int numberOfDataHits = cluster_counter->GetBinContent(histoBin);
173 B2WARNING(label <<
": Number of data hits too small for fitting (" << numberOfDataHits <<
" < " <<
minClusters <<
174 "). Use default value of 0.");
183 B2INFO(
"PXD Cluster Charge Calibration Successful");
195 const string treename = str(format(
"tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
197 vector<double> signals;
200 auto tree = getObjectPtr<TTree>(treename);
201 tree->SetBranchAddress(
"signal", &m_signal);
204 const auto nEntries = tree->GetEntries();
205 for (
int i = 0; i < nEntries; ++i) {
208 double noise = gRandom->Gaus(0.0,
noiseSigma);
209 signals.push_back(m_signal + noise);
214 B2FATAL(
"strategy unavailable, use 0 for medians or 1 for landau fit!");
221 auto size = signals.size();
226 sort(signals.begin(), signals.end());
228 return (signals[size / 2 - 1] + signals[size / 2]) / 2;
230 return signals[size / 2];
237 auto size = signals.size();
238 if (size == 0)
return 0.0;
241 int max = *max_element(signals.begin(), signals.end());
242 int min = *min_element(signals.begin(), signals.end());
246 TH1D* hist_signals =
new TH1D(
"",
"", max - min, min, max);
247 for (
auto it = signals.begin(); it != signals.end(); ++it) {
248 hist_signals->Fill(*it);
252 TF1* landau =
new TF1(
"landau",
"TMath::Landau(x,[0],[1])*[2]", min, max);
253 landau->SetParNames(
"MPV",
"sigma",
"scale");
254 landau->SetParameters(100, 1, 1000);
257 Int_t status = hist_signals->Fit(
"landau",
"Lq",
"", 30, 350);
258 double MPV = landau->GetParameter(
"MPV");
265 if (status == 0)
return MPV;
267 B2WARNING(
"Fit failed! using default value 0.0!");