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