10#include <ecl/calibration/eclTValidationAlgorithm.h>
13#include <ecl/dataobjects/ECLElementNumbers.h>
14#include <ecl/dbobjects/ECLCrystalCalib.h>
15#include <ecl/digitization/EclConfiguration.h>
16#include <ecl/geometry/ECLGeometryPar.h>
17#include <ecl/mapper/ECLChannelMapper.h>
20#include <framework/database/DBObjPtr.h>
21#include <framework/database/DBStore.h>
22#include <framework/dataobjects/EventMetaData.h>
23#include <framework/datastore/DataStore.h>
24#include <framework/datastore/StoreObjPtr.h>
29#include <TGraphAsymmErrors.h>
47 cellIDHi(ECLElementNumbers::c_NCrystals),
48 readPrevCrysPayload(false),
49 meanCleanRebinFactor(1),
50 meanCleanCutMinFactor(0),
51 clusterTimesFractionWindow_maxtime(8),
52 debugFilenameBase(
"eclTValidationAlgorithm")
54 setDescription(
"Fit gaussian function to the cluster times to validate results.");
68 cellIDHi(ECLElementNumbers::c_NCrystals),
69 readPrevCrysPayload(false),
70 meanCleanRebinFactor(1),
71 meanCleanCutMinFactor(0),
72 clusterTimesFractionWindow_maxtime(8),
73 debugFilenameBase(
"eclTValidationAlgorithm")
75 setDescription(
"Fit gaussian function to the cluster times to validate results.");
88 B2INFO(
"eclTValidationAlgorithm parameters:");
98 auto clusterTime = getObjectPtr<TH1F>(
"clusterTime");
99 auto clusterTime_cid = getObjectPtr<TH2F>(
"clusterTime_cid");
100 auto clusterTime_run = getObjectPtr<TH2F>(
"clusterTime_run");
101 auto clusterTimeClusterE = getObjectPtr<TH2F>(
"clusterTimeClusterE");
102 auto dt99_clusterE = getObjectPtr<TH2F>(
"dt99_clusterE");
103 auto eventT0 = getObjectPtr<TH1F>(
"eventT0");
104 auto eventT0Detector = getObjectPtr<TH1F>(
"eventT0Detector");
105 auto clusterTimeE0E1diff = getObjectPtr<TH1F>(
"clusterTimeE0E1diff");
108 auto cutflow = getObjectPtr<TH1F>(
"cutflow");
110 vector <int> binProjectionLeft_Time_vs_E_runDep ;
111 vector <int> binProjectionRight_Time_vs_E_runDep ;
113 for (
int binCounter = 1; binCounter <= clusterTimeClusterE->GetNbinsX(); binCounter++) {
114 binProjectionLeft_Time_vs_E_runDep.push_back(binCounter);
115 binProjectionRight_Time_vs_E_runDep.push_back(binCounter);
130 int numCrysWithNonZeroEntries = 0 ;
131 int numCrysWithGoodFit = 0 ;
133 int minNumEntries = 40;
139 bool minRunNumBool =
false;
140 bool maxRunNumBool =
false;
146 int expNumber = expRun.first;
147 int runNumber = expRun.second;
148 if (!minRunNumBool) {
149 minExpNum = expNumber;
150 minRunNum = runNumber;
151 minRunNumBool =
true;
153 if (!maxRunNumBool) {
154 maxExpNum = expNumber;
155 maxRunNum = runNumber;
156 maxRunNumBool =
true;
158 if (((minRunNum > runNumber) && (minExpNum >= expNumber)) ||
159 (minExpNum > expNumber)) {
160 minExpNum = expNumber;
161 minRunNum = runNumber;
163 if (((maxRunNum < runNumber) && (maxExpNum <= expNumber)) ||
164 (maxExpNum < expNumber))
167 maxExpNum = expNumber;
168 maxRunNum = runNumber;
173 string runNumsString = string(
"_") + to_string(minExpNum) +
"_" + to_string(minRunNum) + string(
"-") +
174 to_string(maxExpNum) +
"_" + to_string(maxRunNum);
182 int eventNumberForCrates = 1;
194 evtPtr.
construct(eventNumberForCrates, minRunNum, minExpNum);
203 B2INFO(
"Uploading payload for exp " << minExpNum <<
", run " << minRunNum <<
", event " << eventNumberForCrates);
206 crystalMapper->initFromDB();
215 B2INFO(
"Dumping payload");
218 std::vector<float> currentValuesCrys = crystalTimeObject->getCalibVector();
219 std::vector<float> currentUncCrys = crystalTimeObject->getCalibUncVector();
222 B2INFO(
"Values read from database. Write out for their values for comparison against those from tcol");
224 B2INFO(
"ts: cellID " << ic + 1 <<
" " << currentValuesCrys[ic] <<
" +/- " << currentUncCrys[ic]);
233 prevValuesCrys = customPrevCrystalTimeObject->getCalibVector();
236 B2INFO(
"Previous values read from database. Write out for their values for comparison against those from tcol");
238 B2INFO(
"ts custom previous payload: cellID " << ic + 1 <<
" " << prevValuesCrys[ic]);
246 B2INFO(
"Debug output rootfile: " << debugFilename);
247 histfile =
new TFile(debugFilename.c_str(),
"recreate");
250 clusterTime ->Write();
251 clusterTime_cid ->Write();
252 clusterTime_run ->Write();
253 clusterTimeClusterE ->Write();
254 dt99_clusterE ->Write();
256 eventT0Detector ->Write();
257 clusterTimeE0E1diff ->Write();
262 double hist_tmin = clusterTime->GetXaxis()->GetXmin();
263 double hist_tmax = clusterTime->GetXaxis()->GetXmax();
264 int hist_nTbins = clusterTime->GetNbinsX();
266 B2INFO(
"hist_tmin = " << hist_tmin);
267 B2INFO(
"hist_tmax = " << hist_tmax);
268 B2INFO(
"hist_nTbins = " << hist_nTbins);
270 double time_fit_min = hist_tmax;
271 double time_fit_max = hist_tmin;
277 auto peakClusterTimes =
new TH1F(
"peakClusterTimes",
278 "-For crystals with at least one hit-;Peak cluster time [ns];Number of crystals",
279 hist_nTbins, hist_tmin, hist_tmax);
280 auto peakClusterTimesGoodFit =
new TH1F(
"peakClusterTimesGoodFit",
281 "-For crystals with a good fit to distribution of hits-;Peak cluster time [ns];Number of crystals",
282 hist_nTbins, hist_tmin, hist_tmax);
284 auto peakClusterTimesGoodFit__cid =
new TH1F(
"peakClusterTimesGoodFit__cid",
285 "-For crystals with a good fit to distribution of hits-;cell id (only crystals with good fit);Peak cluster time [ns]",
290 auto tsNew_MINUS_tsCustomPrev__cid =
new TH1F(
"TsNew_MINUS_TsCustomPrev__cid",
291 ";cell id; ts(new|merged) - ts(old = 'pre-calib'|merged) [ns]",
294 auto tsNew_MINUS_tsCustomPrev =
new TH1F(
"TsNew_MINUS_TsCustomPrev",
295 ";ts(new | merged) - ts(old = 'pre-calib' | merged) [ns];Number of crystals",
296 285, -69.5801, 69.5801);
302 B2INFO(
"timeWindow_maxTime = " << timeWindow_maxTime);
303 int binyLeft = clusterTime_cid->GetYaxis()->FindBin(-timeWindow_maxTime);
304 int binyRight = clusterTime_cid->GetYaxis()->FindBin(timeWindow_maxTime);
305 double windowLowTimeFromBin = clusterTime_cid->GetYaxis()->GetBinLowEdge(binyLeft);
306 double windowHighTimeFromBin = clusterTime_cid->GetYaxis()->GetBinLowEdge(binyRight + 1);
307 std::string s_lowTime = std::to_string(windowLowTimeFromBin);
308 std::string s_highTime = std::to_string(windowHighTimeFromBin);
309 TString fracWindowTitle =
"Fraction of cluster times in window [" + s_lowTime +
", " + s_highTime +
310 "] ns;cell id;Fraction of cluster times in window";
311 B2INFO(
"fracWindowTitle = " << fracWindowTitle);
312 TString fracWindowInGoodECLRingsTitle =
"Fraction of cluster times in window [" + s_lowTime +
", " + s_highTime +
313 "] ns and in good ECL theta rings;cell id;Fraction cluster times in window + good ECL rings";
314 B2INFO(
"fracWindowInGoodECLRingsTitle = " << fracWindowInGoodECLRingsTitle);
315 B2INFO(
"Good ECL rings skip gaps in the acceptance, and includes ECL theta IDs: 3-10, 15-39, 44-56, 61-66.");
317 TString fracWindowHistTitle =
"Fraction of cluster times in window [" + s_lowTime +
", " + s_highTime +
318 "] ns;Fraction of cluster times in window;Number of crystals";
322 auto clusterTimeNumberInWindowInGoodECLRings__cid =
new TH1F(
"clusterTimeNumberInWindowInGoodECLRings__cid", fracWindowTitle,
327 auto clusterTimeFractionInWindow =
new TH1F(
"clusterTimeFractionInWindow", fracWindowHistTitle, 110, 0.0, 1.1);
329 clusterTimeNumberInWindow__cid->Sumw2();
330 clusterTimeNumberInWindowInGoodECLRings__cid->Sumw2();
331 clusterTimeNumber__cid->Sumw2();
341 double clusterTime_mean = 0;
342 double clusterTime_mean_unc = 0;
344 B2INFO(
"Crystal cell id = " << crys_id);
353 TH1D* h_time = clusterTime_cid->ProjectionY((std::string(
"h_time_psi__") + std::to_string(crys_id)).c_str(),
355 TH1D* h_timeMask = (TH1D*)h_time->Clone();
356 TH1D* h_timeMasked = (TH1D*)h_time->Clone((std::string(
"h_time_psi_masked__") + std::to_string(crys_id)).c_str());
357 TH1D* h_timeRebin = (TH1D*)h_time->Clone();
364 h_timeMask->Scale(0.0);
366 time_fit_min = hist_tmax;
367 time_fit_max = hist_tmin;
370 double histRebin_max = h_timeRebin->GetMaximum();
372 bool maskedOutNonZeroBin =
false;
374 for (
int bin = 1; bin <= h_timeRebin->GetNbinsX(); bin++) {
377 if (nonRebinnedBinNumber < h_time->GetNbinsX()) {
379 h_timeMask->SetBinContent(nonRebinnedBinNumber, 1);
382 double x_lower = h_timeRebin->GetXaxis()->GetBinLowEdge(bin);
383 double x_upper = h_timeRebin->GetXaxis()->GetBinUpEdge(bin);
384 if (x_lower < time_fit_min) {
385 time_fit_min = x_lower;
387 if (x_upper > time_fit_max) {
388 time_fit_max = x_upper;
392 if (h_time->GetBinContent(nonRebinnedBinNumber) > 0) {
393 B2DEBUG(22,
"Setting bin " << nonRebinnedBinNumber <<
" from " << h_timeMasked->GetBinContent(nonRebinnedBinNumber) <<
" to 0");
394 maskedOutNonZeroBin =
true;
396 h_timeMasked->SetBinContent(nonRebinnedBinNumber, 0);
401 B2INFO(
"Bins with non-zero values have been masked out: " << maskedOutNonZeroBin);
402 h_timeMasked->ResetStats();
403 h_timeMask->ResetStats();
408 double default_meanMasked = h_timeMasked->GetMean();
410 B2INFO(
"default_meanMasked = " << default_meanMasked);
414 double default_mean = h_time->GetMean();
415 double default_mean_unc = h_time->GetMeanError();
416 double default_sigma = h_time->GetStdDev();
418 B2INFO(
"Fitting crystal between " << time_fit_min <<
" and " << time_fit_max);
421 TF1* gaus =
new TF1(
"func",
"gaus(0)", time_fit_min, time_fit_max);
422 gaus->SetParNames(
"numCrystalHitsNormalization",
"mean",
"sigma");
429 double hist_max = h_time->GetMaximum();
432 double stddev = h_time->GetStdDev();
437 gaus->SetParameter(0, hist_max / 2.);
438 gaus->SetParameter(1, mean);
439 gaus->SetParameter(2, sigma);
446 h_timeMasked->Fit(gaus,
"LQR");
448 double fit_mean = gaus->GetParameter(1);
449 double fit_mean_unc = gaus->GetParError(1);
450 double fit_sigma = gaus->GetParameter(2);
452 double meanDiff = fit_mean - default_mean;
453 double meanUncDiff = fit_mean_unc - default_mean_unc;
454 double sigmaDiff = fit_sigma - default_sigma;
456 bool good_fit =
false;
458 if ((fabs(meanDiff) > 10) ||
459 (fabs(meanUncDiff) > 10) ||
460 (fabs(sigmaDiff) > 10) ||
461 (fit_mean_unc > 3) ||
463 (fit_mean < time_fit_min) ||
464 (fit_mean > time_fit_max)) {
465 B2INFO(
"Crystal cell id = " << crys_id);
466 B2INFO(
"fit mean, default mean = " << fit_mean <<
", " << default_mean);
467 B2INFO(
"fit mean unc, default mean unc = " << fit_mean_unc <<
", " << default_mean_unc);
468 B2INFO(
"fit sigma, default sigma = " << fit_sigma <<
", " << default_sigma);
470 B2INFO(
"crystal fit mean - hist mean = " << meanDiff);
471 B2INFO(
"fit mean unc. - hist mean unc. = " << meanUncDiff);
472 B2INFO(
"fit sigma - hist sigma = " << sigmaDiff);
474 B2INFO(
"fit_mean = " << fit_mean);
475 B2INFO(
"time_fit_min = " << time_fit_min);
476 B2INFO(
"time_fit_max = " << time_fit_max);
478 if (fabs(meanDiff) > 10) B2INFO(
"fit mean diff too large");
479 if (fabs(meanUncDiff) > 10) B2INFO(
"fit mean unc diff too large");
480 if (fabs(sigmaDiff) > 10) B2INFO(
"fit mean sigma diff too large");
481 if (fit_mean_unc > 3) B2INFO(
"fit mean unc too large");
482 if (fit_sigma < 0.1) B2INFO(
"fit sigma too small");
486 numCrysWithGoodFit++;
487 crysHasGoodFit[crys_id - 1] = true ;
491 int numEntries = h_time->GetEntries();
494 if ((numEntries >= minNumEntries) && good_fit) {
495 clusterTime_mean = fit_mean;
496 clusterTime_mean_unc = fit_mean_unc;
497 crysHasGoodFitandStats[crys_id - 1] = true ;
499 clusterTime_mean = default_mean;
500 clusterTime_mean_unc = default_mean_unc;
503 if (numEntries < minNumEntries) B2INFO(
"Number of entries less than minimum");
504 if (numEntries == 0) B2INFO(
"Number of entries == 0");
507 t_offsets[crys_id - 1] = clusterTime_mean ;
508 t_offsets_unc[crys_id - 1] = clusterTime_mean_unc ;
509 numClusterPerCrys[crys_id - 1] = numEntries ;
511 histfile->WriteTObject(h_time, (std::string(
"h_time_psi") + std::to_string(crys_id)).c_str());
512 histfile->WriteTObject(h_timeMasked, (std::string(
"h_time_psi_masked") + std::to_string(crys_id)).c_str());
515 peakClusterTime_cid->SetBinContent(crys_id, t_offsets[crys_id - 1]);
516 peakClusterTime_cid->SetBinError(crys_id, t_offsets_unc[crys_id - 1]);
520 if (numEntries > 0) {
521 peakClusterTimes->Fill(t_offsets[crys_id - 1]);
522 numCrysWithNonZeroEntries++ ;
524 if ((numEntries >= minNumEntries) && good_fit) {
525 peakClusterTimesGoodFit->Fill(t_offsets[crys_id - 1]);
526 peakClusterTimesGoodFit__cid->SetBinContent(crys_id, t_offsets[crys_id - 1]);
527 peakClusterTimesGoodFit__cid->SetBinError(crys_id, t_offsets_unc[crys_id - 1]);
532 double numClusterTimesWithinWindowFraction = h_time->Integral(binyLeft, binyRight) ;
533 double clusterTimesWithinWindowFraction = numClusterTimesWithinWindowFraction;
534 if (numEntries > 0) {
535 clusterTimesWithinWindowFraction /= numEntries;
537 clusterTimesWithinWindowFraction = -0.1;
540 B2INFO(
"Crystal cell id = " << crys_id <<
", theta id = " <<
541 thetaID <<
", clusterTimesWithinWindowFraction = " <<
542 numClusterTimesWithinWindowFraction <<
" / " << numEntries <<
" = " <<
543 clusterTimesWithinWindowFraction);
545 clusterTimeFractionInWindow->Fill(clusterTimesWithinWindowFraction);
546 clusterTimeNumberInWindow__cid->SetBinContent(crys_id, numClusterTimesWithinWindowFraction);
547 clusterTimeNumber__cid->SetBinContent(crys_id, numEntries);
549 if ((thetaID >= 3 && thetaID <= 10) ||
550 (thetaID >= 15 && thetaID <= 39) ||
551 (thetaID >= 44 && thetaID <= 56) ||
552 (thetaID >= 61 && thetaID <= 66)) {
553 clusterTimeNumberInWindowInGoodECLRings__cid->SetBinContent(crys_id, numClusterTimesWithinWindowFraction);
561 auto g_clusterTimeFractionInWindow__cid =
new TGraphAsymmErrors(clusterTimeNumberInWindow__cid, clusterTimeNumber__cid,
"w");
562 auto g_clusterTimeFractionInWindowInGoodECLRings__cid =
new TGraphAsymmErrors(clusterTimeNumberInWindowInGoodECLRings__cid,
563 clusterTimeNumber__cid,
"w");
564 g_clusterTimeFractionInWindow__cid->SetTitle(fracWindowTitle);
565 g_clusterTimeFractionInWindowInGoodECLRings__cid->SetTitle(fracWindowInGoodECLRingsTitle);
568 peakClusterTime_cid->ResetStats();
569 peakClusterTimesGoodFit__cid->ResetStats();
571 histfile->WriteTObject(peakClusterTime_cid,
"peakClusterTime_cid");
572 histfile->WriteTObject(peakClusterTimes,
"peakClusterTimes");
573 histfile->WriteTObject(peakClusterTimesGoodFit__cid,
"peakClusterTimesGoodFit__cid");
574 histfile->WriteTObject(peakClusterTimesGoodFit,
"peakClusterTimesGoodFit");
575 histfile->WriteTObject(g_clusterTimeFractionInWindow__cid,
"g_clusterTimeFractionInWindow__cid");
576 histfile->WriteTObject(g_clusterTimeFractionInWindowInGoodECLRings__cid,
"g_clusterTimeFractionInWindowInGoodECLRings__cid");
577 histfile->WriteTObject(clusterTimeFractionInWindow,
"clusterTimeFractionInWindow");
584 vector <int> binProjectionLeft = binProjectionLeft_Time_vs_E_runDep;
585 vector <int> binProjectionRight = binProjectionRight_Time_vs_E_runDep;
587 auto h2 = clusterTimeClusterE;
590 double max_E = h2->GetXaxis()->GetXmax();
593 vector <double> E_binEdges(binProjectionLeft.size() + 1);
594 for (
long unsigned int x_bin = 0; x_bin < binProjectionLeft.size(); x_bin++) {
595 TH1D* h_E_t_slice = h2->ProjectionX(
"h_E_t_slice", 1, 1) ;
596 E_binEdges[x_bin] = h_E_t_slice->GetXaxis()->GetBinLowEdge(binProjectionLeft[x_bin]) ;
597 B2INFO(
"E_binEdges[" << x_bin <<
"] = " << E_binEdges[x_bin]);
598 if (x_bin == binProjectionLeft.size() - 1) {
599 E_binEdges[x_bin + 1] = max_E ;
600 B2INFO(
"E_binEdges[" << x_bin + 1 <<
"] = " << E_binEdges[x_bin + 1]);
605 auto clusterTimePeak_ClusterEnergy_varBin =
new TH1F(
"clusterTimePeak_ClusterEnergy_varBin",
606 ";ECL cluster energy [GeV];Cluster time fit position [ns]", E_binEdges.size() - 1, &(E_binEdges[0]));
607 auto clusterTimePeakWidth_ClusterEnergy_varBin =
new TH1F(
"clusterTimePeakWidth_ClusterEnergy_varBin",
608 ";ECL cluster energy [GeV];Cluster time fit width [ns]", E_binEdges.size() - 1, &(E_binEdges[0]));
610 int Ebin_counter = 1 ;
613 for (
long unsigned int x_bin = 0; x_bin < binProjectionLeft.size(); x_bin++) {
614 double clusterTime_mean = 0;
615 double clusterTime_mean_unc = 0;
616 double clusterTime_sigma = 0;
618 B2INFO(
"x_bin = " << x_bin);
622 TH1D* h_time = h2->ProjectionY((
"h_time_E_slice_" + std::to_string(x_bin)).c_str(), binProjectionLeft[x_bin],
623 binProjectionRight[x_bin]) ;
626 TH1D* h_E_t_slice = h2->ProjectionX(
"h_E_t_slice", 1, 1) ;
627 double lowE = h_E_t_slice->GetXaxis()->GetBinLowEdge(binProjectionLeft[x_bin]) ;
628 double highE = h_E_t_slice->GetXaxis()->GetBinUpEdge(binProjectionRight[x_bin]) ;
629 double meanE = (lowE + highE) / 2.0 ;
631 B2INFO(
"bin " << Ebin_counter <<
": low E = " << lowE <<
", high E = " << highE <<
" GeV");
633 TH1D* h_timeMask = (TH1D*)h_time->Clone();
634 TH1D* h_timeMasked = (TH1D*)h_time->Clone((std::string(
"h_time_E_slice_masked__") + std::to_string(meanE)).c_str());
635 TH1D* h_timeRebin = (TH1D*)h_time->Clone();
642 h_timeMask->Scale(0.0);
644 time_fit_min = hist_tmax;
645 time_fit_max = hist_tmin;
648 double histRebin_max = h_timeRebin->GetMaximum();
650 bool maskedOutNonZeroBin =
false;
652 for (
int bin = 1; bin <= h_timeRebin->GetNbinsX(); bin++) {
655 if (nonRebinnedBinNumber < h_time->GetNbinsX()) {
657 h_timeMask->SetBinContent(nonRebinnedBinNumber, 1);
660 double x_lower = h_timeRebin->GetXaxis()->GetBinLowEdge(bin);
661 double x_upper = h_timeRebin->GetXaxis()->GetBinUpEdge(bin);
662 if (x_lower < time_fit_min) {
663 time_fit_min = x_lower;
665 if (x_upper > time_fit_max) {
666 time_fit_max = x_upper;
670 if (h_time->GetBinContent(nonRebinnedBinNumber) > 0) {
671 B2DEBUG(22,
"Setting bin " << nonRebinnedBinNumber <<
" from " << h_timeMasked->GetBinContent(nonRebinnedBinNumber) <<
" to 0");
672 maskedOutNonZeroBin =
true;
674 h_timeMasked->SetBinContent(nonRebinnedBinNumber, 0);
679 B2INFO(
"Bins with non-zero values have been masked out: " << maskedOutNonZeroBin);
680 h_timeMasked->ResetStats();
681 h_timeMask->ResetStats();
687 double default_meanMasked = h_timeMasked->GetMean();
689 B2INFO(
"default_meanMasked = " << default_meanMasked);
693 double default_mean = h_time->GetMean();
694 double default_mean_unc = h_time->GetMeanError();
695 double default_sigma = h_time->GetStdDev();
697 B2INFO(
"Fitting crystal between " << time_fit_min <<
" and " << time_fit_max);
700 TF1* gaus =
new TF1(
"func",
"gaus(0)", time_fit_min, time_fit_max);
701 gaus->SetParNames(
"numCrystalHitsNormalization",
"mean",
"sigma");
708 double hist_max = h_time->GetMaximum();
711 double stddev = h_time->GetStdDev();
716 gaus->SetParameter(0, hist_max / 2.);
717 gaus->SetParameter(1, mean);
718 gaus->SetParameter(2, sigma);
725 h_timeMasked->Fit(gaus,
"LQR");
727 double fit_mean = gaus->GetParameter(1);
728 double fit_mean_unc = gaus->GetParError(1);
729 double fit_sigma = gaus->GetParameter(2);
731 double meanDiff = fit_mean - default_mean;
732 double meanUncDiff = fit_mean_unc - default_mean_unc;
733 double sigmaDiff = fit_sigma - default_sigma;
735 bool good_fit =
false;
737 if ((fabs(meanDiff) > 10) ||
738 (fabs(meanUncDiff) > 10) ||
739 (fabs(sigmaDiff) > 10) ||
740 (fit_mean_unc > 3) ||
742 (fit_mean < time_fit_min) ||
743 (fit_mean > time_fit_max)) {
744 B2INFO(
"x_bin = " << x_bin);
745 B2INFO(
"fit mean, default mean = " << fit_mean <<
", " << default_mean);
746 B2INFO(
"fit mean unc, default mean unc = " << fit_mean_unc <<
", " << default_mean_unc);
747 B2INFO(
"fit sigma, default sigma = " << fit_sigma <<
", " << default_sigma);
749 B2INFO(
"crystal fit mean - hist mean = " << meanDiff);
750 B2INFO(
"fit mean unc. - hist mean unc. = " << meanUncDiff);
751 B2INFO(
"fit sigma - hist sigma = " << sigmaDiff);
753 B2INFO(
"fit_mean = " << fit_mean);
754 B2INFO(
"time_fit_min = " << time_fit_min);
755 B2INFO(
"time_fit_max = " << time_fit_max);
757 if (fabs(meanDiff) > 10) B2INFO(
"fit mean diff too large");
758 if (fabs(meanUncDiff) > 10) B2INFO(
"fit mean unc diff too large");
759 if (fabs(sigmaDiff) > 10) B2INFO(
"fit mean sigma diff too large");
760 if (fit_mean_unc > 3) B2INFO(
"fit mean unc too large");
761 if (fit_sigma < 0.1) B2INFO(
"fit sigma too small");
768 int numEntries = h_time->GetEntries();
772 if ((numEntries >= minNumEntries) && good_fit) {
773 clusterTime_mean = fit_mean;
774 clusterTime_mean_unc = fit_mean_unc;
775 clusterTime_sigma = fit_sigma;
777 clusterTime_mean = default_mean;
778 clusterTime_mean_unc = default_mean_unc;
779 clusterTime_sigma = default_sigma;
782 if (numEntries < minNumEntries) B2INFO(
"Number of entries less than minimum");
783 if (numEntries == 0) B2INFO(
"Number of entries == 0");
785 histfile->WriteTObject(h_time, (std::string(
"h_time_E_slice") + std::to_string(meanE)).c_str());
786 histfile->WriteTObject(h_timeMasked, (std::string(
"h_time_E_slice_masked") + std::to_string(meanE)).c_str());
789 clusterTimePeak_ClusterEnergy_varBin->SetBinContent(Ebin_counter, clusterTime_mean);
790 clusterTimePeak_ClusterEnergy_varBin->SetBinError(Ebin_counter, clusterTime_mean_unc);
792 clusterTimePeakWidth_ClusterEnergy_varBin->SetBinContent(Ebin_counter, clusterTime_sigma);
793 clusterTimePeakWidth_ClusterEnergy_varBin->SetBinError(Ebin_counter, 0);
809 vector< pair<double, int> > fitClusterTime__crystalIDBase0__pairs;
813 fitClusterTime__crystalIDBase0__pairs.push_back(make_pair(0.0, cid));
818 fitClusterTime__crystalIDBase0__pairs[crys_id - 1] = make_pair(fabs(t_offsets[crys_id - 1]), crys_id - 1) ;
822 sort(fitClusterTime__crystalIDBase0__pairs.begin(), fitClusterTime__crystalIDBase0__pairs.end());
826 B2INFO(
"-------- List of the (fitted) peak cluster times sorted by their absolute value ----------");
827 B2INFO(
"------------------------------------------------------------------------------------------");
828 B2INFO(
"------------------------------------------------------------------------------------------");
829 B2INFO(
"Quoted # of clusters is before the cutting off of the distribution tails, cellID=1..ECLElementNumbers::c_NCrystals, crysID=0..8735");
831 bool hasHitThresholdBadTimes = false ;
833 int cid = fitClusterTime__crystalIDBase0__pairs[iSortedTimes].second ;
834 if (!hasHitThresholdBadTimes && fitClusterTime__crystalIDBase0__pairs[iSortedTimes].first > 2) {
835 B2INFO(
"======== |t_fit| > Xns threshold ======");
836 hasHitThresholdBadTimes =
true;
839 B2INFO(
"cid = " << cid <<
", peak clust t = " << t_offsets[cid] <<
" +- " << t_offsets_unc[cid] <<
" ns, # clust = " <<
840 numClusterPerCrys[cid] <<
", good fit = " << crysHasGoodFit[cid] <<
", good fit & stats = " << crysHasGoodFitandStats[cid]);
848 B2INFO(
"######## List of poor (fitted) peak cluster times sorted by their absolute value #########");
849 B2INFO(
"##########################################################################################");
850 B2INFO(
"##########################################################################################");
853 int cid = fitClusterTime__crystalIDBase0__pairs[iSortedTimes].second ;
854 if (fitClusterTime__crystalIDBase0__pairs[iSortedTimes].first > 2 && crysHasGoodFitandStats[cid]) {
855 B2INFO(
"WARNING: cid = " << cid <<
", peak clust t = " << t_offsets[cid] <<
" +- " << t_offsets_unc[cid] <<
" ns, # clust = " <<
856 numClusterPerCrys[cid] <<
", good fit = " << crysHasGoodFit[cid] <<
", good fit & stats = " << crysHasGoodFitandStats[cid]);
862 B2INFO(
"Number of crystals with non-zero number of hits = " << numCrysWithNonZeroEntries);
863 B2INFO(
"Number of crystals with good quality fits = " << numCrysWithGoodFit);
866 clusterTimePeak_ClusterEnergy_varBin->ResetStats();
867 clusterTimePeakWidth_ClusterEnergy_varBin->ResetStats();
869 histfile->WriteTObject(clusterTimePeak_ClusterEnergy_varBin,
"clusterTimePeak_ClusterEnergy_varBin");
870 histfile->WriteTObject(clusterTimePeakWidth_ClusterEnergy_varBin,
"clusterTimePeakWidth_ClusterEnergy_varBin");
876 B2INFO(
"Filling histograms for difference in crystal payload values and the pre-calibration values. These older values may be from a previous bucket or older reprocessing of the data.");
878 double tsDiffCustomOld_ns = -999;
880 tsDiffCustomOld_ns = (currentValuesCrys[crys_id - 1] - prevValuesCrys[crys_id - 1]) * TICKS_TO_NS;
881 B2INFO(
"Crystal " << crys_id <<
": ts new merged - 'before 1st iter' merged = (" <<
882 currentValuesCrys[crys_id - 1] <<
" - " << prevValuesCrys[crys_id - 1] <<
883 ") ticks * " << TICKS_TO_NS <<
" ns/tick = " << tsDiffCustomOld_ns <<
" ns");
886 tsNew_MINUS_tsCustomPrev__cid->SetBinContent(crys_id, tsDiffCustomOld_ns);
887 tsNew_MINUS_tsCustomPrev__cid->SetBinError(crys_id, 0);
888 tsNew_MINUS_tsCustomPrev__cid->ResetStats();
890 tsNew_MINUS_tsCustomPrev->Fill(tsDiffCustomOld_ns);
891 tsNew_MINUS_tsCustomPrev->ResetStats();
894 histfile->WriteTObject(tsNew_MINUS_tsCustomPrev__cid,
"tsNew_MINUS_tsCustomPrev__cid");
895 histfile->WriteTObject(tsNew_MINUS_tsCustomPrev,
"tsNew_MINUS_tsCustomPrev");
900 B2INFO(
"Finished validations algorithm");
Base class for calibration algorithms.
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_Failure
Failed =3 in Python.
Class for accessing objects in the database.
Singleton class to cache database objects.
static DataStore & Instance()
Instance of singleton Store.
void setInitializeActive(bool active)
Setter for m_initializeActive.
This class provides access to ECL channel map that is either a) Loaded from the database (see ecl/dbo...
The Class for ECL Geometry Parameters.
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
void Mapping(int cid)
Mapping theta, phi Id.
int GetThetaID()
Get Theta Id.
static double getRF()
See m_rf.
eclTValidationAlgorithm()
..Constructor
int cellIDHi
Fit crystals with cellID0 in the inclusive range [cellIDLo,cellIDHi].
int cellIDLo
Fit crystals with cellID0 in the inclusive range [cellIDLo,cellIDHi].
double meanCleanRebinFactor
Rebinning factor for mean calculation.
double clusterTimesFractionWindow_maxtime
Maximum time for window to calculate cluster time fraction, in ns.
double meanCleanCutMinFactor
After rebinning, create a mask for bins that have values less than meanCleanCutMinFactor times the ma...
bool readPrevCrysPayload
Read the previous crystal payload values for comparison.
EResult calibrate() override
..Run algorithm on events
std::string debugFilenameBase
Name of file with debug output, eclTValidationAlgorithm.root by default.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Type-safe access to single objects in the data store.
bool construct(Args &&... params)
Construct an object of type T in this StoreObjPtr, using the provided constructor arguments.
static DBStore & Instance()
Instance of a singleton DBStore.
void updateEvent()
Updates all intra-run dependent objects.
void update()
Updates all objects that are outside their interval of validity.
const int c_NCrystals
Number of crystals.
Abstract base class for different kinds of events.