Belle II Software development
eclMuMuEAlgorithm.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/eclMuMuEAlgorithm.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 eclNovoConst(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 double qa = tail * sqrt(log(4.));
43 double qb = sinh(qa) / qa;
44 double qx = (x[0] - peak) / width * qb;
45 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
55
57 cellIDHi(ECLElementNumbers::c_NCrystals), minEntries(150),
58 nToRebin(1000), tRatioMin(0.05), tRatioMax(0.40), lowerEdgeThresh(0.10), performFits(true), findExpValues(false), storeConst(0)
59{
61 "Perform energy calibration of ecl crystals by fitting a Novosibirsk function to energy deposited by muons"
62 );
63}
64
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 peakMin(0.4), peakMax(2.2); /*< range for peak of normalized energy distribution */
72 double peakTol = limitTol * (peakMax - peakMin); /*< fit is at limit if it is within peakTol of min or max */
73 double effSigMin(0.02), effSigMax(0.2); /*< range for effective sigma of normalized energy distribution */
74 double effSigTol = limitTol * (effSigMax - effSigMin); /*< fit is at limit if it is within effSigTol of min or max */
75 double etaNom(-0.41); /*< Nominal tail parameter; fixed to this value for low statistics fits */
76 double etaMin(-1.), etaMax(0.); /*< Novosibirsk tail parameter range */
77 double etaTol = limitTol * (etaMax - etaMin); /*< fit is at limit if it is within etaTol of min or max */
78
80 gROOT->SetBatch();
81
84 B2INFO("eclMuMuEAlgorithm parameters:");
85 B2INFO("cellIDLo = " << cellIDLo);
86 B2INFO("cellIDHi = " << cellIDHi);
87 B2INFO("minEntries = " << minEntries);
88 B2INFO("tRatioMin = " << tRatioMin);
89 B2INFO("tRatioMax = " << tRatioMax);
90 B2INFO("lowerEdgeThresh = " << lowerEdgeThresh);
91 B2INFO("performFits = " << performFits);
92 B2INFO("findExpValues = " << findExpValues);
93 B2INFO("storeConst = " << storeConst);
94
97 TH1F* dummy;
98 dummy = (TH1F*)gROOT->FindObject("IntegralVsCrysID");
99 if (dummy) {delete dummy;}
100 dummy = (TH1F*)gROOT->FindObject("AverageExpECrys");
101 if (dummy) {delete dummy;}
102 dummy = (TH1F*)gROOT->FindObject("AverageElecCalib");
103 if (dummy) {delete dummy;}
104 dummy = (TH1F*)gROOT->FindObject("AverageInitCalib");
105 if (dummy) {delete dummy;}
106
109 auto TrkPerCrysID = getObjectPtr<TH1F>("TrkPerCrysID");
110 auto EnVsCrysID = getObjectPtr<TH2F>("EnVsCrysID");
111 auto ExpEvsCrys = getObjectPtr<TH1F>("ExpEvsCrys");
112 auto ElecCalibvsCrys = getObjectPtr<TH1F>("ElecCalibvsCrys");
113 auto InitialCalibvsCrys = getObjectPtr<TH1F>("InitialCalibvsCrys");
114 auto CalibEntriesvsCrys = getObjectPtr<TH1F>("CalibEntriesvsCrys");
115 auto RawDigitAmpvsCrys = getObjectPtr<TH2F>("RawDigitAmpvsCrys");
116 auto RawDigitTimevsCrys = getObjectPtr<TH2F>("RawDigitTimevsCrys");
117 auto hitCrysVsExtrapolatedCrys = getObjectPtr<TH2F>("hitCrysVsExtrapolatedCrys");
118
122 TH1F* IntegralVsCrysID = new TH1F("IntegralVsCrysID", "Integral of EnVsCrysID for each crystal;crystal ID;Entries",
124 TH1F* AverageExpECrys = new TH1F("AverageExpECrys", "Average expected E per crys from collector;Crystal ID;Energy (GeV)",
127 TH1F* AverageElecCalib = new TH1F("AverageElecCalib", "Average electronics calib const vs crystal;Crystal ID;Calibration constant",
129 TH1F* AverageInitCalib = new TH1F("AverageInitCalib", "Average initial calib const vs crystal;Crystal ID;Calibration constant",
131
132 for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
133 TH1D* hEnergy = EnVsCrysID->ProjectionY("hEnergy", crysID + 1, crysID + 1);
134 int Integral = hEnergy->Integral();
135 IntegralVsCrysID->SetBinContent(crysID + 1, Integral);
136
137 double TotEntries = CalibEntriesvsCrys->GetBinContent(crysID + 1);
138
139 double expectedE = 0.;
140 if (TotEntries > 0.) {expectedE = ExpEvsCrys->GetBinContent(crysID + 1) / TotEntries;}
141 AverageExpECrys->SetBinContent(crysID + 1, expectedE);
142
143 double calibconst = 0.;
144 if (TotEntries > 0.) {calibconst = ElecCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
145 AverageElecCalib->SetBinContent(crysID + 1, calibconst);
146
147 calibconst = 0.;
148 if (TotEntries > 0.) {calibconst = InitialCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
149 AverageInitCalib->SetBinContent(crysID + 1, calibconst);
150 }
151
154 TFile* histfile = new TFile("eclMuMuEAlgorithm.root", "recreate");
155 TrkPerCrysID->Write();
156 EnVsCrysID->Write();
157 IntegralVsCrysID->Write();
158 AverageExpECrys->Write();
159 AverageElecCalib->Write();
160 AverageInitCalib->Write();
161 RawDigitAmpvsCrys->Write();
162 RawDigitTimevsCrys->Write();
163 hitCrysVsExtrapolatedCrys->Write();
164
167 if (!performFits) {
168 B2RESULT("eclMuMuEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
169 histfile->Close();
170 return c_NotEnoughData;
171 }
172
175 bool sufficientData = true;
176 for (int crysID = cellIDLo - 1; crysID < cellIDHi; crysID++) {
177 if (IntegralVsCrysID->GetBinContent(crysID + 1) < minEntries) {
178 if (storeConst == 1) {B2RESULT("eclMuMuEAlgorithm: crystal " << crysID << " has insufficient statistics: " << IntegralVsCrysID->GetBinContent(crysID + 1) << ". Requirement is " << minEntries);}
179 sufficientData = false;
180 break;
181 }
182 }
183
185 if (!sufficientData && storeConst == 1) {
186 histfile->Close();
187 return c_NotEnoughData;
188 }
189
194 TH1F* CalibVsCrysID = new TH1F("CalibVsCrysID", "Calibration constant vs crystal ID;crystal ID;counts per GeV",
196 TH1F* ExpEnergyperCrys = new TH1F("ExpEnergyperCrys", "Expected energy per crystal;Crystal ID;Peak energy (GeV)",
198
200 TH1F* PeakVsCrysID = new TH1F("PeakVsCrysID", "Peak of Novo fit vs crystal ID;crystal ID;Peak normalized energy",
203 TH1F* EdgeVsCrysID = new TH1F("EdgeVsCrysID", "Lower edge of Novo fit vs crystal ID;crystal ID", ECLElementNumbers::c_NCrystals, 0,
205 TH1F* effSigVsCrysID = new TH1F("effSigVsCrysID", "effSigma vs crystal ID;crystal ID;sigma)", ECLElementNumbers::c_NCrystals, 0,
207 TH1F* etaVsCrysID = new TH1F("etaVsCrysID", "eta vs crystal ID;crystal ID;Novo eta parameter", ECLElementNumbers::c_NCrystals, 0,
209 TH1F* normVsCrysID = new TH1F("normVsCrysID", "Novosibirsk normalization vs crystal ID;crystal ID;normalization",
211 TH1F* lowerLimitVsCrysID = new TH1F("lowerLimitVsCrysID", "fit range lower limit vs crystal ID;crystal ID;lower fit limit",
214 TH1F* fitLimitVsCrysID = new TH1F("fitLimitVsCrysID", "fit range upper limit vs crystal ID;crystal ID;upper fit limit",
217 TH1F* StatusVsCrysID = new TH1F("StatusVsCrysID", "Fit status vs crystal ID;crystal ID;Fit status", ECLElementNumbers::c_NCrystals,
219
221 TH1F* hStatus = new TH1F("hStatus", "Fit status", 25, -5, 20);
222 TH1F* hPeak = new TH1F("hPeak", "Peaks of normalized energy distributions, successful fits;Peak of Novosibirsk fit", 200, 0.8, 1.2);
223 TH1F* fracPeakUnc = new TH1F("fracPeakUnc", "Fractional uncertainty on peak uncertainty, successful fits;Uncertainty on peak", 100,
224 0, 0.1);
225
226
229 bool allFitsOK = true;
230 for (int crysID = cellIDLo - 1; crysID < cellIDHi; crysID++) {
231
233 TString name = "Enormalized";
234 name += crysID;
235 TH1D* hEnergy = EnVsCrysID->ProjectionY(name, crysID + 1, crysID + 1);
236
238 double histMin = hEnergy->GetXaxis()->GetXmin();
239 double histMax = hEnergy->GetXaxis()->GetXmax();
240 TF1* func = new TF1("eclNovoConst", eclNovoConst, histMin, histMax, 5);
241 func->SetParNames("normalization", "peak", "effSigma", "eta", "const");
242 func->SetParLimits(1, peakMin, peakMax);
243 func->SetParLimits(2, effSigMin, effSigMax);
244 func->SetParLimits(3, etaMin, etaMax);
245
246 //..Currently not using the constant term
247 double constant = 0.;
248 func->FixParameter(4, constant);
249
250 //..If low statistics, rebin, and fix eta
251 if (hEnergy->GetEntries() < nToRebin) {
252 hEnergy->Rebin(2);
253 func->FixParameter(3, etaNom);
254 }
255
257 hEnergy->GetXaxis()->SetRangeUser(peakMin, peakMax);
258 double peak = hEnergy->GetMaximum();
259 int maxBin = hEnergy->GetMaximumBin();
260 double peakE = hEnergy->GetBinLowEdge(maxBin);
261 double peakEUnc = 0.;
262 double normalization = hEnergy->GetMaximum();
263 double normUnc = 0.;
264 double effSigma = hEnergy->GetRMS();
265 double sigmaUnc = 0.;
266 hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
267 double fitProb = 0.;
268
270 double eta = etaNom;
271 double etaUnc = 0.;
272
273 //..Will find the lower edge of normalized energy at the end of the fit
274 double lowerEnEdge = 0.;
275
276 //..Fit range from set of bins with sufficient entries.
277 double targetY = tRatioMin * peak;
278 int iLow = maxBin;
279 while (hEnergy->GetBinContent(iLow) > targetY) {iLow--;}
280 double fitlow = hEnergy->GetBinLowEdge(iLow);
281
282 targetY = tRatioMax * peak;
283 int iHigh = maxBin;
284 while (hEnergy->GetBinContent(iHigh) > targetY) {iHigh++;}
285 double fithigh = hEnergy->GetBinLowEdge(iHigh + 1);
286
289 bool fitHist = IntegralVsCrysID->GetBinContent(crysID + 1) >= minEntries; /* fit only if enough events */
290 if (fitHist) {
291
293 func->SetParameters(normalization, peakE, effSigma, eta, constant);
294
296 hEnergy->Fit(func, "LIQ", "", fitlow, fithigh);
297 normalization = func->GetParameter(0);
298 normUnc = func->GetParError(0);
299 peakE = func->GetParameter(1);
300 peakEUnc = func->GetParError(1);
301 effSigma = func->GetParameter(2);
302 sigmaUnc = func->GetParError(2);
303 eta = func->GetParameter(3);
304 etaUnc = func->GetParError(3);
305 fitProb = func->GetProb();
306
307 //..Lower edge of the fit function is used to find the calibration constant.
308 // Can now use the peak of the fit, instead of the bin content.
309 peak = func->Eval(peakE);
310 targetY = lowerEdgeThresh * peak;
311
313 iHigh = hEnergy->GetXaxis()->FindBin(peakE) + 1;
314 iLow = hEnergy->GetXaxis()->FindBin(fitlow);
315 int iLast = iHigh;
316 for (int ibin = iHigh; ibin > iLow; ibin--) {
317 double xc = hEnergy->GetBinCenter(ibin);
318 if (func->Eval(xc) > targetY) {iLast = ibin;}
319 }
320 double xHigh = hEnergy->GetBinCenter(iLast);
321 double xLow = hEnergy->GetBinCenter(iLast - 1);
322
324 if (func->Eval(xLow) < targetY and func->Eval(xHigh) > targetY) {
325 func->SetNpx(1000);
326 lowerEnEdge = func->GetX(targetY, xLow, xHigh);
327 }
328 }
329
332 int iStatus = fitOK; // success
333
335 if (lowerEnEdge < 0.01) {iStatus = noLowerEdge;}
336
338 if (fitProb <= minFitLimit) {iStatus = poorFit;}
339
341 if ((peakE < peakMin + peakTol) || (peakE > peakMax - peakTol)) {iStatus = atLimit;}
342 if ((effSigma < effSigMin + effSigTol) || (effSigma > effSigMax - effSigTol)) {iStatus = atLimit;}
343 if ((eta < etaMin + etaTol) || (eta > etaMax - etaTol)) {iStatus = atLimit;}
344
345 //** No fit
346 if (!fitHist) {iStatus = notFit;} // not fit
347
349 int histbin = crysID + 1;
350 PeakVsCrysID->SetBinContent(histbin, peakE);
351 PeakVsCrysID->SetBinError(histbin, peakEUnc);
352 EdgeVsCrysID->SetBinContent(histbin, lowerEnEdge);
353 EdgeVsCrysID->SetBinError(histbin, peakEUnc); // approximate
354 effSigVsCrysID->SetBinContent(histbin, effSigma);
355 effSigVsCrysID->SetBinError(histbin, sigmaUnc);
356 etaVsCrysID->SetBinContent(histbin, eta);
357 etaVsCrysID->SetBinError(histbin, etaUnc);
358 normVsCrysID->SetBinContent(histbin, normalization);
359 normVsCrysID->SetBinError(histbin, normUnc);
360 lowerLimitVsCrysID->SetBinContent(histbin, fitlow);
361 lowerLimitVsCrysID->SetBinError(histbin, 0);
362 fitLimitVsCrysID->SetBinContent(histbin, fithigh);
363 fitLimitVsCrysID->SetBinError(histbin, 0);
364 StatusVsCrysID->SetBinContent(histbin, iStatus);
365 StatusVsCrysID->SetBinError(histbin, 0);
366
368 hStatus->Fill(iStatus);
369 if (iStatus >= iterations) {
370 hPeak->Fill(peakE);
371 fracPeakUnc->Fill(peakEUnc / peakE);
372 }
373
375 B2INFO("cellID " << crysID + 1 << " status = " << iStatus << " fit probability = " << fitProb);
376 histfile->cd();
377 hEnergy->Write();
378
379 } /* end of loop over crystals */
380
383 for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
384 int histbin = crysID + 1;
385 double fitstatus = StatusVsCrysID->GetBinContent(histbin);
386 double peakE = PeakVsCrysID->GetBinContent(histbin);
387 double edge = EdgeVsCrysID->GetBinContent(histbin);
388 double fracpeakEUnc = PeakVsCrysID->GetBinError(histbin) / peakE;
389
391 if (fitstatus < iterations) {
392 peakE = -1.;
393 edge = -1.;
394 fracpeakEUnc = 0.;
395 if (histbin >= cellIDLo && histbin <= cellIDHi) {
396 B2RESULT("eclMuMuEAlgorithm: cellID " << histbin << " is not a successful fit. Status = " << fitstatus);
397 allFitsOK = false;
398 }
399 }
400
402 if (findExpValues) {
403 double inputExpE = abs(AverageExpECrys->GetBinContent(histbin));
404 ExpEnergyperCrys->SetBinContent(histbin, inputExpE * edge);
405 ExpEnergyperCrys->SetBinError(histbin, fracpeakEUnc * inputExpE * peakE);
406 } else {
407
409 double inputCalib = abs(AverageInitCalib->GetBinContent(histbin));
410 CalibVsCrysID->SetBinContent(histbin, inputCalib / edge);
411 CalibVsCrysID->SetBinError(histbin, fracpeakEUnc * inputCalib / peakE);
412 }
413 }
414
417 bool DBsuccess = false;
418 if (storeConst == 0 || (storeConst == 1 && allFitsOK)) {
419 DBsuccess = true;
420 if (findExpValues) {
421
423 std::vector<float> tempE;
424 std::vector<float> tempUnc;
425 for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
426 tempE.push_back(ExpEnergyperCrys->GetBinContent(crysID + 1));
427 tempUnc.push_back(ExpEnergyperCrys->GetBinError(crysID + 1));
428 }
429 ECLCrystalCalib* ExpectedE = new ECLCrystalCalib();
430 ExpectedE->setCalibVector(tempE, tempUnc);
431 saveCalibration(ExpectedE, "ECLExpMuMuE");
432 B2RESULT("eclCosmicEAlgorithm: successfully stored expected energies ECLExpMuMuE");
433
434 } else {
435
437 std::vector<float> tempCalib;
438 std::vector<float> tempCalibUnc;
439 for (int crysID = 0; crysID < ECLElementNumbers::c_NCrystals; crysID++) {
440 tempCalib.push_back(CalibVsCrysID->GetBinContent(crysID + 1));
441 tempCalibUnc.push_back(CalibVsCrysID->GetBinError(crysID + 1));
442 }
443 ECLCrystalCalib* MuMuECalib = new ECLCrystalCalib();
444 MuMuECalib->setCalibVector(tempCalib, tempCalibUnc);
445 saveCalibration(MuMuECalib, "ECLCrystalEnergyMuMu");
446 B2RESULT("eclMuMuEAlgorithm: successfully stored ECLCrystalEnergyMuMu calibration constants");
447 }
448 }
449
453 PeakVsCrysID->Write();
454 EdgeVsCrysID->Write();
455 effSigVsCrysID->Write();
456 etaVsCrysID->Write();
457 normVsCrysID->Write();
458 lowerLimitVsCrysID->Write();
459 fitLimitVsCrysID->Write();
460 StatusVsCrysID->Write();
461 hPeak->Write();
462 fracPeakUnc->Write();
463 hStatus->Write();
464
466 if (findExpValues) {
467 ExpEnergyperCrys->Write();
468 } else {
469 CalibVsCrysID->Write();
470 }
471 histfile->Close();
472
475 dummy = (TH1F*)gROOT->FindObject("PeakVsCrysID"); delete dummy;
476 dummy = (TH1F*)gROOT->FindObject("EdgeVsCrysID"); delete dummy;
477 dummy = (TH1F*)gROOT->FindObject("effSigVsCrysID"); delete dummy;
478 dummy = (TH1F*)gROOT->FindObject("etaVsCrysID"); delete dummy;
479 dummy = (TH1F*)gROOT->FindObject("normVsCrysID"); delete dummy;
480 dummy = (TH1F*)gROOT->FindObject("lowerLimitVsCrysID"); delete dummy;
481 dummy = (TH1F*)gROOT->FindObject("fitLimitVsCrysID"); delete dummy;
482 dummy = (TH1F*)gROOT->FindObject("StatusVsCrysID"); delete dummy;
483 dummy = (TH1F*)gROOT->FindObject("fitProbSame"); delete dummy;
484 dummy = (TH1F*)gROOT->FindObject("fracPeakUnc"); delete dummy;
485 dummy = (TH1F*)gROOT->FindObject("hStatus"); delete dummy;
486 dummy = (TH1F*)gROOT->FindObject("ExpEnergyperCrys"); delete dummy;
487 dummy = (TH1F*)gROOT->FindObject("CalibVsCrysID"); delete dummy;
488
489
492 if (storeConst == -1) {
493 B2RESULT("eclMuMuEAlgorithm performed fits but was not asked to store constants");
494 return c_Failure;
495 } else if (!DBsuccess) {
496 if (findExpValues) { B2RESULT("eclMuMuEAlgorithm: failed to store expected values"); }
497 else { B2RESULT("eclMuMuEAlgorithm: failed to store calibration constants"); }
498 return c_Failure;
499 }
500 return c_OK;
501}
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 successfully =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.
double tRatioMin
entries/peak at low edge of fit must be greater than this
int poorFit
low chi square; fit not usable
int iterations
fit reached max number of iterations, but is usable
int noLowerEdge
could not determine lower edge of fit
int storeConst
controls which values are written to the database.
bool performFits
if false, input histograms are copied to output, but no fits are done.
int cellIDHi
Last cellID to be fit.
double lowerEdgeThresh
Lower edge is where the fit = lowerEdgeThresh * peak value.
int cellIDLo
..Parameters to control Novosibirsk fit to energy deposited in each crystal by mu+mu- events
bool findExpValues
if true, fits are used to find expected energy deposit for each crystal instead of the calibration co...
int minEntries
All crystals to be fit must have at least minEntries events in the fit range.
int nToRebin
If fewer entries than this, rebin and fix eta parameter.
virtual EResult calibrate() override
..Run algorithm on events
int atLimit
a parameter is at the limit; fit not usable
double tRatioMax
entries/peak at high edge of fit must be greater than this
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
const int c_NCrystals
Number of crystals.
Abstract base class for different kinds of events.
STL namespace.