69 double limitTol = 0.0005;
70 double minFitLimit = 1e-25;
71 double minFitProbIter = 1e-8;
72 double constRatio = 0.5;
73 double peakMin(0.5), peakMax(1.75);
74 double peakTol = limitTol * (peakMax - peakMin);
75 double effSigMin(0.08), effSigMax(0.5);
76 double effSigTol = limitTol * (effSigMax - effSigMin);
77 double etaMin(-3.), etaMax(1.);
79 double etaTol = limitTol * (etaMax - etaMin);
80 double constTol = 0.1;
88 dummy = (TH1F*)gROOT->FindObject(
"EnvsCrysSameRing");
89 if (dummy) {
delete dummy;}
90 dummy = (TH1F*)gROOT->FindObject(
"EnvsCrysDifferentRing");
91 if (dummy) {
delete dummy;}
92 dummy = (TH1F*)gROOT->FindObject(
"IntegralVsCrysSame");
93 if (dummy) {
delete dummy;}
94 dummy = (TH1F*)gROOT->FindObject(
"IntegralVsCrysDifferent");
95 if (dummy) {
delete dummy;}
96 dummy = (TH1F*)gROOT->FindObject(
"AverageExpECrysSame");
97 if (dummy) {
delete dummy;}
98 dummy = (TH1F*)gROOT->FindObject(
"AverageExpECrysDifferent");
99 if (dummy) {
delete dummy;}
100 dummy = (TH1F*)gROOT->FindObject(
"AverageElecCalibSame");
101 if (dummy) {
delete dummy;}
102 dummy = (TH1F*)gROOT->FindObject(
"AverageElecCalibDifferent");
103 if (dummy) {
delete dummy;}
104 dummy = (TH1F*)gROOT->FindObject(
"AverageInitialCalibSame");
105 if (dummy) {
delete dummy;}
106 dummy = (TH1F*)gROOT->FindObject(
"AverageInitialCalibDifferent");
107 if (dummy) {
delete dummy;}
111 std::vector<std::shared_ptr<TH2F>> EnvsCrys;
115 std::vector<std::shared_ptr<TH1F>> ExpEvsCrys;
119 std::vector<std::shared_ptr<TH1F>> ElecCalibvsCrys;
123 std::vector<std::shared_ptr<TH1F>> InitialCalibvsCrys;
127 std::vector<std::shared_ptr<TH1F>> CalibEntriesvsCrys;
136 TH1F* IntegralVsCrys[2];
137 IntegralVsCrys[0] =
new TH1F(
"IntegralVsCrysSame",
"Integral of EnVsCrys for each crystal, same theta ring;Crystal ID",
140 IntegralVsCrys[1] =
new TH1F(
"IntegralVsCrysDifferent",
"Integral of EnVsCrys for each crystal, different theta rings;Crystal ID",
143 TH1F* AverageExpECrys[2];
144 AverageExpECrys[0] =
new TH1F(
"AverageExpECrysSame",
147 AverageExpECrys[1] =
new TH1F(
"AverageExpECrysDifferent",
151 TH1F* AverageElecCalib[2];
152 AverageElecCalib[0] =
new TH1F(
"AverageElecCalibSame",
155 AverageElecCalib[1] =
new TH1F(
"AverageElecCalibDifferent",
159 TH1F* AverageInitialCalib[2];
160 AverageInitialCalib[0] =
new TH1F(
"AverageInitialCalibSame",
163 AverageInitialCalib[1] =
new TH1F(
"AverageInitialCalibDifferent",
168 int histbin = crysID + 1;
169 for (
int idir = 0; idir < 2; idir++) {
170 TH1D* hEnergy = EnvsCrys[idir]->ProjectionY(
"hEnergy", histbin, histbin);
171 int Integral = hEnergy->Integral();
172 IntegralVsCrys[idir]->SetBinContent(histbin, Integral);
174 double TotEntries = CalibEntriesvsCrys[idir]->GetBinContent(histbin);
176 double expectedE = 0.;
177 if (TotEntries > 0.) {expectedE = ExpEvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
178 AverageExpECrys[idir]->SetBinContent(histbin, expectedE);
179 AverageExpECrys[idir]->SetBinError(histbin, 0.);
181 double calibconst = 0.;
182 if (TotEntries > 0.) {calibconst = ElecCalibvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
183 AverageElecCalib[idir]->SetBinContent(histbin, calibconst);
184 AverageElecCalib[idir]->SetBinError(histbin, 0);
187 if (TotEntries > 0.) {calibconst = InitialCalibvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
188 AverageInitialCalib[idir]->SetBinContent(histbin, calibconst);
189 AverageInitialCalib[idir]->SetBinError(histbin, 0);
195 TFile* histfile =
new TFile(
"eclCosmicEAlgorithm.root",
"recreate");
196 for (
int idir = 0; idir < 2; idir++) {
197 EnvsCrys[idir]->Write();
198 IntegralVsCrys[idir]->Write();
199 AverageExpECrys[idir]->Write();
200 AverageElecCalib[idir]->Write();
201 AverageInitialCalib[idir]->Write();
203 RawDigitAmpvsCrys->Write();
208 B2INFO(
"eclCosmicEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
217 bool sufficientData =
true;
219 int histbin = crysID + 1;
220 bool SameLow = IntegralVsCrys[0]->GetBinContent(histbin) <
minEntries;
221 bool DifferentLow = IntegralVsCrys[1]->GetBinContent(histbin) <
minEntries;
222 if ((SameLow && DifferentLow) || (
findExpValues && (SameLow || DifferentLow))) {
223 if (
storeConst == 1) {B2INFO(
"eclCosmicEAlgorithm: cellID " << histbin <<
" has insufficient statistics: " << IntegralVsCrys[0]->GetBinContent(histbin) <<
" and " << IntegralVsCrys[1]->GetBinContent(histbin) <<
". Requirement is " <<
minEntries);}
224 sufficientData =
false;
238 const TString preName[2] = {
"SameRing",
"DifferentRing"};
240 TH1F* PeakperCrys[2];
241 PeakperCrys[0] =
new TH1F(
"PeakperCrysSame",
"Fit peak per crystal, same theta ring;Crystal ID;Peak normalized energy",
244 PeakperCrys[1] =
new TH1F(
"PeakperCrysDifferent",
"Fit peak per crystal, different theta ring;Crystal ID;Peak normalized energy",
247 TH1F* SigmaperCrys[2];
248 SigmaperCrys[0] =
new TH1F(
"SigmaperCrysSame",
"Fit sigma per crysal, same theta ring;Crystal ID;Sigma (ADC)",
250 SigmaperCrys[1] =
new TH1F(
"SigmaperCrysDifferent",
"Fit sigma per crysal, different theta ring;Crystal ID;Sigma (ADC)",
257 EtaperCrys[1] =
new TH1F(
"EtaperCrysDifferent",
"Fit eta per crysal, different theta ring;Crystal ID;Eta",
260 TH1F* ConstperCrys[2];
261 ConstperCrys[0] =
new TH1F(
"ConstperCrysSame",
"Fit constant per crystal, same theta ring;Crystal ID;Constant",
263 ConstperCrys[1] =
new TH1F(
"ConstperCrysDifferent",
"Fit constant per crystal, different theta ring;Crystal ID;Constant",
267 TH1F* StatusperCrys[2];
270 StatusperCrys[1] =
new TH1F(
"StatusperCrysDifferent",
"Fit status, different theta ring;Crystal ID;Status",
275 hStatus[0] =
new TH1F(
"StatusSame",
"Fit status, same theta ring", 25, -5, 20);
276 hStatus[1] =
new TH1F(
"StatusDifferent",
"Fit status, different theta ring", 25, -5, 20);
278 TH1F* fracPeakUnc[2];
279 fracPeakUnc[0] =
new TH1F(
"fracPeakUncSame",
"Fractional uncertainty on peak location, same theta ring", 100, 0, 0.1);
280 fracPeakUnc[1] =
new TH1F(
"fracPeakUncDifferent",
"Fractional uncertainty on peak location, different theta ring", 100, 0, 0.1);
283 hfitProb[0] =
new TH1F(
"fitProbSame",
"Probability of fit, same theta ring", 200, 0, 0.02);
284 hfitProb[1] =
new TH1F(
"fitProbDifferent",
"Probability of fit, different theta ring", 200, 0, 0.02);
287 TH1F* ExpEnergyperCrys[2];
288 ExpEnergyperCrys[0] =
new TH1F(
"ExpEnergyperCrysSame",
"Expected energy per crystal, same theta ring;Crystal ID;Peak energy (GeV)",
290 ExpEnergyperCrys[1] =
new TH1F(
"ExpEnergyperCrysDifferent",
294 TH1F* CalibvsCrys =
new TH1F(
"CalibvsCrys",
"Energy calibration constant from cosmics;Crystal ID;Calibration constant",
300 bool allFitsOK =
true;
302 int histbin = crysID + 1;
303 for (
int idir = 0; idir < 2; idir++) {
306 TString hname = preName[idir];
307 hname +=
"Enormalized";
309 TH1D* hEnergy = EnvsCrys[idir]->ProjectionY(hname, histbin, histbin);
312 double histMin = hEnergy->GetXaxis()->GetXmin();
313 double histMax = hEnergy->GetXaxis()->GetXmax();
314 TF1* func =
new TF1(
"eclCosmicNovoConst", eclCosmicNovoConst, histMin, histMax, 5);
315 func->SetParNames(
"normalization",
"peak",
"effSigma",
"eta",
"const");
316 func->SetParLimits(1, peakMin, peakMax);
317 func->SetParLimits(2, effSigMin, effSigMax);
318 func->SetParLimits(3, etaMin, etaMax);
321 hEnergy->GetXaxis()->SetRangeUser(peakMin, peakMax);
322 int maxBin = hEnergy->GetMaximumBin();
323 double peakE = hEnergy->GetBinLowEdge(maxBin);
324 double peakEUnc = 0.;
325 double normalization = hEnergy->GetMaximum();
326 double effSigma = hEnergy->GetRMS();
327 double sigmaUnc = 0.;
328 hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
331 double fitlow = 0.25;
332 double fithigh = peakE + 2.5 * effSigma;
335 int il0 = hEnergy->GetXaxis()->FindBin(fitlow);
336 int il1 = hEnergy->GetXaxis()->FindBin(fitlow + 0.1);
337 double constant = hEnergy->Integral(il0, il1) / (1 + il1 - il0);
338 double constUnc = 0.;
345 double dIter = 0.1 * (histMax - histMin) / hEnergy->GetNbinsX();
346 double highold(0.), higholdold(0.);
348 double fitProbDefault(0.);
349 bool fitHist = IntegralVsCrys[idir]->GetBinContent(histbin) >=
minEntries;
350 bool fixConst =
false;
359 func->SetParameters(normalization, peakE, effSigma, eta, constant);
360 if (fixConst) { func->FixParameter(4, 0); }
363 hEnergy->Fit(func,
"LIQ",
"", fitlow, fithigh);
364 normalization = func->GetParameter(0);
365 peakE = func->GetParameter(1);
366 peakEUnc = func->GetParError(1);
367 effSigma = func->GetParameter(2);
368 sigmaUnc = func->GetParError(2);
369 eta = func->GetParameter(3);
370 etaUnc = func->GetParError(3);
371 constant = func->GetParameter(4);
372 constUnc = func->GetParError(4);
373 fitProbDefault = func->GetProb();
377 double peak = func->Eval(peakE) - constant;
378 double tRatio = (func->Eval(fithigh) - constant) / peak;
379 if (tRatio < tRatioMin || tRatio >
tRatioMax) {
381 higholdold = highold;
383 fithigh = func->GetX(targetY, peakE, histMax);
387 if (abs(fithigh - higholdold) < dIter) {fithigh = 0.5 * (highold + higholdold); }
390 if (nIter >
maxIterations - 3) {fithigh = 0.33333 * (fithigh + highold + higholdold); }
394 if (constant < constTol && !fixConst) {
402 B2DEBUG(200,
"cellID = " << histbin <<
" " << nIter <<
" " << preName[idir] <<
" " << peakE <<
" " << constant <<
" " << tRatio <<
411 int lowbin = hEnergy->GetXaxis()->FindBin(fitlow);
412 int highbin = hEnergy->GetXaxis()->FindBin(fithigh);
414 if (fixConst) {npar = 4;}
415 int ndeg = (highbin - lowbin) + 1 - npar;
417 double halfbinwidth = 0.5 * hEnergy->GetBinWidth(1);
418 for (
int ib = lowbin; ib <= highbin; ib++) {
419 double yexp = func->Eval(hEnergy->GetBinLowEdge(ib) + halfbinwidth);
420 double yobs = hEnergy->GetBinContent(ib);
421 double dchi2 = (yexp - yobs) * (yexp - yobs) / yexp;
424 fitProb = 0.5 * (TMath::Prob(chisq, ndeg) + fitProbDefault);
433 if (normalization < constRatio * constant) {iStatus =
noPeak;}
436 if (fitProb <= minFitLimit || (fitProb < minFitProbIter && iStatus ==
iterations)) {iStatus =
poorFit;}
439 if ((peakE < peakMin + peakTol) || (peakE > peakMax - peakTol)) {iStatus =
atLimit;}
440 if ((effSigma < effSigMin + effSigTol) || (effSigma > effSigMax - effSigTol)) {iStatus =
atLimit;}
441 if ((eta < etaMin + etaTol) || (eta > etaMax - etaTol)) {iStatus =
atLimit;}
444 if (nIter == 0) {iStatus =
notFit;}
447 PeakperCrys[idir]->SetBinContent(histbin, peakE);
448 PeakperCrys[idir]->SetBinError(histbin, peakEUnc);
449 SigmaperCrys[idir]->SetBinContent(histbin, effSigma);
450 SigmaperCrys[idir]->SetBinError(histbin, sigmaUnc);
451 EtaperCrys[idir]->SetBinContent(histbin, eta);
452 EtaperCrys[idir]->SetBinError(histbin, etaUnc);
453 ConstperCrys[idir]->SetBinContent(histbin, constant);
454 ConstperCrys[idir]->SetBinError(histbin, constUnc);
455 StatusperCrys[idir]->SetBinContent(histbin, iStatus);
456 hStatus[idir]->Fill(iStatus);
457 fracPeakUnc[idir]->Fill(peakEUnc / peakE);
458 hfitProb[idir]->Fill(fitProb);
461 B2INFO(
"cellID " << histbin <<
" " << preName[idir] <<
" status = " << iStatus <<
" fit probability = " << fitProb);
473 int histbin = crysID + 1;
474 bool atLeastOneOK =
false;
475 for (
int idir = 0; idir < 2; idir++) {
476 double fitstatus = StatusperCrys[idir]->GetBinContent(histbin);
477 double peakE = PeakperCrys[idir]->GetBinContent(histbin);
478 double peakEUnc = PeakperCrys[idir]->GetBinError(histbin);
483 B2INFO(
"eclCosmicEAlgorithm: crystal " << crysID <<
" " << preName[idir] <<
" is not a successful fit. Status = " << fitstatus);
490 double inputExpE = abs(AverageExpECrys[idir]->GetBinContent(histbin));
491 ExpEnergyperCrys[idir]->SetBinContent(histbin, inputExpE * peakE);
492 ExpEnergyperCrys[idir]->SetBinError(histbin, inputExpE * peakEUnc / peakE);
494 if (!atLeastOneOK) {allFitsOK =
false;}
503 int histbin = crysID + 1;
504 double calibConst[2] = {};
505 double calibConstUnc[2] = {999999., 999999.};
506 double weight[2] = {};
507 bool bothFitsBad =
true;
508 for (
int idir = 0; idir < 2; idir++) {
511 double peakE = PeakperCrys[idir]->GetBinContent(histbin);
512 double fracPeakEUnc = PeakperCrys[idir]->GetBinError(histbin) / peakE;
513 double inputConst = AverageInitialCalib[idir]->GetBinContent(histbin);
514 double fitstatus = StatusperCrys[idir]->GetBinContent(histbin);
515 double inputExpE = AverageExpECrys[idir]->GetBinContent(histbin);
516 if (fitstatus >=
iterations && inputConst == 0) {B2FATAL(
"eclCosmicEAlgorithm: input calibration = 0 for idir = " << idir <<
" and crysID = " << crysID);}
519 if (fitstatus >=
iterations && inputExpE > 0.) {
520 calibConst[idir] = abs(inputConst) / peakE;
521 calibConstUnc[idir] = calibConst[idir] * fracPeakEUnc / peakE;
522 weight[idir] = 1. / (calibConstUnc[idir] * calibConstUnc[idir]);
526 B2INFO(
"eclCosmicEAlgorithm: cellID " << histbin <<
" " << preName[idir] <<
" is not a successful fit. Status = " << fitstatus);
528 B2WARNING(
"eclCosmicEAlgorithm: cellID " << histbin <<
" " << preName[idir] <<
" has no expected energy. Status = " << fitstatus);
535 double averageConstUnc;
539 if (histbin >=
cellIDLo && histbin <=
cellIDHi) {B2INFO(
"eclCosmicEAlgorithm: no constant found for cellID = " << histbin <<
" status = " << StatusperCrys[0]->GetBinContent(histbin) <<
" and " << StatusperCrys[1]->GetBinContent(histbin));}
540 averageConst = -1.*abs(AverageInitialCalib[0]->GetBinContent(histbin));
541 averageConstUnc = 0.;
543 averageConst = (calibConst[0] * weight[0] + calibConst[1] * weight[1]) / (weight[0] + weight[1]);
544 averageConstUnc = 1. /
sqrt(weight[0] + weight[1]);
546 CalibvsCrys->SetBinContent(histbin, averageConst);
547 CalibvsCrys->SetBinError(histbin, averageConstUnc);
553 bool DBsuccess =
false;
559 std::vector<std::string> DBname = {
"ECLExpCosmicESame",
"ECLExpCosmicEDifferent"};
560 for (
int idir = 0; idir < 2; idir++) {
561 std::vector<float> tempE;
562 std::vector<float> tempUnc;
564 int histbin = crysID + 1;
565 tempE.push_back(ExpEnergyperCrys[idir]->GetBinContent(histbin));
566 tempUnc.push_back(ExpEnergyperCrys[idir]->GetBinError(histbin));
571 B2INFO(
"eclCosmicEAlgorithm: successfully stored expected values for " << DBname[idir]);
576 std::vector<float> tempCalib;
577 std::vector<float> tempCalibUnc;
579 int histbin = crysID + 1;
580 tempCalib.push_back(CalibvsCrys->GetBinContent(histbin));
581 tempCalibUnc.push_back(CalibvsCrys->GetBinError(histbin));
586 B2INFO(
"eclCosmicEAlgorithm: successfully stored calibration constants");
591 for (
int idir = 0; idir < 2; idir++) {
592 PeakperCrys[idir]->Write();
593 SigmaperCrys[idir]->Write();
594 EtaperCrys[idir]->Write();
595 ConstperCrys[idir]->Write();
596 StatusperCrys[idir]->Write();
597 hStatus[idir]->Write();
598 fracPeakUnc[idir]->Write();
599 hfitProb[idir]->Write();
604 ExpEnergyperCrys[0]->Write();
605 ExpEnergyperCrys[1]->Write();
607 CalibvsCrys->Write();
613 dummy = (TH1F*)gROOT->FindObject(
"PeakperCrysSame");
delete dummy;
614 dummy = (TH1F*)gROOT->FindObject(
"SigmaperCrysSame");
delete dummy;
615 dummy = (TH1F*)gROOT->FindObject(
"EtaperCrysSame");
delete dummy;
616 dummy = (TH1F*)gROOT->FindObject(
"ConstperCrysSame");
delete dummy;
617 dummy = (TH1F*)gROOT->FindObject(
"StatusperCrysSame");
delete dummy;
618 dummy = (TH1F*)gROOT->FindObject(
"StatusSame");
delete dummy;
619 dummy = (TH1F*)gROOT->FindObject(
"fracPeakUncSame");
delete dummy;
620 dummy = (TH1F*)gROOT->FindObject(
"fitProbSame");
delete dummy;
621 dummy = (TH1F*)gROOT->FindObject(
"ExpEnergyperCrysSame");
delete dummy;
622 dummy = (TH1F*)gROOT->FindObject(
"PeakperCrysDifferent");
delete dummy;
623 dummy = (TH1F*)gROOT->FindObject(
"SigmaperCrysDifferent");
delete dummy;
624 dummy = (TH1F*)gROOT->FindObject(
"EtaperCrysDifferent");
delete dummy;
625 dummy = (TH1F*)gROOT->FindObject(
"ConstperCrysDifferent");
delete dummy;
626 dummy = (TH1F*)gROOT->FindObject(
"StatusperCrysDifferent");
delete dummy;
627 dummy = (TH1F*)gROOT->FindObject(
"StatusDifferent");
delete dummy;
628 dummy = (TH1F*)gROOT->FindObject(
"fracPeakUncDifferent");
delete dummy;
629 dummy = (TH1F*)gROOT->FindObject(
"fitProbDifferent");
delete dummy;
630 dummy = (TH1F*)gROOT->FindObject(
"ExpEnergyperCrysDifferent");
delete dummy;
631 dummy = (TH1F*)gROOT->FindObject(
"CalibvsCrys");
delete dummy;
636 B2INFO(
"eclCosmicEAlgorithm performed fits but was not asked to store constants");
638 }
else if (!DBsuccess) {
639 if (
findExpValues) { B2INFO(
"eclCosmicEAlgorithm: failed to store expected values"); }
640 else { B2INFO(
"eclCosmicEAlgorithm: failed to store calibration constants"); }