Belle II Software development
CDCDedxInjectTimeAlgorithm.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
9#include <cdc/calibration/CDCdEdx/CDCDedxInjectTimeAlgorithm.h>
10
11#include <cdc/utilities/CDCDedxMeanPred.h>
12#include <cdc/utilities/CDCDedxSigmaPred.h>
13
14#include <framework/gearbox/Const.h>
15
16#include <TCanvas.h>
17#include <TF1.h>
18#include <TLegend.h>
19
20using namespace Belle2;
21
22//-----------------------------------------------------------------
23// Implementation
24//-----------------------------------------------------------------
26 CalibrationAlgorithm("CDCDedxElectronCollector"),
27 m_sigmaR(2.0),
28 m_dedxBins(250),
29 m_dedxMin(0.0),
30 m_dedxMax(2.5),
31 m_chiBins(250),
32 m_chiMin(-10.0),
33 m_chiMax(10.0),
34 m_countR(0),
35 m_thersE(1000),
36 m_isminStat(false),
37 m_ismakePlots(true),
38 m_isMerge(true),
39 m_prefix("cdcdedx_injcal"),
40 m_suffix("")
41{
42 // Set module properties
43 setDescription("A calibration algorithm for CDC dE/dx injection time gain and reso");
44}
45
46//-----------------------------------------------------------------
47// Run the calibration
48//-----------------------------------------------------------------
50{
51
53
54 //existing inject time payload for merging
55 if (!m_DBInjectTime.isValid())
56 B2FATAL("There is no valid payload for Injection time");
57
58 // Get data objects
59 auto ttree = getObjectPtr<TTree>("tree");
60 if (!ttree) return c_NotEnoughData;
61
62 double dedx = 0.0, injtime = 0.0, injring = 1.0, costh, mom;
63 int nhits;
64 ttree->SetBranchAddress("dedx", &dedx);
65 ttree->SetBranchAddress("injtime", &injtime);
66 ttree->SetBranchAddress("injring", &injring);
67 ttree->SetBranchAddress("costh", &costh);
68 ttree->SetBranchAddress("p", &mom);
69 ttree->SetBranchAddress("nhits", &nhits);
70
71 //way to define time bins/edges only once
72 if (m_countR == 0) {
73 defineTimeBins(); //returns m_vtlocaledges
74 m_tbins = m_vtlocaledges.size() - 1;
76 m_countR++;
77 }
78
79 TH1D* htimes = new TH1D(Form("htimes_%s", m_suffix.data()), "", m_tbins, m_tedges);
80
81 //time bins are changeable from out so vector is used
82 std::array<std::vector<TH1D*>, numdedx::nrings> hdedx, htime, hchi, hdedx_corr;
83 defineHisto(hdedx, "dedx");
84 defineHisto(hdedx_corr, "dedx_corr");
85 defineHisto(htime, "timeinj");
86 defineHisto(hchi, "chi");
87
88 const double tzedges[4] = {0, 2.0e3, 0.1e6, 20e6};
89 std::array<std::array<TH1D*, 3>, numdedx::nrings> hztime;
90 defineTimeHisto(hztime);
91
92 //fill histograms
93 for (int i = 0; i < ttree->GetEntries(); ++i) {
94
95 ttree->GetEvent(i);
96 if (dedx <= 0 || injtime < 0 || injring < 0) continue;
97
98 //add larger times to the last bin
99 if (injtime > m_tedges[m_tbins]) injtime = m_tedges[m_tbins] - 10.0;
100
101 //injection ring
102 int wr = 0;
103 if (injring > 0.5) wr = 1;
104
105 //injection time bin
106 unsigned int tb = htimes->GetXaxis()->FindBin(injtime);
107 if (tb > m_tbins)tb = m_tbins; //overflow
108 tb = tb - 1;
109
110 htimes->Fill(injtime);
111 if (injtime < tzedges[1]) hztime[wr][0]->Fill(injtime);
112 else if (injtime < tzedges[2]) hztime[wr][1]->Fill(injtime);
113 else hztime[wr][2]->Fill(injtime);
114
115 hdedx[wr][tb]->Fill(dedx);
116 htime[wr][tb]->Fill(injtime);
117 }
118
119 //keep merging runs to achieve enough stats
120 m_isminStat = false;
121 checkStatistics(hdedx);
122 if (m_isminStat) {
123 deleteHisto(htime);
124 deleteHisto(hdedx);
125 deleteHisto(hdedx_corr);
126 deleteHisto(hchi);
127 deleteTimeHisto(hztime);
128 delete htimes;
129 return c_NotEnoughData;
130 }
131
132 //clear vector of existing constants
133 std::map<int, std::vector<double>> vmeans, vresos, vtimes;
134
135 // get time vector
136 for (unsigned int ir = 0; ir < c_rings; ir++) {
137 for (unsigned int it = 0; it < m_tbins; it++) {
138 double avgtime = htime[ir][it]->GetMean();
139 double avgtimeerr = htime[ir][it]->GetMeanError();
140 vtimes[ir * 2].push_back(avgtime);
141 vtimes[ir * 2 + 1].push_back(avgtimeerr);
142 }
143 }
144
145 //Fit dedx to get mean and resolution
146 getMeanReso(hdedx, vmeans, vresos);
147
148 //Bin-bias correction to mean
149 std::map<int, std::vector<double>> vmeanscorr;
150 correctBinBias(vmeanscorr, vmeans, vtimes, htimes);
151
152 //scale the mean and merge with old constants
153 std::array<double, numdedx::nrings> scale;
154 std::map<int, std::vector<double>> vmeanscal;
155 createPayload(scale, vmeanscorr, vmeanscal, "mean");
156
157 //................................................
158 // Do the calibration for resolution
159 //................................................
160 CDCDedxMeanPred mbg;
162 mbg.setParameters();
163 sbg.setParameters();
164
165 //fill histograms for the resolution
166 for (int i = 0; i < ttree->GetEntries(); ++i) {
167
168 ttree->GetEvent(i);
169 if (dedx <= 0 || injtime < 0 || injring < 0) continue;
170
171 double corrcetion = getCorrection(injring, injtime, vmeanscal);
172 double old_cor = m_DBInjectTime->getCorrection("mean", injring, injtime);
173
174 dedx = (dedx * old_cor) / corrcetion;
175 //add larger times to the last bin
176 if (injtime > m_tedges[m_tbins]) injtime = m_tedges[m_tbins] - 10.0;
177
178 //injection ring
179 int wr = 0;
180 if (injring > 0.5) wr = 1;
181
182 //injection time bin
183 unsigned int tb = htimes->GetXaxis()->FindBin(injtime);
184 if (tb > m_tbins)tb = m_tbins; //overflow
185 tb = tb - 1;
186
187 double predmean = mbg.getMean(mom / Const::electron.getMass());
188 double predsigma = sbg.nhitPrediction(nhits) * sbg.ionzPrediction(dedx) * sbg.cosPrediction(costh);
189
190 double chi = (dedx - predmean) / predsigma;
191 hdedx_corr[wr][tb]->Fill(dedx);
192 hchi[wr][tb]->Fill(chi);
193 }
194
195 // fit chi to get mean and resolution
196 std::map<int, std::vector<double>> vmeans_chi, vresos_chi;
197 getMeanReso(hchi, vmeans_chi, vresos_chi);
198
199 //bin-bias correction to the resolution
200 std::map<int, std::vector<double>> vresoscorr;
201 correctBinBias(vresoscorr, vresos_chi, vtimes, htimes);
202
203 // scale the resolution
204 std::map<int, std::vector<double>> vresoscal;
205 std::array<double, numdedx::nrings> scale_reso;
206 createPayload(scale_reso, vresoscorr, vresoscal, "reso");
207
208 //Fit the corrected mean to check for consistency
209 std::map<int, std::vector<double>> vmeans_corr, vresos_corr;
210 getMeanReso(hdedx_corr, vmeans_corr, vresos_corr);
211
212 //................................................
213 //preparing final payload
214 //................................................
215 m_vinjPayload.clear();
216 m_vinjPayload.reserve(6);
217 for (int ir = 0; ir < 2; ir++) {
218 m_vinjPayload.push_back(m_vtlocaledges);
219 m_vinjPayload.push_back(vmeanscal[ir * 2]);
220 m_vinjPayload.push_back(vresoscal[ir * 2]);
221 }
222
223 if (m_ismakePlots) {
224
225 //0 plot event track statistics
227
228 //1 plot injection time
229 plotInjectionTime(hztime);
230
231 //2. Draw dedxfits
232 plotBinLevelDist(hdedx, "dedxfits");
233
234 //3. Draw chifits
235 plotBinLevelDist(hchi, "chifits");
236
237 //4. Draw timedist
238 plotBinLevelDist(htime, "timedist");
239
240 //5. plot relative const., bias-bias corrected for dedx
241 plotRelConstants(vmeans, vresos, vmeanscorr, "dedx");
242
243 //6. plot relative const., bias-bias corrected for chi
244 plotRelConstants(vmeans_chi, vresos_chi, vresoscorr, "chi");
245
246 //7. plot mean and resolution of corrected dedx to check for consistency
247 plotRelConstants(vmeans_corr, vresos_corr, vresoscorr, "dedx_corr");
248
249 //8. plot time statistics dist
250 plotTimeStat(htime);
251
252 //9. plot final merged const. and comparison to old
253 plotFinalConstants(vmeanscal, vresoscal, scale, scale_reso);
254 }
255
256 //saving payloads;
258 saveCalibration(gains, "CDCDedxInjectionTime");
259 B2INFO("dE/dx Injection time calibration done");
260
261 //delete all histograms
262 deleteHisto(htime);
263 deleteHisto(hdedx);
264 deleteHisto(hdedx_corr);
265 deleteHisto(hchi);
266 deleteTimeHisto(hztime);
267 delete htimes;
268
269 B2INFO("Saving calibration for: " << m_suffix << "");
270 return c_OK;
271}
272
273//------------------------------------
274void CDCDedxInjectTimeAlgorithm::getMeanReso(std::array<std::vector<TH1D*>, numdedx::nrings>& hvar,
275 std::map<int, std::vector<double>>& vmeans, std::map<int, std::vector<double>>& vresos)
276{
277 for (unsigned int ir = 0; ir < c_rings; ir++) {
278
279 for (unsigned int it = 0; it < m_tbins; it++) {
280 double mean = 1.00, meanerr = 0.0;
281 double reso = 1.00, resoerr = 0.0;
282 if (hvar[ir][it]->Integral() > 250) {
283 fstatus status;
284 fitGaussianWRange(hvar[ir][it], status);
285 if (status != fitOK) {
286 mean = hvar[ir][it]->GetMean();
287 hvar[ir][it]->SetTitle(Form("%s, (%d)", hvar[ir][it]->GetTitle(), status));
288 } else {
289 mean = hvar[ir][it]->GetFunction("gaus")->GetParameter(1);
290 meanerr = hvar[ir][it]->GetFunction("gaus")->GetParError(1);
291 reso = hvar[ir][it]->GetFunction("gaus")->GetParameter(2);
292 resoerr = hvar[ir][it]->GetFunction("gaus")->GetParError(2);
293 std::string title = Form("#mu_{fit}: %0.03f, #sigma_{fit}: %0.03f", mean, reso);
294 hvar[ir][it]->SetTitle(Form("%s, %s", hvar[ir][it]->GetTitle(), title.data()));
295 }
296 }
297
298 vmeans[ir * 2].push_back(mean);
299 vresos[ir * 2].push_back(reso);
300 vmeans[ir * 2 + 1].push_back(meanerr);
301 vresos[ir * 2 + 1].push_back(resoerr);
302 }
303 }
304}
305
306//----------------------------------------
307void CDCDedxInjectTimeAlgorithm::fitGaussianWRange(TH1D*& temphist, fstatus& status)
308{
309 double histmean = temphist->GetMean();
310 double histrms = temphist->GetRMS();
311 temphist->GetXaxis()->SetRangeUser(histmean - 5.0 * histrms, histmean + 5.0 * histrms);
312
313 int fs = temphist->Fit("gaus", "Q0");
314 if (fs != 0) {
315 B2INFO(Form("\tFit (round 1) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
316 status = fitFailed;
317 return;
318 } else {
319 double mean = temphist->GetFunction("gaus")->GetParameter(1);
320 double width = temphist->GetFunction("gaus")->GetParameter(2);
321 temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
322 fs = temphist->Fit("gaus", "QR", "", mean - m_sigmaR * width, mean + m_sigmaR * width);
323 if (fs != 0) {
324 B2INFO(Form("\tFit (round 2) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
325 status = fitFailed;
326 return;
327 } else {
328 temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
329 B2INFO(Form("\tFit for hist (%s) successful (status = %d)", temphist->GetName(), fs));
330 status = fitOK;
331 }
332 }
333}
334
335//------------------------------------
337{
338 int cruns = 0;
339 for (auto expRun : getRunList()) {
340 if (cruns == 0)B2INFO("CDCDedxInjectTimeAlgorithm: start exp " << expRun.first << " and run " << expRun.second << "");
341 cruns++;
342 }
343
344 const auto erStart = getRunList()[0];
345 int estart = erStart.first;
346 int rstart = erStart.second;
347
348 const auto erEnd = getRunList()[cruns - 1];
349 int rend = erEnd.second;
350
351 updateDBObjPtrs(1, erStart.second, erStart.first);
352
353 if (m_isminStat) {
354 m_suffix = Form("e%dr%dto%d_nruns%d", estart, rstart, rend, cruns);
355 B2INFO("\t+ run = " << rend << ", m_suffix = " << m_suffix << "");
356 } else {
357 m_suffix = Form("e%dr%d", estart, rstart);
358 B2INFO("tool run = " << estart << ", exp = " << estart << ", m_suffix = " << m_suffix << "");
359 }
360}
361
362//------------------------------------
364{
365 //empty local vector or find a way to execulate this function
366 //only once
367 if (!m_vtlocaledges.empty()) m_vtlocaledges.clear();
368 if (m_vtedges.empty()) {
369 double fixedges[69];
370 for (int ib = 0; ib < 69; ib++) {
371 fixedges[ib] = ib * 0.5 * 1e3;
372 if (ib > 40 && ib <= 60) fixedges[ib] = fixedges[ib - 1] + 1.0 * 1e3;
373 else if (ib > 60 && ib <= 64) fixedges[ib] = fixedges[ib - 1] + 10.0 * 1e3;
374 else if (ib > 64 && ib <= 65) fixedges[ib] = fixedges[ib - 1] + 420.0 * 1e3;
375 else if (ib > 65 && ib <= 66) fixedges[ib] = fixedges[ib - 1] + 500.0 * 1e3;
376 else if (ib > 66) fixedges[ib] = fixedges[ib - 1] + 2e6;
377 m_vtlocaledges.push_back(fixedges[ib]);
378 }
379 } else {
380 for (unsigned int ib = 0; ib < m_vtedges.size(); ib++)
381 m_vtlocaledges.push_back(m_vtedges.at(ib));
382 }
383}
384
385//------------------------------------
386void CDCDedxInjectTimeAlgorithm::defineHisto(std::array<std::vector<TH1D*>, numdedx::nrings>& htemp, std::string var)
387{
388 for (unsigned int ir = 0; ir < c_rings; ir++) {
389 htemp[ir].resize(m_tbins);
390 for (unsigned int it = 0; it < m_tbins; it++) {
391 std::string label = getTimeBinLabel(m_tedges[it], it);
392 std::string title = Form("%s(%s), time(%s)", m_suffix.data(), m_sring[ir].data(), label.data());
393 std::string hname = Form("h%s_%s_%s_t%d", var.data(), m_sring[ir].data(), m_suffix.data(), it);
394 if (var == "dedx" or var == "dedx_corr") htemp[ir][it] = new TH1D(hname.data(), "", m_dedxBins, m_dedxMin, m_dedxMax);
395 else if (var == "chi") htemp[ir][it] = new TH1D(hname.data(), "", m_chiBins, m_chiMin, m_chiMax);
396 else htemp[ir][it] = new TH1D(hname.data(), "", 50, m_tedges[it], m_tedges[it + 1]);
397 htemp[ir][it]->SetTitle(Form("%s;%s;entries", title.data(), var.data()));
398 }
399 }
400}
401
402//------------------------------------
403void CDCDedxInjectTimeAlgorithm::defineTimeHisto(std::array<std::array<TH1D*, 3>, numdedx::nrings>& htemp)
404{
405 const int tzoom = 3;
406 const int nt[tzoom] = {50, 500, 1000};
407 double tzedges[tzoom + 1] = {0, 2.0e3, 0.1e6, 20e6};
408 std::string stname[tzoom] = {"early", "mid", "later"};
409 std::string stlabel[tzoom] = {"zoom <2ms", "early time <= 100ms", "later time >100ms"};
410 for (unsigned int ir = 0; ir < c_rings; ir++) {
411 for (int wt = 0; wt < tzoom; wt++) {
412 std::string title = Form("inject time (%s), %s (%s)", stlabel[wt].data(), m_sring[ir].data(), m_suffix.data());
413 std::string hname = Form("htimezoom_%s_%s_%s", m_sring[ir].data(), stname[wt].data(), m_suffix.data());
414 htemp[ir][wt] = new TH1D(Form("%s", hname.data()), "", nt[wt], tzedges[wt], tzedges[wt + 1]);
415 htemp[ir][wt]->SetTitle(Form("%s;injection time(#mu-sec);entries", title.data()));
416 }
417 }
418}
419
420//------------------------------------
421void CDCDedxInjectTimeAlgorithm::checkStatistics(std::array<std::vector<TH1D*>, numdedx::nrings>& hvar)
422{
423 for (unsigned int ir = 0; ir < c_rings; ir++) {
424 for (unsigned int it = 3; it < m_tbins; it++) {
425 //check statiscs from 1-40ms
426 if (m_tedges[it] < 4e4 && hvar[ir][it]->Integral() < m_thersE) {
427 m_isminStat = true;
428 break;
429 } else continue;
430 }
431 }
432}
433
434//------------------------------------
435void CDCDedxInjectTimeAlgorithm::correctBinBias(std::map<int, std::vector<double>>& varcorr,
436 std::map<int, std::vector<double>>& var,
437 std::map<int, std::vector<double>>& time, TH1D*& htimes)
438{
439 //Deep copy OK
440 varcorr = var;
441
442 for (int ir = 0; ir < 2; ir++) {
443
444 for (int ix = varcorr[ir * 2].size(); ix -- > 0;) {
445
446 double var_thisbin = 1.0;
447 var_thisbin = var[ir * 2].at(ix);
448
449 double atime_thisbin = time[ir * 2].at(ix);
450 double ctime_thisbin = htimes->GetBinCenter(ix + 1);
451
452 if (atime_thisbin > 0 && atime_thisbin < 4e4 * 0.99) {
453
454 double var_nextbin = 1.0;
455 var_nextbin = var[ir * 2].at(ix + 1);
456 double var_diff = var_nextbin - var_thisbin;
457
458 double atime_nextbin = time[ir * 2].at(ix + 1);
459 double atime_diff = atime_nextbin - atime_thisbin;
460
461 double slope = (atime_diff > 0) ? var_diff / atime_diff : -1.0;
462
463 //extrapolation after veto only
464 if (var_diff > 0 && slope > 0)varcorr[ir * 2].at(ix) = var_thisbin + (ctime_thisbin - atime_thisbin) * (slope);
465 printf("\t %s ix = %d, center = %0.2f(%0.3f), var = %0.5f(%0.5f) \n", m_sring[ir].data(), ix, ctime_thisbin, atime_thisbin,
466 var_thisbin, varcorr[ir * 2].at(ix));
467 } else {
468 printf("\t %s --> ix = %d, center = %0.2f(%0.3f), var = %0.5f(%0.5f) \n", m_sring[ir].data(), ix, ctime_thisbin, atime_thisbin,
469 var_thisbin, varcorr[ir * 2].at(ix));
470 }
471 }
472 }
473}
474
475//-------------------------------------
476void CDCDedxInjectTimeAlgorithm::createPayload(std::array<double, numdedx::nrings>& scale,
477 std::map<int, std::vector<double>>& var,
478 std::map<int, std::vector<double>>& varscal, std::string svar)
479{
480 varscal = var;
481
482 B2INFO("CDCDedxInjectTimeAlgorithm: normalising constants with plateau");
483 for (unsigned int ir = 0; ir < c_rings; ir++) {
484 //scaling means with time >40ms
485 unsigned int msize = varscal[ir * 2].size();
486 int countsum = 0;
487 scale[ir] = 0;
488 for (unsigned int im = 0; im < msize; im++) {
489 double time = m_vtlocaledges.at(im);
490 double mean = varscal[ir * 2].at(im);
491 if (time > 4e4 && mean > 0) {
492 scale[ir] += mean;
493 countsum++;
494 }
495 }
496 if (countsum > 0 && scale[ir] > 0) {
497 scale[ir] /= countsum;
498 for (unsigned int im = 0; im < msize; im++) {
499 varscal[ir * 2].at(im) /= scale[ir];
500 }
501 }
502 }
503 if (m_isMerge && svar == "mean") {
504 //merge only no change in payload structure
505 bool incomp_bin = false;
506 std::vector<std::vector<double>> oldvectors;
507 if (m_DBInjectTime) oldvectors = m_DBInjectTime->getConstVector();
508 int vsize = oldvectors.size();
509 if (vsize != 6) incomp_bin = true;
510 else {
511 for (int iv = 0; iv < 2; iv++) {
512 if (oldvectors[iv * 3 + 1].size() != varscal[iv * 2].size()) incomp_bin = true;
513 }
514 }
515 if (!incomp_bin) {
516 B2INFO("CDCDedxInjectTimeAlgorithm: started merging relative constants");
517 for (int ir = 0; ir < 2; ir++) {//merging only means
518 unsigned int msize = varscal[ir * 2].size();
519 for (unsigned int im = 0; im < msize; im++) {
520 double relvalue = varscal[ir * 2].at(im);
521 double oldvalue = oldvectors[ir * 3 + 1].at(im);
522 double merged = oldvalue * relvalue;
523 printf("%s: rel %0.03f, old %0.03f, merged %0.03f\n", m_suffix.data(), relvalue, oldvalue, merged);
524 varscal[ir * 2].at(im) *= oldvectors[ir * 3 + 1].at(im) ;
525 }
526 }
527 } else B2ERROR("CDCDedxInjectTimeAlgorithm: found incompatible bins for merging");
528 } else B2INFO("CDCDedxInjectTimeAlgorithm: saving final (abs) calibration");
529}
530
531//------------------------------------
532void CDCDedxInjectTimeAlgorithm::plotBinLevelDist(std::array<std::vector<TH1D*>, numdedx::nrings>& hvar, std::string var)
533{
534 TCanvas cfit("cfit", "cfit", 1000, 500);
535 cfit.Divide(2, 1);
536 std::stringstream psname_fit;
537 psname_fit << Form("%s_%s_%s.pdf[", m_prefix.data(), var.data(), m_suffix.data());
538 cfit.Print(psname_fit.str().c_str());
539 psname_fit.str("");
540 psname_fit << Form("%s_%s_%s.pdf", m_prefix.data(), var.data(), m_suffix.data());
541 for (unsigned int it = 0; it < m_tbins; it++) {
542 for (unsigned int ir = 0; ir < c_rings; ir++) {
543 cfit.cd(ir + 1);
544 hvar[ir][it]->SetFillColorAlpha(ir + 5, 0.25);
545 hvar[ir][it]->Draw();
546 }
547 cfit.Print(psname_fit.str().c_str());
548 }
549 psname_fit.str("");
550 psname_fit << Form("%s_%s_%s.pdf]", m_prefix.data(), var.data(), m_suffix.data());
551 cfit.Print(psname_fit.str().c_str());
552}
553
554//----------------------------------------
556{
557 // draw event and track statistics
558 TCanvas cestat("cestat", "cestat", 1000, 500);
559 cestat.SetBatch(kTRUE);
560 cestat.Divide(2, 1);
561
562 cestat.cd(1);
563 auto hestats = getObjectPtr<TH1I>("hestats");
564 if (hestats) {
565 hestats->SetName(Form("hestats_%s", m_suffix.data()));
566 hestats->SetStats(0);
567 hestats->Draw("hist text");
568 }
569 cestat.cd(2);
570 auto htstats = getObjectPtr<TH1I>("htstats");
571 if (htstats) {
572 htstats->SetName(Form("htstats_%s", m_suffix.data()));
573 htstats->SetStats(0);
574 htstats->Draw("hist text");
575 }
576 cestat.Print(Form("%s_eventstat_%s.pdf", m_prefix.data(), m_suffix.data()));
577}
578
579//------------------------------------
580void CDCDedxInjectTimeAlgorithm::plotInjectionTime(std::array<std::array<TH1D*, 3>, numdedx::nrings>& hvar)
581{
582 TCanvas ctzoom("ctzoom", "ctzoom", 1500, 450);
583 ctzoom.SetBatch(kTRUE);
584 ctzoom.Divide(3, 1);
585 for (int wt = 0; wt < 3; wt++) {
586 ctzoom.cd(wt + 1);
587 if (wt == 2) gPad->SetLogy();
588 for (unsigned int ir = 0; ir < c_rings; ir++) {
589 hvar[ir][wt]->SetStats(0);
590 hvar[ir][wt]->SetFillColorAlpha(5 + ir, 0.20);
591 if (ir == 0) {
592 double max1 = hvar[ir][wt]->GetMaximum();
593 double max2 = hvar[c_rings - 1][wt]->GetMaximum();
594 if (max2 > max1) hvar[ir][wt]->SetMaximum(max2 * 1.05);
595 hvar[ir][wt]->Draw("");
596 } else hvar[ir][wt]->Draw("same");
597 }
598 }
599 ctzoom.Print(Form("%s_timezoom_%s.pdf]", m_prefix.data(), m_suffix.data()));
600}
601
602//------------------------------------
603void CDCDedxInjectTimeAlgorithm::plotRelConstants(std::map<int, std::vector<double>>& vmeans,
604 std::map<int, std::vector<double>>& vresos, std::map<int, std::vector<double>>& corr, std::string svar)
605{
606 std::string sname[3] = {"mean", "reso"};
607 const int lcolors[c_rings] = {2, 4};
608
609 TH1D* hmean[c_rings], *hcorr[c_rings];
610 TH1D* hreso[c_rings];
611
612 TCanvas* cconst[2];
613 for (int ic = 0; ic < 2; ic++) {
614 cconst[ic] = new TCanvas(Form("c%sconst", sname[ic].data()), "", 900, 500);
615 cconst[ic]->SetGridy(1);
616 }
617
618 TLegend* mleg = new TLegend(0.50, 0.54, 0.80, 0.75, NULL, "brNDC");
619 mleg->SetBorderSize(0);
620 mleg->SetFillStyle(0);
621
622 for (unsigned int ir = 0; ir < c_rings; ir++) {
623
624 std::string mtitle = Form("#mu(dedx), relative const. compare, (%s)", m_suffix.data());
625 hmean[ir] = new TH1D(Form("hmean_%s_%s", m_suffix.data(), m_sring[ir].data()), "", m_tbins, 0, m_tbins);
626
627 std::string rtitle = Form("#sigma(#chi), relative const. compare, (%s)", m_suffix.data());
628 hreso[ir] = new TH1D(Form("hreso_%s_%s", m_suffix.data(), m_sring[ir].data()), "", m_tbins, 0, m_tbins);
629
630 hcorr[ir] = new TH1D(Form("hcorr_%s_%s", m_suffix.data(), m_sring[ir].data()), "", m_tbins, 0, m_tbins);
631 if (svar == "chi") {
632 hmean[ir]->SetTitle(Form("%s;injection time(#mu-second);#mu (#chi-fit)", rtitle.data()));
633 hreso[ir]->SetTitle(Form("%s;injection time(#mu-second);#sigma (#chi-fit)", rtitle.data()));
634 hcorr[ir]->SetTitle(Form("%s;injection time(#mu-second);#sigma (#chi-fit bin-bais-corr)", rtitle.data()));
635
636 } else {
637 hmean[ir]->SetTitle(Form("%s;injection time(#mu-second);#mu (dedx-fit)", mtitle.data()));
638 hreso[ir]->SetTitle(Form("%s;injection time(#mu-second);#sigma (dedx-fit)", mtitle.data()));
639 hcorr[ir]->SetTitle(Form("%s;injection time(#mu-second);#mu (dedx-fit, bin-bais-corr)", mtitle.data()));
640 }
641
642
643 for (unsigned int it = 0; it < m_tbins; it++) {
644
645 std::string label = getTimeBinLabel(m_tedges[it], it);
646
647 hmean[ir]->SetBinContent(it + 1, vmeans[ir * 2].at(it));
648 hmean[ir]->SetBinError(it + 1, vmeans[ir * 2 + 1].at(it));
649 hmean[ir]->GetXaxis()->SetBinLabel(it + 1, label.data());
650
651 hreso[ir]->SetBinContent(it + 1, vresos[ir * 2].at(it));
652 hreso[ir]->SetBinError(it + 1, vresos[ir * 2 + 1].at(it));
653 hreso[ir]->GetXaxis()->SetBinLabel(it + 1, label.data());
654
655 hcorr[ir]->SetBinContent(it + 1, corr[ir * 2].at(it));
656 hcorr[ir]->SetBinError(it + 1, corr[ir * 2 + 1].at(it));
657 hcorr[ir]->GetXaxis()->SetBinLabel(it + 1, label.data());
658 }
659
660 mleg->AddEntry(hmean[ir], Form("%s", m_sring[ir].data()), "lep");
661 cconst[0]->cd();
662 if (svar == "chi") setHistStyle(hmean[ir], lcolors[ir], ir + 24, -0.60, 0.60);
663 else if (svar == "dedx") setHistStyle(hmean[ir], lcolors[ir], ir + 24, 0.60, 1.10);
664 else setHistStyle(hmean[ir], lcolors[ir], ir + 24, 0.9, 1.10);
665
666 if (ir == 0)hmean[ir]->Draw("");
667 else hmean[ir]->Draw("same");
668 if (svar == "dedx") {
669 mleg->AddEntry(hcorr[ir], Form("%s (bin-bias-corr)", m_sring[ir].data()), "lep");
670 setHistStyle(hcorr[ir], lcolors[ir], ir + 20, 0.60, 1.10);
671 hcorr[ir]->Draw("same");
672 }
673 if (ir == 1)mleg->Draw("same");
674
675 cconst[1]->cd();
676 if (svar == "chi") setHistStyle(hreso[ir], lcolors[ir], ir + 24, 0.5, 1.50);
677 else setHistStyle(hreso[ir], lcolors[ir], ir + 24, 0.0, 0.15);
678 if (ir == 0)hreso[ir]->Draw("");
679 else hreso[ir]->Draw("same");
680 if (svar == "chi") {
681 mleg->AddEntry(hcorr[ir], Form("%s (bin-bias-corr)", m_sring[ir].data()), "lep");
682 setHistStyle(hcorr[ir], lcolors[ir], ir + 20, 0.5, 1.50);
683 hcorr[ir]->Draw("same");
684 }
685 if (ir == 1)mleg->Draw("same");
686 }
687
688 for (int ic = 0; ic < 2; ic++) {
689 cconst[ic]->SaveAs(Form("%s_relconst_%s_%s_%s.pdf", m_prefix.data(), sname[ic].data(), svar.data(), m_suffix.data()));
690 cconst[ic]->SaveAs(Form("%s_relconst_%s_%s_%s.root", m_prefix.data(), sname[ic].data(), svar.data(), m_suffix.data()));
691 delete cconst[ic];
692 }
693 for (int ic = 0; ic < 2; ic++) {
694
695 delete hmean[ic];
696 delete hreso[ic];
697 delete hcorr[ic];
698 }
699}
700
701//------------------------------------
702void CDCDedxInjectTimeAlgorithm::plotTimeStat(std::array<std::vector<TH1D*>, numdedx::nrings>& htime)
703{
704 const int lcolors[c_rings] = {2, 4};
705
706 TH1D* htimestat[c_rings];
707
708 TCanvas* cconst = new TCanvas("ctimestatconst", "", 900, 500);
709 cconst->SetGridy(1);
710
711 TLegend* rleg = new TLegend(0.40, 0.60, 0.80, 0.72, NULL, "brNDC");
712 rleg->SetBorderSize(0);
713 rleg->SetFillStyle(0);
714
715 for (unsigned int ir = 0; ir < c_rings; ir++) {
716
717 std::string title = Form("injection time, her-ler comparison, (%s)", m_suffix.data());
718 htimestat[ir] = new TH1D(Form("htimestat_%s_%s", m_suffix.data(), m_sring[ir].data()), "", m_tbins, 0, m_tbins);
719 htimestat[ir]->SetTitle(Form("%s;injection time(#mu-second);norm. entries", title.data()));
720
721 for (unsigned int it = 0; it < m_tbins; it++) {
722 std::string label = getTimeBinLabel(m_tedges[it], it);
723 htimestat[ir]->SetBinContent(it + 1, htime[ir][it]->Integral());
724 htimestat[ir]->SetBinError(it + 1, 0);
725 htimestat[ir]->GetXaxis()->SetBinLabel(it + 1, label.data());
726 }
727
728 cconst->cd();
729 double norm = htimestat[ir]->GetMaximum();
730 rleg->AddEntry(htimestat[ir], Form("%s (scaled with %0.02f)", m_sring[ir].data(), norm), "lep");
731 htimestat[ir]->Scale(1.0 / norm);
732 setHistStyle(htimestat[ir], lcolors[ir], ir + 20, 0.0, 1.10);
733 htimestat[ir]->SetFillColorAlpha(lcolors[ir], 0.30);
734 if (ir == 0) htimestat[ir]->Draw("hist");
735 else htimestat[ir]->Draw("hist same");
736 if (ir == 1)rleg->Draw("same");
737 }
738
739 cconst->SaveAs(Form("%s_relconst_timestat_%s.pdf", m_prefix.data(), m_suffix.data()));
740 cconst->SaveAs(Form("%s_relconst_timestat_%s.root", m_prefix.data(), m_suffix.data()));
741 delete cconst;
742}
743
744//------------------------------------
745void CDCDedxInjectTimeAlgorithm::plotFinalConstants(std::map<int, std::vector<double>>& vmeans,
746 std::map<int, std::vector<double>>& vresos,
747 std::array<double, numdedx::nrings>& scale, std::array<double, numdedx::nrings>& scale_reso)
748{
749
750 std::vector<std::vector<double>> oldvectors;
751 if (m_DBInjectTime)oldvectors = m_DBInjectTime->getConstVector();
752
753 const int c_type = 2; //old and new
754 std::string sname[c_rings] = {"mean", "reso"};
755 std::string stype[c_type] = {"new", "old"};
756 const int lcolors[c_rings] = {2, 4};
757 const int lmarker[c_type] = {20, 24}; //+2 for different rings
758
759 TH1D* hmean[c_rings][c_type], *hreso[c_rings][c_type];
760
761 TCanvas* cconst[c_rings];
762 for (int ic = 0; ic < 2; ic++) {
763 cconst[ic] = new TCanvas(Form("c%sconst", sname[ic].data()), "", 900, 500);
764 cconst[ic]->SetGridy(1);
765 }
766
767 TLegend* mleg = new TLegend(0.50, 0.54, 0.80, 0.75, NULL, "brNDC");
768 mleg->SetBorderSize(0);
769 mleg->SetFillStyle(0);
770
771 TLegend* rleg = new TLegend(0.50, 0.54, 0.80, 0.75, NULL, "brNDC");
772 rleg->SetBorderSize(0);
773 rleg->SetFillStyle(0);
774
775 for (unsigned int ip = 0; ip < c_type; ip++) {
776
777 for (unsigned int ir = 0; ir < c_rings; ir++) {
778
779 std::string hname = Form("hfmean_%s_%s_%s", m_sring[ir].data(), stype[ip].data(), m_suffix.data());
780 hmean[ir][ip] = new TH1D(hname.data(), "", m_tbins, 0, m_tbins);
781 std::string title = Form("#mu(dedx), final-mean-compare (%s)", m_suffix.data());
782 hmean[ir][ip]->SetTitle(Form("%s;injection time(#mu-second);#mu (dedx-fit)", title.data()));
783
784 hname = Form("hfreso_%s_%s_%s", m_sring[ir].data(), stype[ip].data(), m_suffix.data());
785 hreso[ir][ip] = new TH1D(hname.data(), "", m_tbins, 0, m_tbins);
786 title = Form("#sigma(#chi), final-reso-compare (%s)", m_suffix.data());
787 hreso[ir][ip]->SetTitle(Form("%s;injection time(#mu-second);#sigma (#chi-fit)", title.data()));
788
789 for (unsigned int it = 0; it < m_tbins; it++) {
790
791 std::string label = getTimeBinLabel(m_tedges[it], it);
792 double mean = 0.0, reso = 0.0;
793 if (ip == 0) {
794 mean = m_vinjPayload[ir * 3 + 1].at(it);
795 //reso is reso/mu (reso is relative so mean needs to be relative)
796 reso = m_vinjPayload[ir * 3 + 2].at(it);
797
798 } else {
799 mean = oldvectors[ir * 3 + 1].at(it);
800 reso = oldvectors[ir * 3 + 2].at(it);
801
802 }
803
804 //old payloads
805 hmean[ir][ip]->SetBinContent(it + 1, mean);
806 hmean[ir][ip]->SetBinError(it + 1, vmeans[ir * 2 + 1].at(it));
807 hmean[ir][ip]->GetXaxis()->SetBinLabel(it + 1, label.data());
808
809 hreso[ir][ip]->SetBinContent(it + 1, reso);
810 hreso[ir][ip]->SetBinError(it + 1, vresos[ir * 2 + 1].at(it));
811 hreso[ir][ip]->GetXaxis()->SetBinLabel(it + 1, label.data());
812 }
813
814 cconst[0]->cd();
815 if (ip == 1)mleg->AddEntry(hmean[ir][ip], Form("%s, %s", m_sring[ir].data(), stype[ip].data()), "lep");
816 else mleg->AddEntry(hmean[ir][ip], Form("%s, %s (scaled by %0.03f)", m_sring[ir].data(), stype[ip].data(), scale[ir]), "lep");
817 setHistStyle(hmean[ir][ip], lcolors[ir], lmarker[ip] + ir * 2, 0.60, 1.05);
818 if (ir == 0 && ip == 0)hmean[ir][ip]->Draw("");
819 else hmean[ir][ip]->Draw("same");
820 if (ir == 1 && ip == 1)mleg->Draw("same");
821
822 cconst[1]->cd();
823 if (ip == 1)rleg->AddEntry(hreso[ir][ip], Form("%s, %s", m_sring[ir].data(), stype[ip].data()), "lep");
824 else rleg->AddEntry(hreso[ir][ip], Form("%s, %s (scaled by %0.03f)", m_sring[ir].data(), stype[ip].data(), scale_reso[ir]), "lep");
825 setHistStyle(hreso[ir][ip], lcolors[ir], lmarker[ip] + ir * 2, 0.6, 1.9);
826 if (ir == 0 && ip == 0)hreso[ir][ip]->Draw("");
827 else hreso[ir][ip]->Draw("same");
828 if (ir == 1 && ip == 1)rleg->Draw("same");
829 }
830 }
831
832 for (int ic = 0; ic < 2; ic++) {
833 cconst[ic]->SaveAs(Form("%s_finalconst_%s_%s.pdf", m_prefix.data(), sname[ic].data(), m_suffix.data()));
834 cconst[ic]->SaveAs(Form("%s_finalconst_%s_%s.root", m_prefix.data(), sname[ic].data(), m_suffix.data()));
835 delete cconst[ic];
836 }
837
838 for (unsigned int ir = 0; ir < c_rings; ir++) {
839 for (unsigned int ip = 0; ip < c_type; ip++) {
840 delete hmean[ir][ip];
841 delete hreso[ir][ip];
842 }
843 }
844}
845
846//------------------------------------
847double CDCDedxInjectTimeAlgorithm::getCorrection(unsigned int ring, unsigned int time,
848 std::map<int, std::vector<double>>& vmeans)
849{
850
851 unsigned int iv = ring * 2;
852
853 unsigned int sizet = m_vtlocaledges.size(); //time
854
855 std::vector<unsigned int> tedges(sizet); //time edges array
856 std::copy(m_vtlocaledges.begin(), m_vtlocaledges.end(), tedges.begin());
857
858 if (time >= 5e6) time = 5e6 - 10;
859 unsigned int it = m_DBInjectTime->getTimeBin(tedges, time);
860
861 double center = 0.5 * (m_vtlocaledges.at(it) + m_vtlocaledges.at(it + 1));
862
863 //no corr before veto bin (usually one or two starting bin)
864 //intrapolation for entire range except
865 //--extrapolation (for first half and last half of intended bin)
866 int thisbin = it, nextbin = it;
867 if (center != time && it > 0) {
868
869 if (time < center) {
870 thisbin = it - 1;
871 } else {
872 if (it < sizet - 2)nextbin = it + 1;
873 else thisbin = it - 1;
874 }
875
876 if (it <= 2) {
877 double diff = vmeans[iv].at(2) - vmeans[iv].at(1) ;
878 if (diff < -0.015) { //difference above 1.0%
879 thisbin = it;
880 if (it == 1) nextbin = it;
881 else nextbin = it + 1;
882 } else {
883 if (it == 1) {
884 thisbin = it;
885 nextbin = it + 1;
886 }
887 }
888 }
889 }
890
891 double thisdedx = vmeans[iv].at(thisbin);
892 double nextdedx = vmeans[iv].at(nextbin);
893
894 double thistime = 0.5 * (m_vtlocaledges.at(thisbin) + m_vtlocaledges.at(thisbin + 1));
895 double nexttime = 0.5 * (m_vtlocaledges.at(nextbin) + m_vtlocaledges.at(nextbin + 1));
896
897 double newdedx = vmeans[iv].at(it);
898 if (thisbin != nextbin)
899 newdedx = thisdedx + ((nextdedx - thisdedx) / (nexttime - thistime)) * (time - thistime);
900
901 return newdedx;
902}
bool m_isMerge
merge payload when rel constant
void plotBinLevelDist(std::array< std::vector< TH1D * >, numdedx::nrings > &hvar, std::string var)
function to draw dedx, chi and time dist.
std::vector< double > m_vtlocaledges
internal time vector
double m_sigmaR
fit dedx dist in sigma range
void correctBinBias(std::map< int, std::vector< double > > &varcorr, std::map< int, std::vector< double > > &var, std::map< int, std::vector< double > > &time, TH1D *&htimes)
function to correct dedx mean/reso and return corrected vector map
void defineHisto(std::array< std::vector< TH1D * >, numdedx::nrings > &htemp, std::string var)
function to define histograms for dedx and time dist.
double getCorrection(unsigned int ring, unsigned int time, std::map< int, std::vector< double > > &vmeans)
function to get the correction factor of mean
std::string m_prefix
string prefix for plot names
void getExpRunInfo()
function to get exp/run information (payload object, plotting)
void setHistStyle(TH1D *&htemp, const int ic, const int is, const double min, const double max)
function to set histogram cosmetics
void plotFinalConstants(std::map< int, std::vector< double > > &vmeans, std::map< int, std::vector< double > > &vresos, std::array< double, numdedx::nrings > &scale, std::array< double, numdedx::nrings > &scale_reso)
function to final constant from merging or abs fits
int m_countR
a hack for running functions once
double * m_tedges
internal time array (copy of vtlocaledges)
bool m_ismakePlots
produce plots for monitoring
std::string m_suffix
string suffix for object names
static const int c_rings
injection ring constants
void createPayload(std::array< double, numdedx::nrings > &scale, std::map< int, std::vector< double > > &vmeans, std::map< int, std::vector< double > > &varscal, std::string svar)
function to store payloads after full calibration
void plotInjectionTime(std::array< std::array< TH1D *, 3 >, numdedx::nrings > &hvar)
function to injection time distributions (HER/LER in three bins)
void plotTimeStat(std::array< std::vector< TH1D * >, numdedx::nrings > &htime)
function to draw time stats
void fitGaussianWRange(TH1D *&temphist, fstatus &status)
function to perform gauss fit for input histogram
void defineTimeBins()
function to set/reset time bins
void plotRelConstants(std::map< int, std::vector< double > > &vmeans, std::map< int, std::vector< double > > &vresos, std::map< int, std::vector< double > > &corr, std::string svar)
function to relative constant from dedx fit mean and chi fit reso
void defineTimeHisto(std::array< std::array< TH1D *, 3 >, numdedx::nrings > &htemp)
function to define injection time bins histograms (monitoring only)
void deleteHisto(std::array< std::vector< TH1D * >, numdedx::nrings > &htemp)
function to delete histograms for dedx and time dist.
void deleteTimeHisto(std::array< std::array< TH1D *, 3 >, numdedx::nrings > &htemp)
function to define injection time bins histograms (monitoring only)
void plotEventStats()
function to draw event/track statistics plots
virtual EResult calibrate() override
CDC dE/dx Injection time algorithm.
void checkStatistics(std::array< std::vector< TH1D * >, numdedx::nrings > &hvar)
check statistics for obtaining calibration const.
std::vector< double > m_vtedges
external time vector
std::vector< std::vector< double > > m_vinjPayload
vector to store payload values
bool m_isminStat
flag to merge runs for statistics thershold
CDCDedxInjectTimeAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
std::array< std::string, numdedx::nrings > m_sring
injection ring name
DBObjPtr< CDCDedxInjectionTime > m_DBInjectTime
Injection time DB object.
int m_thersE
min tracks to start calibration
void getMeanReso(std::array< std::vector< TH1D * >, numdedx::nrings > &hvar, std::map< int, std::vector< double > > &vmeans, std::map< int, std::vector< double > > &vresos)
function to get mean and reso of histogram
std::string getTimeBinLabel(const double &tedges, const int &it)
function to return time label for histograms labeling
dE/dx injection time calibration constants
Class to hold the prediction of mean as a function of beta-gamma (bg)
double getMean(double bg)
Return the predicted mean value as a function of beta-gamma (bg)
void setParameters(std::string infile)
set the parameters from file
Class to hold the prediction of resolution depending dE/dx, nhit, and cos(theta)
double ionzPrediction(double dedx)
Return sigma from the ionization parameterization.
double cosPrediction(double cos)
Return sigma from the cos parameterization.
double nhitPrediction(double nhit)
Return sigma from the nhit parameterization.
void setParameters(std::string infile)
set the parameters from file
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(....
static const ChargedStable electron
electron particle
Definition Const.h:659
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.