Belle II Software  release-08-01-10
eclGammaGammaEAlgorithm.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/eclGammaGammaEAlgorithm.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 
29 // cppcheck-suppress constParameter ; TF1 fit functions cannot have const parameters
30 double eclGammaGammaNovoConst(double* x, double* par)
31 {
32  double qc = 0;
33 
34  double peak = par[1];
35  double width = par[2];
36  double sln4 = sqrt(log(4));
37  double y = par[3] * sln4;
38  double tail = -log(y + sqrt(1 + y * y)) / sln4;
39 
40  if (TMath::Abs(tail) < 1.e-7) {
41  qc = 0.5 * TMath::Power(((x[0] - peak) / width), 2);
42  } else {
43  double qa = tail * sqrt(log(4.));
44  double qb = sinh(qa) / qa;
45  double qx = (x[0] - peak) / width * qb;
46  double qy = 1. + tail * qx;
47 
48  if (qy > 1.E-7)
49  qc = 0.5 * (TMath::Power((log(qy) / tail), 2) + tail * tail);
50  else
51  qc = 15.0;
52  }
53  return par[0] * exp(-qc) + par[4];
54 }
55 
56 
58 eclGammaGammaEAlgorithm::eclGammaGammaEAlgorithm(): CalibrationAlgorithm("eclGammaGammaECollector")
59 {
60  setDescription(
61  "Perform energy calibration of ecl crystals by fitting a Novosibirsk function to energy deposited by photons in e+e- --> gamma gamma"
62  );
63 }
64 
65 CalibrationAlgorithm::EResult eclGammaGammaEAlgorithm::calibrate()
66 {
69  const double limitTol = 0.0005; /*< tolerance for checking if a parameter is at the limit */
70  const double minFitLimit = 1e-25; /*< cut off for labeling a fit as poor */
71  const double minFitProbIter = 1e-8; /*< cut off for labeling a fit as poor if it also has many iterations */
72  const double constRatio = 0.5; /*< Novosibirsk normalization must be greater than constRatio x constant term */
73  const double peakMin(0.4), peakMax(2.2); /*< range for peak of normalized energy distribution */
74  const double peakTol = limitTol * (peakMax - peakMin); /*< fit is at limit if it is within peakTol of min or max */
75  const double effSigMin(0.02), effSigMax(0.4); /*< range for effective sigma of normalized energy distribution */
76  const double effSigTol = limitTol * (effSigMax - effSigMin); /*< fit is at limit if it is within effSigTol of min or max */
77  const double etaNom(0.41); /*< Nominal tail parameter */
78  const double etaMin(0.), etaMax(5.0); /*< Novosibirsk tail parameter range */
79  const double etaTol = limitTol * (etaMax - etaMin); /*< fit is at limit if it is within etaTol of min or max */
80  const double constTol = 0.001; /*< if constant is less than constTol, it will be fixed to 0 */
81 
83  gROOT->SetBatch();
84 
87  B2INFO("eclGammaGammaAlgorithm parameters:");
88  B2INFO("outputName = " << m_outputName);
89  B2INFO("cellIDLo = " << m_cellIDLo);
90  B2INFO("cellIDHi = " << m_cellIDHi);
91  B2INFO("minEntries = " << m_minEntries);
92  B2INFO("highStatEntries = " << m_highStatEntries);
93  B2INFO("maxIterations = " << m_maxIterations);
94  B2INFO("tRatioMinNom = " << m_tRatioMinNom);
95  B2INFO("tRatioMaxNom = " << m_tRatioMaxNom);
96  B2INFO("tRatioMinHiStat = " << m_tRatioMinHiStat);
97  B2INFO("tRatioMaxHiStat = " << m_tRatioMaxHiStat);
98  B2INFO("upperEdgeThresh = " << m_upperEdgeThresh);
99  B2INFO("performFits = " << m_performFits);
100  B2INFO("findExpValues = " << m_findExpValues);
101  B2INFO("storeConst = " << m_storeConst);
102 
105  TH1F* dummy;
106  dummy = (TH1F*)gROOT->FindObject("IntegralVsCrysID");
107  if (dummy) {delete dummy;}
108  dummy = (TH1F*)gROOT->FindObject("AverageExpECrys");
109  if (dummy) {delete dummy;}
110  dummy = (TH1F*)gROOT->FindObject("AverageElecCalib");
111  if (dummy) {delete dummy;}
112  dummy = (TH1F*)gROOT->FindObject("AverageInitCalib");
113  if (dummy) {delete dummy;}
114 
117  auto EnVsCrysID = getObjectPtr<TH2F>("EnVsCrysID");
118  auto ExpEvsCrys = getObjectPtr<TH1F>("ExpEvsCrys");
119  auto ElecCalibvsCrys = getObjectPtr<TH1F>("ElecCalibvsCrys");
120  auto InitialCalibvsCrys = getObjectPtr<TH1F>("InitialCalibvsCrys");
121  auto CalibEntriesvsCrys = getObjectPtr<TH1F>("CalibEntriesvsCrys");
122 
126  TH1F* IntegralVsCrysID = new TH1F("IntegralVsCrysID", "Integral of EnVsCrysID for each crystal;crystal ID;Entries",
127  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
128  TH1F* AverageExpECrys = new TH1F("AverageExpECrys", "Average expected E per crys from collector;Crystal ID;Energy (GeV)",
129  ECLElementNumbers::c_NCrystals, 0,
130  ECLElementNumbers::c_NCrystals);
131  TH1F* AverageElecCalib = new TH1F("AverageElecCalib", "Average electronics calib const vs crystal;Crystal ID;Calibration constant",
132  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
133  TH1F* AverageInitCalib = new TH1F("AverageInitCalib", "Average initial calib const vs crystal;Crystal ID;Calibration constant",
134  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
135 
136  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
137  TH1D* hEnergy = EnVsCrysID->ProjectionY("hEnergy", crysID + 1, crysID + 1);
138  int Integral = hEnergy->Integral();
139  IntegralVsCrysID->SetBinContent(crysID + 1, Integral);
140 
141  double TotEntries = CalibEntriesvsCrys->GetBinContent(crysID + 1);
142 
143  double expectedE = 0.;
144  if (TotEntries > 0.) {expectedE = ExpEvsCrys->GetBinContent(crysID + 1) / TotEntries;}
145  AverageExpECrys->SetBinContent(crysID + 1, expectedE);
146 
147  double calibconst = 0.;
148  if (TotEntries > 0.) {calibconst = ElecCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
149  AverageElecCalib->SetBinContent(crysID + 1, calibconst);
150 
151  calibconst = 0.;
152  if (TotEntries > 0.) {calibconst = InitialCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
153  AverageInitCalib->SetBinContent(crysID + 1, calibconst);
154  }
155 
158  TString fName = m_outputName;
159  TFile* histfile = new TFile(fName, "recreate");
160  EnVsCrysID->Write();
161  IntegralVsCrysID->Write();
162  AverageExpECrys->Write();
163  AverageElecCalib->Write();
164  AverageInitCalib->Write();
165 
168  if (!m_performFits) {
169  B2INFO("eclGammaGammaEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
170  histfile->Close();
171  return c_NotEnoughData;
172  }
173 
176  bool sufficientData = true;
177  for (int crysID = m_cellIDLo - 1; crysID < m_cellIDHi; crysID++) {
178  if (IntegralVsCrysID->GetBinContent(crysID + 1) < m_minEntries) {
179  if (m_storeConst == 1) {B2INFO("eclGammaGammaEAlgorithm: crystal " << crysID << " has insufficient statistics: " << IntegralVsCrysID->GetBinContent(crysID + 1) << ". Requirement is " << m_minEntries);}
180  sufficientData = false;
181  break;
182  }
183  }
184 
186  if (!sufficientData && m_storeConst == 1) {
187  histfile->Close();
188  return c_NotEnoughData;
189  }
190 
195  TH1F* CalibVsCrysID = new TH1F("CalibVsCrysID", "Calibration constant vs crystal ID;crystal ID;counts per GeV",
196  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
197  TH1F* ExpEnergyperCrys = new TH1F("ExpEnergyperCrys", "Expected energy per crystal;Crystal ID;Peak energy (GeV)",
198  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
199 
201  TH1F* PeakVsCrysID = new TH1F("PeakVsCrysID", "Peak of Novo fit vs crystal ID;crystal ID;Peak normalized energy",
202  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
203  TH1F* EdgeVsCrysID = new TH1F("EdgeVsCrysID", "Upper edge of Novo fit vs crystal ID;crystal ID;Maximum normalized energy",
204  ECLElementNumbers::c_NCrystals, 0,
205  ECLElementNumbers::c_NCrystals);
206  TH1F* effSigVsCrysID = new TH1F("effSigVsCrysID", "effSigma vs crystal ID;crystal ID;sigma)", ECLElementNumbers::c_NCrystals, 0,
207  ECLElementNumbers::c_NCrystals);
208  TH1F* etaVsCrysID = new TH1F("etaVsCrysID", "eta vs crystal ID;crystal ID;Novo eta parameter", ECLElementNumbers::c_NCrystals, 0,
209  ECLElementNumbers::c_NCrystals);
210  TH1F* constVsCrysID = new TH1F("constVsCrysID", "fit constant vs crystal ID;crystal ID;fit constant",
211  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
212  TH1F* normVsCrysID = new TH1F("normVsCrysID", "Novosibirsk normalization vs crystal ID;crystal ID;normalization",
213  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
214  TH1F* fitLimitVsCrysID = new TH1F("fitLimitVsCrysID", "fit range lower limit vs crystal ID;crystal ID;upper fit limit",
215  ECLElementNumbers::c_NCrystals, 0,
216  ECLElementNumbers::c_NCrystals);
217  TH1F* StatusVsCrysID = new TH1F("StatusVsCrysID", "Fit status vs crystal ID;crystal ID;Fit status", ECLElementNumbers::c_NCrystals,
218  0, ECLElementNumbers::c_NCrystals);
219  TH1F* FitProbVsCrysID = new TH1F("FitProbVsCrysID", "Fit probability vs crystal id;crystal ID;Fit probability",
220  ECLElementNumbers::c_NCrystals, 0, ECLElementNumbers::c_NCrystals);
221 
223  TH1F* hStatus = new TH1F("hStatus", "Fit status", 25, -5, 20);
224  TH1F* hPeak = new TH1F("hPeak", "Peaks of normalized energy distributions, successful fits;Peak of Novosibirsk fit", 200, 0.8, 1.2);
225  TH1F* fracPeakUnc = new TH1F("fracPeakUnc", "Fractional uncertainty on peak uncertainty, successful fits;Uncertainty on peak", 100,
226  0, 0.1);
227  TH1F* nIterations = new TH1F("nIterations", "Number of times histogram was fit;Number of iterations", 20, -0.5, 19.5);
228 
229 
232  bool allFitsOK = true;
233  for (int crysID = m_cellIDLo - 1; crysID < m_cellIDHi; crysID++) {
234 
236  TString name = "Enormalized";
237  name += crysID;
238  TH1D* hEnergy = EnVsCrysID->ProjectionY(name, crysID + 1, crysID + 1);
239 
241  double histMin = hEnergy->GetXaxis()->GetXmin();
242  double histMax = hEnergy->GetXaxis()->GetXmax();
243  TF1* func = new TF1("eclGammaGammaNovoConst", eclGammaGammaNovoConst, histMin, histMax, 5);
244  func->SetParNames("normalization", "peak", "effSigma", "eta", "const");
245  func->SetParLimits(1, peakMin, peakMax);
246  func->SetParLimits(2, effSigMin, effSigMax);
247  func->SetParLimits(3, etaMin, etaMax);
248  func->SetParLimits(4, 0., hEnergy->GetMaximum());
249 
251  hEnergy->GetXaxis()->SetRangeUser(peakMin, peakMax);
252  int maxBin = hEnergy->GetMaximumBin();
253  double peakE = hEnergy->GetBinLowEdge(maxBin);
254  double peakEUnc = 0.;
255  double normalization = hEnergy->GetMaximum();
256  double normUnc = 0.;
257  double effSigma = hEnergy->GetRMS();
258  double sigmaUnc = 0.;
259  hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
260 
262  double fitlow = peakE - effSigma;
263  double fithigh = histMax;
264 
266  double eta = etaNom;
267  double etaUnc = 0.;
268  double constant = 0.01 * normalization;
269  double constUnc = 0.;
270 
272  double dIter = 0.1 * (histMax - histMin) / hEnergy->GetNbinsX();
273  double fitProb(0.);
274  double fitProbDefault(0.);
275  double lowold(0.), lowoldold(0.);
276  bool fixConst = false;
277  int nIter = 0;
278  double histIntegral = IntegralVsCrysID->GetBinContent(crysID + 1);
279  bool fitHist = histIntegral >= m_minEntries; /* fit only if enough events */
280 
282  double m_tRatioMin = m_tRatioMinNom;
283  double m_tRatioMax = m_tRatioMaxNom;
284  if (histIntegral > m_highStatEntries) {
285  m_tRatioMin = m_tRatioMinHiStat;
286  m_tRatioMax = m_tRatioMaxHiStat;
287  }
288 
291  while (fitHist) {
292 
294  func->SetParameters(normalization, peakE, effSigma, eta, constant);
295  if (fixConst) { func->FixParameter(4, 0); }
296 
298  hEnergy->Fit(func, "LIQ", "", fitlow, fithigh);
299  nIter++;
300  fitHist = false;
301  normalization = func->GetParameter(0);
302  normUnc = func->GetParError(0);
303  peakE = func->GetParameter(1);
304  peakEUnc = func->GetParError(1);
305  effSigma = func->GetParameter(2);
306  sigmaUnc = func->GetParError(2);
307  eta = func->GetParameter(3);
308  etaUnc = func->GetParError(3);
309  constant = func->GetParameter(4);
310  constUnc = func->GetParError(4);
311  fitProbDefault = func->GetProb();
312 
314  double peak = func->Eval(peakE) - constant;
315  double tRatio = (func->Eval(fitlow) - constant) / peak;
316  if (tRatio < m_tRatioMin || tRatio > m_tRatioMax) {
317  double targetY = constant + 0.5 * (m_tRatioMin + m_tRatioMax) * peak;
318  lowoldold = lowold;
319  lowold = fitlow;
320  fitlow = func->GetX(targetY, histMin, peakE);
321  fitHist = true;
322 
324  if (abs(fitlow - lowoldold) < dIter) {fitlow = 0.5 * (lowold + lowoldold); }
325 
327  if (nIter > m_maxIterations - 3) {fitlow = 0.33333 * (fitlow + lowold + lowoldold); }
328  }
329 
331  if (constant < constTol && !fixConst) {
332  constant = 0;
333  fixConst = true;
334  fitHist = true;
335  }
336 
338  if (nIter == m_maxIterations) {fitHist = false;}
339  B2DEBUG(10, crysID << " " << nIter << " " << peakE << " " << constant << " " << tRatio << " " << fitlow);
340  }
341 
344  fitProb = 0.;
345  if (nIter > 0) {
346  int lowbin = hEnergy->GetXaxis()->FindBin(fitlow);
347  int highbin = hEnergy->GetXaxis()->FindBin(fithigh);
348  int npar = 5;
349  if (fixConst) {npar = 4;}
350  int ndeg = -npar;
351  double chisq = 0.;
352  double binwidth = hEnergy->GetBinWidth(1);
353  for (int ib = lowbin; ib <= highbin; ib++) {
354  double xlow = hEnergy->GetBinLowEdge(ib);
355  double yexp = func->Integral(xlow, xlow + binwidth) / binwidth;
356 
358  if (yexp > constTol) {
359  double yobs = hEnergy->GetBinContent(ib);
360  double dnom = yexp;
361  if (yexp < 0.9999 && yobs > yexp) {dnom = yobs;}
362  double dchi2 = (yexp - yobs) * (yexp - yobs) / dnom;
363  chisq += dchi2;
364  ndeg++;
365  }
366  }
367  fitProb = TMath::Prob(chisq, ndeg);
368  }
369 
372  int iStatus = fitOK; // success
373  if (nIter == m_maxIterations) {iStatus = iterations;} // too many iterations
374 
376  if (normalization < constRatio * constant) {iStatus = noPeak;}
377 
379  if (fitProb <= minFitLimit || (fitProb < minFitProbIter && iStatus == iterations)) {iStatus = poorFit;}
380 
382  if ((peakE < peakMin + peakTol) || (peakE > peakMax - peakTol)) {iStatus = atLimit;}
383  if ((effSigma < effSigMin + effSigTol) || (effSigma > effSigMax - effSigTol)) {iStatus = atLimit;}
384  if ((eta < etaMin + etaTol) || (eta > etaMax - etaTol)) {iStatus = atLimit;}
385 
386  //** No fit
387  if (nIter == 0) {iStatus = notFit;} // not fit
388 
391  double upperEdge = peakE;
392  double edgeUnc = peakEUnc;
393 
394  if (iStatus >= iterations) {
395 
397  double targetY = constant + m_upperEdgeThresh * (func->Eval(peakE) - constant);
398 
400  int iLow = hEnergy->GetXaxis()->FindBin(peakE) + 1;
401  int iHigh = hEnergy->GetNbinsX();
402  int iLast = iLow;
403  for (int ibin = iLow; ibin < iHigh; ibin++) {
404  double xc = hEnergy->GetBinCenter(ibin);
405  if (func->Eval(xc) > targetY) {iLast = ibin;}
406  }
407  double xLow = hEnergy->GetBinCenter(iLast);
408  double xHigh = hEnergy->GetBinCenter(iLast + 1);
409 
411  func->SetNpx(1000);
412  upperEdge = func->GetX(targetY, xLow, xHigh);
413 
414 
415  } else if (iStatus > notFit) {
416 
418  int iLast = -1;
419  int thisBin = hEnergy->GetBinContent(1);
420  for (int ibin = 2; ibin < hEnergy->GetNbinsX(); ibin++) {
421  int prevBin = thisBin;
422  thisBin = hEnergy->GetBinContent(ibin);
423  if (thisBin > 0 && thisBin + prevBin >= 2) {iLast = ibin;}
424  }
425  if (iLast > 0) {upperEdge = hEnergy->GetBinLowEdge(iLast);} // lower edge is a better estimate than center
426  }
427 
430  int histbin = crysID + 1;
431  PeakVsCrysID->SetBinContent(histbin, peakE);
432  PeakVsCrysID->SetBinError(histbin, peakEUnc);
433  EdgeVsCrysID->SetBinContent(histbin, upperEdge);
434  EdgeVsCrysID->SetBinError(histbin, edgeUnc);
435  effSigVsCrysID->SetBinContent(histbin, effSigma);
436  effSigVsCrysID->SetBinError(histbin, sigmaUnc);
437  etaVsCrysID->SetBinContent(histbin, eta);
438  etaVsCrysID->SetBinError(histbin, etaUnc);
439  constVsCrysID->SetBinContent(histbin, constant);
440  constVsCrysID->SetBinError(histbin, constUnc);
441  normVsCrysID->SetBinContent(histbin, normalization);
442  normVsCrysID->SetBinError(histbin, normUnc);
443  fitLimitVsCrysID->SetBinContent(histbin, fitlow);
444  fitLimitVsCrysID->SetBinError(histbin, 0);
445  StatusVsCrysID->SetBinContent(histbin, iStatus);
446  StatusVsCrysID->SetBinError(histbin, 0);
447  FitProbVsCrysID->SetBinContent(histbin, fitProb);
448  FitProbVsCrysID->SetBinError(histbin, 0);
449 
451  hStatus->Fill(iStatus);
452  nIterations->Fill(nIter);
453  if (iStatus >= iterations) {
454  hPeak->Fill(peakE);
455  fracPeakUnc->Fill(peakEUnc / peakE);
456  }
457 
459  B2INFO("cellID " << crysID + 1 << " status = " << iStatus << " fit probability = " << fitProb << " default prob = " <<
460  fitProbDefault);
461  histfile->cd();
462  hEnergy->Write();
463 
464  } /* end of loop over crystals */
465 
468  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
469  int histbin = crysID + 1;
470  double fitstatus = StatusVsCrysID->GetBinContent(histbin);
471  double upperEdge = EdgeVsCrysID->GetBinContent(histbin);
472  double fracEdgeUnc = EdgeVsCrysID->GetBinError(histbin) / upperEdge;
473 
475  if (fitstatus < 0) {
476  upperEdge = -1.;
477  fracEdgeUnc = 0.;
478  if (histbin >= m_cellIDLo && histbin <= m_cellIDHi) {
479  B2INFO("eclGammaGammaEAlgorithm: cellID " << histbin << " is not a successful fit. Status = " << fitstatus);
480  allFitsOK = false;
481  }
482  }
483 
485  if (m_findExpValues) {
486  double inputExpE = abs(AverageExpECrys->GetBinContent(histbin));
487  ExpEnergyperCrys->SetBinContent(histbin, inputExpE * upperEdge);
488  ExpEnergyperCrys->SetBinError(histbin, fracEdgeUnc * inputExpE * upperEdge);
489  } else {
490 
492  double inputCalib = abs(AverageInitCalib->GetBinContent(histbin));
493  CalibVsCrysID->SetBinContent(histbin, inputCalib / upperEdge);
494  CalibVsCrysID->SetBinError(histbin, fracEdgeUnc * inputCalib / upperEdge);
495  }
496  }
497 
500  bool DBsuccess = false;
501  if (m_storeConst == 0 || (m_storeConst == 1 && allFitsOK)) {
502  DBsuccess = true;
503  if (m_findExpValues) {
504 
506  std::vector<float> tempE;
507  std::vector<float> tempUnc;
508  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
509  tempE.push_back(ExpEnergyperCrys->GetBinContent(crysID + 1));
510  tempUnc.push_back(ExpEnergyperCrys->GetBinError(crysID + 1));
511  }
512  ECLCrystalCalib* ExpectedE = new ECLCrystalCalib();
513  ExpectedE->setCalibVector(tempE, tempUnc);
514  saveCalibration(ExpectedE, "ECLExpGammaGammaE");
515  B2INFO("eclCosmicEAlgorithm: successfully stored expected energies ECLExpGammaGammaE");
516 
517  } else {
518 
520  std::vector<float> tempCalib;
521  std::vector<float> tempCalibUnc;
522  for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
523  tempCalib.push_back(CalibVsCrysID->GetBinContent(crysID + 1));
524  tempCalibUnc.push_back(CalibVsCrysID->GetBinError(crysID + 1));
525  }
526  ECLCrystalCalib* GammaGammaECalib = new ECLCrystalCalib();
527  GammaGammaECalib->setCalibVector(tempCalib, tempCalibUnc);
528  saveCalibration(GammaGammaECalib, "ECLCrystalEnergyGammaGamma");
529  B2INFO("eclGammaGammaEAlgorithm: successfully stored ECLCrystalEnergyGammaGamma calibration constants");
530  }
531  }
532 
536  PeakVsCrysID->Write();
537  EdgeVsCrysID->Write();
538  effSigVsCrysID->Write();
539  etaVsCrysID->Write();
540  constVsCrysID->Write();
541  normVsCrysID->Write();
542  fitLimitVsCrysID->Write();
543  StatusVsCrysID->Write();
544  FitProbVsCrysID->Write();
545  hPeak->Write();
546  fracPeakUnc->Write();
547  nIterations->Write();
548  hStatus->Write();
549 
551  if (m_findExpValues) {
552  ExpEnergyperCrys->Write();
553  } else {
554  CalibVsCrysID->Write();
555  }
556  histfile->Close();
557 
560  dummy = (TH1F*)gROOT->FindObject("PeakVsCrysID"); delete dummy;
561  dummy = (TH1F*)gROOT->FindObject("EdgeVsCrysID"); delete dummy;
562  dummy = (TH1F*)gROOT->FindObject("effSigVsCrysID"); delete dummy;
563  dummy = (TH1F*)gROOT->FindObject("etaVsCrysID"); delete dummy;
564  dummy = (TH1F*)gROOT->FindObject("constVsCrysID"); delete dummy;
565  dummy = (TH1F*)gROOT->FindObject("normVsCrysID"); delete dummy;
566  dummy = (TH1F*)gROOT->FindObject("fitLimitVsCrysID"); delete dummy;
567  dummy = (TH1F*)gROOT->FindObject("StatusVsCrysID"); delete dummy;
568  dummy = (TH1F*)gROOT->FindObject("FitProbVsCrysID"); delete dummy;
569  dummy = (TH1F*)gROOT->FindObject("fracPeakUnc"); delete dummy;
570  dummy = (TH1F*)gROOT->FindObject("nIterations"); delete dummy;
571  dummy = (TH1F*)gROOT->FindObject("hStatus"); delete dummy;
572  dummy = (TH1F*)gROOT->FindObject("ExpEnergyperCrys"); delete dummy;
573  dummy = (TH1F*)gROOT->FindObject("CalibVsCrysID"); delete dummy;
574 
575 
578  if (m_storeConst == -1) {
579  B2INFO("eclGammaGammaEAlgorithm performed fits but was not asked to store contants");
580  return c_Failure;
581  } else if (!DBsuccess) {
582  if (m_findExpValues) { B2INFO("eclGammaGammaEAlgorithm: failed to store expected values"); }
583  else { B2INFO("eclGammaGammaEAlgorithm: failed to store calibration constants"); }
584  return c_Failure;
585  }
586  return c_OK;
587 }
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::eclGammaGammaEAlgorithm::m_minEntries
int m_minEntries
Minimum entries to fit a crystal.
Definition: eclGammaGammaEAlgorithm.h:124
Belle2::ECL::eclGammaGammaEAlgorithm::poorFit
int poorFit
low chi square; upper edge is found from histogram, not fit
Definition: eclGammaGammaEAlgorithm.h:147
Belle2::ECL::eclGammaGammaEAlgorithm::m_maxIterations
int m_maxIterations
no more than maxIteration iterations
Definition: eclGammaGammaEAlgorithm.h:126
Belle2::ECL::eclGammaGammaEAlgorithm::m_tRatioMaxHiStat
double m_tRatioMaxHiStat
Fit range is adjusted so that fit at lower endpoint is between tRatioMin and tRatioMax of peak.
Definition: eclGammaGammaEAlgorithm.h:132
Belle2::ECL::eclGammaGammaEAlgorithm::iterations
int iterations
fit reached max number of iterations, but is useable
Definition: eclGammaGammaEAlgorithm.h:145
Belle2::ECL::eclGammaGammaEAlgorithm::m_tRatioMaxNom
double m_tRatioMaxNom
Fit range is adjusted so that fit at lower endpoint is between tRatioMin and tRatioMax of peak.
Definition: eclGammaGammaEAlgorithm.h:129
Belle2::ECL::eclGammaGammaEAlgorithm::m_tRatioMinNom
double m_tRatioMinNom
Fit range is adjusted so that fit at lower endpoint is between tRatioMin and tRatioMax of peak.
Definition: eclGammaGammaEAlgorithm.h:127
Belle2::ECL::eclGammaGammaEAlgorithm::m_highStatEntries
int m_highStatEntries
Adjust fit range above this many entries.
Definition: eclGammaGammaEAlgorithm.h:125
Belle2::ECL::eclGammaGammaEAlgorithm::m_outputName
std::string m_outputName
..Parameters to control Novosibirsk fit to energy deposited in each crystal by mu+mu- events
Definition: eclGammaGammaEAlgorithm.h:121
Belle2::ECL::eclGammaGammaEAlgorithm::m_storeConst
int m_storeConst
controls which values are written to the database.
Definition: eclGammaGammaEAlgorithm.h:138
Belle2::ECL::eclGammaGammaEAlgorithm::notFit
int notFit
no fit performed; no constants found for this crystal
Definition: eclGammaGammaEAlgorithm.h:149
Belle2::ECL::eclGammaGammaEAlgorithm::m_findExpValues
bool m_findExpValues
if true, fits are used to find expected energy deposit for each crystal instead of the calibration co...
Definition: eclGammaGammaEAlgorithm.h:136
Belle2::ECL::eclGammaGammaEAlgorithm::calibrate
virtual EResult calibrate() override
..Run algorithm on events
Definition: eclGammaGammaEAlgorithm.cc:65
Belle2::ECL::eclGammaGammaEAlgorithm::m_performFits
bool m_performFits
if false, input histograms are copied to output, but no fits are done
Definition: eclGammaGammaEAlgorithm.h:135
Belle2::ECL::eclGammaGammaEAlgorithm::m_upperEdgeThresh
double m_upperEdgeThresh
Upper edge is where the fit = upperEdgeThresh * peak value.
Definition: eclGammaGammaEAlgorithm.h:134
Belle2::ECL::eclGammaGammaEAlgorithm::noPeak
int noPeak
Novosibirsk component of fit is negligible; upper edge is found from histogram, not fit.
Definition: eclGammaGammaEAlgorithm.h:148
Belle2::ECL::eclGammaGammaEAlgorithm::m_tRatioMinHiStat
double m_tRatioMinHiStat
Fit range is adjusted so that fit at lower endpoint is between tRatioMin and tRatioMax of peak.
Definition: eclGammaGammaEAlgorithm.h:130
Belle2::ECL::eclGammaGammaEAlgorithm::atLimit
int atLimit
a parameter is at the limit; upper edge is found from histogram, not fit
Definition: eclGammaGammaEAlgorithm.h:146
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