Belle II Software development
CDCDedx1DCellAlgorithm.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/CDCDedx1DCellAlgorithm.h>
10
11#include <TCanvas.h>
12#include <TLegend.h>
13#include <TMath.h>
14#include <TPad.h>
15#include <TRandom.h>
16#include <TStyle.h>
17
18#include <cmath>
19
20using namespace Belle2;
21
22//-----------------------------------------------------------------
23// Implementation
24//-----------------------------------------------------------------
26 CalibrationAlgorithm("CDCDedxElectronCollector"),
27 m_eaMin(-TMath::Pi() / 2),
28 m_eaMax(+TMath::Pi() / 2),
29 m_eaB(316),
30 m_dedxMin(0.0),
31 m_dedxMax(5.0),
32 m_dedxBin(250),
33 m_ptMax(8.0),
34 m_cosMax(1.0),
35 m_truncMin(0.05),
36 m_truncMax(0.75),
37 m_binSplit(3),
38 m_chargeType(0),
39 m_adjustFac(1.00),
40 isFixTrunc(false),
41 isVarBins(true),
42 isRotSymm(false),
43 isMakePlots(true),
44 isPrintLog(false),
45 isMerge(true),
46 m_suffix("")
47{
48 // Set module properties
49 setDescription("A calibration algorithm for the CDC dE/dx entrance angle cleanup correction");
50}
51
52//-----------------------------------------------------------------
53// Run the calibration
54//-----------------------------------------------------------------
56{
57
59
60 if (!m_DBOneDCell.isValid())
61 B2FATAL("There is no valid previous payload for CDCDedx1DCell");
62
63 //reading radiative electron collector TREE
64 auto ttree = getObjectPtr<TTree>("tree");
65 if (!ttree) return c_NotEnoughData;
66
67 std::vector<double>* dedxhit = 0, *enta = 0;
68 std::vector<int>* layer = 0;
69 double pt = 0, costh = 0;
70
71 ttree->SetBranchAddress("dedxhit", &dedxhit);
72 ttree->SetBranchAddress("entaRS", &enta);
73 ttree->SetBranchAddress("layer", &layer);
74 ttree->SetBranchAddress("pt", &pt);
75 ttree->SetBranchAddress("costh", &costh);
76
77 //repair nbins if they are not divisible accordingly
80
81 //Settings of variables bins
83
84 if (isPrintLog) {
85 B2INFO("inner layers bins: " << m_eaBinLocal[0]);
86 B2INFO("outer layers bins: " << m_eaBinLocal[1]);
87 }
88
89 // dedxhit vector to store dE/dx values for each enta bin
90 std::vector<TH1D*> hdedxhit[2];
91 TH1D* hdedxlay[2];
92 TH1D* hentalay[2];
93
94 TH2D* hptcosth = new TH2D("hptcosth", "pt vs costh dist;pt;costh", 1000, -8.0, 8.0, 1000, -1.0, 1.0);
95
96 defineHisto(hdedxhit, hdedxlay, hentalay);
97
98 //Star filling histogram defined above
99 for (int i = 0; i < ttree->GetEntries(); ++i) {
100
101 ttree->GetEvent(i);
102
103 if (std::abs(costh) > m_cosMax) continue;
104
105 // remove wide angle bhabha tracks
106 // double mom = pt/sqrt(1-costh*costh);
107 // if(abs(pt)<2.4 && abs(mom)>3.6)continue;
108
109 if (std::abs(pt) > m_ptMax) continue;
110
111 //change to random 10%
112 int rand = gRandom->Integer(100);
113 if (rand < 10) hptcosth->Fill(pt, costh);
114
115 for (unsigned int j = 0; j < dedxhit->size(); ++j) {
116
117 if (dedxhit->at(j) == 0) continue;
118
119 double entaval = enta->at(j);
120 //Mapped bin corresponds to entaval
121 int ibin = std::floor((entaval - m_eaMin) / m_eaBW);
122 if (ibin < 0 || ibin > m_eaBin) continue;
123
124 int mL = -1;
125 if (layer->at(j) < 8)mL = 0;
126 else mL = 1;
127
128 hdedxlay[mL]->Fill(dedxhit->at(j));
129 if (rand < 10) hentalay[mL]->Fill(entaval);
130
131 int jbinea = ibin;
132 if (isVarBins) jbinea = m_binIndex[mL].at(ibin);
133 hdedxhit[mL][jbinea]->Fill(dedxhit->at(j));
134 }
135 }
136
137 for (int il = 0; il < 2; il++) {
138
139 int minlay = 0, maxlay = 0;
140
141 if (isFixTrunc) {
142 getTruncatedBins(hdedxlay[il], minlay, maxlay);
143 hdedxlay[il]->SetTitle(Form("%s;%d;%d", hdedxlay[il]->GetTitle(), minlay, maxlay));
144 }
145
146 std::vector<double>tempconst;
147 tempconst.reserve(m_eaBinLocal[il]);
148
149 for (int iea = 0; iea < m_eaBinLocal[il]; iea++) {
150
151 int jea = iea;
152
153 // rotation symmtery for 1<->3 and 4<->2 but only symmetric bin
154 if (!isVarBins && isRotSymm) jea = rotationalBin(m_eaBinLocal[il], jea);
155
156 TH1D* htemp = (TH1D*)hdedxhit[il][jea]->Clone(Form("h_%s_b%d_c", m_label[il].data(), jea));
157
158 int minbin = 1, maxbin = 1;
159 if (isFixTrunc) {
160 minbin = minlay;
161 maxbin = maxlay;
162 } else {
163 //extract truncation window per bin
164 getTruncatedBins(htemp, minbin, maxbin);
165 }
166
167 double dedxmean;
168 dedxmean = getTruncationMean(htemp, minbin, maxbin);
169 tempconst.push_back(dedxmean);
170
171 hdedxhit[il][iea]->SetTitle(Form("%s, #mu_{truc} = %0.5f;%d;%d", hdedxhit[il][iea]->GetTitle(), dedxmean, minbin, maxbin));
172 }
173
174 //Expending constants
175 std::vector<double>layerconst;
176 layerconst.reserve(m_eaBin);
177
178 for (int iea = 0; iea < m_eaBin; iea++) {
179 int jea = iea;
180 if (isVarBins) jea = m_binIndex[il].at(iea);
181 layerconst.push_back(tempconst.at(jea));
182 }
183
184 // plot the rel constants var/sym bins
185 if (isMakePlots) plotRelConst(tempconst, layerconst, il);
186 m_onedcors.push_back(layerconst);
187
188 layerconst.clear();
189 tempconst.clear();
190 }
191
192 //Saving final constants
194
195 if (isMakePlots) {
196
197 //1. dE/dx dist. for entrance angle bins
198 plotdedxHist(hdedxhit);
199
200 //3. Inner and Outer layer dE/dx distributions
201 plotLayerDist(hdedxlay);
202
203 //4. entrance angle distribution sym/var bins
204 plotQaPars(hentalay, hptcosth);
205
206 //6. draw the final constants
208
209 //7. plot statistics related plots here
211 }
212
213 for (int il = 0; il < 2; il++) {
214 delete hentalay[il];
215 delete hdedxlay[il];
216 for (int iea = 0; iea < m_eaBinLocal[il]; iea++)
217 delete hdedxhit[il][iea];
218 }
219
220 delete hptcosth;
221 m_eaBinLocal.clear();
222 for (int il = 0; il < 2; il++) {
223 m_binValue[il].clear();
224 m_binIndex[il].clear();
225 }
226 return c_OK;
227}
228
229//--------------------------------------------------
231{
232
233 int cruns = 0;
234 for (auto expRun : getRunList()) {
235 if (cruns == 0) B2INFO("CDCDedxBadWires: start exp " << expRun.first << " and run " << expRun.second << "");
236 cruns++;
237 }
238
239 const auto erStart = getRunList()[0];
240 int estart = erStart.first;
241 int rstart = erStart.second;
242
243 const auto erEnd = getRunList()[cruns - 1];
244 int eend = erEnd.first;
245 int rend = erEnd.second;
246
247 updateDBObjPtrs(1, rstart, estart);
248
249 m_runExp = Form("Range (%d:%d,%d:%d)", estart, rstart, eend, rend);
250 if (m_suffix.length() > 0) m_suffix = Form("%s_e%d_r%dr%d", m_suffix.data(), estart, rstart, rend);
251 else m_suffix = Form("e%d_r%dr%d", estart, rstart, rend);
252}
253
254//--------------------------------------------------
256{
257
258 std::map<int, std::vector<double>> bounds;
259 std::map<int, std::vector<int>> steps;
260
261 const std::array<int, 2> nDev{8, 4};
262 bounds[0] = {0, 108, 123, 133, 158, 183, 193, 208, 316}; //il boundaries
263 steps[0] = {9, 3, 2, 1, 1, 2, 3, 9}; //il steps
264 bounds[1] = {0, 38, 158, 278, 316}; //OL boundaries
265 steps[1] = {2, 1, 1, 2}; //OL steps
266
267 for (int il = 0; il < 2; il++) {
268
269 for (int ibin = 0; ibin <= nDev[il]; ibin++) bounds[il][ibin] = bounds[il][ibin] * m_binSplit;
270
271 int ieaprime = -1, temp = -99, ibin = 0;
272
273 double pastbin = m_eaMin;
274 m_binValue[il].push_back(pastbin);
275
276 for (int iea = 0; iea < m_eaBin; iea++) {
277
278 if (isVarBins) {
279 if (iea % int(bounds[il][ibin + 1]) == 0 && iea > 0) ibin++;
280 int diff = iea - int(bounds[il][ibin]);
281 if (diff % steps[il][ibin] == 0) ieaprime++;
282 } else ieaprime = iea;
283
284 m_binIndex[il].push_back(ieaprime);
285
286 if (ieaprime != temp) {
287 double binwidth = m_eaBW;
288 if (isVarBins) binwidth = m_eaBW * steps[il][ibin];
289 double binvalue = pastbin + binwidth;
290 pastbin = binvalue;
291 if (std::abs(binvalue) < 1e-5)binvalue = 0;
292 m_binValue[il].push_back(binvalue);
293 }
294 temp = ieaprime;
295 }
296 m_eaBinLocal.push_back(int(m_binValue[il].size()) - 1) ;
297 }
298 if (isMakePlots) plotMergeFactor(bounds, nDev, steps);
299}
300
301//--------------------------------------------------
302void CDCDedx1DCellAlgorithm::defineHisto(std::vector<TH1D*> hdedxhit[2], TH1D* hdedxlay[2], TH1D* hentalay[2])
303{
304 for (int il = 0; il < 2; il++) {
305
306 std::string title = Form("dedxhit dist (%s): %s ; dedxhit;entries", m_label[il].data(), m_runExp.data());
307 hdedxlay[il] = new TH1D(Form("hdedxlay%s", m_label[il].data()), "", m_dedxBin, m_dedxMin, m_dedxMax);
308 hdedxlay[il]->SetTitle(Form("%s", title.data()));
309
310 Double_t* nvarBins;
311 nvarBins = &m_binValue[il][0];
312
313 if (isVarBins) title = Form("entaRS dist (variable bins): %s: (%s); entaRS (#alpha);entries", m_label[il].data(), m_runExp.data());
314 else title = Form("entaRS dist (sym. bins): %s: (%s); entaRS (#alpha);entries", m_label[il].data(), m_runExp.data());
315
316 hentalay[il] = new TH1D(Form("hentalay%s", m_label[il].data()), "", m_eaBinLocal[il], nvarBins);
317 hentalay[il]->SetTitle(Form("%s", title.data()));
318
319 for (int iea = 0; iea < m_eaBinLocal[il]; iea++) {
320
321 double min = m_binValue[il].at(iea);
322 double max = m_binValue[il].at(iea + 1);
323 double width = max - min;
324
325 if (isPrintLog) B2INFO("bin: " << iea << " ], min:" << min << " , max: " << max << " , width: " << width);
326
327 title = Form("%s: entaRS = (%0.03f to %0.03f)", m_label[il].data(), min, max);
328 hdedxhit[il].push_back(new TH1D(Form("hdedxhit_%s_bin%d", m_label[il].data(), iea), "", m_dedxBin, m_dedxMin, m_dedxMax));
329 hdedxhit[il][iea]->SetTitle(Form("%s", title.data()));
330 }
331 }
332}
333
334//--------------------------------------------------
335void CDCDedx1DCellAlgorithm::getTruncatedBins(TH1D* hist, int& binlow, int& binhigh)
336{
337
338 //calculating truncation average
339 double sum = hist->Integral();
340 if (sum <= 0 || hist->GetNbinsX() <= 0) {
341 binlow = 1; binhigh = 1;
342 return ;
343 }
344
345 binlow = 1.0; binhigh = 1.0;
346 double sumPer5 = 0.0, sumPer75 = 0.0;
347 for (int ibin = 1; ibin <= hist->GetNbinsX(); ibin++) {
348 double bcdedx = hist->GetBinContent(ibin);
349 if (sumPer5 <= m_truncMin * sum) {
350 sumPer5 += bcdedx;
351 binlow = ibin;
352 }
353 if (sumPer75 <= m_truncMax * sum) {
354 sumPer75 += bcdedx;
355 binhigh = ibin;
356 }
357 }
358 return;
359}
360
361//--------------------------
362double CDCDedx1DCellAlgorithm::getTruncationMean(TH1D* hist, int binlow, int binhigh)
363{
364
365 //calculating truncation average
366 if (hist->Integral() < 100) return 1.0;
367
368 if (binlow <= 0 || binhigh > hist->GetNbinsX())return 1.0;
369
370 double binweights = 0., sumofbc = 0.;
371 for (int ibin = binlow; ibin <= binhigh; ibin++) {
372 double bcdedx = hist->GetBinContent(ibin);
373 if (bcdedx > 0) {
374 binweights += (bcdedx * hist->GetBinCenter(ibin));
375 sumofbc += bcdedx;
376 }
377 }
378 if (sumofbc > 0) return binweights / sumofbc;
379 else return 1.0;
380}
381
382//--------------------------------------------------
384{
385
386 B2INFO("dE/dx one cell calibration: Generating payloads");
387
388 for (unsigned int il = 0; il < 2; il++) {
389
390 if (isMerge) {
391 unsigned int nbins = m_DBOneDCell->getNBins(il);
392
393 if (int(nbins) != m_eaBin)
394 B2ERROR("merging failed because of unmatch bins (old " << m_eaBin << " new " << nbins << ")");
395
396 for (unsigned int iea = 0; iea < nbins; iea++) {
397 double prev = m_DBOneDCell->getMean(8 * il + 1, iea);
398 m_onedcors[il][iea] *= prev;
399 }
400 }
401
402 double binsize = TMath::Pi() / m_onedcors[il].size();
403
404 auto computeAverages = [&](double angLow, double angHigh, const std::string & label) {
405 unsigned int binLow = std::floor((angLow + TMath::Pi() / 2.0) / binsize);
406 unsigned int binHigh = std::floor((angHigh + TMath::Pi() / 2.0) / binsize);
407 double sum_new = 0.0, sum_prev = 0.0;
408 int count = 0;
409
410 for (unsigned int iea = binLow; iea < binHigh; ++iea) {
411 sum_new += m_onedcors[il][iea];
412 sum_prev += m_DBOneDCell->getMean(8 * il + 1, iea);
413 ++count;
414 }
415
416 double avg_new = (count > 0) ? sum_new / count : 1.0;
417 double avg_prev = (count > 0) ? sum_prev / count : 1.0;
418 return std::make_pair(avg_new, avg_prev);
419 };
420
421 double negLow = -0.75, negHigh = -0.25;
422 double posLow = 0.25, posHigh = 0.75;
423
424 if (il == 0) {
425 negLow = -0.5; negHigh = -0.2;
426 posLow = 0.2; posHigh = 0.5;
427 }
428
429 auto [avgNewNeg, avgPrevNeg] = computeAverages(negLow, negHigh, "Negative");
430 auto [avgNewPos, avgPrevPos] = computeAverages(posLow, posHigh, "Positive");
431
432 double avgNew = (avgNewNeg + avgNewPos) / 2.0;
433 double avgPrev = (avgPrevNeg + avgPrevPos) / 2.0;
434 double scaleFactor = avgPrev / avgNew;
435
436 for (unsigned int iea = 0; iea < m_eaBin; iea++) {
437 m_onedcors[il][iea] *= scaleFactor;
438 }
439
440 if (m_chargeType > 0)
441 for (int ie = 0; ie < m_eaBin / 2; ie++) m_onedcors[il][ie] *= m_adjustFac;
442 if (m_chargeType < 0)
443 for (int ie = m_eaBin / 2; ie < m_eaBin; ie++) m_onedcors[il][ie] *= m_adjustFac;
444
445 }
446 //Saving constants
447 B2INFO("dE/dx Calibration done for CDCDedx1DCell");
448 CDCDedx1DCell* gain = new CDCDedx1DCell(0, m_onedcors);
449 saveCalibration(gain, "CDCDedx1DCell");
450}
451
452//--------------------------------------------------
453void CDCDedx1DCellAlgorithm::plotMergeFactor(std::map<int, std::vector<double>> bounds, const std::array<int, 2> nDev,
454 std::map<int, std::vector<int>> steps)
455{
456
457 TCanvas cmfactor("cmfactor", "Merging factors", 800, 400);
458 cmfactor.Divide(2, 1);
459
460 for (int il = 0; il < 2; il++) {
461 Double_t* nvarBins;
462 nvarBins = &bounds[il][0];
463
464 TH1I* hist = new TH1I(Form("hist_%s", m_label[il].data()), "", nDev[il], nvarBins);
465 hist->SetTitle(Form("Merging factor for %s bins;binindex;merge-factors", m_label[il].data()));
466
467 for (int ibin = 0; ibin < nDev[il]; ibin++) hist->SetBinContent(ibin + 1, steps[il][ibin]);
468
469 cmfactor.cd(il + 1);
470 hist->SetFillColor(kYellow);
471 hist->Draw("hist");
472 delete hist;
473 }
474
475 cmfactor.SaveAs(Form("cdcdedx_1dcell_mergefactor%s.pdf", m_suffix.data()));
476 cmfactor.SaveAs(Form("cdcdedx_1dcell_mergefactor%s.root", m_suffix.data()));
477}
478
479//--------------------------------------------------
480void CDCDedx1DCellAlgorithm::plotdedxHist(std::vector<TH1D*> hdedxhit[2])
481{
482
483 TCanvas ctmp("tmp", "tmp", 1200, 1200);
484 ctmp.Divide(4, 4);
485 std::stringstream psname;
486
487 psname << Form("cdcdedx_1dcell_dedxhit%s.pdf[", m_suffix.data());
488 ctmp.Print(psname.str().c_str());
489 psname.str("");
490 psname << Form("cdcdedx_1dcell_dedxhit%s.pdf", m_suffix.data());
491
492 for (int il = 0; il < 2; il++) {
493
494 for (int jea = 0; jea < m_eaBinLocal[il]; jea++) {
495
496 int minbin = std::stoi(hdedxhit[il][jea]->GetXaxis()->GetTitle());
497 int maxbin = std::stoi(hdedxhit[il][jea]->GetYaxis()->GetTitle());
498
499 ctmp.cd(jea % 16 + 1);
500 hdedxhit[il][jea]->SetFillColor(4 + il);
501
502 hdedxhit[il][jea]->SetTitle(Form("%s;dedxhit;entries", hdedxhit[il][jea]->GetTitle()));
503 hdedxhit[il][jea]->DrawClone("hist");
504 TH1D* htempC = (TH1D*)hdedxhit[il][jea]->Clone(Form("%sc2", hdedxhit[il][jea]->GetName()));
505 htempC->GetXaxis()->SetRange(minbin, maxbin);
506 htempC->SetFillColor(kGray);
507 htempC->DrawClone("same hist");
508
509 if (jea % 16 == 15 || (jea == m_eaBinLocal[il] - 1)) {
510 ctmp.Print(psname.str().c_str());
511 gPad->Clear("D");
512 ctmp.Clear("D");
513 }
514 delete htempC;
515 }
516 }
517 psname.str("");
518 psname << Form("cdcdedx_1dcell_dedxhit%s.pdf]", m_suffix.data());
519 ctmp.Print(psname.str().c_str());
520}
521
522//--------------------------------------------------
524{
525
526 TCanvas cdedxlayer("layerdedxhit", "Inner and Outer Layer dedxhit dist", 900, 400);
527 cdedxlayer.Divide(2, 1);
528
529 for (int il = 0; il < 2; il++) {
530 int minlay = 0, maxlay = 0;
531 if (isFixTrunc) {
532 minlay = std::stoi(hdedxlay[il]->GetXaxis()->GetTitle());
533 maxlay = std::stoi(hdedxlay[il]->GetYaxis()->GetTitle());
534 double lowedge = hdedxlay[il]->GetXaxis()->GetBinLowEdge(minlay);
535 double upedge = hdedxlay[il]->GetXaxis()->GetBinUpEdge(maxlay);
536 hdedxlay[il]->SetTitle(Form("%s, trunc #rightarrow: %0.02f - %0.02f;dedxhit;entries", hdedxlay[il]->GetTitle(), lowedge, upedge));
537 }
538
539 cdedxlayer.cd(il + 1);
540 hdedxlay[il]->SetFillColor(kYellow);
541 hdedxlay[il]->Draw("histo");
542
543 if (isFixTrunc) {
544 TH1D* hdedxlayC = (TH1D*)hdedxlay[il]->Clone(Form("hdedxlayC%d", il));
545 hdedxlayC->GetXaxis()->SetRange(minlay, maxlay);
546 hdedxlayC->SetFillColor(kAzure + 1);
547 hdedxlayC->Draw("same histo");
548 }
549 }
550
551 cdedxlayer.SaveAs(Form("cdcdedx_1dcell_dedxlay%s.pdf", m_suffix.data()));
552 cdedxlayer.SaveAs(Form("cdcdedx_1dcell_dedxlay%s.root", m_suffix.data()));
553}
554
555//--------------------------------------------------
556void CDCDedx1DCellAlgorithm::plotQaPars(TH1D* hentalay[2], TH2D* hptcosth)
557{
558
559 TCanvas ceadist("ceadist", "Enta distributions", 800, 400);
560 ceadist.Divide(2, 1);
561
562 for (int il = 0; il < 2; il++) {
563
564 ceadist.cd(il + 1);
565 gPad->SetLogy();
566 hentalay[il]->SetFillColor(kYellow);
567 hentalay[il]->Draw("hist");
568 }
569
570 TCanvas cptcos("cptcos", "pt vs costh dist.", 400, 400);
571 cptcos.cd();
572 hptcosth->Draw("colz");
573
574 cptcos.SaveAs(Form("cdcdedx_ptcosth_%s.pdf", m_suffix.data()));
575 ceadist.SaveAs(Form("cdcdedx_1dcell_enta%s.pdf", m_suffix.data()));
576 ceadist.SaveAs(Form("cdcdedx_1dcell_enta%s.root", m_suffix.data()));
577}
578
579//--------------------------------------------------
580void CDCDedx1DCellAlgorithm::plotRelConst(std::vector<double>tempconst, std::vector<double>layerconst, int il)
581{
582
583 TH1D* hconst, *hconstvar;
584
585 Double_t* nvarBins;
586 nvarBins = &m_binValue[il][0];
587
588 hconst = new TH1D(Form("hconst%s", m_label[il].data()), "", m_eaBin, m_eaMin, m_eaMax);
589 std::string title = Form("calibration const dist: %s: (%s); entaRS (#alpha); entries", m_label[il].data(), m_runExp.data());
590 hconst->SetTitle(Form("%s", title.data()));
591
592 hconstvar = new TH1D(Form("hconstvar%s", m_label[il].data()), "", m_eaBinLocal[il], nvarBins);
593 title = Form("calibration const dist (var bins): %s: (%s); entaRS (#alpha);entries", m_label[il].data(), m_runExp.data());
594 hconstvar->SetTitle(Form("%s", title.data()));
595
596 if (isVarBins) {
597 for (int iea = 0; iea < m_eaBinLocal[il]; iea++)
598 hconstvar->SetBinContent(iea + 1, tempconst.at(iea));
599 }
600
601 for (int jea = 0; jea < m_eaBin; jea++) hconst->SetBinContent(jea + 1, layerconst.at(jea));
602
603 gStyle->SetOptStat("ne");
604 TCanvas cconst("cconst", "Calirbation Constants", 800, 400);
605 if (isVarBins) {
606 cconst.Divide(2, 1);
607 cconst.SetWindowSize(1000, 800);
608 }
609
610 cconst.cd(1);
611 hconst->SetFillColor(kYellow);
612 hconst->Draw("histo");
613 if (isVarBins) {
614 cconst.cd(2);
615 hconstvar->SetFillColor(kBlue);
616 hconstvar->Draw("hist");
617 }
618 cconst.SaveAs(Form("cdcdedx_1dcell_relconst%s_%s.pdf", m_label[il].data(), m_suffix.data()));
619 cconst.SaveAs(Form("cdcdedx_1dcell_relconst%s_%s.root", m_label[il].data(), m_suffix.data()));
620
621 delete hconst;
622 delete hconstvar;
623}
624
625//--------------------------------------------------
627{
628
629 //Draw New/Old final constants
630 TH1D* hnewconst[2], *holdconst[2];
631 double min[2], max[2];
632
633 for (unsigned int il = 0; il < 2; il++) {
634 unsigned int nbins = m_DBOneDCell->getNBins(il);
635
636 std::string title = Form("final calibration const dist (%s): %s; entaRS (#alpha); entries", m_label[il].data(), m_runExp.data());
637 hnewconst[il] = new TH1D(Form("hnewconst_%s", m_label[il].data()), "", m_eaBin, m_eaMin, m_eaMax);
638 hnewconst[il]->SetTitle(Form("%s", title.data()));
639
640 title = Form("old calibration const dist (%s): %s; entaRS (#alpha); entries", m_label[il].data(), m_runExp.data());
641 holdconst[il] = new TH1D(Form("holdconst_%s", m_label[il].data()), "", m_eaBin, m_eaMin, m_eaMax);
642 holdconst[il]->SetTitle(Form("%s", title.data()));
643
644 for (unsigned int iea = 0; iea < nbins; iea++) {
645 double prev = m_DBOneDCell->getMean(8 * il + 1, iea);
646 holdconst[il]->SetBinContent(iea + 1, prev);
647 hnewconst[il]->SetBinContent(iea + 1, m_onedcors[il][iea]);
648 }
649 min[il] = hnewconst[il]->GetMinimum();
650 max[il] = hnewconst[il]->GetMaximum();
651 }
652
653 //Plotting final constants
654 if (max[1] < max[0])max[1] = max[0];
655 if (min[1] > min[0])min[1] = min[0];
656
657 gStyle->SetOptStat("ne");
658 TCanvas cfconst("cfconst", "Final calirbation constants", 800, 400);
659 cfconst.Divide(2, 1);
660
661 for (int il = 0; il < 2; il++) {
662 cfconst.cd(il + 1);
663 hnewconst[il]->GetYaxis()->SetRangeUser(min[1] * 0.95, max[1] * 1.05);
664 hnewconst[il]->SetLineColor(kBlack);
665 hnewconst[il]->Draw("histo");
666 holdconst[il]->SetLineColor(kRed);
667 holdconst[il]->Draw("histo same");
668
669 auto legend = new TLegend(0.4, 0.75, 0.56, 0.85);
670 legend->AddEntry(holdconst[il], "Old", "lep");
671 legend->AddEntry(hnewconst[il], "New", "lep");
672 legend->Draw();
673 }
674
675 cfconst.SaveAs(Form("cdcdedx_1dcell_fconsts%s.pdf", m_suffix.data()));
676 cfconst.SaveAs(Form("cdcdedx_1dcell_fconsts%s.root", m_suffix.data()));
677
678 for (int il = 0; il < 2; il++) {
679 delete hnewconst[il];
680 delete holdconst[il];
681 }
682}
683
684//------------------------------------
686{
687
688 TCanvas cstats("cstats", "cstats", 1000, 500);
689 cstats.SetBatch(kTRUE);
690 cstats.Divide(2, 1);
691
692 cstats.cd(1);
693 auto hestats = getObjectPtr<TH1I>("hestats");
694 if (hestats) {
695 hestats->SetName(Form("hestats_%s", m_suffix.data()));
696 hestats->SetStats(0);
697 hestats->DrawCopy("");
698 }
699
700 cstats.cd(2);
701 auto htstats = getObjectPtr<TH1I>("htstats");
702 if (htstats) {
703 hestats->SetName(Form("htstats_%s", m_suffix.data()));
704 htstats->DrawCopy("");
705 hestats->SetStats(0);
706 }
707 cstats.Print(Form("cdcdedx_1dcell_stats_%s.pdf", m_suffix.data()));
708}
CDCDedx1DCellAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
std::string m_label[2]
add inner/outer layer label
double m_eaMax
upper edge of entrance angle
void plotMergeFactor(std::map< int, std::vector< double > > bounds, const std::array< int, 2 > nDev, std::map< int, std::vector< int > > steps)
function to plot merging factor
double m_truncMax
upper threshold on truncation
int m_binSplit
multiply nbins by this factor in full range
std::array< std::vector< int >, 2 > m_binIndex
symm/Var bin numbers
double m_truncMin
lower threshold on truncation
double m_adjustFac
factor with that one what to adjust baseline
void getTruncatedBins(TH1D *hist, int &binlow, int &binhigh)
function to get bins of truncation from histogram
void CreateBinMapping()
class function to create vectors for bin mapping (Var->symm)
double m_chargeType
charge type for baseline adj
void getExpRunInfo()
function to get extract calibration run/exp
DBObjPtr< CDCDedx1DCell > m_DBOneDCell
One cell correction DB object.
double m_cosMax
a limit on cos theta
bool isPrintLog
print more debug information
std::array< std::vector< double >, 2 > m_binValue
enta Var bin values
std::string m_suffix
add suffix to all plot name
int m_eaB
reset # of bins for entrance angle for each experiment
double getTruncationMean(TH1D *hist, int binlow, int binhigh)
function to get truncated mean
double m_ptMax
a limit on transverse momentum
void plotConstants()
function to draw the old/new final constants
bool isFixTrunc
true = fix window for all out/inner layers
bool isVarBins
true: if variable bin size is requested
void plotdedxHist(std::vector< TH1D * > hdedxhit[2])
function to draw the dE/dx histogram in enta bins
void defineHisto(std::vector< TH1D * > hdedxhit[2], TH1D *hdedxlay[2], TH1D *hentalay[2])
function to define histograms
double m_eaBW
binwdith of entrance angle bin
bool isRotSymm
if rotation symmetry requested
std::string m_runExp
add suffix to all plot name
void plotEventStats()
function to draw the stats plots
int rotationalBin(int nbin, int ibin)
class function to set rotation symmetry
virtual EResult calibrate() override
1D cell algorithm
void plotQaPars(TH1D *hentalay[2], TH2D *hptcosth)
function to draw pt vs costh and entrance angle distribution for Inner/Outer layer
double m_dedxMax
upper edge of dedxhit
void createPayload()
function to generate final constants
bool isMakePlots
produce plots for status
void plotRelConst(std::vector< double >tempconst, std::vector< double >layerconst, int il)
function to draw symm/Var layer constant
std::vector< std::vector< double > > m_onedcors
final vectors of calibration
bool isMerge
print more debug information
double m_dedxMin
lower edge of dedxhit
void plotLayerDist(TH1D *hdedxL[2])
function to draw dedx dist.
double m_eaMin
lower edge of entrance angle
dE/dx wire gain calibration constants
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(....
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.