10#include <ecl/calibration/eclMuMuEAlgorithm.h>
13#include <ecl/dataobjects/ECLElementNumbers.h>
14#include <ecl/dbobjects/ECLCrystalCalib.h>
29double eclNovoConst(
double* x,
double* par)
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;
39 if (TMath::Abs(tail) < 1.e-7) {
40 qc = 0.5 * TMath::Power(((x[0] - peak) / width), 2);
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;
48 qc = 0.5 * (TMath::Power((log(qy) / tail), 2) + tail * tail);
52 return par[0] * exp(-qc) + par[4];
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)
61 "Perform energy calibration of ecl crystals by fitting a Novosibirsk function to energy deposited by muons"
69 double limitTol = 0.0005;
70 double minFitLimit = 1e-25;
71 double peakMin(0.4), peakMax(2.2);
72 double peakTol = limitTol * (peakMax - peakMin);
73 double effSigMin(0.02), effSigMax(0.2);
74 double effSigTol = limitTol * (effSigMax - effSigMin);
76 double etaMin(-1.), etaMax(0.);
77 double etaTol = limitTol * (etaMax - etaMin);
84 B2INFO(
"eclMuMuEAlgorithm parameters:");
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;}
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");
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",
133 TH1D* hEnergy = EnVsCrysID->ProjectionY(
"hEnergy", crysID + 1, crysID + 1);
134 int Integral = hEnergy->Integral();
135 IntegralVsCrysID->SetBinContent(crysID + 1, Integral);
137 double TotEntries = CalibEntriesvsCrys->GetBinContent(crysID + 1);
139 double expectedE = 0.;
140 if (TotEntries > 0.) {expectedE = ExpEvsCrys->GetBinContent(crysID + 1) / TotEntries;}
141 AverageExpECrys->SetBinContent(crysID + 1, expectedE);
143 double calibconst = 0.;
144 if (TotEntries > 0.) {calibconst = ElecCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
145 AverageElecCalib->SetBinContent(crysID + 1, calibconst);
148 if (TotEntries > 0.) {calibconst = InitialCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
149 AverageInitCalib->SetBinContent(crysID + 1, calibconst);
154 TFile* histfile =
new TFile(
"eclMuMuEAlgorithm.root",
"recreate");
155 TrkPerCrysID->Write();
157 IntegralVsCrysID->Write();
158 AverageExpECrys->Write();
159 AverageElecCalib->Write();
160 AverageInitCalib->Write();
161 RawDigitAmpvsCrys->Write();
162 RawDigitTimevsCrys->Write();
163 hitCrysVsExtrapolatedCrys->Write();
168 B2RESULT(
"eclMuMuEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
175 bool sufficientData =
true;
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;
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)",
200 TH1F* PeakVsCrysID =
new TH1F(
"PeakVsCrysID",
"Peak of Novo fit vs crystal ID;crystal ID;Peak normalized energy",
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",
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,
229 bool allFitsOK =
true;
233 TString name =
"Enormalized";
235 TH1D* hEnergy = EnVsCrysID->ProjectionY(name, crysID + 1, crysID + 1);
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);
247 double constant = 0.;
248 func->FixParameter(4, constant);
251 if (hEnergy->GetEntries() <
nToRebin) {
253 func->FixParameter(3, etaNom);
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();
264 double effSigma = hEnergy->GetRMS();
265 double sigmaUnc = 0.;
266 hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
274 double lowerEnEdge = 0.;
279 while (hEnergy->GetBinContent(iLow) > targetY) {iLow--;}
280 double fitlow = hEnergy->GetBinLowEdge(iLow);
284 while (hEnergy->GetBinContent(iHigh) > targetY) {iHigh++;}
285 double fithigh = hEnergy->GetBinLowEdge(iHigh + 1);
289 bool fitHist = IntegralVsCrysID->GetBinContent(crysID + 1) >=
minEntries;
293 func->SetParameters(normalization, peakE, effSigma, eta, constant);
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();
309 peak = func->Eval(peakE);
313 iHigh = hEnergy->GetXaxis()->FindBin(peakE) + 1;
314 iLow = hEnergy->GetXaxis()->FindBin(fitlow);
316 for (
int ibin = iHigh; ibin > iLow; ibin--) {
317 double xc = hEnergy->GetBinCenter(ibin);
318 if (func->Eval(xc) > targetY) {iLast = ibin;}
320 double xHigh = hEnergy->GetBinCenter(iLast);
321 double xLow = hEnergy->GetBinCenter(iLast - 1);
324 if (func->Eval(xLow) < targetY and func->Eval(xHigh) > targetY) {
326 lowerEnEdge = func->GetX(targetY, xLow, xHigh);
338 if (fitProb <= minFitLimit) {iStatus =
poorFit;}
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;}
346 if (!fitHist) {iStatus =
notFit;}
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);
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);
368 hStatus->Fill(iStatus);
371 fracPeakUnc->Fill(peakEUnc / peakE);
375 B2INFO(
"cellID " << crysID + 1 <<
" status = " << iStatus <<
" fit probability = " << fitProb);
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;
396 B2RESULT(
"eclMuMuEAlgorithm: cellID " << histbin <<
" is not a successful fit. Status = " << fitstatus);
403 double inputExpE = abs(AverageExpECrys->GetBinContent(histbin));
404 ExpEnergyperCrys->SetBinContent(histbin, inputExpE * edge);
405 ExpEnergyperCrys->SetBinError(histbin, fracpeakEUnc * inputExpE * peakE);
409 double inputCalib = abs(AverageInitCalib->GetBinContent(histbin));
410 CalibVsCrysID->SetBinContent(histbin, inputCalib / edge);
411 CalibVsCrysID->SetBinError(histbin, fracpeakEUnc * inputCalib / peakE);
417 bool DBsuccess =
false;
423 std::vector<float> tempE;
424 std::vector<float> tempUnc;
426 tempE.push_back(ExpEnergyperCrys->GetBinContent(crysID + 1));
427 tempUnc.push_back(ExpEnergyperCrys->GetBinError(crysID + 1));
432 B2RESULT(
"eclCosmicEAlgorithm: successfully stored expected energies ECLExpMuMuE");
437 std::vector<float> tempCalib;
438 std::vector<float> tempCalibUnc;
440 tempCalib.push_back(CalibVsCrysID->GetBinContent(crysID + 1));
441 tempCalibUnc.push_back(CalibVsCrysID->GetBinError(crysID + 1));
446 B2RESULT(
"eclMuMuEAlgorithm: successfully stored ECLCrystalEnergyMuMu calibration constants");
453 PeakVsCrysID->Write();
454 EdgeVsCrysID->Write();
455 effSigVsCrysID->Write();
456 etaVsCrysID->Write();
457 normVsCrysID->Write();
458 lowerLimitVsCrysID->Write();
459 fitLimitVsCrysID->Write();
460 StatusVsCrysID->Write();
462 fracPeakUnc->Write();
467 ExpEnergyperCrys->Write();
469 CalibVsCrysID->Write();
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;
493 B2RESULT(
"eclMuMuEAlgorithm performed fits but was not asked to store constants");
495 }
else if (!DBsuccess) {
496 if (
findExpValues) { B2RESULT(
"eclMuMuEAlgorithm: failed to store expected values"); }
497 else { B2RESULT(
"eclMuMuEAlgorithm: 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 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.
eclMuMuEAlgorithm()
..Constructor
int notFit
no fit performed
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
const int c_NCrystals
Number of crystals.
Abstract base class for different kinds of events.