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