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