Belle II Software development
XTCalibrationAlgorithm.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 <cdc/calibration/XTCalibrationAlgorithm.h>
9#include <cdc/calibration/XTFunction.h>
10#include <cdc/geometry/CDCGeometryPar.h>
11#include <cdc/dbobjects/CDCXtRelations.h>
12
13#include <TError.h>
14#include <TStopwatch.h>
15#include <TROOT.h>
16#include <TProfile.h>
17#include <TF1.h>
18#include <TFile.h>
19#include <TTree.h>
20
21#include <memory>
22
23#include <framework/database/DBObjPtr.h>
24
25using namespace std;
26using namespace Belle2;
27using namespace CDC;
28
29typedef std::array<float, 3> array3;
31{
33 " -------------------------- XT Calibration Algorithm -------------------------\n"
34 );
35}
36
38{
39
40 B2INFO("create and fill histo");
41 /*Create histogram*/
42 for (int i = 0; i < 56; ++i) {
43 for (int lr = 0; lr < 2; ++lr) {
44 for (int al = 0; al < m_nAlphaBins; ++al) {
45 for (int th = 0; th < m_nThetaBins; ++th) {
46 m_hProf[i][lr][al][th] = new TProfile(Form("m_hProf%d_%d_%d_%d", i, lr, al, th),
47 Form("(L=%d)-(lr=%d)-(#alpha=%3.0f)-(#theta=%3.0f); Drift time (ns);Drift Length (cm)",
48 i, lr, m_iAlpha[al], m_iTheta[th]), 210, -20, 600, 0, 1.2, "i");
49 m_hist2d[i][lr][al][th] = new TH2F(Form("h%d_%d_%d_%d", i, lr, al, th),
50 Form("(L=%d)-(lr=%d)-(#alpha=%3.0f)-(#theta=%3.0f); Drift time (ns);Drift Length (cm)",
51 i, lr, m_iAlpha[al], m_iTheta[th]), 210, -20, 600, 110, 0, 1.2);
52 if (lr == 1)
53 m_hist2dDraw[i][al][th] = new TH2F(Form("h_draw%d_%d_%d", i, al, th),
54 Form("(L=%d)-(#alpha=%3.0f)-(#theta=%3.0f); Drift time (ns);Drift Length (cm)",
55 i, m_iAlpha[al], m_iTheta[th]), 210, -20, 600, 220, -1.2, 1.2);
56 }
57 }
58 }
59 }
60
61 /* Read data and make histo*/
62
63 auto tree = getObjectPtr<TTree>("tree");
64
65 UChar_t lay;
66 Float_t dt;
67 Float_t dx;
68 Float_t Pval, alpha, theta;
69 Float_t ndf;
70
71 tree->SetBranchAddress("lay", &lay);
72 tree->SetBranchAddress("t", &dt);
73 tree->SetBranchAddress("x_u", &dx);
74 tree->SetBranchAddress("alpha", &alpha);
75 tree->SetBranchAddress("theta", &theta);
76 tree->SetBranchAddress("Pval", &Pval);
77 tree->SetBranchAddress("ndf", &ndf);
78
79 /* Disable unused branch */
80 std::vector<TString> list_vars = {"lay", "t", "x_u", "alpha", "theta", "Pval", "ndf"};
81 tree->SetBranchStatus("*", 0);
82
83 for (TString brname : list_vars) {
84 tree->SetBranchStatus(brname, 1);
85 }
86
87
88 int al = 0;
89 int th = 0;
90 TStopwatch time;
91 time.Start();
92 const Long64_t nEntries = tree->GetEntries();
93 B2INFO("Number of entries " << nEntries);
94 for (Long64_t i = 0; i < nEntries; ++i) {
95 tree->GetEntry(i);
96
97 if (Pval < m_minPval || ndf < m_minNdf) continue;
98
99 for (int k = 0; k < m_nAlphaBins; ++k) {
100 if (alpha < m_upperAlpha[k]) {
101 al = k;
102 break;
103 }
104 }
105 for (int j = 0; j < m_nThetaBins; ++j) {
106 if (theta < m_upperTheta[j]) {
107 th = j;
108 break;
109 }
110 }
111 int lr = dx > 0 ? c_Right : c_Left;
112 if (m_LRseparate) {
113 m_hProf[lay][lr][al][th]->Fill(dt, abs(dx));
114 m_hist2d[lay][lr][al][th]->Fill(dt, abs(dx));
115 } else {
116 m_hProf[lay][0][al][th]->Fill(dt, abs(dx));
117 m_hist2d[lay][0][al][th]->Fill(dt, abs(dx));
118 m_hProf[lay][1][al][th]->Fill(dt, abs(dx));
119 m_hist2d[lay][1][al][th]->Fill(dt, abs(dx));
120 }
121 m_hist2dDraw[lay][al][th]->Fill(dt, dx);
122 }
123 time.Stop();
124 B2INFO("Time to fill histograms: " << time.RealTime() << "s");
125 // time.Print();
126}
127
129{
130 gROOT->SetBatch(1);
131 gPrintViaErrorHandler = true; // Suppress huge log output from TMinuit
132 gErrorIgnoreLevel = 3001;
133 B2INFO("Start calibration");
134
135
136 const auto exprun = getRunList()[0];
137 B2INFO("ExpRun used for DB Geometry : " << exprun.first << " " << exprun.second);
138 updateDBObjPtrs(1, exprun.second, exprun.first);
139 B2INFO("Creating CDCGeometryPar object");
141
142 prepare();
143 createHisto();
144
145 B2INFO("Start Fitting");
146 std::unique_ptr<XTFunction> xt;
147 for (int l = 0; l < 56; ++l) {
148 for (int lr = 0; lr < 2; ++lr) {
149 for (int al = 0; al < m_nAlphaBins; ++al) {
150 for (int th = 0; th < m_nThetaBins; ++th) {
151 if (m_hist2d[l][lr][al][th]->GetEntries() < m_minEntriesRequired) {
152 m_fitStatus[l][lr][al][th] = FitStatus::c_lowStat;
153 continue;
154 }
155 double p0, p1, tmin;
156 TF1* fpol1;
157 if (m_useSliceFit) {
158 m_hist2d[l][lr][al][th]->FitSlicesY(0, 0, -1, 5);
159 m_hist2d_1[l][lr][al][th] = (TH1F*)gDirectory->Get(Form("h%d_%d_%d_%d_1", l, lr, al, th));
160 if (!m_hist2d_1[l][lr][al][th]) {
161 m_fitStatus[l][lr][al][th] = FitStatus::c_lowStat;
162 B2WARNING("Error, not found results of slices fit");
163 continue;
164 }
165 m_hist2d_1[l][lr][al][th]->Fit("pol1", "Q", "", 30, 60);
166 fpol1 = (TF1*)m_hProf[l][lr][al][th]->GetFunction("pol1");
167 } else {
168 /*Set Error for low statistic bin*/
169 for (int n = 0; n < m_hProf[l][lr][al][th]->GetNbinsX(); ++n) {
170 if (m_hProf[l][lr][al][th]->GetBinEntries(n) < 5 && m_hProf[l][lr][al][th]->GetBinEntries(n) > 1) {
171 m_hProf[l][lr][al][th]->SetBinError(n, 0.3 / m_hProf[l][lr][al][th]->GetBinEntries(n));
172 }
173 }
174 m_hProf[l][lr][al][th]->Fit("pol1", "Q", "", 30, 60);
175 fpol1 = (TF1*)m_hProf[l][lr][al][th]->GetFunction("pol1");
176 }
177
178 if (fpol1) {
179 //determine tmin in fitting
180 p0 = fpol1->GetParameter(0);
181 p1 = fpol1->GetParameter(1);
182 tmin = -1 * p0 / p1 + 15;
183 } else {
184 p0 = 0;
185 p1 = 0.005;
186 tmin = 12;
187 }
188
189 // B2INFO("layer " << l << ", lr " << lr << ", alpha " << m_iAlpha[al] << ", theta " << m_iTheta[th]);
190 if (m_useSliceFit) { // if slice fit results exist.
191 xt.reset(new XTFunction(m_hist2d_1[l][lr][al][th], m_xtMode));
192 } else { // from TProfile.
193 xt.reset(new XTFunction(m_hProf[l][lr][al][th], m_xtMode));
194 }
195
196 if (m_bField) {
197 int ial_old = 0;
198 int ith_old = 0;
199 for (int k = 0; k < m_nAlphaBins; ++k) {
200 if (m_iAlpha[al] < m_upperAlpha[k]) {
201 ial_old = k;
202 break;
203 }
204 }
205 for (int j = 0; j < m_nThetaBins; ++j) {
206 if (m_iTheta[th] < m_upperTheta[j]) {
207 ith_old = j;
208 break;
209 }
210 }
211
212 double p6 = m_xtPrior[l][lr][ial_old][ith_old][6];
213 if (p6 > 400) {
214 p6 = 400;
215 }
216
217 if (m_xtMode == m_xtModePrior) {
218 xt->setXTParams(m_xtPrior[l][lr][ial_old][ith_old]);
219 xt->setP6(p6);
220 } else {
221 xt->setXTParams(p0, p1, 0., 0., 0., 0., p6, m_xtPrior[l][lr][ial_old][ith_old][7]);
222 }
223 xt->setFitRange(tmin, p6 + 100);
224 } else {
225 xt->setXTParams(p0, p1, 0., 0., 0., 0., m_par6[l], 0.0001);
226 xt->setFitRange(tmin, m_par6[l] + 100);
227 }
228 xt->setDebug(m_debug);
229 xt->setBField(m_bField);
230 xt->fitXT();
231 if (xt->isValid() == false) {
232 B2WARNING("Empty xt");
233 m_fitStatus[l][lr][al][th] = c_fitFailure;
234 continue;
235 }
236 if (xt->getFitStatus() != 1) {
237 B2WARNING("Fit failed");
238 m_fitStatus[l][lr][al][th] = c_fitFailure;
239 continue;
240 }
241 if (xt->validate() == true) {
242 m_fitStatus[l][lr][al][th] = xt->getFitStatus();
243 m_xtFunc[l][lr][al][th] = (TF1*)xt->getXTFunction();
244
245 if (m_useSliceFit) {
246 m_hist2d_1[l][lr][al][th] = (TH1F*)xt->getFittedHisto();
247 } else {
248 m_hProf[l][lr][al][th] = (TProfile*)xt->getFittedHisto();
249 }
250 } else {
251 m_fitStatus[l][lr][al][th] = c_fitFailure;
252 }
253 }
254 }
255 }
256 }
258 write();
259 storeHisto();
260 return checkConvergence();
261}
262
264{
265 const double tMax = 500; // max drift time (nsec)
266 for (int l = 0; l < 56; ++l) {
267 for (int lr = 0; lr < 2; ++lr) {
268 for (int al = 0; al < m_nAlphaBins; ++al) {
269 for (int th = 0; th < m_nThetaBins; ++th) {
270 if (m_fitStatus[l][lr][al][th] == FitStatus::c_OK) {
271 TF1* fun = m_xtFunc[l][lr][al][th];
272 double y = fun->Eval(tMax);
273 if (y < 0) {
274 B2INFO("Strange XT function l " << l << " lr " << lr << " alpha " << al << " theta " << th
275 << ", replaced by initial one");
276 fun->SetParameters(m_xtPrior[l][lr][al][th]);
277 }
278 }
279 }
280 }
281 }
282 }
283}
285{
286
287 const int nTotal = 56 * 2 * m_nAlphaBins * m_nThetaBins;
288 int nFitCompleted = 0;
289 for (int l = 0; l < 56; ++l) {
290 for (int lr = 0; lr < 2; ++lr) {
291 for (int al = 0; al < m_nAlphaBins; ++al) {
292 for (int th = 0; th < m_nThetaBins; ++th) {
293 if (m_fitStatus[l][lr][al][th] == FitStatus::c_OK) {
294 nFitCompleted++;
295 }
296 }
297 }
298 }
299 }
300
301 if (static_cast<double>(nFitCompleted) / nTotal < m_threshold) {
302 B2WARNING("Less than " << m_threshold * 100 << " % of XTs were fitted.");
303 return c_NotEnoughData;
304 }
305 return c_OK;
306}
307
309{
310 B2INFO("Prepare calibration of XT");
311 const double rad2deg = 180 / M_PI;
312
314 m_nAlphaBins = dbXT->getNoOfAlphaBins();
315 m_nThetaBins = dbXT->getNoOfThetaBins();
316 for (unsigned short i = 0; i < m_nAlphaBins; ++i) {
317 array3 alpha = dbXT->getAlphaBin(i);
318 m_lowerAlpha[i] = alpha[0] * rad2deg;
319 m_upperAlpha[i] = alpha[1] * rad2deg;
320 m_iAlpha[i] = alpha[2] * rad2deg;
321 }
322
323 for (unsigned short i = 0; i < m_nThetaBins; ++i) {
324 array3 theta = dbXT->getThetaBin(i);
325 m_lowerTheta[i] = theta[0] * rad2deg;
326 m_upperTheta[i] = theta[1] * rad2deg;
327 m_iTheta[i] = theta[2] * rad2deg;
328 }
329
330 m_xtModePrior = dbXT->getXtParamMode();
331 if (!(m_xtModePrior == c_Chebyshev || m_xtModePrior == c_Polynomial)) {
332 B2FATAL("Function type before calibration is wrong " << m_xtModePrior);
333 }
334
335 B2INFO("Number of alpha bins " << m_nAlphaBins);
336 B2INFO("Number of theta bins " << m_nThetaBins);
337 B2INFO("Function type " << m_xtMode);
338
339 for (unsigned short iCL = 0; iCL < 56; ++iCL) {
340 for (unsigned short iLR = 0; iLR < 2; ++iLR) {
341 for (unsigned short iA = 0; iA < m_nAlphaBins; ++iA) {
342 for (unsigned short iT = 0; iT < m_nThetaBins; ++iT) {
343 const std::vector<float> params = dbXT->getXtParams(iCL, iLR, iA, iT);
344 unsigned short np = params.size();
345 for (unsigned short i = 0; i < np; ++i) {
346 m_xtPrior[iCL][iLR][iA][iT][i] = params[i];
347 }
348 }
349 }
350 }
351 }
352}
353
355{
356 B2INFO("write calibrated XT");
357 double par[8];
358
359
360 int nfitted = 0;
361 int nfailure = 0;
362
363 //
364 // Save to the localDB
365 //
366
367 CDCXtRelations* xtRel = new CDCXtRelations();
368 const float deg2rad = static_cast<float>(Unit::deg);
369
370 for (int i = 0; i < m_nAlphaBins; ++i) {
371 std::array<float, 3> alpha3 = {m_lowerAlpha[i]* deg2rad,
372 m_upperAlpha[i]* deg2rad,
373 m_iAlpha[i]* deg2rad
374 };
375 xtRel->setAlphaBin(alpha3);
376 }
377
378 for (int i = 0; i < m_nThetaBins; ++i) {
379 std::array<float, 3> theta3 = {m_lowerTheta[i]* deg2rad,
380 m_upperTheta[i]* deg2rad,
381 m_iTheta[i]* deg2rad
382 };
383 xtRel->setThetaBin(theta3);
384 }
385
386 xtRel->setXtParamMode(m_xtMode);
387
388 for (int th = 0; th < m_nThetaBins; ++th) {
389 for (int al = 0; al < m_nAlphaBins; ++al) {
390 for (int l = 0; l < 56; ++l) {
391 for (int lr = 0; lr < 2; ++lr) {
392 if (m_fitStatus[l][lr][al][th] != FitStatus::c_OK) {
393 nfailure += 1;
394 B2DEBUG(21, "fit failure status = " << m_fitStatus[l][lr][al][th]);
395 B2DEBUG(21, "layer " << l << ", r " << lr << ", alpha " << m_iAlpha[al] << ", theta " << m_iTheta[th]);
396 B2DEBUG(21, "number of event: " << m_hProf[l][lr][al][th]->GetEntries());
397 if (m_fitStatus[l][lr][al][th] != FitStatus::c_lowStat) {
398 if (m_xtFunc[l][lr][al][th]) {
399 B2DEBUG(21, "Probability of fit: " << m_xtFunc[l][lr][al][th]->GetProb());
400 }
401 }
402 // If fit is failed
403 // and mode of input xt (prior) is same as output, previous xt is used.
404 if (m_xtMode == m_xtModePrior) {
405 for (int i = 0; i < 8; ++i) {
406 par[i] = m_xtPrior[l][lr][al][th][i];
407 }
408 } else {
409 B2FATAL("XT mode before/after calibration is different!");
410 }
411
412 } else {
413 if (par[1] < 0) { // if negative c1, previous xt is kept.
414 for (int i = 0; i < 8; ++i) {
415 par[i] = m_xtPrior[l][lr][al][th][i];
416 }
417 } else {
418 m_xtFunc[l][lr][al][th]->GetParameters(par);
419 nfitted += 1;
420 }
421 }
422 std::vector<float> xtbuff;
423 for (int i = 0; i < 8; ++i) {
424 xtbuff.push_back(par[i]);
425 }
426 xtRel->setXtParams(l, lr, al, th, xtbuff);
427 }//lr
428 }//layer
429 }//alpha
430 }//theta
431
432 if (m_textOutput == true) {
434 }
435
436 saveCalibration(xtRel, "CDCXtRelations");
437
438 B2RESULT("Total number of xt fit: " << m_nAlphaBins * m_nThetaBins * 2 * 56);
439 B2RESULT("Successfully Fitted: " << nfitted);
440 B2RESULT("Failure Fit: " << nfailure);
441
442}
443
445{
446
447 auto hNDF = getObjectPtr<TH1F>("hNDF");
448 auto hPval = getObjectPtr<TH1F>("hPval");
449 auto hEvtT0 = getObjectPtr<TH1F>("hEventT0");
450 B2INFO("saving histograms");
451 TFile* fout = new TFile(m_histName.c_str(), "RECREATE");
452 TDirectory* top = gDirectory;
453 //store NDF, P-val. EventT0 histogram for monitoring during calibration
454 if (hNDF && hPval && hEvtT0) {
455 hEvtT0->Write();
456 hPval->Write();
457 hNDF->Write();
458 }
459 // for each layer
460
461 TDirectory* Direct[56];
462 int nhisto = 0;
463 for (int l = 0; l < 56; ++l) {
464 top->cd();
465 Direct[l] = gDirectory->mkdir(Form("lay_%d", l));
466 Direct[l]->cd();
467 for (int th = 0; th < m_nThetaBins; ++th) {
468 for (int al = 0; al < m_nAlphaBins; ++al) {
469 m_hist2dDraw[l][al][th]->Write();
470 for (int lr = 0; lr < 2; ++lr) {
471 m_hist2d[l][lr][al][th]->Write();
472 if (m_fitStatus[l][lr][al][th] != 1) continue;
473 if (m_useSliceFit) {
474 if (m_hist2d_1[l][lr][al][th]) {
475 m_hist2d_1[l][lr][al][th]->Write();
476 nhisto += 1;
477 }
478 } else {
479 m_hProf[l][lr][al][th]->Write();
480 nhisto += 1;
481 }
482 }
483 }
484 }
485 }
486 top->cd();
487
488 fout->Close();
489 B2RESULT(" " << nhisto << " histograms was stored.");
490}
491
Database object for xt-relations.
void setThetaBin(const array3 &theta)
Set theta-angle bin (rad)
void outputToFile(std::string fileName) const
Output the contents in test file format.
void setXtParamMode(unsigned short mode)
Set xt parameterization mode.
void setXtParams(const XtID xtID, const std::vector< float > &params)
Set xt parameters for the specified id.
void setAlphaBin(const array3 &alpha)
Set alpha-angle bin (rad)
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
void storeHisto()
Store histogram to file.
void prepare()
Prepare the calibration of XT.
void sanitaryCheck()
Check if there are any wrong xt functions.
double m_par6[56]
boundary parameter for fitting, semi-experiment number
float m_lowerTheta[7]
Lower boundaries of theta bins.
double m_xtPrior[56][2][18][7][8]
parameters of XT before calibration
TF1 * m_xtFunc[56][2][20][10]
XTFunction.
double m_threshold
minimal requirement for the fraction of fitted results
bool m_LRseparate
Separate LR in calibration or mix.
void createHisto()
Create histogram for calibration.
int m_minEntriesRequired
minimum number of hit per hitosgram.
TH2F * m_hist2d[56][2][20][10]
2D histo of xt
double m_minPval
minimum pvalue required
TProfile * m_hProf[56][2][20][10]
Profile xt histo.
float m_iAlpha[18]
Represented alpha in alpha bins.
float m_lowerAlpha[18]
Lower boundaries of alpha bins.
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
bool m_useSliceFit
Use slice fit or profile.
void write()
Store calibrated constant.
float m_upperAlpha[18]
Upper boundaries of alpha bins.
bool m_textOutput
output text file if true
EResult calibrate() override
Run algo on data.
EResult checkConvergence()
Check the convergence of XT fit.
int m_xtModePrior
Mode of xt before calibration; 0 is polynomial;1 is Chebyshev.
TH2F * m_hist2dDraw[56][20][10]
2d histo for draw
float m_upperTheta[7]
Upper boundaries of theta bins.
int m_xtMode
Mode of xt; 0 is polynomial;1 is Chebyshev.
TH1F * m_hist2d_1[56][2][20][10]
1D xt histo, results of slice fit
std::string m_outputFileName
Output xt filename.
float m_iTheta[7]
Represented theta in theta bins.
Class to perform fitting for each xt function.
Definition XTFunction.h:69
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
Class for accessing objects in the database.
Definition DBObjPtr.h:21
static const double deg
degree to radians
Definition Unit.h:109
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...
Abstract base class for different kinds of events.
STL namespace.