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