Belle II Software prerelease-11-00-00a
SVDClusterAbsoluteTimeShifterAlgorithm.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include <svd/calibration/SVDClusterAbsoluteTimeShifterAlgorithm.h>
9
10#include <svd/dbobjects/SVDAbsoluteClusterTimeShift.h>
11
12#include <TF1.h>
13#include <TH1F.h>
14#include <TH2F.h>
15#include <TCanvas.h>
16#include <TStyle.h>
17#include <TPaveStats.h>
18#include <framework/logging/Logger.h>
19#include <iostream>
20#include <TString.h>
21#include <TFitResult.h>
22
23using namespace Belle2;
24
25SVDClusterAbsoluteTimeShifterAlgorithm::SVDClusterAbsoluteTimeShifterAlgorithm(const std::string& str) :
26 CalibrationAlgorithm("SVDClusterAbsoluteTimeShifterCollector")
27 , m_id(str)
28{
29 setDescription("SVDClusterAbsoluteTimeShifter calibration algorithm");
30}
31
32CalibrationAlgorithm::EResult SVDClusterAbsoluteTimeShifterAlgorithm::calibrate()
33{
34
35 gROOT->SetBatch(true);
36
37 gStyle->SetFillColor(0);
38 gStyle->SetFillStyle(0);
39 gStyle->SetOptStat(0);
40 gStyle->SetOptFit(0);
41
42 FileStat_t info;
43 int cal_rev = 1;
44 while (gSystem->GetPathInfo(Form("algorithm_svdClusterAbsoluteTimeShifter_%s_output_rev_%d.root", m_timeAlgorithms[0].Data(),
45 cal_rev),
46 info) == 0)
47 cal_rev++;
48
49 auto payload = new Belle2::SVDAbsoluteClusterTimeShift(Form("SVDClusterAbsoluteTimeShifter_%s_rev_%d", m_id.data(), cal_rev));
50
51 // single gauss fit function
52 auto fn_singleGaus = std::make_unique<TF1>("fn_singleGaus", AbssingleGaus, -50., 50., 4);
53 fn_singleGaus->SetLineColor(kGreen + 2);
54 fn_singleGaus->SetParName(0, "N");
55 fn_singleGaus->SetParName(1, "#mu");
56 fn_singleGaus->SetParName(2, "#sigma");
57 fn_singleGaus->SetParName(3, "C");
58
59 // double gauss fit function
60 auto fn_doubleGaus = std::make_unique<TF1>("fn_doubleGaus", AbsdoubleGaus, -50., 50., 7);
61 fn_doubleGaus->SetLineColor(kRed + 2);
62 fn_doubleGaus->SetParName(0, "N");
63 fn_doubleGaus->SetParName(1, "f");
64 fn_doubleGaus->SetParName(2, "#mu_{1}");
65 fn_doubleGaus->SetParName(3, "#sigma_{1}");
66 fn_doubleGaus->SetParName(4, "#mu_{2}");
67 fn_doubleGaus->SetParName(5, "#sigma_{2}");
68 fn_doubleGaus->SetParName(6, "C");
69 fn_doubleGaus->SetParLimits(1, 0.01, 0.99);
70 fn_doubleGaus->SetParLimits(3, 0.5, 25.);
71 fn_doubleGaus->SetParLimits(5, 0.5, 25.);
72
73 int nBins = (m_outerLayer - m_innerLayer + 1) * 2;
74
75 for (auto alg : m_timeAlgorithms) {
76
77 B2INFO("Calculating shift for algorithm " << alg);
78
79 auto h_ClustersOnTrack = getObjectPtr<TH2F>(("hClsTimeOnTracks_" + alg).Data()); // __hClsOnTrack__
80
81 // auto __CDCEventT0__ = getObjectPtr<TH1F>("hCDCEventT0_");
82 // map : shift values
83 std::map< TString, Double_t > shiftValues;
84
85 std::unique_ptr<TFile> f(new TFile(Form("algorithm_svdClusterAbsoluteTimeShifter_%s_output_rev_%d.root", alg.Data(), cal_rev),
86 "RECREATE"));
87
88 B2INFO("ROOT file created at: " << gSystem->WorkingDirectory() << "/" << f->GetName());
89
90 TH1D* hShiftMean = new TH1D("hShiftMean",
91 Form("Fitted shift mean (%s);Bin;Mean (ns)", alg.Data()), nBins, 0.5, nBins + 0.5);
92 TH1D* hShiftSigma = new TH1D("hShiftSigma",
93 Form("Fitted shift sigma (%s);Bin;Sigma (ns)", alg.Data()), nBins, 0.5, nBins + 0.5);
94 for (int l = m_innerLayer; l <= m_outerLayer; l++)
95 for (int s = 0; s < 2; s++) {
96 int b = 2 * (l - m_innerLayer) + s + 1;
97 hShiftMean ->GetXaxis()->SetBinLabel(b, TString::Format("L%dS%c", l, s == 0 ? 'U' : 'V'));
98 hShiftSigma->GetXaxis()->SetBinLabel(b, TString::Format("L%dS%c", l, s == 0 ? 'U' : 'V'));
99 }
100
101 TString outPDF = Form("algorithm_svdClusterAbsoluteTimeShifter_%s_output_rev_%d.pdf", alg.Data(), cal_rev);
102 TCanvas c1("c1", "c1", 640, 480);
103 c1.Print(outPDF + "[");
104 TPad onePad("onePad", "onePad", 0, 0, 1, 1, kWhite);
105 onePad.SetFillColor(0);
106 onePad.SetBorderMode(0);
107 onePad.SetBorderSize(2);
108 onePad.SetRightMargin(0.1339713);
109 onePad.SetBottomMargin(0.15);
110 onePad.SetFrameBorderMode(0);
111 onePad.SetFrameBorderMode(0);
112 onePad.SetNumber(1);
113 onePad.Draw();
114
115
116 for (int layer = m_innerLayer; layer <= m_outerLayer; layer++) {
117 for (int side = 0; side < 2; side++) { // 0 for U, 1 for V
118 int LayerSensorID = 2 * layer - side; // 1-based index
119
120
121
122 TString binLabel = TString::Format("L%iS%c", layer, (side == 0 ? 'U' : 'V'));
123 TH1F* hist = (TH1F*)h_ClustersOnTrack->ProjectionX(Form("hClsTimeOnTracks_L%dS%c", layer, (side == 0 ? 'U' : 'V')), LayerSensorID,
124 LayerSensorID, "");
125 hist->SetTitle(Form("Cluster Time in L%dS%c", layer, (side == 0 ? 'U' : 'V')));
126 hist->SetDirectory(0);
127
128 B2INFO("Histogram: " << hist->GetName() <<
129 " Entries (n. clusters): " << hist->GetEntries());
130 if (hist->GetEntries() < m_minEntries) {
131 B2INFO("Histogram: " << hist->GetName() <<
132 " Entries (n. clusters): " << hist->GetEntries() <<
133 " Entries required: " << m_minEntries);
134 B2WARNING("Not enough data, adding one run to the collector");
135 delete hist;
136 c1.Print(outPDF + "]");
137 f->Close();
138 gSystem->Unlink(Form("algorithm_svdClusterAbsoluteTimeShifter_%s_output_rev_%d.root", m_timeAlgorithms[0].Data(), cal_rev));
139 return c_NotEnoughData;
140 }
141 if (hist->GetEntries() == 0) {
142 B2INFO("Histogram: " << hist->GetName() <<
143 " Entries (n. clusters): " << hist->GetEntries() <<
144 " No entries, setting shift to 0.");
145 shiftValues[binLabel] = 0.;
146 delete hist;
147 continue;
148 }
149 onePad.Clear();
150 onePad.cd();
151 hist->Draw();
152 if (hist->GetMaximum() < 200.) {
153 int rebinValue = 200. / hist->GetMaximum() + 1.;
154 while ((hist->GetNbinsX()) % rebinValue != 0) rebinValue++;
155 hist->Rebin(rebinValue);
156 }
157
158 int fitCount = 0;
159 bool isSingleGausFitValid, isDoubleGausFitValid;
160 while (fitCount++ < 5) {
161 double histMean = hist->GetMean();
162 double histStd = hist->GetStdDev();
163 fn_singleGaus->SetParameter(0, hist->GetSumOfWeights() * hist->GetBinWidth(1));
164 fn_singleGaus->SetParameter(1, histMean);
165 fn_singleGaus->SetParameter(2, histStd * 0.75);
166 fn_singleGaus->SetParameter(3, 1.);
167 fn_singleGaus->SetParLimits(1, histMean - histStd, histMean + histStd);
168 auto singleGausFitStatus = hist->Fit("fn_singleGaus", "SQ", "",
169 histMean - 2. * histStd, histMean + 2. * histStd);
170 isSingleGausFitValid = singleGausFitStatus->IsValid();
171
172 fn_doubleGaus->SetParameter(0, hist->GetSumOfWeights() * hist->GetBinWidth(1));
173 fn_doubleGaus->SetParameter(1, 0.95);
174 fn_doubleGaus->SetParameter(2, fn_singleGaus->GetParameter(1));
175 fn_doubleGaus->SetParameter(3, std::fabs(fn_singleGaus->GetParameter(2)));
176 fn_doubleGaus->SetParameter(4, fn_singleGaus->GetParameter(1) - 3.);
177 fn_doubleGaus->SetParameter(5, std::fabs(fn_singleGaus->GetParameter(2)) + 5.);
178 fn_doubleGaus->SetParameter(6, 10.);
179 fn_doubleGaus->SetParLimits(2,
180 fn_doubleGaus->GetParameter(2) - m_maximumAllowedShift,
181 fn_doubleGaus->GetParameter(2) + m_maximumAllowedShift);
182 fn_doubleGaus->SetParLimits(4,
183 fn_doubleGaus->GetParameter(4) - m_maximumAllowedShift,
184 fn_doubleGaus->GetParameter(4) + m_maximumAllowedShift);
185
186 auto doubleGausFitStatus = hist->Fit("fn_doubleGaus", "SQ+");
187 isDoubleGausFitValid = doubleGausFitStatus->IsValid();
188 if (isDoubleGausFitValid) break;
189 int rebinValue = 2;
190 while ((hist->GetNbinsX()) % rebinValue != 0) rebinValue++;
191 hist->Rebin(rebinValue);
192 }
193
194 TPaveStats* ptstats = new TPaveStats(0.55, 0.73, 0.85, 0.88, "brNDC");
195 ptstats->SetName("stats1");
196 ptstats->SetBorderSize(1);
197 ptstats->SetFillColor(0);
198 ptstats->SetTextAlign(12);
199 ptstats->SetTextFont(42);
200 ptstats->SetTextColor(kGreen + 2);
201 ptstats->SetOptStat(11);
202 ptstats->SetParent(hist);
203 ptstats->Draw();
204 ptstats->AddText("Single Gauss");
205 for (int npar = 0; npar < (fn_singleGaus->GetNpar()); npar++)
206 ptstats->AddText(TString::Format("%s = %.3f #pm %.4f",
207 fn_singleGaus->GetParName(npar),
208 fn_singleGaus->GetParameter(npar),
209 fn_singleGaus->GetParError(npar)));
210 ptstats = new TPaveStats(0.55, 0.49, 0.85, 0.73, "brNDC");
211 ptstats->SetName("stats2");
212 ptstats->SetBorderSize(1);
213 ptstats->SetFillColor(0);
214 ptstats->SetTextAlign(12);
215 ptstats->SetTextFont(42);
216 ptstats->SetTextColor(kRed + 2);
217 ptstats->SetOptStat(11);
218 ptstats->SetParent(hist);
219 ptstats->Draw();
220 ptstats->AddText("Double Gauss");
221 for (int npar = 0; npar < (fn_doubleGaus->GetNpar()); npar++)
222 ptstats->AddText(TString::Format("%s = %.3f #pm %.4f",
223 fn_doubleGaus->GetParName(npar),
224 fn_doubleGaus->GetParameter(npar),
225 fn_doubleGaus->GetParError(npar)));
226 ptstats = new TPaveStats(0.55, 0.43, 0.85, 0.49, "brNDC");
227 ptstats->SetName("stats3");
228 ptstats->SetBorderSize(1);
229 ptstats->SetFillColor(0);
230 ptstats->SetTextAlign(12);
231 ptstats->SetTextFont(42);
232 ptstats->SetTextColor(kBlue + 2);
233 ptstats->SetOptStat(11);
234 ptstats->SetParent(hist);
235 ptstats->Draw();
236
237 // setting `fillShiftVal` a higher than `m_maximumAllowedShift`. If both fit statuses are invalid,
238 // then this condition will set shift value to the histogram mean.
239 double fillShiftVal = m_maximumAllowedShift + 1.;
240
241 if (isDoubleGausFitValid) {
242 fillShiftVal = (fn_doubleGaus->GetParameter(1) > 0.5 ?
243 fn_doubleGaus->GetParameter(2) : fn_doubleGaus->GetParameter(4));
244 ptstats->AddText(TString::Format("#splitline{Shift Value from Double Gauss}{%.3f}", fillShiftVal));
245 } else if (isSingleGausFitValid) {
246 fillShiftVal = fn_singleGaus->GetParameter(1);
247 B2WARNING("Fit failed for " << hist->GetName() <<
248 "; using mean from single gauss fit. ");
249 ptstats->AddText(TString::Format("#splitline{Shift Value from Single Gauss}{%.3f}", fillShiftVal));
250 }
251
252 if (std::fabs(fillShiftVal) > m_maximumAllowedShift) {
253 B2WARNING("Shift value is more than allowed or fit failed in " <<
254 hist->GetName() << " : " <<
255 shiftValues[binLabel] <<
256 "; using mean of the histogram.");
257 fillShiftVal = hist->GetMean();
258 ptstats->AddText(TString::Format("#splitline{Shift Value from Histogram Mean}{%.3f}", fillShiftVal));
259 }
260
261 shiftValues[binLabel] = 1. * int(1000. * fillShiftVal) / 1000.;
262
263 double fillSigmaVal;
264 if (isDoubleGausFitValid)
265 fillSigmaVal = (fn_doubleGaus->GetParameter(1) > 0.5 ?
266 fn_doubleGaus->GetParameter(3) : fn_doubleGaus->GetParameter(5));
267 else if (isSingleGausFitValid)
268 fillSigmaVal = fn_singleGaus->GetParameter(2);
269 else
270 fillSigmaVal = hist->GetStdDev();
271
272 int binIdx = 2 * (layer - 3) + side + 1;
273 hShiftMean ->SetBinContent(binIdx, fillShiftVal);
274 hShiftSigma->SetBinContent(binIdx, fillSigmaVal);
275
276 // Ensure the plot is visually restricted to [-50, 50] after all fitting and drawing
277 hist->GetXaxis()->SetRangeUser(-50, 50);
278 gPad->Range(-50, 0, 50, hist->GetMaximum() * 1.2);
279 gPad->Update();
280
281 c1.Print(outPDF, TString("Title:") + hist->GetName());
282 f->cd();
283 hist->Write();
284
285 delete hist;
286 } // loop over cluster size
287
288 f->cd();
289
290 } // loop over sensors
291
292 for (auto item : shiftValues)
293 payload->setClusterTimeShift(alg, item.first, item.second);
294 B2INFO("Shift values for algorithm " + alg + " : ");
295 for (const auto& item : shiftValues) {
296 B2INFO(" " << item.first << " : " << item.second);
297
298 }
299
300 onePad.cd();
301
302 f->cd();
303
304 c1.Print(outPDF + "]");
305 f->cd();
306 hShiftMean->Write();
307 hShiftSigma->Write();
308 delete hShiftMean;
309 delete hShiftSigma;
310 f->Close();
311 } // loop over algorithms
312
313 saveCalibration(payload, "SVDAbsoluteClusterTimeShift");
314
315 return c_OK;
316}
317
318bool SVDClusterAbsoluteTimeShifterAlgorithm::isBoundaryRequired(const Calibration::ExpRun& currentRun)
319{
320
321 TString alg = "CoG3";
322 if (std::find(m_timeAlgorithms.begin(), m_timeAlgorithms.end(), "CoG3") == m_timeAlgorithms.end()) {
323 alg = m_timeAlgorithms[0];
324 }
325 float meanTimeL3V = 0;
326 auto h_ClustersOnTrack = getObjectPtr<TH2F>(("hClsTimeOnTracks_" + alg).Data());
327 int layer = 3;
328 int side = 0;
329 int LayerSensorID = 2 * layer - side;
330 auto timeL3V = (TH1F*)h_ClustersOnTrack->ProjectionX(Form("hClsTimeOnTracks_L%dS%c", layer, (side == 0 ? 'U' : 'V')), LayerSensorID,
331 LayerSensorID, "");
332 if (!timeL3V) {
333 if (m_previousTimeMeanL3V)
334 meanTimeL3V = m_previousTimeMeanL3V.value();
335 } else {
336 if (timeL3V->GetEntries() > m_minEntries)
337 meanTimeL3V = timeL3V->GetMean();
338 else {
339 if (m_previousTimeMeanL3V)
340 meanTimeL3V = m_previousTimeMeanL3V.value();
341 }
342 }
343 if (!m_previousTimeMeanL3V) {
344 B2INFO("Setting start payload boundary to be the first run ("
345 << currentRun.first << "," << currentRun.second << ")");
346 m_previousTimeMeanL3V.emplace(meanTimeL3V);
347
348 return true;
349 } else if (abs(meanTimeL3V - m_previousTimeMeanL3V.value()) > m_allowedTimeShift) {
350 B2INFO("Histogram mean has shifted from " << m_previousTimeMeanL3V.value()
351 << " to " << meanTimeL3V << ". We are requesting a new payload boundary for ("
352 << currentRun.first << "," << currentRun.second << ")");
353 m_previousTimeMeanL3V.emplace(meanTimeL3V);
354 return true;
355 } else {
356 return false;
357 }
358}
359
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition CalibrationAlgorithm.cc:297
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition CalibrationAlgorithm.h:321
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition CalibrationAlgorithm.h:40
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfully =0 in Python.
Definition CalibrationAlgorithm.h:41
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition CalibrationAlgorithm.h:43
Belle2::CalibrationAlgorithm::CalibrationAlgorithm
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
Definition CalibrationAlgorithm.h:140
Belle2::SVDAbsoluteClusterTimeShift
This class store the shift of the cluster time w.r.t a full calibration.
Definition SVDAbsoluteClusterTimeShift.h:29
Belle2::SVDClusterAbsoluteTimeShifterAlgorithm::m_minEntries
int m_minEntries
Set the minimum number of entries required in the histograms.
Definition SVDClusterAbsoluteTimeShifterAlgorithm.h:80
Belle2::SVDClusterAbsoluteTimeShifterAlgorithm::SVDClusterAbsoluteTimeShifterAlgorithm
SVDClusterAbsoluteTimeShifterAlgorithm(const std::string &str)
Constructor set the prefix to SVDTimeCalibrationCollector.
Definition SVDClusterAbsoluteTimeShifterAlgorithm.cc:25
Belle2::SVDClusterAbsoluteTimeShifterAlgorithm::m_timeAlgorithms
std::vector< TString > m_timeAlgorithms
Hardcoded list of time algorithms to calibrate.
Definition SVDClusterAbsoluteTimeShifterAlgorithm.h:81
Belle2::SVDClusterAbsoluteTimeShifterAlgorithm::m_id
std::string m_id
Parameter given to set the UniqueID of the payload.
Definition SVDClusterAbsoluteTimeShifterAlgorithm.h:76
Belle2::SVDClusterAbsoluteTimeShifterAlgorithm::m_maximumAllowedShift
float m_maximumAllowedShift
Allowed deviation of clsOnTracks histo wrt EventT0 histo.
Definition SVDClusterAbsoluteTimeShifterAlgorithm.h:79
Belle2::SVDClusterAbsoluteTimeShifterAlgorithm::m_previousTimeMeanL3V
std::optional< float > m_previousTimeMeanL3V
Cluster time of the previous run.
Definition SVDClusterAbsoluteTimeShifterAlgorithm.h:77
Belle2::SVDClusterAbsoluteTimeShifterAlgorithm::isBoundaryRequired
virtual bool isBoundaryRequired(const Calibration::ExpRun &currentRun) override
If the event T0 changes significantly return true.
Definition SVDClusterAbsoluteTimeShifterAlgorithm.cc:318
Belle2::CalibrationAlgorithm::getObjectPtr
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
Definition CalibrationAlgorithm.h:414
Belle2::AbsdoubleGaus
double AbsdoubleGaus(const double *x, const double *par)
Double gauss function, where N is a normalization constant, f is the fractional contribution of the f...
Definition SVDClusterAbsoluteTimeShifterAlgorithm.h:118
Belle2::AbssingleGaus
double AbssingleGaus(const double *x, const double *par)
Single gauss function, where N is a normalization constant, a is the mean of the Gauss distribution,...
Definition SVDClusterAbsoluteTimeShifterAlgorithm.h:97
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17