Belle II Software development
CDCDedxCosineAlgorithm.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 <reconstruction/calibration/CDCDedxCosineAlgorithm.h>
10
11#include <TF1.h>
12#include <TLine.h>
13#include <TCanvas.h>
14#include <TH1I.h>
15#include <vector>
16
17using namespace Belle2;
18//-----------------------------------------------------------------
19// Implementation
20//-----------------------------------------------------------------
21CDCDedxCosineAlgorithm::CDCDedxCosineAlgorithm() :
22 CalibrationAlgorithm("CDCDedxElectronCollector"),
23 isMethodSep(true),
24 isMakePlots(true),
25 isMergePayload(true),
26 fSigLim(2.5),
27 fCosbins(100),
28 fCosMin(-1.0),
29 fCosMax(1.0),
30 fHistbins(600),
31 fdEdxMin(0.0),
32 fdEdxMax(3.0),
33 fStartRun(0)
34
35{
36 // Set module properties
37 setDescription("A calibration algorithm for CDC dE/dx electron cos(theta) dependence");
38
39}
40
41//-----------------------------------------------------------------
42// Run the calibration
43//-----------------------------------------------------------------
44CalibrationAlgorithm::EResult CDCDedxCosineAlgorithm::calibrate()
45{
46 B2INFO("Preparing dE/dx calibration for CDC dE/dx electron saturation");
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 const auto expRun = getRunList()[0];
58 updateDBObjPtrs(1, expRun.second, expRun.first);
59 fStartRun = expRun.second;
60
61 // make histograms to store dE/dx values in bins of cos(theta)
62 // bin size can be arbitrary, but for now just make uniform bins
63 std::vector<TH1D*> hdEdx_elCosbin(fCosbins), hdEdx_poCosbin(fCosbins), hdEdx_epCosbin(fCosbins);
64
65 const double binW = (fCosMax - fCosMin) / fCosbins;
66
67 for (unsigned int i = 0; i < fCosbins; ++i) {
68
69 double coslow = i * binW + fCosMin, coshigh = coslow + binW;
70
71 hdEdx_elCosbin[i] = new TH1D(Form("hdEdx_elCosbin%d_fRun%d", i, fStartRun), "", fHistbins, fdEdxMin, fdEdxMax);
72 hdEdx_elCosbin[i]->SetTitle(Form("dE/dx dist (e-) in costh (%0.02f, %0.02f), run start: %d", coslow, coshigh, fStartRun));
73 hdEdx_elCosbin[i]->GetXaxis()->SetTitle("dE/dx (no had sat, for e-)");
74 hdEdx_elCosbin[i]->GetYaxis()->SetTitle("Entries");
75
76 hdEdx_poCosbin[i] = new TH1D(Form("hdEdx_poCosbin%d_fRun%d", i, fStartRun), "", fHistbins, fdEdxMin, fdEdxMax);
77 hdEdx_poCosbin[i]->SetTitle(Form("dE/dx dist (e+) in costh (%0.02f, %0.02f), run start: %d", coslow, coshigh, fStartRun));
78 hdEdx_poCosbin[i]->GetXaxis()->SetTitle("dE/dx (no had sat, for e+)");
79 hdEdx_poCosbin[i]->GetYaxis()->SetTitle("Entries");
80
81 hdEdx_epCosbin[i] = new TH1D(Form("hdEdx_epCosbin%d_fRun%d", i, fStartRun), "", fHistbins, fdEdxMin, fdEdxMax);
82 hdEdx_epCosbin[i]->SetTitle(Form("dE/dx dist (e-,e+) in costh (%0.02f, %0.02f), run start: %d", coslow, coshigh, fStartRun));
83 hdEdx_epCosbin[i]->GetXaxis()->SetTitle("dE/dx (no had sat, for e-,e+)");
84 hdEdx_epCosbin[i]->GetYaxis()->SetTitle("Entries");
85 }
86
87 // fill histograms, bin size may be arbitrary
88 TH1D* hCosth_el = new TH1D(Form("hCosth_el_fRun%d", fStartRun),
89 Form("cos(#theta) dist (e- and e+), start run: %d; cos(#theta); Entries", fStartRun), fCosbins, fCosMin, fCosMax);
90 TH1D* hCosth_po = new TH1D(Form("hCosth_po_fRun%d", fStartRun), Form("cos(#theta) dist (e+), start run: %d; cos(#theta); Entries",
91 fStartRun), fCosbins, fCosMin, fCosMax);
92 TH1D* hCosth_ep = new TH1D(Form("hCosth_ep_fRun%d", fStartRun),
93 Form("cos(#theta) dist (e- and e+), start run: %d; cos(#theta); Entries", fStartRun), fCosbins, fCosMin, fCosMax);
94
95 for (int i = 0; i < ttree->GetEntries(); ++i) {
96
97 ttree->GetEvent(i);
98
99 //if track is a junk
100 if (dedx <= 0 || charge == 0) continue;
101
102 //if track is in CDC acceptance (though it is inbuilt in collector module)
103 if (costh < TMath::Cos(150 * TMath::DegToRad()) || costh > TMath::Cos(17 * TMath::DegToRad())) continue;
104
105 int bin = int((costh - fCosMin) / binW);
106 if (bin < 0 || bin >= int(fCosbins)) continue;
107
108 if (isMethodSep) {
109 if (charge < 0) {
110 hCosth_el->Fill(costh);
111 hdEdx_elCosbin[bin]->Fill(dedx);
112 } else if (charge > 0) {
113 hCosth_po->Fill(costh);
114 hdEdx_poCosbin[bin]->Fill(dedx);
115 }
116 } else {
117 hCosth_ep->Fill(costh);
118 hdEdx_epCosbin[bin]->Fill(dedx);
119 }
120 }
121
122 //Plot constants
123 TH1D* hdEdxMeanvsCos_po = new TH1D(Form("hdEdxMeanvsCos_po_fRun%d", fStartRun),
124 Form("dE/dx(e+) rel means, start run: %d; cos(#theta); dE/dx (#mu_{fit})", fStartRun), fCosbins, fCosMin, fCosMax);
125 TH1D* hdEdxSigmavsCos_po = new TH1D(Form("hdEdxSigmavsCos_po_fRun%d", fStartRun),
126 Form("dE/dx(e+) rel means, start run: %d; cos(#theta); dE/dx (#mu_{fit})", fStartRun), fCosbins, fCosMin, fCosMax);
127
128 TH1D* hdEdxMeanvsCos_el = new TH1D(Form("hdEdxMeanvsCos_el_fRun%d", fStartRun),
129 Form("dE/dx(e-) rel means, start run: %d; cos(#theta); dE/dx (#mu_{fit})", fStartRun), fCosbins, fCosMin, fCosMax);
130 TH1D* hdEdxSigmavsCos_el = new TH1D(Form("hdEdxSigmavsCos_el_fRun%d", fStartRun),
131 Form("dE/dx(e-) rel means, start run: %d; cos(#theta); dE/dx (#mu_{fit})", fStartRun), fCosbins, fCosMin, fCosMax);
132
133 TH1D* hdEdxMeanvsCos_ep = new TH1D(Form("hdEdxMeanvsCos_ep_fRun%d", fStartRun),
134 Form("dE/dx(e+, e-) rel means, start run: %d; cos(#theta); dE/dx (#mu_{fit})", fStartRun), fCosbins, fCosMin, fCosMax);
135 TH1D* hdEdxSigmavsCos_ep = new TH1D(Form("hdEdxSigmavsCos_ep_fRun%d", fStartRun),
136 Form("dE/dx(e+, e-) rel means, start run: %d; cos(#theta); dE/dx (#mu_{fit})", fStartRun), fCosbins, fCosMin, fCosMax);
137
138 // more validation plots
139 TCanvas* ctmp_ep = new TCanvas("ctmp_ep", "ctmp_ep", 800, 400);
140 if (isMethodSep)ctmp_ep->Divide(2, 1);
141 else {
142 ctmp_ep->Divide(2, 2);
143 ctmp_ep->SetCanvasSize(800, 800);
144 }
145 std::stringstream psname_ep;
146
147 //validation plots: individual bin dedx dist and fits
148 if (isMakePlots) {
149 psname_ep << Form("cdcdedx_coscal_fits_frun%d.pdf[", fStartRun);
150 ctmp_ep->Print(psname_ep.str().c_str());
151 psname_ep.str("");
152 psname_ep << Form("cdcdedx_coscal_fits_frun%d.pdf", fStartRun);
153 }
154
155 // fit histograms to get gains in bins of cos(theta)
156 std::vector<double> cosine;
157 for (unsigned int i = 0; i < fCosbins; ++i) {
158
159
160 TLine* tl = new TLine();
161 tl->SetLineColor(kBlack);
162
163 double fdEdxMean = 1.0; //This is what we need for calibration
164 double fdEdxMeanErr = 0.0;
165
166 if (!isMethodSep) {
167
168 TString status = "";
169
170 double fdEdxSigma = 0.0, fdEdxSigmaErr = 0.0;
171 FitGaussianWRange(hdEdx_epCosbin[i], status);
172
173 if (status != "FitOK") {
174 fdEdxMean = 1.0;
175 hdEdx_epCosbin[i]->SetTitle(Form("%s, Fit(%s)", hdEdx_epCosbin[i]->GetTitle(), status.Data()));
176 } else {
177 fdEdxMean = hdEdx_epCosbin[i]->GetFunction("gaus")->GetParameter(1);
178 fdEdxMeanErr = hdEdx_epCosbin[i]->GetFunction("gaus")->GetParError(1);
179 fdEdxSigma = hdEdx_epCosbin[i]->GetFunction("gaus")->GetParameter(2);
180 fdEdxSigmaErr = hdEdx_epCosbin[i]->GetFunction("gaus")->GetParError(2);
181 hdEdx_epCosbin[i]->SetTitle(Form("%s, Fit (%s), #mu_{fit}: %0.04f#pm%0.04f,, #sigma_{fit}: %0.04f", hdEdx_epCosbin[i]->GetTitle(),
182 status.Data(), fdEdxMean, fdEdxMeanErr, fdEdxSigma));
183 }
184
185 hdEdxMeanvsCos_ep->SetBinContent(i + 1, fdEdxMean);
186 hdEdxMeanvsCos_ep->SetBinError(i + 1, fdEdxMeanErr);
187 hdEdxSigmavsCos_ep->SetBinContent(i + 1, fdEdxSigma);
188 hdEdxSigmavsCos_ep->SetBinError(i + 1, fdEdxSigmaErr);
189
190 if (isMakePlots) {
191 ctmp_ep->cd(i % 4 + 1); // each canvas is 2x2
192 hdEdx_epCosbin[i]->SetFillColorAlpha(kYellow, 0.25);
193 hdEdx_epCosbin[i]->DrawCopy("hist");
194
195 tl->SetX1(fdEdxMean); tl->SetX2(fdEdxMean);
196 tl->SetY1(0); tl->SetY2(hdEdx_epCosbin[i]->GetMaximum());
197 tl->DrawClone("same");
198 if ((i + 1) % 4 == 0 || (i + 1 == fCosbins))ctmp_ep->Print(psname_ep.str().c_str());
199 }
200 } else {
201
202 double fdEdxMean_el = 1.0, fdEdxMean_elErr = 0.0;
203 double fdEdxSigma_el = 0.0, fdEdxSigma_elErr = 0.0;
204 double fdEdxMean_po = 1.0, fdEdxMean_poErr = 0.0;
205 double fdEdxSigma_po = 0.0, fdEdxSigma_poErr = 0.0;
206 TString status_el = "", status_po = "";
207
208 //Fit _eltrons in cos bins
209 FitGaussianWRange(hdEdx_elCosbin[i], status_el);
210 if (status_el != "FitOK") {
211 fdEdxMean_el = 1.0;
212 hdEdx_elCosbin[i]->SetTitle(Form("%s, Fit(%s)", hdEdx_elCosbin[i]->GetTitle(), status_el.Data()));
213 } else {
214 fdEdxMean_el = hdEdx_elCosbin[i]->GetFunction("gaus")->GetParameter(1);
215 fdEdxMean_elErr = hdEdx_elCosbin[i]->GetFunction("gaus")->GetParError(1);
216 fdEdxSigma_el = hdEdx_elCosbin[i]->GetFunction("gaus")->GetParameter(2);
217 fdEdxSigma_elErr = hdEdx_elCosbin[i]->GetFunction("gaus")->GetParError(2);
218 hdEdx_elCosbin[i]->SetTitle(Form("%s, Fit (%s), #mu_{fit}: %0.04f#pm%0.04f,, #sigma_{fit}: %0.04f", hdEdx_elCosbin[i]->GetTitle(),
219 status_el.Data(), fdEdxMean_el, fdEdxMean_elErr, fdEdxSigma_el));
220 }
221
222 hdEdxMeanvsCos_el->SetBinContent(i + 1, fdEdxMean_el);
223 hdEdxMeanvsCos_el->SetBinError(i + 1, fdEdxMean_elErr);
224 hdEdxSigmavsCos_el->SetBinContent(i + 1, fdEdxSigma_el);
225 hdEdxSigmavsCos_el->SetBinError(i + 1, fdEdxSigma_elErr);
226
227 //Fit _potron in cos bins
228 FitGaussianWRange(hdEdx_poCosbin[i], status_po);
229 if (status_po != "FitOK") {
230 fdEdxMean_po = 1.0;
231 hdEdx_poCosbin[i]->SetTitle(Form("%s, Fit(%s)", hdEdx_poCosbin[i]->GetTitle(), status_po.Data()));
232 } else {
233 fdEdxMean_po = hdEdx_poCosbin[i]->GetFunction("gaus")->GetParameter(1);
234 fdEdxMean_poErr = hdEdx_poCosbin[i]->GetFunction("gaus")->GetParError(1);
235 fdEdxSigma_po = hdEdx_poCosbin[i]->GetFunction("gaus")->GetParameter(2);
236 fdEdxSigma_poErr = hdEdx_poCosbin[i]->GetFunction("gaus")->GetParError(2);
237 hdEdx_poCosbin[i]->SetTitle(Form("%s, Fit (%s), #mu_{fit}: %0.04f#pm%0.04f,, #sigma_{fit}: %0.04f", hdEdx_poCosbin[i]->GetTitle(),
238 status_po.Data(), fdEdxMean_po, fdEdxMean_poErr, fdEdxSigma_po));
239 }
240
241 if (status_po != "FitOK" && status_el == "FitOK") {
242 fdEdxMean_po = fdEdxMean_el;
243 hdEdx_poCosbin[i]->SetTitle(Form("%s, mean (manual) = elec left", hdEdx_poCosbin[i]->GetTitle()));
244 } else if (status_el != "FitOK" && status_po == "FitOK") {
245 fdEdxMean_el = fdEdxMean_po;
246 hdEdx_elCosbin[i]->SetTitle(Form("%s, mean (manual) = posi right", hdEdx_elCosbin[i]->GetTitle()));
247 } else if (status_el != "FitOK" && status_po != "FitOK") {
248 fdEdxMean_po = 1.0; fdEdxMean_el = 1.0;
249 }
250
251 hdEdxMeanvsCos_po->SetBinContent(i + 1, fdEdxMean_po);
252 hdEdxMeanvsCos_po->SetBinError(i + 1, fdEdxMean_poErr);
253 hdEdxSigmavsCos_po->SetBinContent(i + 1, fdEdxSigma_po);
254 hdEdxSigmavsCos_po->SetBinError(i + 1, fdEdxSigma_poErr);
255
256 //for validation purpose
257 if (isMakePlots) {
258
259 ctmp_ep->cd(1); // each canvas is 2x2
260 hdEdx_elCosbin[i]->SetFillColorAlpha(kYellow, 0.25);
261 hdEdx_elCosbin[i]->DrawCopy("");
262 tl->SetX1(fdEdxMean_el); tl->SetX2(fdEdxMean_el);
263 tl->SetY1(0); tl->SetY2(hdEdx_elCosbin[i]->GetMaximum());
264 tl->DrawClone("same");
265
266 ctmp_ep->cd(2); // each canvas is 2x2
267 hdEdx_poCosbin[i]->SetFillColorAlpha(kBlue, 0.25);
268 hdEdx_poCosbin[i]->DrawCopy("");
269 tl->SetX1(fdEdxMean_po); tl->SetX2(fdEdxMean_po);
270 tl->SetY1(0); tl->SetY2(hdEdx_poCosbin[i]->GetMaximum());
271 tl->DrawClone("same");
272 ctmp_ep->Print(psname_ep.str().c_str());
273 }
274
275 //avg of both e+ and e- fdEdxMean
276 fdEdxMean = 0.5 * (fdEdxMean_po + fdEdxMean_el);
277 if (fdEdxMean <= 0)fdEdxMean = 1.0; //protection only
278 fdEdxMeanErr = 0.5 * TMath::Sqrt(fdEdxMean_elErr * fdEdxMean_elErr + fdEdxMean_poErr * fdEdxMean_poErr);
279 hdEdxMeanvsCos_ep->SetBinContent(i + 1, fdEdxMean);
280 hdEdxMeanvsCos_ep->SetBinError(i + 1, fdEdxMeanErr);
281 }
282
283 cosine.push_back(fdEdxMean);
284 delete tl;
285 }
286
287 //more validation plots for debugging
288 if (isMakePlots) {
289
290 psname_ep.str("");
291 psname_ep << Form("cdcdedx_coscal_fits_frun%d.pdf]", fStartRun);
292 ctmp_ep->Print(psname_ep.str().c_str());
293 delete ctmp_ep;
294
295 TCanvas* cstats = new TCanvas("cstats", "cstats", 1000, 500);
296 cstats->SetBatch(kTRUE);
297 cstats->Divide(2, 1);
298 cstats->cd(1);
299 auto hestats = getObjectPtr<TH1I>("hestats");
300 if (hestats) {
301 hestats->SetName(Form("hestats_fRun%d", fStartRun));
302 hestats->SetStats(0);
303 hestats->DrawCopy("");
304 }
305 cstats->cd(2);
306 auto htstats = getObjectPtr<TH1I>("htstats");
307 if (htstats) {
308 hestats->SetName(Form("htstats_fRun%d", fStartRun));
309 htstats->DrawCopy("");
310 hestats->SetStats(0);
311 }
312 cstats->Print(Form("cdcdedx_coscal_stats_frun%d.pdf", fStartRun));
313 delete cstats;
314
315 TCanvas* ctmp_epConst = new TCanvas("ctmp_epConst", "ctmp_epConst", 800, 400);
316 ctmp_epConst->Divide(2, 1);
317
318 TCanvas* ctmp_epCosth = new TCanvas("ctmp_epCosth", "ctmp_epCosth", 600, 500);
319
320 if (isMethodSep) {
321
322 ctmp_epConst->cd(1);
323 gPad->SetGridy(1);
324 hdEdxMeanvsCos_el->SetMarkerStyle(20);
325 hdEdxMeanvsCos_el->SetMarkerSize(0.60);
326 hdEdxMeanvsCos_el->SetMarkerColor(kRed);
327 hdEdxMeanvsCos_el->SetStats(0);
328 hdEdxMeanvsCos_el->SetTitle("comparison of dedx #mu_{fit}^{rel}: (e-=red, e+=blue, avg=black)");
329 hdEdxMeanvsCos_el->GetYaxis()->SetRangeUser(0.96, 1.04);
330 hdEdxMeanvsCos_el->DrawCopy("");
331
332 hdEdxMeanvsCos_po->SetMarkerStyle(20);
333 hdEdxMeanvsCos_po->SetMarkerSize(0.60);
334 hdEdxMeanvsCos_po->SetMarkerColor(kBlue);
335 hdEdxMeanvsCos_po->SetStats(0);
336 hdEdxMeanvsCos_po->DrawCopy("same");
337
338 hdEdxMeanvsCos_ep->SetMarkerStyle(20);
339 hdEdxMeanvsCos_ep->SetMarkerSize(0.60);
340 hdEdxMeanvsCos_ep->SetMarkerColor(kBlack);
341 hdEdxMeanvsCos_ep->SetStats(0);
342 hdEdxMeanvsCos_ep->DrawCopy("same");
343
344 ctmp_epConst->cd(2);
345 gPad->SetGridy(1);
346 hdEdxSigmavsCos_el->SetMarkerStyle(4);
347 hdEdxSigmavsCos_el->SetMarkerColor(kRed);
348 hdEdxSigmavsCos_el->SetMarkerSize(0.90);
349 hdEdxSigmavsCos_el->SetTitle("comparison of dedx #mu_{fit}^{rel}: (e-=open, e+=closed)");
350 hdEdxSigmavsCos_el->GetYaxis()->SetRangeUser(0.4, 0.12);
351 hdEdxSigmavsCos_el->SetStats(0);
352 hdEdxSigmavsCos_el->DrawCopy("");
353
354 hdEdxSigmavsCos_po->SetMarkerStyle(8);
355 hdEdxSigmavsCos_po->SetMarkerSize(0.80);
356 hdEdxSigmavsCos_po->SetMarkerColor(kBlue);
357 hdEdxSigmavsCos_po->SetStats(0);
358 hdEdxSigmavsCos_po->DrawCopy("same");
359
360 ctmp_epCosth->cd();
361 hCosth_el->SetStats(0);
362 hCosth_el->SetLineColor(kRed);
363 hCosth_el->SetFillColorAlpha(kYellow, 0.55);
364 hCosth_el->DrawCopy("");
365 hCosth_po->SetStats(0);
366 hCosth_po->SetLineColor(kBlue);
367 hCosth_po->SetFillColorAlpha(kGray, 0.35);
368 hCosth_po->DrawCopy("same");
369
370 } else {
371
372 ctmp_epConst->cd(1);
373 gPad->SetGridy(1);
374 hdEdxMeanvsCos_ep->SetMarkerStyle(20);
375 hdEdxMeanvsCos_ep->SetMarkerSize(0.60);
376 hdEdxMeanvsCos_ep->SetMarkerColor(kBlack);
377 hdEdxMeanvsCos_ep->SetStats(0);
378 hdEdxMeanvsCos_ep->SetTitle("dedx rel(#mu_{fit}) for e- and e+ combined");
379 hdEdxMeanvsCos_ep->GetYaxis()->SetRangeUser(0.97, 1.04);
380 hdEdxMeanvsCos_ep->DrawCopy("");
381
382 ctmp_epConst->cd(2);
383 gPad->SetGridy(1);
384 hdEdxSigmavsCos_ep->SetMarkerStyle(20);
385 hdEdxSigmavsCos_ep->SetMarkerColor(kRed);
386 hdEdxSigmavsCos_ep->SetMarkerSize(1.1);
387 hdEdxSigmavsCos_ep->SetTitle("dedx rel(#sigma_{fit}) for e- and e+ combined");
388 hdEdxSigmavsCos_ep->GetYaxis()->SetRangeUser(0.4, 0.12);
389 hdEdxSigmavsCos_ep->SetStats(0);
390 hdEdxSigmavsCos_ep->DrawCopy("");
391
392 ctmp_epCosth->cd();
393 hCosth_ep->SetStats(0);
394 hCosth_ep->SetLineColor(kGray);
395 hCosth_ep->SetFillColorAlpha(kGray, 0.25);
396 hCosth_ep->DrawCopy("same");
397 }
398
399 ctmp_epCosth->SaveAs(Form("cdcdedx_coscal_costhdist_frun%d.pdf", fStartRun));
400 delete ctmp_epCosth;
401
402 ctmp_epConst->SaveAs(Form("cdcdedx_coscal_relmeans_frun%d.pdf", fStartRun));
403 ctmp_epConst->SaveAs(Form("cdcdedx_coscal_relmeans_frun%d.root", fStartRun));
404 delete ctmp_epConst;
405 }
406
407 generateNewPayloads(cosine);
408 return c_OK;
409}
410
411void CDCDedxCosineAlgorithm::generateNewPayloads(std::vector<double> cosine)
412{
413
414 TH1D* hCosCorrOld = new TH1D(Form("hCosCorrOld_fRun%d", fStartRun),
415 Form("cos corr const comparison (red=old, blue=new), start run: %d;cos(#theta);dE/dx #mu_{fit}", fStartRun), fCosbins, fCosMin,
416 fCosMax);
417 TH1D* hCosCorrNew = new TH1D(Form("hCosCorrNew_fRun%d", fStartRun), Form("coss corr, start run: %d;cos(#theta);dE/dx #mu_{fit}",
418 fStartRun), fCosbins, fCosMin, fCosMax);
419 TH1D* hCosCorrRel = new TH1D(Form("hCosCorrRel_fRun%d", fStartRun),
420 Form("new relative cos corr, start run: %d;cos(#theta);dE/dx #mu_{fit}", fStartRun), fCosbins, fCosMin, fCosMax);
421
422 if (isMergePayload) {
423 const auto expRun = getRunList()[0];
424 updateDBObjPtrs(1, expRun.second, expRun.first);
425 // bool refchange = m_DBCosineCor.hasChanged(); //Add this feature for major processing
426 B2INFO("Saving new rung for (Exp, Run) : (" << expRun.first << "," << expRun.second << ")");
427 for (unsigned int ibin = 0; ibin < m_DBCosineCor->getSize(); ibin++) {
428 hCosCorrOld->SetBinContent(ibin + 1, (double)m_DBCosineCor->getMean(ibin));
429 hCosCorrRel->SetBinContent(ibin + 1, cosine.at(ibin));
430 B2INFO("Cosine Corr for Bin # " << ibin << ", Previous = " << m_DBCosineCor->getMean(ibin) << ", Relative = " << cosine.at(
431 ibin) << ", Merged = " << m_DBCosineCor->getMean(ibin)*cosine.at(ibin));
432 cosine.at(ibin) *= (double)m_DBCosineCor->getMean(ibin);
433 hCosCorrNew->SetBinContent(ibin + 1, cosine.at(ibin));
434 }
435 }
436
437 if (isMakePlots) {
438 TCanvas* ctmp_const = new TCanvas("ctmp_const", "ctmp_const", 900, 450);
439 ctmp_const->Divide(2, 1);
440
441 ctmp_const->cd(1);
442 gPad->SetGridy(1);
443 gPad->SetGridx(1);
444 hCosCorrOld->SetStats(0);
445 hCosCorrOld->SetLineColor(kRed);
446 hCosCorrOld->GetYaxis()->SetRangeUser(0.64, 1.20);
447 hCosCorrOld->DrawCopy("");
448 hCosCorrNew->SetStats(0);
449 hCosCorrNew->SetLineColor(kBlue);
450 hCosCorrNew->DrawCopy("same");
451
452 ctmp_const->cd(2);
453 gPad->SetGridy(1);
454 hCosCorrRel->SetStats(0);
455 hCosCorrRel->GetYaxis()->SetRangeUser(0.97, 1.03);
456 hCosCorrRel->SetLineColor(kBlack);
457 hCosCorrRel->DrawCopy("");
458
459 ctmp_const->SaveAs(Form("cdcdedx_coscal_constants_frun%d.pdf", fStartRun));
460 ctmp_const->SaveAs(Form("cdcdedx_coscal_constants_frun%d.root", fStartRun));
461 delete ctmp_const;
462 }
463
464 B2INFO("dE/dx calibration done for CDC dE/dx _eltron saturation");
465 CDCDedxCosineCor* gain = new CDCDedxCosineCor(cosine);
466 saveCalibration(gain, "CDCDedxCosineCor");
467}
468
469void CDCDedxCosineAlgorithm::FitGaussianWRange(TH1D*& temphist, TString& status)
470{
471 if (temphist->Integral() < 2000) { //at least 1k bhabha events
472 B2INFO(Form("\tThis hist (%s) have insufficient entries to perform fit (%0.03f)", temphist->GetName(), temphist->Integral()));
473 status = "LowStats";
474 return;
475 } else {
476 temphist->GetXaxis()->SetRange(temphist->FindFirstBinAbove(0, 1), temphist->FindLastBinAbove(0, 1));
477 int fs = temphist->Fit("gaus", "QR");
478 if (fs != 0) {
479 B2INFO(Form("\tFit (round 1) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
480 status = "FitFailed";
481 return;
482 } else {
483 double fdEdxMean = temphist->GetFunction("gaus")->GetParameter(1);
484 double width = temphist->GetFunction("gaus")->GetParameter(2);
485 temphist->GetXaxis()->SetRangeUser(fdEdxMean - 5.0 * width, fdEdxMean + 5.0 * width);
486 fs = temphist->Fit("gaus", "QR", "", fdEdxMean - fSigLim * width, fdEdxMean + fSigLim * width);
487 if (fs != 0) {
488 B2INFO(Form("\tFit (round 2) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
489 status = "FitFailed";
490 return;
491 } else {
492 temphist->GetXaxis()->SetRangeUser(fdEdxMean - 5.0 * width, fdEdxMean + 5.0 * width);
493 B2INFO(Form("\tFit for hist (%s) successful (status = %d)", temphist->GetName(), fs));
494 status = "FitOK";
495 }
496 }
497 }
498}
Belle2::CDCDedxCosineAlgorithm::generateNewPayloads
void generateNewPayloads(std::vector< double > cosine)
function to store new payload after full calibration
Definition: CDCDedxCosineAlgorithm.cc:411
Belle2::CDCDedxCosineAlgorithm::fStartRun
int fStartRun
boundary start at this run
Definition: CDCDedxCosineAlgorithm.h:107
Belle2::CDCDedxCosineAlgorithm::fdEdxMax
double fdEdxMax
max dedx range for gain cal
Definition: CDCDedxCosineAlgorithm.h:106
Belle2::CDCDedxCosineAlgorithm::isMergePayload
bool isMergePayload
merge payload at the of calibration
Definition: CDCDedxCosineAlgorithm.h:99
Belle2::CDCDedxCosineAlgorithm::fCosMax
double fCosMax
max cosine angle for cal
Definition: CDCDedxCosineAlgorithm.h:103
Belle2::CDCDedxCosineAlgorithm::m_DBCosineCor
DBObjPtr< CDCDedxCosineCor > m_DBCosineCor
Electron saturation correction DB object.
Definition: CDCDedxCosineAlgorithm.h:108
Belle2::CDCDedxCosineAlgorithm::fCosMin
double fCosMin
min cosine angle for cal
Definition: CDCDedxCosineAlgorithm.h:102
Belle2::CDCDedxCosineAlgorithm::CDCDedxCosineAlgorithm
CDCDedxCosineAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
Definition: CDCDedxCosineAlgorithm.cc:21
Belle2::CDCDedxCosineAlgorithm::calibrate
virtual EResult calibrate() override
Cosine algorithm.
Definition: CDCDedxCosineAlgorithm.cc:44
Belle2::CDCDedxCosineAlgorithm::fdEdxMin
double fdEdxMin
min dedx range for gain cal
Definition: CDCDedxCosineAlgorithm.h:105
Belle2::CDCDedxCosineAlgorithm::FitGaussianWRange
void FitGaussianWRange(TH1D *&temphist, TString &status)
function to fit histogram in each cosine bin
Definition: CDCDedxCosineAlgorithm.cc:469
Belle2::CDCDedxCosineAlgorithm::fHistbins
int fHistbins
number of bins for dedx histogram
Definition: CDCDedxCosineAlgorithm.h:104
Belle2::CDCDedxCosineAlgorithm::isMethodSep
bool isMethodSep
if e+e- need to be consider sep
Definition: CDCDedxCosineAlgorithm.h:97
Belle2::CDCDedxCosineAlgorithm::fCosbins
unsigned int fCosbins
number of bins across cosine range
Definition: CDCDedxCosineAlgorithm.h:101
Belle2::CDCDedxCosineAlgorithm::isMakePlots
bool isMakePlots
produce plots for status
Definition: CDCDedxCosineAlgorithm.h:98
Belle2::CDCDedxCosineAlgorithm::fSigLim
double fSigLim
gaussian fit sigma limit
Definition: CDCDedxCosineAlgorithm.h:100
Belle2::CDCDedxCosineCor
dE/dx wire gain calibration constants
Definition: CDCDedxCosineCor.h:26
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:37
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:297
Belle2::CalibrationAlgorithm::updateDBObjPtrs
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
Definition: CalibrationAlgorithm.cc:404
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:321
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition: CalibrationAlgorithm.h:266
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:40
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:41
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:43
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17