Belle II Software  release-08-01-10
eclCosmicEAlgorithm.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 /* Own header. */
10 #include <ecl/calibration/eclCosmicEAlgorithm.h>
11 
12 /* ECL headers. */
13 #include <ecl/dataobjects/ECLElementNumbers.h>
14 #include <ecl/dbobjects/ECLCrystalCalib.h>
15 
16 /* ROOT headers. */
17 #include <TF1.h>
18 #include <TFile.h>
19 #include <TH2F.h>
20 #include <TMath.h>
21 #include <TROOT.h>
22 
23 using namespace std;
24 using namespace Belle2;
25 using namespace ECL;
26 
28 // cppcheck-suppress constParameter ; TF1 fit functions cannot have const parameters
29 double eclCosmicNovoConst(double* x, double* par)
30 {
31  double qc = 0.;
32 
33  double peak = par[1];
34  double width = par[2];
35  double sln4 = sqrt(log(4));
36  double y = par[3] * sln4;
37  double tail = -log(y + sqrt(1 + y * y)) / sln4;
38 
39  if (TMath::Abs(tail) < 1.e-7) {
40  qc = 0.5 * TMath::Power(((x[0] - peak) / width), 2);
41  } else {
42  const double qa = tail * sqrt(log(4.));
43  const double qb = sinh(qa) / qa;
44  const double qx = (x[0] - peak) / width * qb;
45  const double qy = 1. + tail * qx;
46 
47  if (qy > 1.E-7)
48  qc = 0.5 * (TMath::Power((log(qy) / tail), 2) + tail * tail);
49  else
50  qc = 15.0;
51  }
52  return par[0] * exp(-qc) + par[4];
53 }
54 
55 eclCosmicEAlgorithm::eclCosmicEAlgorithm(): CalibrationAlgorithm("eclCosmicECollector"), cellIDLo(1),
56  cellIDHi(ECLElementNumbers::c_NCrystals),
57  minEntries(150), maxIterations(10), tRatioMin(0.2), tRatioMax(0.25), performFits(true), findExpValues(false), storeConst(0)
58 {
59  setDescription(
60  "Perform energy calibration of ecl crystals by fitting a Novosibirsk function to energy deposited by cosmic rays"
61  );
62 }
63 
64 CalibrationAlgorithm::EResult eclCosmicEAlgorithm::calibrate()
65 {
66 
69  double limitTol = 0.0005; /*< tolerance for checking if a parameter is at the limit */
70  double minFitLimit = 1e-25; /*< cut off for labeling a fit as poor */
71  double minFitProbIter = 1e-8; /*< cut off for labeling a fit as poor if it also has many iterations */
72  double constRatio = 0.5; /*< Novosibirsk normalization must be greater than constRatio x constant term */
73  double peakMin(0.5), peakMax(1.75); /*< range for peak of measured energy distribution */
74  double peakTol = limitTol * (peakMax - peakMin); /*< fit is at limit if it is within peakTol of min or max */
75  double effSigMin(0.08), effSigMax(0.5); /*< range for effective sigma of measured energy distribution */
76  double effSigTol = limitTol * (effSigMax - effSigMin);
77  double etaMin(-3.), etaMax(1.); /*< Novosibirsk tail parameter range */
78  double etaNom(-0.41); /*< Nominal tail parameter */
79  double etaTol = limitTol * (etaMax - etaMin);
80  double constTol = 0.1; /*< if constant is less than constTol, it will be fixed to 0 */
81 
83  gROOT->SetBatch();
84 
87  TH1F* dummy;
88  dummy = (TH1F*)gROOT->FindObject("EnvsCrysSameRing");
89  if (dummy) {delete dummy;}
90  dummy = (TH1F*)gROOT->FindObject("EnvsCrysDifferentRing");
91  if (dummy) {delete dummy;}
92  dummy = (TH1F*)gROOT->FindObject("IntegralVsCrysSame");
93  if (dummy) {delete dummy;}
94  dummy = (TH1F*)gROOT->FindObject("IntegralVsCrysDifferent");
95  if (dummy) {delete dummy;}
96  dummy = (TH1F*)gROOT->FindObject("AverageExpECrysSame");
97  if (dummy) {delete dummy;}
98  dummy = (TH1F*)gROOT->FindObject("AverageExpECrysDifferent");
99  if (dummy) {delete dummy;}
100  dummy = (TH1F*)gROOT->FindObject("AverageElecCalibSame");
101  if (dummy) {delete dummy;}
102  dummy = (TH1F*)gROOT->FindObject("AverageElecCalibDifferent");
103  if (dummy) {delete dummy;}
104  dummy = (TH1F*)gROOT->FindObject("AverageInitialCalibSame");
105  if (dummy) {delete dummy;}
106  dummy = (TH1F*)gROOT->FindObject("AverageInitialCalibDifferent");
107  if (dummy) {delete dummy;}
108 
111  std::vector<std::shared_ptr<TH2F>> EnvsCrys;
112  EnvsCrys.push_back(getObjectPtr<TH2F>("EnvsCrysSameRing"));
113  EnvsCrys.push_back(getObjectPtr<TH2F>("EnvsCrysDifferentRing"));
114 
115  std::vector<std::shared_ptr<TH1F>> ExpEvsCrys;
116  ExpEvsCrys.push_back(getObjectPtr<TH1F>("ExpEvsCrysSameRing"));
117  ExpEvsCrys.push_back(getObjectPtr<TH1F>("ExpEvsCrysDifferentRing"));
118 
119  std::vector<std::shared_ptr<TH1F>> ElecCalibvsCrys;
120  ElecCalibvsCrys.push_back(getObjectPtr<TH1F>("ElecCalibvsCrysSameRing"));
121  ElecCalibvsCrys.push_back(getObjectPtr<TH1F>("ElecCalibvsCrysDifferentRing"));
122 
123  std::vector<std::shared_ptr<TH1F>> InitialCalibvsCrys;
124  InitialCalibvsCrys.push_back(getObjectPtr<TH1F>("InitialCalibvsCrysSameRing"));
125  InitialCalibvsCrys.push_back(getObjectPtr<TH1F>("InitialCalibvsCrysDifferentRing"));
126 
127  std::vector<std::shared_ptr<TH1F>> CalibEntriesvsCrys;
128  CalibEntriesvsCrys.push_back(getObjectPtr<TH1F>("CalibEntriesvsCrysSameRing"));
129  CalibEntriesvsCrys.push_back(getObjectPtr<TH1F>("CalibEntriesvsCrysDifferentRing"));
130 
131  auto RawDigitAmpvsCrys = getObjectPtr<TH2F>("RawDigitAmpvsCrys");
132 
136  TH1F* IntegralVsCrys[2];
137  IntegralVsCrys[0] = new TH1F("IntegralVsCrysSame", "Integral of EnVsCrys for each crystal, same theta ring;Crystal ID",
138  ECLElementNumbers::c_NCrystals, 0,
139  ECLElementNumbers::c_NCrystals);
140  IntegralVsCrys[1] = new TH1F("IntegralVsCrysDifferent", "Integral of EnVsCrys for each crystal, different theta rings;Crystal ID",
141  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
142 
143  TH1F* AverageExpECrys[2];
144  AverageExpECrys[0] = new TH1F("AverageExpECrysSame",
145  "Average expected E per crys from collector, same theta ring;Crystal ID;Energy (GeV)", ECLElementNumbers::c_NCrystals, 0,
146  ECLElementNumbers::c_NCrystals);
147  AverageExpECrys[1] = new TH1F("AverageExpECrysDifferent",
148  "Average expected E per crys from collector, different theta ring;Crystal ID;Energy (GeV)", ECLElementNumbers::c_NCrystals, 0,
149  ECLElementNumbers::c_NCrystals);
150 
151  TH1F* AverageElecCalib[2];
152  AverageElecCalib[0] = new TH1F("AverageElecCalibSame",
153  "Average electronics calib const vs crys, same theta ring;Crystal ID;Calibration constant", ECLElementNumbers::c_NCrystals, 0,
154  ECLElementNumbers::c_NCrystals);
155  AverageElecCalib[1] = new TH1F("AverageElecCalibDifferent",
156  "Average electronics calib const vs crys, different theta rings;Crystal ID;Calibration constant", ECLElementNumbers::c_NCrystals, 0,
157  ECLElementNumbers::c_NCrystals);
158 
159  TH1F* AverageInitialCalib[2];
160  AverageInitialCalib[0] = new TH1F("AverageInitialCalibSame",
161  "Average initial cosmic calib const vs crys, same theta ring;Crystal ID;Calibration constant", ECLElementNumbers::c_NCrystals, 0,
162  ECLElementNumbers::c_NCrystals);
163  AverageInitialCalib[1] = new TH1F("AverageInitialCalibDifferent",
164  "Average initial cosmic calib const vs crys, different theta rings;Crystal ID;Calibration constant", ECLElementNumbers::c_NCrystals,
165  0, ECLElementNumbers::c_NCrystals);
166 
167  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
168  int histbin = crysID + 1;
169  for (int idir = 0; idir < 2; idir++) {
170  TH1D* hEnergy = EnvsCrys[idir]->ProjectionY("hEnergy", histbin, histbin);
171  int Integral = hEnergy->Integral();
172  IntegralVsCrys[idir]->SetBinContent(histbin, Integral);
173 
174  double TotEntries = CalibEntriesvsCrys[idir]->GetBinContent(histbin);
175 
176  double expectedE = 0.;
177  if (TotEntries > 0.) {expectedE = ExpEvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
178  AverageExpECrys[idir]->SetBinContent(histbin, expectedE);
179  AverageExpECrys[idir]->SetBinError(histbin, 0.);
180 
181  double calibconst = 0.;
182  if (TotEntries > 0.) {calibconst = ElecCalibvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
183  AverageElecCalib[idir]->SetBinContent(histbin, calibconst);
184  AverageElecCalib[idir]->SetBinError(histbin, 0);
185 
186  calibconst = 0.;
187  if (TotEntries > 0.) {calibconst = InitialCalibvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
188  AverageInitialCalib[idir]->SetBinContent(histbin, calibconst);
189  AverageInitialCalib[idir]->SetBinError(histbin, 0);
190  }
191  }
192 
195  TFile* histfile = new TFile("eclCosmicEAlgorithm.root", "recreate");
196  for (int idir = 0; idir < 2; idir++) {
197  EnvsCrys[idir]->Write();
198  IntegralVsCrys[idir]->Write();
199  AverageExpECrys[idir]->Write();
200  AverageElecCalib[idir]->Write();
201  AverageInitialCalib[idir]->Write();
202  }
203  RawDigitAmpvsCrys->Write();
204 
207  if (!performFits) {
208  B2INFO("eclCosmicEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
209  histfile->Close();
210  return c_NotEnoughData;
211  }
212 
217  bool sufficientData = true;
218  for (int crysID = cellIDLo - 1; crysID < cellIDHi; crysID++) {
219  int histbin = crysID + 1;
220  bool SameLow = IntegralVsCrys[0]->GetBinContent(histbin) < minEntries;
221  bool DifferentLow = IntegralVsCrys[1]->GetBinContent(histbin) < minEntries;
222  if ((SameLow && DifferentLow) || (findExpValues && (SameLow || DifferentLow))) {
223  if (storeConst == 1) {B2INFO("eclCosmicEAlgorithm: cellID " << histbin << " has insufficient statistics: " << IntegralVsCrys[0]->GetBinContent(histbin) << " and " << IntegralVsCrys[1]->GetBinContent(histbin) << ". Requirement is " << minEntries);}
224  sufficientData = false;
225  break;
226  }
227  }
228 
231  if (!sufficientData && storeConst == 1) {
232  histfile->Close();
233  return c_NotEnoughData;
234  }
235 
238  const TString preName[2] = {"SameRing", "DifferentRing"};
239 
240  TH1F* PeakperCrys[2];
241  PeakperCrys[0] = new TH1F("PeakperCrysSame", "Fit peak per crystal, same theta ring;Crystal ID;Peak normalized energy",
242  ECLElementNumbers::c_NCrystals, 0,
243  ECLElementNumbers::c_NCrystals);
244  PeakperCrys[1] = new TH1F("PeakperCrysDifferent", "Fit peak per crystal, different theta ring;Crystal ID;Peak normalized energy",
245  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
246 
247  TH1F* SigmaperCrys[2];
248  SigmaperCrys[0] = new TH1F("SigmaperCrysSame", "Fit sigma per crysal, same theta ring;Crystal ID;Sigma (ADC)",
249  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
250  SigmaperCrys[1] = new TH1F("SigmaperCrysDifferent", "Fit sigma per crysal, different theta ring;Crystal ID;Sigma (ADC)",
251  ECLElementNumbers::c_NCrystals, 0,
252  ECLElementNumbers::c_NCrystals);
253 
254  TH1F* EtaperCrys[2];
255  EtaperCrys[0] = new TH1F("EtaperCrysSame", "Fit eta per crysal, same theta ring;Crystal ID;Eta", ECLElementNumbers::c_NCrystals, 0,
256  ECLElementNumbers::c_NCrystals);
257  EtaperCrys[1] = new TH1F("EtaperCrysDifferent", "Fit eta per crysal, different theta ring;Crystal ID;Eta",
258  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
259 
260  TH1F* ConstperCrys[2];
261  ConstperCrys[0] = new TH1F("ConstperCrysSame", "Fit constant per crystal, same theta ring;Crystal ID;Constant",
262  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
263  ConstperCrys[1] = new TH1F("ConstperCrysDifferent", "Fit constant per crystal, different theta ring;Crystal ID;Constant",
264  ECLElementNumbers::c_NCrystals, 0,
265  ECLElementNumbers::c_NCrystals);
266 
267  TH1F* StatusperCrys[2];
268  StatusperCrys[0] = new TH1F("StatusperCrysSame", "Fit status, same theta ring;Crystal ID;Status", ECLElementNumbers::c_NCrystals, 0,
269  ECLElementNumbers::c_NCrystals);
270  StatusperCrys[1] = new TH1F("StatusperCrysDifferent", "Fit status, different theta ring;Crystal ID;Status",
271  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
272 
274  TH1F* hStatus[2];
275  hStatus[0] = new TH1F("StatusSame", "Fit status, same theta ring", 25, -5, 20);
276  hStatus[1] = new TH1F("StatusDifferent", "Fit status, different theta ring", 25, -5, 20);
277 
278  TH1F* fracPeakUnc[2];
279  fracPeakUnc[0] = new TH1F("fracPeakUncSame", "Fractional uncertainty on peak location, same theta ring", 100, 0, 0.1);
280  fracPeakUnc[1] = new TH1F("fracPeakUncDifferent", "Fractional uncertainty on peak location, different theta ring", 100, 0, 0.1);
281 
282  TH1F* hfitProb[2];
283  hfitProb[0] = new TH1F("fitProbSame", "Probability of fit, same theta ring", 200, 0, 0.02);
284  hfitProb[1] = new TH1F("fitProbDifferent", "Probability of fit, different theta ring", 200, 0, 0.02);
285 
287  TH1F* ExpEnergyperCrys[2];
288  ExpEnergyperCrys[0] = new TH1F("ExpEnergyperCrysSame", "Expected energy per crystal, same theta ring;Crystal ID;Peak energy (GeV)",
289  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
290  ExpEnergyperCrys[1] = new TH1F("ExpEnergyperCrysDifferent",
291  "Expected energy per crystal, different theta ring;Crystal ID;Peak energy (GeV)", ECLElementNumbers::c_NCrystals, 0,
292  ECLElementNumbers::c_NCrystals);
293 
294  TH1F* CalibvsCrys = new TH1F("CalibvsCrys", "Energy calibration constant from cosmics;Crystal ID;Calibration constant",
295  ECLElementNumbers::c_NCrystals, 0,
296  ECLElementNumbers::c_NCrystals);
297 
300  bool allFitsOK = true;
301  for (int crysID = cellIDLo - 1; crysID < cellIDHi; crysID++) {
302  int histbin = crysID + 1;
303  for (int idir = 0; idir < 2; idir++) {
304 
306  TString hname = preName[idir];
307  hname += "Enormalized";
308  hname += crysID;
309  TH1D* hEnergy = EnvsCrys[idir]->ProjectionY(hname, histbin, histbin);
310 
312  double histMin = hEnergy->GetXaxis()->GetXmin();
313  double histMax = hEnergy->GetXaxis()->GetXmax();
314  TF1* func = new TF1("eclCosmicNovoConst", eclCosmicNovoConst, histMin, histMax, 5);
315  func->SetParNames("normalization", "peak", "effSigma", "eta", "const");
316  func->SetParLimits(1, peakMin, peakMax);
317  func->SetParLimits(2, effSigMin, effSigMax);
318  func->SetParLimits(3, etaMin, etaMax);
319 
321  hEnergy->GetXaxis()->SetRangeUser(peakMin, peakMax);
322  int maxBin = hEnergy->GetMaximumBin();
323  double peakE = hEnergy->GetBinLowEdge(maxBin);
324  double peakEUnc = 0.;
325  double normalization = hEnergy->GetMaximum();
326  double effSigma = hEnergy->GetRMS();
327  double sigmaUnc = 0.;
328  hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
329 
331  double fitlow = 0.25;
332  double fithigh = peakE + 2.5 * effSigma;
333 
335  int il0 = hEnergy->GetXaxis()->FindBin(fitlow);
336  int il1 = hEnergy->GetXaxis()->FindBin(fitlow + 0.1);
337  double constant = hEnergy->Integral(il0, il1) / (1 + il1 - il0);
338  double constUnc = 0.;
339 
341  double eta = etaNom;
342  double etaUnc = 0.;
343 
345  double dIter = 0.1 * (histMax - histMin) / hEnergy->GetNbinsX();
346  double highold(0.), higholdold(0.);
347  double fitProb(0.);
348  double fitProbDefault(0.);
349  bool fitHist = IntegralVsCrys[idir]->GetBinContent(histbin) >= minEntries; /* fit only if enough events */
350  bool fixConst = false;
351  int nIter = 0;
352 
355  while (fitHist) {
356  nIter++;
357 
359  func->SetParameters(normalization, peakE, effSigma, eta, constant);
360  if (fixConst) { func->FixParameter(4, 0); }
361 
363  hEnergy->Fit(func, "LIQ", "", fitlow, fithigh);
364  normalization = func->GetParameter(0);
365  peakE = func->GetParameter(1);
366  peakEUnc = func->GetParError(1);
367  effSigma = func->GetParameter(2);
368  sigmaUnc = func->GetParError(2);
369  eta = func->GetParameter(3);
370  etaUnc = func->GetParError(3);
371  constant = func->GetParameter(4);
372  constUnc = func->GetParError(4);
373  fitProbDefault = func->GetProb();
374 
376  fitHist = false;
377  double peak = func->Eval(peakE) - constant;
378  double tRatio = (func->Eval(fithigh) - constant) / peak;
379  if (tRatio < tRatioMin || tRatio > tRatioMax) {
380  double targetY = constant + 0.5 * (tRatioMin + tRatioMax) * peak;
381  higholdold = highold;
382  highold = fithigh;
383  fithigh = func->GetX(targetY, peakE, histMax);
384  fitHist = true;
385 
387  if (abs(fithigh - higholdold) < dIter) {fithigh = 0.5 * (highold + higholdold); }
388 
390  if (nIter > maxIterations - 3) {fithigh = 0.33333 * (fithigh + highold + higholdold); }
391  }
392 
394  if (constant < constTol && !fixConst) {
395  constant = 0;
396  fixConst = true;
397  fitHist = true;
398  }
399 
401  if (nIter == maxIterations) {fitHist = false;}
402  B2DEBUG(200, "cellID = " << histbin << " " << nIter << " " << preName[idir] << " " << peakE << " " << constant << " " << tRatio <<
403  " " << fithigh);
404 
405  }
406 
409  fitProb = 0.;
410  if (nIter > 0) {
411  int lowbin = hEnergy->GetXaxis()->FindBin(fitlow);
412  int highbin = hEnergy->GetXaxis()->FindBin(fithigh);
413  int npar = 5;
414  if (fixConst) {npar = 4;}
415  int ndeg = (highbin - lowbin) + 1 - npar;
416  double chisq = 0.;
417  double halfbinwidth = 0.5 * hEnergy->GetBinWidth(1);
418  for (int ib = lowbin; ib <= highbin; ib++) {
419  double yexp = func->Eval(hEnergy->GetBinLowEdge(ib) + halfbinwidth);
420  double yobs = hEnergy->GetBinContent(ib);
421  double dchi2 = (yexp - yobs) * (yexp - yobs) / yexp;
422  chisq += dchi2;
423  }
424  fitProb = 0.5 * (TMath::Prob(chisq, ndeg) + fitProbDefault);
425  }
426 
429  int iStatus = fitOK; // success
430  if (nIter == maxIterations) {iStatus = iterations;} // too many iterations
431 
433  if (normalization < constRatio * constant) {iStatus = noPeak;}
434 
436  if (fitProb <= minFitLimit || (fitProb < minFitProbIter && iStatus == iterations)) {iStatus = poorFit;}
437 
439  if ((peakE < peakMin + peakTol) || (peakE > peakMax - peakTol)) {iStatus = atLimit;}
440  if ((effSigma < effSigMin + effSigTol) || (effSigma > effSigMax - effSigTol)) {iStatus = atLimit;}
441  if ((eta < etaMin + etaTol) || (eta > etaMax - etaTol)) {iStatus = atLimit;}
442 
443  //**..No fit
444  if (nIter == 0) {iStatus = notFit;} // not fit
445 
447  PeakperCrys[idir]->SetBinContent(histbin, peakE);
448  PeakperCrys[idir]->SetBinError(histbin, peakEUnc);
449  SigmaperCrys[idir]->SetBinContent(histbin, effSigma);
450  SigmaperCrys[idir]->SetBinError(histbin, sigmaUnc);
451  EtaperCrys[idir]->SetBinContent(histbin, eta);
452  EtaperCrys[idir]->SetBinError(histbin, etaUnc);
453  ConstperCrys[idir]->SetBinContent(histbin, constant);
454  ConstperCrys[idir]->SetBinError(histbin, constUnc);
455  StatusperCrys[idir]->SetBinContent(histbin, iStatus);
456  hStatus[idir]->Fill(iStatus);
457  fracPeakUnc[idir]->Fill(peakEUnc / peakE);
458  hfitProb[idir]->Fill(fitProb);
459 
461  B2INFO("cellID " << histbin << " " << preName[idir] << " status = " << iStatus << " fit probability = " << fitProb);
462  histfile->cd();
463  hEnergy->Write();
464  }
465  }
466 
469  if (findExpValues) {
470 
472  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
473  int histbin = crysID + 1;
474  bool atLeastOneOK = false;
475  for (int idir = 0; idir < 2; idir++) {
476  double fitstatus = StatusperCrys[idir]->GetBinContent(histbin);
477  double peakE = PeakperCrys[idir]->GetBinContent(histbin);
478  double peakEUnc = PeakperCrys[idir]->GetBinError(histbin);
479 
480  //**..For failed fits, store the negative of the input expected energy */
481  if (fitstatus < iterations) {
482  if (histbin >= cellIDLo && histbin <= cellIDHi) {
483  B2INFO("eclCosmicEAlgorithm: crystal " << crysID << " " << preName[idir] << " is not a successful fit. Status = " << fitstatus);
484  }
485  peakE = -1.;
486  peakEUnc = 0.;
487  } else {
488  atLeastOneOK = true;
489  }
490  double inputExpE = abs(AverageExpECrys[idir]->GetBinContent(histbin));
491  ExpEnergyperCrys[idir]->SetBinContent(histbin, inputExpE * peakE);
492  ExpEnergyperCrys[idir]->SetBinError(histbin, inputExpE * peakEUnc / peakE);
493  }
494  if (!atLeastOneOK) {allFitsOK = false;}
495  }
496 
499  } else {
500 
502  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
503  int histbin = crysID + 1;
504  double calibConst[2] = {};
505  double calibConstUnc[2] = {999999., 999999.};
506  double weight[2] = {};
507  bool bothFitsBad = true;
508  for (int idir = 0; idir < 2; idir++) {
509 
511  double peakE = PeakperCrys[idir]->GetBinContent(histbin);
512  double fracPeakEUnc = PeakperCrys[idir]->GetBinError(histbin) / peakE;
513  double inputConst = AverageInitialCalib[idir]->GetBinContent(histbin);
514  double fitstatus = StatusperCrys[idir]->GetBinContent(histbin);
515  double inputExpE = AverageExpECrys[idir]->GetBinContent(histbin);
516  if (fitstatus >= iterations && inputConst == 0) {B2FATAL("eclCosmicEAlgorithm: input calibration = 0 for idir = " << idir << " and crysID = " << crysID);}
517 
518  //** Find constant only if fit was successful and we have a value for the expected energy */
519  if (fitstatus >= iterations && inputExpE > 0.) {
520  calibConst[idir] = abs(inputConst) / peakE;
521  calibConstUnc[idir] = calibConst[idir] * fracPeakEUnc / peakE;
522  weight[idir] = 1. / (calibConstUnc[idir] * calibConstUnc[idir]);
523  bothFitsBad = false;
524  }
525  if (fitstatus < iterations && histbin >= cellIDLo && histbin <= cellIDHi) {
526  B2INFO("eclCosmicEAlgorithm: cellID " << histbin << " " << preName[idir] << " is not a successful fit. Status = " << fitstatus);
527  } else if (inputExpE < 0. && histbin >= cellIDLo && histbin <= cellIDHi) {
528  B2WARNING("eclCosmicEAlgorithm: cellID " << histbin << " " << preName[idir] << " has no expected energy. Status = " << fitstatus);
529  }
530  }
531 
532 
534  double averageConst;
535  double averageConstUnc;
536 
538  if (bothFitsBad) {
539  if (histbin >= cellIDLo && histbin <= cellIDHi) {B2INFO("eclCosmicEAlgorithm: no constant found for cellID = " << histbin << " status = " << StatusperCrys[0]->GetBinContent(histbin) << " and " << StatusperCrys[1]->GetBinContent(histbin));}
540  averageConst = -1.*abs(AverageInitialCalib[0]->GetBinContent(histbin));
541  averageConstUnc = 0.;
542  } else {
543  averageConst = (calibConst[0] * weight[0] + calibConst[1] * weight[1]) / (weight[0] + weight[1]);
544  averageConstUnc = 1. / sqrt(weight[0] + weight[1]);
545  }
546  CalibvsCrys->SetBinContent(histbin, averageConst);
547  CalibvsCrys->SetBinError(histbin, averageConstUnc);
548  }
549  }
550 
553  bool DBsuccess = false;
554  if (storeConst == 0 || (storeConst == 1 && allFitsOK)) {
555  DBsuccess = true;
556 
558  if (findExpValues) {
559  std::vector<std::string> DBname = {"ECLExpCosmicESame", "ECLExpCosmicEDifferent"};
560  for (int idir = 0; idir < 2; idir++) {
561  std::vector<float> tempE;
562  std::vector<float> tempUnc;
563  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
564  int histbin = crysID + 1;
565  tempE.push_back(ExpEnergyperCrys[idir]->GetBinContent(histbin));
566  tempUnc.push_back(ExpEnergyperCrys[idir]->GetBinError(histbin));
567  }
568  ECLCrystalCalib* ExpectedE = new ECLCrystalCalib();
569  ExpectedE->setCalibVector(tempE, tempUnc);
570  saveCalibration(ExpectedE, DBname[idir]);
571  B2INFO("eclCosmicEAlgorithm: successfully stored expected values for " << DBname[idir]);
572  }
573 
575  } else {
576  std::vector<float> tempCalib;
577  std::vector<float> tempCalibUnc;
578  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
579  int histbin = crysID + 1;
580  tempCalib.push_back(CalibvsCrys->GetBinContent(histbin));
581  tempCalibUnc.push_back(CalibvsCrys->GetBinError(histbin));
582  }
583  ECLCrystalCalib* CosmicECalib = new ECLCrystalCalib();
584  CosmicECalib->setCalibVector(tempCalib, tempCalibUnc);
585  saveCalibration(CosmicECalib, "ECLCrystalEnergyCosmic");
586  B2INFO("eclCosmicEAlgorithm: successfully stored calibration constants");
587  }
588  }
589 
591  for (int idir = 0; idir < 2; idir++) {
592  PeakperCrys[idir]->Write();
593  SigmaperCrys[idir]->Write();
594  EtaperCrys[idir]->Write();
595  ConstperCrys[idir]->Write();
596  StatusperCrys[idir]->Write();
597  hStatus[idir]->Write();
598  fracPeakUnc[idir]->Write();
599  hfitProb[idir]->Write();
600  }
601 
603  if (findExpValues) {
604  ExpEnergyperCrys[0]->Write();
605  ExpEnergyperCrys[1]->Write();
606  } else {
607  CalibvsCrys->Write();
608  }
609  histfile->Close();
610 
613  dummy = (TH1F*)gROOT->FindObject("PeakperCrysSame"); delete dummy;
614  dummy = (TH1F*)gROOT->FindObject("SigmaperCrysSame"); delete dummy;
615  dummy = (TH1F*)gROOT->FindObject("EtaperCrysSame"); delete dummy;
616  dummy = (TH1F*)gROOT->FindObject("ConstperCrysSame"); delete dummy;
617  dummy = (TH1F*)gROOT->FindObject("StatusperCrysSame"); delete dummy;
618  dummy = (TH1F*)gROOT->FindObject("StatusSame"); delete dummy;
619  dummy = (TH1F*)gROOT->FindObject("fracPeakUncSame"); delete dummy;
620  dummy = (TH1F*)gROOT->FindObject("fitProbSame"); delete dummy;
621  dummy = (TH1F*)gROOT->FindObject("ExpEnergyperCrysSame"); delete dummy;
622  dummy = (TH1F*)gROOT->FindObject("PeakperCrysDifferent"); delete dummy;
623  dummy = (TH1F*)gROOT->FindObject("SigmaperCrysDifferent"); delete dummy;
624  dummy = (TH1F*)gROOT->FindObject("EtaperCrysDifferent"); delete dummy;
625  dummy = (TH1F*)gROOT->FindObject("ConstperCrysDifferent"); delete dummy;
626  dummy = (TH1F*)gROOT->FindObject("StatusperCrysDifferent"); delete dummy;
627  dummy = (TH1F*)gROOT->FindObject("StatusDifferent"); delete dummy;
628  dummy = (TH1F*)gROOT->FindObject("fracPeakUncDifferent"); delete dummy;
629  dummy = (TH1F*)gROOT->FindObject("fitProbDifferent"); delete dummy;
630  dummy = (TH1F*)gROOT->FindObject("ExpEnergyperCrysDifferent"); delete dummy;
631  dummy = (TH1F*)gROOT->FindObject("CalibvsCrys"); delete dummy;
632 
635  if (storeConst == -1) {
636  B2INFO("eclCosmicEAlgorithm performed fits but was not asked to store contants");
637  return c_Failure;
638  } else if (!DBsuccess) {
639  if (findExpValues) { B2INFO("eclCosmicEAlgorithm: failed to store expected values"); }
640  else { B2INFO("eclCosmicEAlgorithm: failed to store calibration constants"); }
641  return c_Failure;
642  }
643  return c_OK;
644 }
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::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::c_Failure
@ c_Failure
Failed =3 in Python.
Definition: CalibrationAlgorithm.h:44
Belle2::ECLCrystalCalib
General DB object to store one calibration number per ECL crystal.
Definition: ECLCrystalCalib.h:27
Belle2::ECLCrystalCalib::setCalibVector
void setCalibVector(const std::vector< float > &CalibConst, const std::vector< float > &CalibConstUnc)
Set vector of constants with uncertainties.
Definition: ECLCrystalCalib.h:41
Belle2::ECL::eclCosmicEAlgorithm::tRatioMin
double tRatioMin
Fit range is adjusted so that fit at upper endpoint is between tRatioMin and tRatioMax of peak.
Definition: eclCosmicEAlgorithm.h:38
Belle2::ECL::eclCosmicEAlgorithm::poorFit
int poorFit
low chi square; fit not useable
Definition: eclCosmicEAlgorithm.h:56
Belle2::ECL::eclCosmicEAlgorithm::iterations
int iterations
fit reached max number of iterations, but is useable
Definition: eclCosmicEAlgorithm.h:54
Belle2::ECL::eclCosmicEAlgorithm::storeConst
int storeConst
controls which values are written to the database.
Definition: eclCosmicEAlgorithm.h:42
Belle2::ECL::eclCosmicEAlgorithm::performFits
bool performFits
if false, input histograms are copied to output, but no fits are done.
Definition: eclCosmicEAlgorithm.h:40
Belle2::ECL::eclCosmicEAlgorithm::cellIDHi
int cellIDHi
Last cellID to be fit.
Definition: eclCosmicEAlgorithm.h:35
Belle2::ECL::eclCosmicEAlgorithm::cellIDLo
int cellIDLo
Parameters to control Novosibirsk fit to signal measured in each crystal.
Definition: eclCosmicEAlgorithm.h:34
Belle2::ECL::eclCosmicEAlgorithm::findExpValues
bool findExpValues
if true, fits are used to find expected energy deposit for each crystal instead of the calibration co...
Definition: eclCosmicEAlgorithm.h:41
Belle2::ECL::eclCosmicEAlgorithm::minEntries
int minEntries
All crystals to be fit must have at least minEntries events in the fit range.
Definition: eclCosmicEAlgorithm.h:36
Belle2::ECL::eclCosmicEAlgorithm::calibrate
virtual EResult calibrate() override
Run algorithm on events.
Definition: eclCosmicEAlgorithm.cc:64
Belle2::ECL::eclCosmicEAlgorithm::noPeak
int noPeak
Novosibirsk component of fit is negligible; fit not useable.
Definition: eclCosmicEAlgorithm.h:57
Belle2::ECL::eclCosmicEAlgorithm::atLimit
int atLimit
a parameter is at the limit; fit not useable
Definition: eclCosmicEAlgorithm.h:55
Belle2::ECL::eclCosmicEAlgorithm::maxIterations
int maxIterations
no more than maxIteration iterations
Definition: eclCosmicEAlgorithm.h:37
Belle2::ECL::eclCosmicEAlgorithm::tRatioMax
double tRatioMax
Fit range is adjusted so that fit at upper endpoint is between tRatioMin and tRatioMax of peak.
Definition: eclCosmicEAlgorithm.h:39
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::ECLElementNumbers::c_NCrystals
const int c_NCrystals
Number of crystals.
Definition: ECLElementNumbers.h:23
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17