69 const double limitTol = 0.0005;
70 const double minFitLimit = 1e-25;
71 const double minFitProbIter = 1e-8;
72 const double constRatio = 0.5;
73 const double peakMin(0.4), peakMax(2.2);
74 const double peakTol = limitTol * (peakMax - peakMin);
75 const double effSigMin(0.02), effSigMax(0.4);
76 const double effSigTol = limitTol * (effSigMax - effSigMin);
77 const double etaNom(0.41);
78 const double etaMin(0.), etaMax(5.0);
79 const double etaTol = limitTol * (etaMax - etaMin);
80 const double constTol = 0.001;
87 B2INFO(
"eclGammaGammaAlgorithm parameters:");
106 dummy = (TH1F*)gROOT->FindObject(
"IntegralVsCrysID");
107 if (dummy) {
delete dummy;}
108 dummy = (TH1F*)gROOT->FindObject(
"AverageExpECrys");
109 if (dummy) {
delete dummy;}
110 dummy = (TH1F*)gROOT->FindObject(
"AverageElecCalib");
111 if (dummy) {
delete dummy;}
112 dummy = (TH1F*)gROOT->FindObject(
"AverageInitCalib");
113 if (dummy) {
delete dummy;}
126 TH1F* IntegralVsCrysID =
new TH1F(
"IntegralVsCrysID",
"Integral of EnVsCrysID for each crystal;crystal ID;Entries",
128 TH1F* AverageExpECrys =
new TH1F(
"AverageExpECrys",
"Average expected E per crys from collector;Crystal ID;Energy (GeV)",
131 TH1F* AverageElecCalib =
new TH1F(
"AverageElecCalib",
"Average electronics calib const vs crystal;Crystal ID;Calibration constant",
133 TH1F* AverageInitCalib =
new TH1F(
"AverageInitCalib",
"Average initial calib const vs crystal;Crystal ID;Calibration constant",
137 TH1D* hEnergy = EnVsCrysID->ProjectionY(
"hEnergy", crysID + 1, crysID + 1);
138 int Integral = hEnergy->Integral();
139 IntegralVsCrysID->SetBinContent(crysID + 1, Integral);
141 double TotEntries = CalibEntriesvsCrys->GetBinContent(crysID + 1);
143 double expectedE = 0.;
144 if (TotEntries > 0.) {expectedE = ExpEvsCrys->GetBinContent(crysID + 1) / TotEntries;}
145 AverageExpECrys->SetBinContent(crysID + 1, expectedE);
147 double calibconst = 0.;
148 if (TotEntries > 0.) {calibconst = ElecCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
149 AverageElecCalib->SetBinContent(crysID + 1, calibconst);
152 if (TotEntries > 0.) {calibconst = InitialCalibvsCrys->GetBinContent(crysID + 1) / TotEntries;}
153 AverageInitCalib->SetBinContent(crysID + 1, calibconst);
159 TFile* histfile =
new TFile(fName,
"recreate");
161 IntegralVsCrysID->Write();
162 AverageExpECrys->Write();
163 AverageElecCalib->Write();
164 AverageInitCalib->Write();
169 B2INFO(
"eclGammaGammaEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
176 bool sufficientData =
true;
178 if (IntegralVsCrysID->GetBinContent(crysID + 1) <
m_minEntries) {
179 if (
m_storeConst == 1) {B2INFO(
"eclGammaGammaEAlgorithm: crystal " << crysID <<
" has insufficient statistics: " << IntegralVsCrysID->GetBinContent(crysID + 1) <<
". Requirement is " <<
m_minEntries);}
180 sufficientData =
false;
195 TH1F* CalibVsCrysID =
new TH1F(
"CalibVsCrysID",
"Calibration constant vs crystal ID;crystal ID;counts per GeV",
197 TH1F* ExpEnergyperCrys =
new TH1F(
"ExpEnergyperCrys",
"Expected energy per crystal;Crystal ID;Peak energy (GeV)",
201 TH1F* PeakVsCrysID =
new TH1F(
"PeakVsCrysID",
"Peak of Novo fit vs crystal ID;crystal ID;Peak normalized energy",
203 TH1F* EdgeVsCrysID =
new TH1F(
"EdgeVsCrysID",
"Upper edge of Novo fit vs crystal ID;crystal ID;Maximum normalized energy",
210 TH1F* constVsCrysID =
new TH1F(
"constVsCrysID",
"fit constant vs crystal ID;crystal ID;fit constant",
212 TH1F* normVsCrysID =
new TH1F(
"normVsCrysID",
"Novosibirsk normalization vs crystal ID;crystal ID;normalization",
214 TH1F* fitLimitVsCrysID =
new TH1F(
"fitLimitVsCrysID",
"fit range lower limit vs crystal ID;crystal ID;upper fit limit",
219 TH1F* FitProbVsCrysID =
new TH1F(
"FitProbVsCrysID",
"Fit probability vs crystal id;crystal ID;Fit probability",
223 TH1F* hStatus =
new TH1F(
"hStatus",
"Fit status", 25, -5, 20);
224 TH1F* hPeak =
new TH1F(
"hPeak",
"Peaks of normalized energy distributions, successful fits;Peak of Novosibirsk fit", 200, 0.8, 1.2);
225 TH1F* fracPeakUnc =
new TH1F(
"fracPeakUnc",
"Fractional uncertainty on peak uncertainty, successful fits;Uncertainty on peak", 100,
227 TH1F* nIterations =
new TH1F(
"nIterations",
"Number of times histogram was fit;Number of iterations", 20, -0.5, 19.5);
232 bool allFitsOK =
true;
236 TString name =
"Enormalized";
238 TH1D* hEnergy = EnVsCrysID->ProjectionY(name, crysID + 1, crysID + 1);
241 double histMin = hEnergy->GetXaxis()->GetXmin();
242 double histMax = hEnergy->GetXaxis()->GetXmax();
243 TF1* func =
new TF1(
"eclGammaGammaNovoConst", eclGammaGammaNovoConst, histMin, histMax, 5);
244 func->SetParNames(
"normalization",
"peak",
"effSigma",
"eta",
"const");
245 func->SetParLimits(1, peakMin, peakMax);
246 func->SetParLimits(2, effSigMin, effSigMax);
247 func->SetParLimits(3, etaMin, etaMax);
248 func->SetParLimits(4, 0., hEnergy->GetMaximum());
251 hEnergy->GetXaxis()->SetRangeUser(peakMin, peakMax);
252 int maxBin = hEnergy->GetMaximumBin();
253 double peakE = hEnergy->GetBinLowEdge(maxBin);
254 double peakEUnc = 0.;
255 double normalization = hEnergy->GetMaximum();
257 double effSigma = hEnergy->GetRMS();
258 double sigmaUnc = 0.;
259 hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
262 double fitlow = peakE - effSigma;
263 double fithigh = histMax;
268 double constant = 0.01 * normalization;
269 double constUnc = 0.;
272 double dIter = 0.1 * (histMax - histMin) / hEnergy->GetNbinsX();
274 double fitProbDefault(0.);
275 double lowold(0.), lowoldold(0.);
276 bool fixConst =
false;
278 double histIntegral = IntegralVsCrysID->GetBinContent(crysID + 1);
294 func->SetParameters(normalization, peakE, effSigma, eta, constant);
295 if (fixConst) { func->FixParameter(4, 0); }
298 hEnergy->Fit(func,
"LIQ",
"", fitlow, fithigh);
301 normalization = func->GetParameter(0);
302 normUnc = func->GetParError(0);
303 peakE = func->GetParameter(1);
304 peakEUnc = func->GetParError(1);
305 effSigma = func->GetParameter(2);
306 sigmaUnc = func->GetParError(2);
307 eta = func->GetParameter(3);
308 etaUnc = func->GetParError(3);
309 constant = func->GetParameter(4);
310 constUnc = func->GetParError(4);
311 fitProbDefault = func->GetProb();
314 double peak = func->Eval(peakE) - constant;
315 double tRatio = (func->Eval(fitlow) - constant) / peak;
316 if (tRatio < m_tRatioMin || tRatio > m_tRatioMax) {
317 double targetY = constant + 0.5 * (m_tRatioMin + m_tRatioMax) * peak;
320 fitlow = func->GetX(targetY, histMin, peakE);
324 if (abs(fitlow - lowoldold) < dIter) {fitlow = 0.5 * (lowold + lowoldold); }
327 if (nIter >
m_maxIterations - 3) {fitlow = 0.33333 * (fitlow + lowold + lowoldold); }
331 if (constant < constTol && !fixConst) {
339 B2DEBUG(10, crysID <<
" " << nIter <<
" " << peakE <<
" " << constant <<
" " << tRatio <<
" " << fitlow);
346 int lowbin = hEnergy->GetXaxis()->FindBin(fitlow);
347 int highbin = hEnergy->GetXaxis()->FindBin(fithigh);
349 if (fixConst) {npar = 4;}
352 double binwidth = hEnergy->GetBinWidth(1);
353 for (
int ib = lowbin; ib <= highbin; ib++) {
354 double xlow = hEnergy->GetBinLowEdge(ib);
355 double yexp = func->Integral(xlow, xlow + binwidth) / binwidth;
358 if (yexp > constTol) {
359 double yobs = hEnergy->GetBinContent(ib);
361 if (yexp < 0.9999 && yobs > yexp) {dnom = yobs;}
362 double dchi2 = (yexp - yobs) * (yexp - yobs) / dnom;
367 fitProb = TMath::Prob(chisq, ndeg);
376 if (normalization < constRatio * constant) {iStatus =
noPeak;}
379 if (fitProb <= minFitLimit || (fitProb < minFitProbIter && iStatus ==
iterations)) {iStatus =
poorFit;}
382 if ((peakE < peakMin + peakTol) || (peakE > peakMax - peakTol)) {iStatus =
atLimit;}
383 if ((effSigma < effSigMin + effSigTol) || (effSigma > effSigMax - effSigTol)) {iStatus =
atLimit;}
384 if ((eta < etaMin + etaTol) || (eta > etaMax - etaTol)) {iStatus =
atLimit;}
387 if (nIter == 0) {iStatus =
notFit;}
391 double upperEdge = peakE;
392 double edgeUnc = peakEUnc;
400 int iLow = hEnergy->GetXaxis()->FindBin(peakE) + 1;
401 int iHigh = hEnergy->GetNbinsX();
403 for (
int ibin = iLow; ibin < iHigh; ibin++) {
404 double xc = hEnergy->GetBinCenter(ibin);
405 if (func->Eval(xc) > targetY) {iLast = ibin;}
407 double xLow = hEnergy->GetBinCenter(iLast);
408 double xHigh = hEnergy->GetBinCenter(iLast + 1);
412 upperEdge = func->GetX(targetY, xLow, xHigh);
415 }
else if (iStatus >
notFit) {
419 int thisBin = hEnergy->GetBinContent(1);
420 for (
int ibin = 2; ibin < hEnergy->GetNbinsX(); ibin++) {
421 int prevBin = thisBin;
422 thisBin = hEnergy->GetBinContent(ibin);
423 if (thisBin > 0 && thisBin + prevBin >= 2) {iLast = ibin;}
425 if (iLast > 0) {upperEdge = hEnergy->GetBinLowEdge(iLast);}
430 int histbin = crysID + 1;
431 PeakVsCrysID->SetBinContent(histbin, peakE);
432 PeakVsCrysID->SetBinError(histbin, peakEUnc);
433 EdgeVsCrysID->SetBinContent(histbin, upperEdge);
434 EdgeVsCrysID->SetBinError(histbin, edgeUnc);
435 effSigVsCrysID->SetBinContent(histbin, effSigma);
436 effSigVsCrysID->SetBinError(histbin, sigmaUnc);
437 etaVsCrysID->SetBinContent(histbin, eta);
438 etaVsCrysID->SetBinError(histbin, etaUnc);
439 constVsCrysID->SetBinContent(histbin, constant);
440 constVsCrysID->SetBinError(histbin, constUnc);
441 normVsCrysID->SetBinContent(histbin, normalization);
442 normVsCrysID->SetBinError(histbin, normUnc);
443 fitLimitVsCrysID->SetBinContent(histbin, fitlow);
444 fitLimitVsCrysID->SetBinError(histbin, 0);
445 StatusVsCrysID->SetBinContent(histbin, iStatus);
446 StatusVsCrysID->SetBinError(histbin, 0);
447 FitProbVsCrysID->SetBinContent(histbin, fitProb);
448 FitProbVsCrysID->SetBinError(histbin, 0);
451 hStatus->Fill(iStatus);
452 nIterations->Fill(nIter);
455 fracPeakUnc->Fill(peakEUnc / peakE);
459 B2INFO(
"cellID " << crysID + 1 <<
" status = " << iStatus <<
" fit probability = " << fitProb <<
" default prob = " <<
469 int histbin = crysID + 1;
470 double fitstatus = StatusVsCrysID->GetBinContent(histbin);
471 double upperEdge = EdgeVsCrysID->GetBinContent(histbin);
472 double fracEdgeUnc = EdgeVsCrysID->GetBinError(histbin) / upperEdge;
479 B2INFO(
"eclGammaGammaEAlgorithm: cellID " << histbin <<
" is not a successful fit. Status = " << fitstatus);
486 double inputExpE = abs(AverageExpECrys->GetBinContent(histbin));
487 ExpEnergyperCrys->SetBinContent(histbin, inputExpE * upperEdge);
488 ExpEnergyperCrys->SetBinError(histbin, fracEdgeUnc * inputExpE * upperEdge);
492 double inputCalib = abs(AverageInitCalib->GetBinContent(histbin));
493 CalibVsCrysID->SetBinContent(histbin, inputCalib / upperEdge);
494 CalibVsCrysID->SetBinError(histbin, fracEdgeUnc * inputCalib / upperEdge);
500 bool DBsuccess =
false;
506 std::vector<float> tempE;
507 std::vector<float> tempUnc;
509 tempE.push_back(ExpEnergyperCrys->GetBinContent(crysID + 1));
510 tempUnc.push_back(ExpEnergyperCrys->GetBinError(crysID + 1));
515 B2INFO(
"eclCosmicEAlgorithm: successfully stored expected energies ECLExpGammaGammaE");
520 std::vector<float> tempCalib;
521 std::vector<float> tempCalibUnc;
523 tempCalib.push_back(CalibVsCrysID->GetBinContent(crysID + 1));
524 tempCalibUnc.push_back(CalibVsCrysID->GetBinError(crysID + 1));
529 B2INFO(
"eclGammaGammaEAlgorithm: successfully stored ECLCrystalEnergyGammaGamma calibration constants");
536 PeakVsCrysID->Write();
537 EdgeVsCrysID->Write();
538 effSigVsCrysID->Write();
539 etaVsCrysID->Write();
540 constVsCrysID->Write();
541 normVsCrysID->Write();
542 fitLimitVsCrysID->Write();
543 StatusVsCrysID->Write();
544 FitProbVsCrysID->Write();
546 fracPeakUnc->Write();
547 nIterations->Write();
552 ExpEnergyperCrys->Write();
554 CalibVsCrysID->Write();
560 dummy = (TH1F*)gROOT->FindObject(
"PeakVsCrysID");
delete dummy;
561 dummy = (TH1F*)gROOT->FindObject(
"EdgeVsCrysID");
delete dummy;
562 dummy = (TH1F*)gROOT->FindObject(
"effSigVsCrysID");
delete dummy;
563 dummy = (TH1F*)gROOT->FindObject(
"etaVsCrysID");
delete dummy;
564 dummy = (TH1F*)gROOT->FindObject(
"constVsCrysID");
delete dummy;
565 dummy = (TH1F*)gROOT->FindObject(
"normVsCrysID");
delete dummy;
566 dummy = (TH1F*)gROOT->FindObject(
"fitLimitVsCrysID");
delete dummy;
567 dummy = (TH1F*)gROOT->FindObject(
"StatusVsCrysID");
delete dummy;
568 dummy = (TH1F*)gROOT->FindObject(
"FitProbVsCrysID");
delete dummy;
569 dummy = (TH1F*)gROOT->FindObject(
"fracPeakUnc");
delete dummy;
570 dummy = (TH1F*)gROOT->FindObject(
"nIterations");
delete dummy;
571 dummy = (TH1F*)gROOT->FindObject(
"hStatus");
delete dummy;
572 dummy = (TH1F*)gROOT->FindObject(
"ExpEnergyperCrys");
delete dummy;
573 dummy = (TH1F*)gROOT->FindObject(
"CalibVsCrysID");
delete dummy;
579 B2INFO(
"eclGammaGammaEAlgorithm performed fits but was not asked to store constants");
581 }
else if (!DBsuccess) {
582 if (
m_findExpValues) { B2INFO(
"eclGammaGammaEAlgorithm: failed to store expected values"); }
583 else { B2INFO(
"eclGammaGammaEAlgorithm: failed to store calibration constants"); }