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