9 #include <ecl/calibration/eclGammaGammaEAlgorithm.h>
10 #include <ecl/dbobjects/ECLCrystalCalib.h>
24 double eclGammaGammaNovoConst(
double* x,
double* par)
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;
34 if (TMath::Abs(tail) < 1.e-7) {
35 qc = 0.5 * TMath::Power(((x[0] - peak) / width), 2);
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;
43 qc = 0.5 * (TMath::Power((log(qy) / tail), 2) + tail * tail);
47 return par[0] * exp(-qc) + par[4];
55 "Perform energy calibration of ecl crystals by fitting a Novosibirsk function to energy deposited by photons in e+e- --> gamma gamma"
63 const double limitTol = 0.0005;
64 const double minFitLimit = 1e-25;
65 const double minFitProbIter = 1e-8;
66 const double constRatio = 0.5;
67 const double peakMin(0.4), peakMax(2.2);
68 const double peakTol = limitTol * (peakMax - peakMin);
69 const double effSigMin(0.02), effSigMax(0.4);
70 const double effSigTol = limitTol * (effSigMax - effSigMin);
71 const double etaNom(0.41);
72 const double etaMin(0.), etaMax(5.0);
73 const double etaTol = limitTol * (etaMax - etaMin);
74 const double constTol = 0.001;
81 B2INFO(
"eclGammaGammaAlgorithm parameters:");
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;}
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");
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,
123 TH1F* AverageElecCalib =
new TH1F(
"AverageElecCalib",
"Average electronics calib const vs crystal;Crystal ID;Calibration constant",
125 TH1F* AverageInitCalib =
new TH1F(
"AverageInitCalib",
"Average initial calib const vs crystal;Crystal ID;Calibration constant",
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);
133 double TotEntries = CalibEntriesvsCrys->GetBinContent(crysID + 1);
135 double expectedE = 0.;
136 if (TotEntries > 0.) {expectedE = ExpEvsCrys->GetBinContent(crysID + 1) / TotEntries;}
137 AverageExpECrys->SetBinContent(crysID + 1, expectedE);
139 double calibconst = 0.;
140 if (TotEntries > 0.) {calibconst = ElecCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
141 AverageElecCalib->SetBinContent(crysID + 1, calibconst);
144 if (TotEntries > 0.) {calibconst = InitialCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
145 AverageInitCalib->SetBinContent(crysID + 1, calibconst);
151 TFile* histfile =
new TFile(fName,
"recreate");
153 IntegralVsCrysID->Write();
154 AverageExpECrys->Write();
155 AverageElecCalib->Write();
156 AverageInitCalib->Write();
161 B2INFO(
"eclGammaGammaEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
168 bool sufficientData =
true;
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;
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);
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,
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,
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);
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,
208 TH1F* nIterations =
new TH1F(
"nIterations",
"Number of times histogram was fit;Number of iterations", 20, -0.5, 19.5);
213 bool allFitsOK =
true;
217 TString name =
"Enormalized";
219 TH1D* hEnergy = EnVsCrysID->ProjectionY(name, crysID + 1, crysID + 1);
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());
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();
238 double effSigma = hEnergy->GetRMS();
239 double sigmaUnc = 0.;
240 hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
243 double fitlow = peakE - effSigma;
244 double fithigh = histMax;
249 double constant = 0.01 * normalization;
250 double constUnc = 0.;
253 double dIter = 0.1 * (histMax - histMin) / hEnergy->GetNbinsX();
255 double fitProbDefault(0.);
256 double lowold(0.), lowoldold(0.);
257 bool fixConst =
false;
259 double histIntegral = IntegralVsCrysID->GetBinContent(crysID + 1);
275 func->SetParameters(normalization, peakE, effSigma, eta, constant);
276 if (fixConst) { func->FixParameter(4, 0); }
279 hEnergy->Fit(func,
"LIQ",
"", fitlow, fithigh);
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();
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;
301 fitlow = func->GetX(targetY, histMin, peakE);
305 if (abs(fitlow - lowoldold) < dIter) {fitlow = 0.5 * (lowold + lowoldold); }
308 if (nIter >
m_maxIterations - 3) {fitlow = 0.33333 * (fitlow + lowold + lowoldold); }
312 if (constant < constTol && !fixConst) {
320 B2DEBUG(10, crysID <<
" " << nIter <<
" " << peakE <<
" " << constant <<
" " << tRatio <<
" " << fitlow);
327 int lowbin = hEnergy->GetXaxis()->FindBin(fitlow);
328 int highbin = hEnergy->GetXaxis()->FindBin(fithigh);
330 if (fixConst) {npar = 4;}
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;
339 if (yexp > constTol) {
340 double yobs = hEnergy->GetBinContent(ib);
342 if (yexp < 0.9999 && yobs > yexp) {dnom = yobs;}
343 double dchi2 = (yexp - yobs) * (yexp - yobs) / dnom;
348 fitProb = TMath::Prob(chisq, ndeg);
357 if (normalization < constRatio * constant) {iStatus =
noPeak;}
360 if (fitProb <= minFitLimit || (fitProb < minFitProbIter && iStatus ==
iterations)) {iStatus =
poorFit;}
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;}
368 if (nIter == 0) {iStatus =
notFit;}
372 double upperEdge = peakE;
373 double edgeUnc = peakEUnc;
381 int iLow = hEnergy->GetXaxis()->FindBin(peakE) + 1;
382 int iHigh = hEnergy->GetNbinsX();
384 for (
int ibin = iLow; ibin < iHigh; ibin++) {
385 double xc = hEnergy->GetBinCenter(ibin);
386 if (func->Eval(xc) > targetY) {iLast = ibin;}
388 double xLow = hEnergy->GetBinCenter(iLast);
389 double xHigh = hEnergy->GetBinCenter(iLast + 1);
393 upperEdge = func->GetX(targetY, xLow, xHigh);
396 }
else if (iStatus >
notFit) {
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;}
406 if (iLast > 0) {upperEdge = hEnergy->GetBinLowEdge(iLast);}
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);
432 hStatus->Fill(iStatus);
433 nIterations->Fill(nIter);
436 fracPeakUnc->Fill(peakEUnc / peakE);
440 B2INFO(
"cellID " << crysID + 1 <<
" status = " << iStatus <<
" fit probability = " << fitProb <<
" default prob = " <<
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;
460 B2INFO(
"eclGammaGammaEAlgorithm: cellID " << histbin <<
" is not a successful fit. Status = " << fitstatus);
467 double inputExpE = abs(AverageExpECrys->GetBinContent(histbin));
468 ExpEnergyperCrys->SetBinContent(histbin, inputExpE * upperEdge);
469 ExpEnergyperCrys->SetBinError(histbin, fracEdgeUnc * inputExpE * upperEdge);
473 double inputCalib = abs(AverageInitCalib->GetBinContent(histbin));
474 CalibVsCrysID->SetBinContent(histbin, inputCalib / upperEdge);
475 CalibVsCrysID->SetBinError(histbin, fracEdgeUnc * inputCalib / upperEdge);
481 bool DBsuccess =
false;
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));
496 B2INFO(
"eclCosmicEAlgorithm: successfully stored expected energies ECLExpGammaGammaE");
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));
510 B2INFO(
"eclGammaGammaEAlgorithm: successfully stored ECLCrystalEnergyGammaGamma calibration constants");
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();
527 fracPeakUnc->Write();
528 nIterations->Write();
533 ExpEnergyperCrys->Write();
535 CalibVsCrysID->Write();
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;
560 B2INFO(
"eclGammaGammaEAlgorithm performed fits but was not asked to store contants");
562 }
else if (!DBsuccess) {
563 if (
m_findExpValues) { B2INFO(
"eclGammaGammaEAlgorithm: failed to store expected values"); }
564 else { B2INFO(
"eclGammaGammaEAlgorithm: failed to store calibration constants"); }
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.
int m_cellIDLo
First cellID to be fit.
int m_cellIDHi
Last cellID to be fit.
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...
int fitOK
Characterize fit status.
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.