Belle II Software development
CDCDedxCosEdgeAlgorithm.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 <TF1.h>
10#include <TCanvas.h>
11#include <TH1I.h>
12
13#include <cdc/calibration/CDCdEdx/CDCDedxCosEdgeAlgorithm.h>
14
15using namespace Belle2;
16using namespace std;
17
18//-----------------------------------------------------------------
19// Implementation
20//-----------------------------------------------------------------
22 CalibrationAlgorithm("CDCDedxElectronCollector"),
23 m_isMakePlots(true),
24 m_isMerge(true),
25 m_sigLim(2.5),
26 m_npBins(20),
27 m_negMin(-0.870),
28 m_negMax(-0.850),
29 m_posMin(0.950),
30 m_posMax(0.960),
31 m_dedxBins(600),
32 m_dedxMin(0.0),
33 m_dedxMax(3.0),
34 m_suffix("")
35{
36 // Set module properties
37 setDescription("A calibration algorithm for CDC dE/dx electron cos(theta) dependence");
38}
39
40//-----------------------------------------------------------------
41// Run the calibration
42//-----------------------------------------------------------------
44{
45
47
48 // Get data objects
49 auto ttree = getObjectPtr<TTree>("tree");
50 if (ttree->GetEntries() < 100) return c_NotEnoughData;
51
52 double dedx, costh; int charge;
53 ttree->SetBranchAddress("dedx", &dedx);
54 ttree->SetBranchAddress("costh", &costh);
55 ttree->SetBranchAddress("charge", &charge);
56
57 // make histograms to store dE/dx values in bins of cos(theta)
58 vector<TH1D*> hdedx_negi, hdedx_posi;
59
60 const double bwnegi = (m_negMax - m_negMin) / m_npBins;
61 const double bwposi = (m_posMax - m_posMin) / m_npBins;
62
63 for (unsigned int i = 0; i < m_npBins; ++i) {
64
65 double mincos = i * bwposi + m_posMin;
66 double maxcos = mincos + bwposi;
67 string title = Form("costh: %0.04f, %0.04f(%s)", mincos, maxcos, m_suffix.data());
68 hdedx_posi.push_back(new TH1D(Form("hdedx_posi%d_%s", i, m_suffix.data()), "", m_dedxBins, m_dedxMin, m_dedxMax));
69 hdedx_posi[i]->SetTitle(Form("%s;dedx;entries", title.data()));
70
71 mincos = i * bwnegi + m_negMin;
72 maxcos = mincos + bwnegi;
73 title = Form("costh: %0.04f, %0.04f(%s)", mincos, maxcos, m_suffix.data());
74 hdedx_negi.push_back(new TH1D(Form("hdedx_negi%d_%s", i, m_suffix.data()), "", m_dedxBins, m_dedxMin, m_dedxMax));
75 hdedx_negi[i]->SetTitle(Form("%s;dedx;entries", title.data()));
76 }
77
78 int icosbin = -99.;
79 for (int i = 0; i < ttree->GetEntries(); ++i) {
80
81 ttree->GetEvent(i);
82
83 //if track is a junk
84 if (dedx <= 0 || charge == 0) continue;
85 if (costh > -0.850 && costh < 0.950) continue;
86
87 if (costh > 0) {
88 icosbin = int((costh - m_posMin) / bwposi);
89 hdedx_posi[icosbin]->Fill(dedx);
90 } else if (costh < 0) {
91 icosbin = int((costh - m_negMin) / bwnegi);
92 hdedx_negi[icosbin]->Fill(dedx);
93 }
94 }
95
96 map<int, vector<double>> vneg_fitpars;
97 map<int, vector<double>> vpos_fitpars;
98
99 vector<double> vneg_const, vpos_const;
100 vector<vector<double>> vfinal_const;
101
102 for (unsigned int i = 0; i < m_npBins; ++i) {
103
104 fitstatus status;
105
106 //Fit dedx in negative cos bins
107 fitGaussianWRange(hdedx_negi[i], status);
108 if (status != FitOK) {
109 vneg_fitpars[0].push_back(1.0);
110 vneg_fitpars[1].push_back(0.0);
111 vneg_fitpars[2].push_back(0.0);
112 vneg_fitpars[3].push_back(0.0);
113 hdedx_negi[i]->SetTitle(Form("%s, Fit(%d)", hdedx_negi[i]->GetTitle(), status));
114 } else {
115 vneg_fitpars[0].push_back(hdedx_negi[i]->GetFunction("gaus")->GetParameter(1));
116 vneg_fitpars[1].push_back(hdedx_negi[i]->GetFunction("gaus")->GetParError(1));
117 vneg_fitpars[2].push_back(hdedx_negi[i]->GetFunction("gaus")->GetParameter(2));
118 vneg_fitpars[3].push_back(hdedx_negi[i]->GetFunction("gaus")->GetParError(2));
119 }
120
121 vneg_const.push_back(vneg_fitpars[0][i]);
122
123 //Fit dedx in positive cos bins
124 fitGaussianWRange(hdedx_posi[i], status);
125 if (status != FitOK) {
126 vpos_fitpars[0].push_back(1.0);
127 vpos_fitpars[1].push_back(0.0);
128 vpos_fitpars[2].push_back(0.0);
129 vpos_fitpars[3].push_back(0.0);
130 hdedx_posi[i]->SetTitle(Form("%s, Fit(%d)", hdedx_posi[i]->GetTitle(), status));
131 } else {
132 vpos_fitpars[0].push_back(hdedx_posi[i]->GetFunction("gaus")->GetParameter(1));
133 vpos_fitpars[1].push_back(hdedx_posi[i]->GetFunction("gaus")->GetParError(1));
134 vpos_fitpars[2].push_back(hdedx_posi[i]->GetFunction("gaus")->GetParameter(2));
135 vpos_fitpars[3].push_back(hdedx_posi[i]->GetFunction("gaus")->GetParError(2));
136 }
137
138 vpos_const.push_back(vpos_fitpars[0][i]);
139 }
140
141 vfinal_const.push_back(vneg_const);
142 vfinal_const.push_back(vpos_const);
143
144 createPayload(vfinal_const);
145
146 if (m_isMakePlots) {
147
148 //1. draw dedx dist of individual bins
149 plotHist(hdedx_posi, vpos_fitpars, "pos");
150 plotHist(hdedx_negi, vneg_fitpars, "neg");
151
152 //2. plot relative const or fit parameters
153 plotFitPars(vneg_fitpars, vpos_fitpars);
154
155 //3. compare new/old calibration constants
156 plotConstants(vfinal_const);
157
158 //4. plot statistics related plots here
159 plotStats();
160 }
161
162 m_suffix.clear();
163 return c_OK;
164}
165
166
167//------------------------------------
169{
170 int cruns = 0;
171 for (auto expRun : getRunList()) {
172 if (cruns == 0)B2INFO("CDCDedxCosEdge: start exp " << expRun.first << " and run " << expRun.second << "");
173 cruns++;
174 }
175
176 const auto erStart = getRunList()[0];
177 int estart = erStart.first;
178 int rstart = erStart.second;
179
180 const auto erEnd = getRunList()[cruns - 1];
181 int rend = erEnd.second;
182
183 updateDBObjPtrs(1, rstart, estart);
184
185 if (m_suffix.length() > 0) m_suffix = Form("%s_e%d_r%dr%d", m_suffix.data(), estart, rstart, rend);
186 else m_suffix = Form("e%d_r%dr%d", estart, rstart, rend);
187}
188
189//------------------------------------
190void CDCDedxCosEdgeAlgorithm::createPayload(vector<vector<double>>& vfinalconst)
191{
192 if (m_isMerge) {
193 if (m_DBCosineCor->getSize(-1) != int(m_npBins) || m_DBCosineCor->getSize(1) != int(m_npBins))
194 B2FATAL("CDCDedxCosEdgeAlgorithm: Can't merge paylaods with different size");
195
196 for (unsigned int ibin = 0; ibin < m_npBins; ibin++) {
197
198 //costh < 0
199 double prevg = m_DBCosineCor->getMean(-1, ibin);
200 double relg = vfinalconst[0].at(ibin);
201 double newg = prevg * relg;
202 B2INFO("CosEdge Const (<0), bin# " << ibin << ", rel " << relg << ", previous " << prevg << ", merged " << newg);
203 vfinalconst[0].at(ibin) *= (double)m_DBCosineCor->getMean(-1, ibin);
204 vfinalconst[0].at(ibin) /= (0.5 * (vfinalconst[0].at(m_npBins - 1) + vfinalconst[0].at(m_npBins - 2)));
205
206 //costh > 0
207 prevg = m_DBCosineCor->getMean(1, ibin);
208 relg = vfinalconst[1].at(ibin);
209 newg = prevg * relg;
210 B2INFO("CosEdge Const (>0), bin# " << ibin << ", rel " << relg << ", previous " << prevg << ", merged " << newg);
211 vfinalconst[1].at(ibin) *= (double)m_DBCosineCor->getMean(1, ibin);
212 vfinalconst[1].at(ibin) /= (0.5 * (vfinalconst[1].at(0) + vfinalconst[1].at(1)));
213 }
214 }
215
216 B2INFO("CDCDedxCosineEdge calibration done");
217 CDCDedxCosineEdge* gain = new CDCDedxCosineEdge(vfinalconst);
218 saveCalibration(gain, "CDCDedxCosineEdge");
219}
220
221//------------------------------------
222void CDCDedxCosEdgeAlgorithm::fitGaussianWRange(TH1D*& temphist, fitstatus& status)
223{
224 if (temphist->Integral() < 500) { //at least 1k bhabha events
225 B2INFO(Form("\t insufficient fit stats (%0.00f) for (%s)", temphist->Integral(), temphist->GetName()));
226 status = LowStats;
227 return;
228 } else {
229 temphist->GetXaxis()->SetRange(temphist->FindFirstBinAbove(0, 1), temphist->FindLastBinAbove(0, 1));
230 int fs = temphist->Fit("gaus", "QR");
231 if (fs != 0) {
232 B2INFO(Form("\tFit (round 1) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
233 status = FitFailed;
234 return;
235 } else {
236 double fdEdxMean = temphist->GetFunction("gaus")->GetParameter(1);
237 double width = temphist->GetFunction("gaus")->GetParameter(2);
238 temphist->GetXaxis()->SetRangeUser(fdEdxMean - 5.0 * width, fdEdxMean + 5.0 * width);
239 fs = temphist->Fit("gaus", "QR", "", fdEdxMean - m_sigLim * width, fdEdxMean + m_sigLim * width);
240 if (fs != 0) {
241 B2INFO(Form("\tFit (round 2) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
242 status = FitFailed;
243 return;
244 } else {
245 temphist->GetXaxis()->SetRangeUser(fdEdxMean - 5.0 * width, fdEdxMean + 5.0 * width);
246 B2INFO(Form("\tFit for hist (%s) successful (status = %d)", temphist->GetName(), fs));
247 status = FitOK;
248 }
249 }
250 }
251}
252
253//------------------------------------
254void CDCDedxCosEdgeAlgorithm::plotHist(vector<TH1D*>& hdedx, map<int, vector<double>>& vpars,
255 string type)
256{
257 TCanvas ctmp("ctmp", "ctmp", 1200, 1200);
258 ctmp.Divide(5, 4);
259
260 stringstream psname;
261 psname << Form("cdcdedx_cosedgecal_fits_%s_%s.pdf[", type.data(), m_suffix.data());
262 ctmp.Print(psname.str().c_str());
263 psname.str("");
264 psname << Form("cdcdedx_cosedgecal_fits_%s_%s.pdf", type.data(), m_suffix.data());
265
266 for (unsigned int i = 0; i < m_npBins; ++i) {
267
268 ctmp.cd(i % 20 + 1); // each canvas is 2x2
269 hdedx[i]->SetStats(0);
270 hdedx[i]->SetFillColor(kAzure - 9);
271 hdedx[i]->DrawCopy();
272
273 TPaveText* pt = new TPaveText(0.5, 0.73, 0.8, 0.89, "NBNDC");
274 setTextCosmetics(pt, 0.04258064);
275 pt->AddText(Form("#mu_{fit}: %0.03f#pm%0.03f", vpars[0][i], vpars[1][i]));
276 pt->AddText(Form("#sigma_{fit}: %0.03f#pm%0.03f", vpars[2][i], vpars[3][i]));
277 pt->Draw("same");
278
279 if ((i + 1) % 20 == 0 || (i + 1) == m_npBins) {
280 ctmp.Print(psname.str().c_str());
281 ctmp.Clear("D");
282 }
283
284 delete hdedx[i];
285 }
286
287 psname.str("");
288 psname << Form("cdcdedx_cosedgecal_fits_%s_%s.pdf]", type.data(), m_suffix.data());
289 ctmp.Print(psname.str().c_str());
290}
291
292//------------------------------------
293void CDCDedxCosEdgeAlgorithm::plotFitPars(map<int, vector<double>>& vneg_fitpars,
294 map<int, vector<double>>& vpos_fitpars)
295{
296 // For qa pars
297 TCanvas cQa("cQa", "cQa", 1200, 1200);
298 cQa.Divide(2, 2);
299
300 double min[2] = {0.85, 0.04};
301 double max[2] = {1.05, 0.3};
302
303 string vars[2] = {"#mu_{fit}", "#sigma_{fit}"};
304 string side[2] = {"pcos", "ncos"};
305
306 for (int is = 0; is < 2; is++) {
307
308 double minp = m_negMin, maxp = m_negMax;
309 if (is == 0) {
310 minp = m_posMin;
311 maxp = m_posMax;
312 }
313
314 for (int iv = 0; iv < 2; iv++) {
315
316 string hname = Form("hpar_%s_%s_%s", vars[iv].data(), side[is].data(), m_suffix.data());
317
318 TH1D htemp(Form("%s", hname.data()), "", m_npBins, minp, maxp);
319 htemp.SetTitle(Form("Constant (%s), cos#theta:(%0.02f, %0.02f);cos(#theta);const", vars[iv].data(), minp, maxp));
320
321 for (unsigned int ib = 0; ib < m_npBins; ib++) {
322 if (is == 0) {
323 htemp.SetBinContent(ib + 1, vpos_fitpars[2 * iv][ib]);
324 htemp.SetBinError(ib + 1, vpos_fitpars[2 * iv + 1][ib]);
325 } else {
326 htemp.SetBinContent(ib + 1, vneg_fitpars[2 * iv][ib]);
327 htemp.SetBinError(ib + 1, vneg_fitpars[2 * iv + 1][ib]);
328 }
329 }
330
331 setHistCosmetics(htemp, iv * 2 + 2, min[iv], max[iv], 0.025);
332
333 cQa.cd(2 * iv + 1 + is); //1,3,2,4
334 gPad->SetGridx(1);
335 gPad->SetGridy(1);
336 htemp.DrawCopy("");
337 }
338 }
339 cQa.SaveAs(Form("cdcdedx_cosedgecal_relconst_%s.pdf", m_suffix.data()));
340 cQa.SaveAs(Form("cdcdedx_cosedgecal_relconst_%s.root", m_suffix.data()));
341}
342
343//------------------------------------
344void CDCDedxCosEdgeAlgorithm::plotConstants(vector<vector<double>>& vfinalconst)
345{
346
347 TCanvas ctmp_const("ctmp_const", "ctmp_const", 900, 450);
348 ctmp_const.Divide(2, 1);
349
350 for (int i = 0; i < 2; i++) {
351
352 double min = m_negMin, max = m_negMax;
353 if (i == 1) {
354 min = m_posMin;
355 max = m_posMax;
356 }
357
358 TH1D holdconst(Form("holdconst%d_%s", i, m_suffix.data()), "", m_npBins, min, max);
359 holdconst.SetTitle(Form("constant comparison, cos#theta:(%0.02f, %0.02f);cos(#theta);const", min, max));
360
361 TH1D hnewconst(Form("hnewconst%d_%s", i, m_suffix.data()), "", m_npBins, min, max);
362
363 int iside = 2 * i - 1; //-1 or +1 for neg and pos cosine side
364 for (int ibin = 0; ibin < m_DBCosineCor->getSize(iside); ibin++) {
365 holdconst.SetBinContent(ibin + 1, (double)m_DBCosineCor->getMean(iside, ibin));
366 hnewconst.SetBinContent(ibin + 1, vfinalconst[i].at(ibin));
367 }
368
369 ctmp_const.cd(i + 1);
370 gPad->SetGridx(1);
371 gPad->SetGridy(1);
372
373 setHistCosmetics(holdconst, kBlack, 0.0, 1.10, 0.025);
374 holdconst.DrawCopy("");
375
376 setHistCosmetics(hnewconst, kRed, 0.0, 1.10, 0.025);
377 hnewconst.DrawCopy("same");
378
379 TPaveText* pt = new TPaveText(0.47, 0.73, 0.77, 0.89, "NBNDC");
380 setTextCosmetics(pt, 0.02258064);
381
382 TText* told = pt->AddText("old const");
383 told->SetTextColor(kBlack);
384 TText* tnew = pt->AddText("new const");
385 tnew->SetTextColor(kRed);
386 pt->Draw("same");
387
388 }
389
390 ctmp_const.SaveAs(Form("cdcdedx_cosedgecal_constcomp_%s.pdf", m_suffix.data()));
391 ctmp_const.SaveAs(Form("cdcdedx_cosedgecal_constcomp_%s.root", m_suffix.data()));
392}
393
394//------------------------------------
396{
397
398 TCanvas cstats("cstats", "cstats", 1000, 500);
399 cstats.SetBatch(kTRUE);
400 cstats.Divide(2, 1);
401
402 cstats.cd(1);
403 auto hestats = getObjectPtr<TH1I>("hestats");
404 if (hestats) {
405 hestats->SetName(Form("hestats_%s", m_suffix.data()));
406 hestats->SetStats(0);
407 hestats->DrawCopy("");
408 }
409
410 cstats.cd(2);
411 auto htstats = getObjectPtr<TH1I>("htstats");
412 if (htstats) {
413 htstats->SetName(Form("htstats_%s", m_suffix.data()));
414 htstats->SetStats(0);
415 htstats->DrawCopy("");
416 }
417 cstats.Print(Form("cdcdedx_cosedgecal_stats_%s.pdf", m_suffix.data()));
418}
void plotFitPars(std::map< int, std::vector< double > > &fPars_Neg, std::map< int, std::vector< double > > &fPars_Pos)
function to draw the fit parameters (relative gains and resolutions)
CDCDedxCosEdgeAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
bool m_isMerge
merge payload if calculated relative
void plotStats()
function to draw the stats plots
void setHistCosmetics(TH1D &hist, Color_t color, double min, double max, double size)
function to change histogram styles
void plotHist(std::vector< TH1D * > &hdedx, std::map< int, std::vector< double > > &fPars, std::string type)
function to draw dedx histograms for each bin
void getExpRunInfo()
function to get info about current exp and run
DBObjPtr< CDCDedxCosineEdge > m_DBCosineCor
CoseEdge correction DB object.
unsigned int m_npBins
number of bins across cosine range
void createPayload(std::vector< std::vector< double > > &vfinalconst)
function to store new payload after full calibration
std::string m_suffix
suffix for better plot naming
void setTextCosmetics(TPaveText *pt, double size)
function to change text styles
int m_dedxBins
number of bins for dedx histogram
void plotConstants(std::vector< std::vector< double > > &vfinalconst)
function to draw the final calibration constants and comparison with old constants
virtual EResult calibrate() override
Cosine edge algorithm.
double m_sigLim
gaussian fit sigma limit
void fitGaussianWRange(TH1D *&temphist, fitstatus &status)
function to perform gauss fit for given histogram
dE/dx special large cosine calibration to fix bending shoulder at large costh
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.
STL namespace.