Belle II Software  release-08-01-10
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 
17 using namespace Belle2;
18 //-----------------------------------------------------------------
19 // Implementation
20 //-----------------------------------------------------------------
21 CDCDedxCosineAlgorithm::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 //-----------------------------------------------------------------
44 CalibrationAlgorithm::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 accpetance (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 
411 void 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 
469 void CDCDedxCosineAlgorithm::FitGaussianWRange(TH1D*& temphist, TString& status)
470 {
471  if (temphist->Integral() < 2000) { //atleast 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) sucessfull (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::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