2 #include <ecl/calibration/eclMuMuEAlgorithm.h>
3 #include <ecl/dbobjects/ECLCrystalCalib.h>
16 double eclNovoConst(
double* x,
double* par)
21 double width = par[2];
22 double sln4 = sqrt(log(4));
23 double y = par[3] * sln4;
24 double tail = -log(y + sqrt(1 + y * y)) / sln4;
26 if (TMath::Abs(tail) < 1.e-7) {
27 qc = 0.5 * TMath::Power(((x[0] - peak) / width), 2);
29 double qa = tail * sqrt(log(4.));
30 double qb = sinh(qa) / qa;
31 double qx = (x[0] - peak) / width * qb;
32 double qy = 1. + tail * qx;
35 qc = 0.5 * (TMath::Power((log(qy) / tail), 2) + tail * tail);
39 return par[0] * exp(-qc) + par[4];
43 eclMuMuEAlgorithm::eclMuMuEAlgorithm():
CalibrationAlgorithm(
"eclMuMuECollector"), cellIDLo(1), cellIDHi(8736), minEntries(150),
44 maxIterations(10), tRatioMin(0.2), tRatioMax(0.25), performFits(true), findExpValues(false), storeConst(0)
47 "Perform energy calibration of ecl crystals by fitting a Novosibirsk function to energy deposited by muons"
55 double limitTol = 0.0005;
56 double minFitLimit = 1e-25;
57 double minFitProbIter = 1e-8;
58 double constRatio = 0.5;
59 double peakMin(0.4), peakMax(2.2);
60 double peakTol = limitTol * (peakMax - peakMin);
61 double effSigMin(0.02), effSigMax(0.2);
62 double effSigTol = limitTol * (effSigMax - effSigMin);
64 double etaMin(-0.7), etaMax(-0.2);
65 double etaTol = limitTol * (etaMax - etaMin);
66 double constMin(0.), constMax(10.);
67 double constTol = limitTol * constMax;
75 dummy = (TH1F*)gROOT->FindObject(
"IntegralVsCrysID");
76 if (dummy) {
delete dummy;}
77 dummy = (TH1F*)gROOT->FindObject(
"AverageExpECrys");
78 if (dummy) {
delete dummy;}
79 dummy = (TH1F*)gROOT->FindObject(
"AverageElecCalib");
80 if (dummy) {
delete dummy;}
81 dummy = (TH1F*)gROOT->FindObject(
"AverageInitCalib");
82 if (dummy) {
delete dummy;}
86 auto TrkPerCrysID = getObjectPtr<TH1F>(
"TrkPerCrysID");
87 auto EnVsCrysID = getObjectPtr<TH2F>(
"EnVsCrysID");
88 auto ExpEvsCrys = getObjectPtr<TH1F>(
"ExpEvsCrys");
89 auto ElecCalibvsCrys = getObjectPtr<TH1F>(
"ElecCalibvsCrys");
90 auto InitialCalibvsCrys = getObjectPtr<TH1F>(
"InitialCalibvsCrys");
91 auto CalibEntriesvsCrys = getObjectPtr<TH1F>(
"CalibEntriesvsCrys");
92 auto RawDigitAmpvsCrys = getObjectPtr<TH2F>(
"RawDigitAmpvsCrys");
93 auto RawDigitTimevsCrys = getObjectPtr<TH2F>(
"RawDigitTimevsCrys");
94 auto hitCrysVsExtrapolatedCrys = getObjectPtr<TH2F>(
"hitCrysVsExtrapolatedCrys");
99 TH1F* IntegralVsCrysID =
new TH1F(
"IntegralVsCrysID",
"Integral of EnVsCrysID for each crystal;crystal ID;Entries", 8736, 0, 8736);
100 TH1F* AverageExpECrys =
new TH1F(
"AverageExpECrys",
"Average expected E per crys from collector;Crystal ID;Energy (GeV)", 8736, 0,
102 TH1F* AverageElecCalib =
new TH1F(
"AverageElecCalib",
"Average electronics calib const vs crystal;Crystal ID;Calibration constant",
104 TH1F* AverageInitCalib =
new TH1F(
"AverageInitCalib",
"Average initial calib const vs crystal;Crystal ID;Calibration constant",
107 for (
int crysID = 0; crysID < 8736; crysID++) {
108 TH1D* hEnergy = EnVsCrysID->ProjectionY(
"hEnergy", crysID + 1, crysID + 1);
109 int Integral = hEnergy->Integral();
110 IntegralVsCrysID->SetBinContent(crysID + 1, Integral);
112 double TotEntries = CalibEntriesvsCrys->GetBinContent(crysID + 1);
114 double expectedE = 0.;
115 if (TotEntries > 0.) {expectedE = ExpEvsCrys->GetBinContent(crysID + 1) / TotEntries;}
116 AverageExpECrys->SetBinContent(crysID + 1, expectedE);
118 double calibconst = 0.;
119 if (TotEntries > 0.) {calibconst = ElecCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
120 AverageElecCalib->SetBinContent(crysID + 1, calibconst);
123 if (TotEntries > 0.) {calibconst = InitialCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
124 AverageInitCalib->SetBinContent(crysID + 1, calibconst);
129 TFile* histfile =
new TFile(
"eclMuMuEAlgorithm.root",
"recreate");
130 TrkPerCrysID->Write();
132 IntegralVsCrysID->Write();
133 AverageExpECrys->Write();
134 AverageElecCalib->Write();
135 AverageInitCalib->Write();
136 RawDigitAmpvsCrys->Write();
137 RawDigitTimevsCrys->Write();
138 hitCrysVsExtrapolatedCrys->Write();
143 B2RESULT(
"eclMuMuEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
150 bool sufficientData =
true;
152 if (IntegralVsCrysID->GetBinContent(crysID + 1) <
minEntries) {
153 if (
storeConst == 1) {B2RESULT(
"eclMuMuEAlgorithm: crystal " << crysID <<
" has insufficient statistics: " << IntegralVsCrysID->GetBinContent(crysID + 1) <<
". Requirement is " <<
minEntries);}
154 sufficientData =
false;
169 TH1F* CalibVsCrysID =
new TH1F(
"CalibVsCrysID",
"Calibration constant vs crystal ID;crystal ID;counts per GeV", 8736, 0, 8736);
170 TH1F* ExpEnergyperCrys =
new TH1F(
"ExpEnergyperCrys",
"Expected energy per crystal;Crystal ID;Peak energy (GeV)", 8736, 0, 8736);
173 TH1F* PeakVsCrysID =
new TH1F(
"PeakVsCrysID",
"Peak of Novo fit vs crystal ID;crystal ID;Peak normalized energy", 8736, 0,
175 TH1F* effSigVsCrysID =
new TH1F(
"effSigVsCrysID",
"effSigma vs crystal ID;crystal ID;sigma)", 8736, 0, 8736);
176 TH1F* etaVsCrysID =
new TH1F(
"etaVsCrysID",
"eta vs crystal ID;crystal ID;Novo eta parameter", 8736, 0, 8736);
177 TH1F* constVsCrysID =
new TH1F(
"constVsCrysID",
"fit constant vs crystal ID;crystal ID;fit constant", 8736, 0, 8736);
178 TH1F* normVsCrysID =
new TH1F(
"normVsCrysID",
"Novosibirsk normalization vs crystal ID;crystal ID;normalization", 8736, 0, 8736);
179 TH1F* fitLimitVsCrysID =
new TH1F(
"fitLimitVsCrysID",
"fit range upper limit vs crystal ID;crystal ID;upper fit limit", 8736, 0,
181 TH1F* StatusVsCrysID =
new TH1F(
"StatusVsCrysID",
"Fit status vs crystal ID;crystal ID;Fit status", 8736, 0, 8736);
184 TH1F* hStatus =
new TH1F(
"hStatus",
"Fit status", 25, -5, 20);
185 TH1F* hPeak =
new TH1F(
"hPeak",
"Peaks of normalized energy distributions, successful fits;Peak of Novosibirsk fit", 200, 0.8, 1.2);
186 TH1F* fracPeakUnc =
new TH1F(
"fracPeakUnc",
"Fractional uncertainty on peak uncertainty, successful fits;Uncertainty on peak", 100,
188 TH1F* nIterations =
new TH1F(
"nIterations",
"Number of times histogram was fit;Number of iterations", 20, -0.5, 19.5);
193 bool allFitsOK =
true;
197 TString name =
"Enormalized";
199 TH1D* hEnergy = EnVsCrysID->ProjectionY(name, crysID + 1, crysID + 1);
202 double histMin = hEnergy->GetXaxis()->GetXmin();
203 double histMax = hEnergy->GetXaxis()->GetXmax();
204 TF1* func =
new TF1(
"eclNovoConst", eclNovoConst, histMin, histMax, 5);
205 func->SetParNames(
"normalization",
"peak",
"effSigma",
"eta",
"const");
206 func->SetParLimits(1, peakMin, peakMax);
207 func->SetParLimits(2, effSigMin, effSigMax);
208 func->SetParLimits(3, etaMin, etaMax);
209 func->SetParLimits(4, constMin, constMax);
212 hEnergy->GetXaxis()->SetRangeUser(peakMin, peakMax);
213 int maxBin = hEnergy->GetMaximumBin();
214 double peakE = hEnergy->GetBinLowEdge(maxBin);
215 double peakEUnc = 0.;
216 double normalization = hEnergy->GetMaximum();
218 double effSigma = hEnergy->GetRMS();
219 double sigmaUnc = 0.;
220 hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
223 double fitlow = 0.25;
224 double fithigh = peakE + 2.5 * effSigma;
231 int il0 = hEnergy->GetXaxis()->FindBin(fitlow);
232 int il1 = hEnergy->GetXaxis()->FindBin(fitlow + 0.1);
233 double constant = hEnergy->Integral(il0, il1) / (1 + il1 - il0);
234 if (constant < 0.01 * normalization) {constant = 0.01 * normalization;}
235 double constUnc = 0.;
238 double dIter = 0.1 * (histMax - histMin) / hEnergy->GetNbinsX();
240 double highold(0.), higholdold(0.);
241 bool fixConst =
false;
243 bool fitHist = IntegralVsCrysID->GetBinContent(crysID + 1) >=
minEntries;
250 func->SetParameters(normalization, peakE, effSigma, eta, constant);
251 if (fixConst) { func->FixParameter(4, 0); }
254 if (crysID < 1152 || crysID > 7775 || !
findExpValues) {func->FixParameter(3, etaNom);}
257 hEnergy->Fit(func,
"LIQ",
"", fitlow, fithigh);
260 normalization = func->GetParameter(0);
261 normUnc = func->GetParError(0);
262 peakE = func->GetParameter(1);
263 peakEUnc = func->GetParError(1);
264 effSigma = func->GetParameter(2);
265 sigmaUnc = func->GetParError(2);
266 eta = func->GetParameter(3);
267 etaUnc = func->GetParError(3);
268 constant = func->GetParameter(4);
269 constUnc = func->GetParError(4);
270 fitProb = func->GetProb();
273 double peak = func->Eval(peakE) - constant;
274 double tRatio = (func->Eval(fithigh) - constant) / peak;
275 if (tRatio < tRatioMin || tRatio >
tRatioMax) {
277 higholdold = highold;
279 fithigh = func->GetX(targetY, peakE, histMax);
283 if (abs(fithigh - higholdold) < dIter) {fithigh = 0.5 * (highold + higholdold); }
286 if (nIter >
maxIterations - 3) {fithigh = 0.33333 * (fithigh + highold + higholdold); }
290 if (constant < constTol && !fixConst) {
298 B2DEBUG(200, crysID <<
" " << nIter <<
" " << peakE <<
" " << constant <<
" " << tRatio <<
" " << fithigh);
307 if (normalization < constRatio * constant) {iStatus =
noPeak;}
310 if (fitProb <= minFitLimit || (fitProb < minFitProbIter && iStatus ==
iterations)) {iStatus =
poorFit;}
313 if ((peakE < peakMin + peakTol) || (peakE > peakMax - peakTol)) {iStatus =
atLimit;}
314 if ((effSigma < effSigMin + effSigTol) || (effSigma > effSigMax - effSigTol)) {iStatus =
atLimit;}
315 if ((eta < etaMin + etaTol) || (eta > etaMax - etaTol)) {iStatus =
atLimit;}
316 if (constant > constMax - constTol) {iStatus =
atLimit;}
319 if (nIter == 0) {iStatus =
notFit;}
322 int histbin = crysID + 1;
323 PeakVsCrysID->SetBinContent(histbin, peakE);
324 PeakVsCrysID->SetBinError(histbin, peakEUnc);
325 effSigVsCrysID->SetBinContent(histbin, effSigma);
326 effSigVsCrysID->SetBinError(histbin, sigmaUnc);
327 etaVsCrysID->SetBinContent(histbin, eta);
328 etaVsCrysID->SetBinError(histbin, etaUnc);
329 constVsCrysID->SetBinContent(histbin, constant);
330 constVsCrysID->SetBinError(histbin, constUnc);
331 normVsCrysID->SetBinContent(histbin, normalization);
332 normVsCrysID->SetBinError(histbin, normUnc);
333 fitLimitVsCrysID->SetBinContent(histbin, fithigh);
334 fitLimitVsCrysID->SetBinError(histbin, 0);
335 StatusVsCrysID->SetBinContent(histbin, iStatus);
336 StatusVsCrysID->SetBinError(histbin, 0);
339 hStatus->Fill(iStatus);
340 nIterations->Fill(nIter);
343 fracPeakUnc->Fill(peakEUnc / peakE);
347 B2INFO(
"cellID " << crysID + 1 <<
" status = " << iStatus <<
" fit probability = " << fitProb);
355 for (
int crysID = 0; crysID < 8736; crysID++) {
356 int histbin = crysID + 1;
357 double fitstatus = StatusVsCrysID->GetBinContent(histbin);
358 double peakE = PeakVsCrysID->GetBinContent(histbin);
359 double fracpeakEUnc = PeakVsCrysID->GetBinError(histbin) / peakE;
366 B2RESULT(
"eclMuMuEAlgorithm: cellID " << histbin <<
" is not a successful fit. Status = " << fitstatus);
373 double inputExpE = abs(AverageExpECrys->GetBinContent(histbin));
374 ExpEnergyperCrys->SetBinContent(histbin, inputExpE * peakE);
375 ExpEnergyperCrys->SetBinError(histbin, fracpeakEUnc * inputExpE * peakE);
379 double inputCalib = abs(AverageInitCalib->GetBinContent(histbin));
380 CalibVsCrysID->SetBinContent(histbin, inputCalib / peakE);
381 CalibVsCrysID->SetBinError(histbin, fracpeakEUnc * inputCalib / peakE);
387 bool DBsuccess =
false;
393 std::vector<float> tempE;
394 std::vector<float> tempUnc;
395 for (
int crysID = 0; crysID < 8736; crysID++) {
396 tempE.push_back(ExpEnergyperCrys->GetBinContent(crysID + 1));
397 tempUnc.push_back(ExpEnergyperCrys->GetBinError(crysID + 1));
402 B2RESULT(
"eclCosmicEAlgorithm: successfully stored expected energies ECLExpMuMuE");
407 std::vector<float> tempCalib;
408 std::vector<float> tempCalibUnc;
409 for (
int crysID = 0; crysID < 8736; crysID++) {
410 tempCalib.push_back(CalibVsCrysID->GetBinContent(crysID + 1));
411 tempCalibUnc.push_back(CalibVsCrysID->GetBinError(crysID + 1));
416 B2RESULT(
"eclMuMuEAlgorithm: successfully stored ECLCrystalEnergyMuMu calibration constants");
423 PeakVsCrysID->Write();
424 effSigVsCrysID->Write();
425 etaVsCrysID->Write();
426 constVsCrysID->Write();
427 normVsCrysID->Write();
428 fitLimitVsCrysID->Write();
429 StatusVsCrysID->Write();
431 fracPeakUnc->Write();
432 nIterations->Write();
437 ExpEnergyperCrys->Write();
439 CalibVsCrysID->Write();
445 dummy = (TH1F*)gROOT->FindObject(
"PeakVsCrysID");
delete dummy;
446 dummy = (TH1F*)gROOT->FindObject(
"effSigVsCrysID");
delete dummy;
447 dummy = (TH1F*)gROOT->FindObject(
"etaVsCrysID");
delete dummy;
448 dummy = (TH1F*)gROOT->FindObject(
"constVsCrysID");
delete dummy;
449 dummy = (TH1F*)gROOT->FindObject(
"normVsCrysID");
delete dummy;
450 dummy = (TH1F*)gROOT->FindObject(
"fitLimitVsCrysID");
delete dummy;
451 dummy = (TH1F*)gROOT->FindObject(
"StatusVsCrysID");
delete dummy;
452 dummy = (TH1F*)gROOT->FindObject(
"fitProbSame");
delete dummy;
453 dummy = (TH1F*)gROOT->FindObject(
"fracPeakUnc");
delete dummy;
454 dummy = (TH1F*)gROOT->FindObject(
"nIterations");
delete dummy;
455 dummy = (TH1F*)gROOT->FindObject(
"hStatus");
delete dummy;
456 dummy = (TH1F*)gROOT->FindObject(
"ExpEnergyperCrys");
delete dummy;
457 dummy = (TH1F*)gROOT->FindObject(
"CalibVsCrysID");
delete dummy;
463 B2RESULT(
"eclMuMuEAlgorithm performed fits but was not asked to store contants");
465 }
else if (!DBsuccess) {
466 if (
findExpValues) { B2RESULT(
"eclMuMuEAlgorithm: failed to store expected values"); }
467 else { B2RESULT(
"eclMuMuEAlgorithm: failed to store calibration constants"); }