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 clusterTimeE0E1diff = getObjectPtr<TH1F>(
"clusterTimeE0E1diff");
107 auto cutflow = getObjectPtr<TH1F>(
"cutflow");
109 vector <int> binProjectionLeft_Time_vs_E_runDep ;
110 vector <int> binProjectionRight_Time_vs_E_runDep ;
112 for (
int binCounter = 1; binCounter <= clusterTimeClusterE->GetNbinsX(); binCounter++) {
113 binProjectionLeft_Time_vs_E_runDep.push_back(binCounter);
114 binProjectionRight_Time_vs_E_runDep.push_back(binCounter);
129 int numCrysWithNonZeroEntries = 0 ;
130 int numCrysWithGoodFit = 0 ;
132 int minNumEntries = 40;
138 bool minRunNumBool =
false;
139 bool maxRunNumBool =
false;
145 int expNumber = expRun.first;
146 int runNumber = expRun.second;
147 if (!minRunNumBool) {
148 minExpNum = expNumber;
149 minRunNum = runNumber;
150 minRunNumBool =
true;
152 if (!maxRunNumBool) {
153 maxExpNum = expNumber;
154 maxRunNum = runNumber;
155 maxRunNumBool =
true;
157 if (((minRunNum > runNumber) && (minExpNum >= expNumber)) ||
158 (minExpNum > expNumber)) {
159 minExpNum = expNumber;
160 minRunNum = runNumber;
162 if (((maxRunNum < runNumber) && (maxExpNum <= expNumber)) ||
163 (maxExpNum < expNumber))
166 maxExpNum = expNumber;
167 maxRunNum = runNumber;
172 string runNumsString = string(
"_") + to_string(minExpNum) +
"_" + to_string(minRunNum) + string(
"-") +
173 to_string(maxExpNum) +
"_" + to_string(maxRunNum);
181 int eventNumberForCrates = 1;
193 evtPtr.
construct(eventNumberForCrates, minRunNum, minExpNum);
202 B2INFO(
"Uploading payload for exp " << minExpNum <<
", run " << minRunNum <<
", event " << eventNumberForCrates);
205 crystalMapper->initFromDB();
214 B2INFO(
"Dumping payload");
217 std::vector<float> currentValuesCrys = crystalTimeObject->getCalibVector();
218 std::vector<float> currentUncCrys = crystalTimeObject->getCalibUncVector();
221 B2INFO(
"Values read from database. Write out for their values for comparison against those from tcol");
223 B2INFO(
"ts: cellID " << ic + 1 <<
" " << currentValuesCrys[ic] <<
" +/- " << currentUncCrys[ic]);
232 prevValuesCrys = customPrevCrystalTimeObject->getCalibVector();
235 B2INFO(
"Previous values read from database. Write out for their values for comparison against those from tcol");
237 B2INFO(
"ts custom previous payload: cellID " << ic + 1 <<
" " << prevValuesCrys[ic]);
245 B2INFO(
"Debug output rootfile: " << debugFilename);
246 histfile =
new TFile(debugFilename.c_str(),
"recreate");
249 clusterTime ->Write();
250 clusterTime_cid ->Write();
251 clusterTime_run ->Write();
252 clusterTimeClusterE ->Write();
253 dt99_clusterE ->Write();
255 clusterTimeE0E1diff ->Write();
260 double hist_tmin = clusterTime->GetXaxis()->GetXmin();
261 double hist_tmax = clusterTime->GetXaxis()->GetXmax();
262 int hist_nTbins = clusterTime->GetNbinsX();
264 B2INFO(
"hist_tmin = " << hist_tmin);
265 B2INFO(
"hist_tmax = " << hist_tmax);
266 B2INFO(
"hist_nTbins = " << hist_nTbins);
268 double time_fit_min = hist_tmax;
269 double time_fit_max = hist_tmin;
275 auto peakClusterTimes =
new TH1F(
"peakClusterTimes",
276 "-For crystals with at least one hit-;Peak cluster time [ns];Number of crystals",
277 hist_nTbins, hist_tmin, hist_tmax);
278 auto peakClusterTimesGoodFit =
new TH1F(
"peakClusterTimesGoodFit",
279 "-For crystals with a good fit to distribution of hits-;Peak cluster time [ns];Number of crystals",
280 hist_nTbins, hist_tmin, hist_tmax);
282 auto peakClusterTimesGoodFit__cid =
new TH1F(
"peakClusterTimesGoodFit__cid",
283 "-For crystals with a good fit to distribution of hits-;cell id (only crystals with good fit);Peak cluster time [ns]",
288 auto tsNew_MINUS_tsCustomPrev__cid =
new TH1F(
"TsNew_MINUS_TsCustomPrev__cid",
289 ";cell id; ts(new|merged) - ts(old = 'pre-calib'|merged) [ns]",
292 auto tsNew_MINUS_tsCustomPrev =
new TH1F(
"TsNew_MINUS_TsCustomPrev",
293 ";ts(new | merged) - ts(old = 'pre-calib' | merged) [ns];Number of crystals",
294 285, -69.5801, 69.5801);
300 B2INFO(
"timeWindow_maxTime = " << timeWindow_maxTime);
301 int binyLeft = clusterTime_cid->GetYaxis()->FindBin(-timeWindow_maxTime);
302 int binyRight = clusterTime_cid->GetYaxis()->FindBin(timeWindow_maxTime);
303 double windowLowTimeFromBin = clusterTime_cid->GetYaxis()->GetBinLowEdge(binyLeft);
304 double windowHighTimeFromBin = clusterTime_cid->GetYaxis()->GetBinLowEdge(binyRight + 1);
305 std::string s_lowTime = std::to_string(windowLowTimeFromBin);
306 std::string s_highTime = std::to_string(windowHighTimeFromBin);
307 TString fracWindowTitle =
"Fraction of cluster times in window [" + s_lowTime +
", " + s_highTime +
308 "] ns;cell id;Fraction of cluster times in window";
309 B2INFO(
"fracWindowTitle = " << fracWindowTitle);
310 TString fracWindowInGoodECLRingsTitle =
"Fraction of cluster times in window [" + s_lowTime +
", " + s_highTime +
311 "] ns and in good ECL theta rings;cell id;Fraction cluster times in window + good ECL rings";
312 B2INFO(
"fracWindowInGoodECLRingsTitle = " << fracWindowInGoodECLRingsTitle);
313 B2INFO(
"Good ECL rings skip gaps in the acceptance, and includes ECL theta IDs: 3-10, 15-39, 44-56, 61-66.");
315 TString fracWindowHistTitle =
"Fraction of cluster times in window [" + s_lowTime +
", " + s_highTime +
316 "] ns;Fraction of cluster times in window;Number of crystals";
320 auto clusterTimeNumberInWindowInGoodECLRings__cid =
new TH1F(
"clusterTimeNumberInWindowInGoodECLRings__cid", fracWindowTitle,
325 auto clusterTimeFractionInWindow =
new TH1F(
"clusterTimeFractionInWindow", fracWindowHistTitle, 110, 0.0, 1.1);
327 clusterTimeNumberInWindow__cid->Sumw2();
328 clusterTimeNumberInWindowInGoodECLRings__cid->Sumw2();
329 clusterTimeNumber__cid->Sumw2();
339 double clusterTime_mean = 0;
340 double clusterTime_mean_unc = 0;
342 B2INFO(
"Crystal cell id = " << crys_id);
351 TH1D* h_time = clusterTime_cid->ProjectionY((std::string(
"h_time_psi__") + std::to_string(crys_id)).c_str(),
353 TH1D* h_timeMask = (TH1D*)h_time->Clone();
354 TH1D* h_timeMasked = (TH1D*)h_time->Clone((std::string(
"h_time_psi_masked__") + std::to_string(crys_id)).c_str());
355 TH1D* h_timeRebin = (TH1D*)h_time->Clone();
362 h_timeMask->Scale(0.0);
364 time_fit_min = hist_tmax;
365 time_fit_max = hist_tmin;
368 double histRebin_max = h_timeRebin->GetMaximum();
370 bool maskedOutNonZeroBin =
false;
372 for (
int bin = 1; bin <= h_timeRebin->GetNbinsX(); bin++) {
375 if (nonRebinnedBinNumber < h_time->GetNbinsX()) {
377 h_timeMask->SetBinContent(nonRebinnedBinNumber, 1);
380 double x_lower = h_timeRebin->GetXaxis()->GetBinLowEdge(bin);
381 double x_upper = h_timeRebin->GetXaxis()->GetBinUpEdge(bin);
382 if (x_lower < time_fit_min) {
383 time_fit_min = x_lower;
385 if (x_upper > time_fit_max) {
386 time_fit_max = x_upper;
390 if (h_time->GetBinContent(nonRebinnedBinNumber) > 0) {
391 B2DEBUG(22,
"Setting bin " << nonRebinnedBinNumber <<
" from " << h_timeMasked->GetBinContent(nonRebinnedBinNumber) <<
" to 0");
392 maskedOutNonZeroBin =
true;
394 h_timeMasked->SetBinContent(nonRebinnedBinNumber, 0);
399 B2INFO(
"Bins with non-zero values have been masked out: " << maskedOutNonZeroBin);
400 h_timeMasked->ResetStats();
401 h_timeMask->ResetStats();
406 double default_meanMasked = h_timeMasked->GetMean();
408 B2INFO(
"default_meanMasked = " << default_meanMasked);
412 double default_mean = h_time->GetMean();
413 double default_mean_unc = h_time->GetMeanError();
414 double default_sigma = h_time->GetStdDev();
416 B2INFO(
"Fitting crystal between " << time_fit_min <<
" and " << time_fit_max);
419 TF1* gaus =
new TF1(
"func",
"gaus(0)", time_fit_min, time_fit_max);
420 gaus->SetParNames(
"numCrystalHitsNormalization",
"mean",
"sigma");
427 double hist_max = h_time->GetMaximum();
430 double stddev = h_time->GetStdDev();
435 gaus->SetParameter(0, hist_max / 2.);
436 gaus->SetParameter(1, mean);
437 gaus->SetParameter(2, sigma);
444 h_timeMasked->Fit(gaus,
"LQR");
446 double fit_mean = gaus->GetParameter(1);
447 double fit_mean_unc = gaus->GetParError(1);
448 double fit_sigma = gaus->GetParameter(2);
450 double meanDiff = fit_mean - default_mean;
451 double meanUncDiff = fit_mean_unc - default_mean_unc;
452 double sigmaDiff = fit_sigma - default_sigma;
454 bool good_fit =
false;
456 if ((fabs(meanDiff) > 10) ||
457 (fabs(meanUncDiff) > 10) ||
458 (fabs(sigmaDiff) > 10) ||
459 (fit_mean_unc > 3) ||
461 (fit_mean < time_fit_min) ||
462 (fit_mean > time_fit_max)) {
463 B2INFO(
"Crystal cell id = " << crys_id);
464 B2INFO(
"fit mean, default mean = " << fit_mean <<
", " << default_mean);
465 B2INFO(
"fit mean unc, default mean unc = " << fit_mean_unc <<
", " << default_mean_unc);
466 B2INFO(
"fit sigma, default sigma = " << fit_sigma <<
", " << default_sigma);
468 B2INFO(
"crystal fit mean - hist mean = " << meanDiff);
469 B2INFO(
"fit mean unc. - hist mean unc. = " << meanUncDiff);
470 B2INFO(
"fit sigma - hist sigma = " << sigmaDiff);
472 B2INFO(
"fit_mean = " << fit_mean);
473 B2INFO(
"time_fit_min = " << time_fit_min);
474 B2INFO(
"time_fit_max = " << time_fit_max);
476 if (fabs(meanDiff) > 10) B2INFO(
"fit mean diff too large");
477 if (fabs(meanUncDiff) > 10) B2INFO(
"fit mean unc diff too large");
478 if (fabs(sigmaDiff) > 10) B2INFO(
"fit mean sigma diff too large");
479 if (fit_mean_unc > 3) B2INFO(
"fit mean unc too large");
480 if (fit_sigma < 0.1) B2INFO(
"fit sigma too small");
484 numCrysWithGoodFit++;
485 crysHasGoodFit[crys_id - 1] = true ;
489 int numEntries = h_time->GetEntries();
492 if ((numEntries >= minNumEntries) && good_fit) {
493 clusterTime_mean = fit_mean;
494 clusterTime_mean_unc = fit_mean_unc;
495 crysHasGoodFitandStats[crys_id - 1] = true ;
497 clusterTime_mean = default_mean;
498 clusterTime_mean_unc = default_mean_unc;
501 if (numEntries < minNumEntries) B2INFO(
"Number of entries less than minimum");
502 if (numEntries == 0) B2INFO(
"Number of entries == 0");
505 t_offsets[crys_id - 1] = clusterTime_mean ;
506 t_offsets_unc[crys_id - 1] = clusterTime_mean_unc ;
507 numClusterPerCrys[crys_id - 1] = numEntries ;
509 histfile->WriteTObject(h_time, (std::string(
"h_time_psi") + std::to_string(crys_id)).c_str());
510 histfile->WriteTObject(h_timeMasked, (std::string(
"h_time_psi_masked") + std::to_string(crys_id)).c_str());
513 peakClusterTime_cid->SetBinContent(crys_id, t_offsets[crys_id - 1]);
514 peakClusterTime_cid->SetBinError(crys_id, t_offsets_unc[crys_id - 1]);
518 if (numEntries > 0) {
519 peakClusterTimes->Fill(t_offsets[crys_id - 1]);
520 numCrysWithNonZeroEntries++ ;
522 if ((numEntries >= minNumEntries) && good_fit) {
523 peakClusterTimesGoodFit->Fill(t_offsets[crys_id - 1]);
524 peakClusterTimesGoodFit__cid->SetBinContent(crys_id, t_offsets[crys_id - 1]);
525 peakClusterTimesGoodFit__cid->SetBinError(crys_id, t_offsets_unc[crys_id - 1]);
530 double numClusterTimesWithinWindowFraction = h_time->Integral(binyLeft, binyRight) ;
531 double clusterTimesWithinWindowFraction = numClusterTimesWithinWindowFraction;
532 if (numEntries > 0) {
533 clusterTimesWithinWindowFraction /= numEntries;
535 clusterTimesWithinWindowFraction = -0.1;
538 B2INFO(
"Crystal cell id = " << crys_id <<
", theta id = " <<
539 thetaID <<
", clusterTimesWithinWindowFraction = " <<
540 numClusterTimesWithinWindowFraction <<
" / " << numEntries <<
" = " <<
541 clusterTimesWithinWindowFraction);
543 clusterTimeFractionInWindow->Fill(clusterTimesWithinWindowFraction);
544 clusterTimeNumberInWindow__cid->SetBinContent(crys_id, numClusterTimesWithinWindowFraction);
545 clusterTimeNumber__cid->SetBinContent(crys_id, numEntries);
547 if ((thetaID >= 3 && thetaID <= 10) ||
548 (thetaID >= 15 && thetaID <= 39) ||
549 (thetaID >= 44 && thetaID <= 56) ||
550 (thetaID >= 61 && thetaID <= 66)) {
551 clusterTimeNumberInWindowInGoodECLRings__cid->SetBinContent(crys_id, numClusterTimesWithinWindowFraction);
559 auto g_clusterTimeFractionInWindow__cid =
new TGraphAsymmErrors(clusterTimeNumberInWindow__cid, clusterTimeNumber__cid,
"w");
560 auto g_clusterTimeFractionInWindowInGoodECLRings__cid =
new TGraphAsymmErrors(clusterTimeNumberInWindowInGoodECLRings__cid,
561 clusterTimeNumber__cid,
"w");
562 g_clusterTimeFractionInWindow__cid->SetTitle(fracWindowTitle);
563 g_clusterTimeFractionInWindowInGoodECLRings__cid->SetTitle(fracWindowInGoodECLRingsTitle);
566 peakClusterTime_cid->ResetStats();
567 peakClusterTimesGoodFit__cid->ResetStats();
569 histfile->WriteTObject(peakClusterTime_cid,
"peakClusterTime_cid");
570 histfile->WriteTObject(peakClusterTimes,
"peakClusterTimes");
571 histfile->WriteTObject(peakClusterTimesGoodFit__cid,
"peakClusterTimesGoodFit__cid");
572 histfile->WriteTObject(peakClusterTimesGoodFit,
"peakClusterTimesGoodFit");
573 histfile->WriteTObject(g_clusterTimeFractionInWindow__cid,
"g_clusterTimeFractionInWindow__cid");
574 histfile->WriteTObject(g_clusterTimeFractionInWindowInGoodECLRings__cid,
"g_clusterTimeFractionInWindowInGoodECLRings__cid");
575 histfile->WriteTObject(clusterTimeFractionInWindow,
"clusterTimeFractionInWindow");
582 vector <int> binProjectionLeft = binProjectionLeft_Time_vs_E_runDep;
583 vector <int> binProjectionRight = binProjectionRight_Time_vs_E_runDep;
585 auto h2 = clusterTimeClusterE;
588 double max_E = h2->GetXaxis()->GetXmax();
591 vector <double> E_binEdges(binProjectionLeft.size() + 1);
592 for (
long unsigned int x_bin = 0; x_bin < binProjectionLeft.size(); x_bin++) {
593 TH1D* h_E_t_slice = h2->ProjectionX(
"h_E_t_slice", 1, 1) ;
594 E_binEdges[x_bin] = h_E_t_slice->GetXaxis()->GetBinLowEdge(binProjectionLeft[x_bin]) ;
595 B2INFO(
"E_binEdges[" << x_bin <<
"] = " << E_binEdges[x_bin]);
596 if (x_bin == binProjectionLeft.size() - 1) {
597 E_binEdges[x_bin + 1] = max_E ;
598 B2INFO(
"E_binEdges[" << x_bin + 1 <<
"] = " << E_binEdges[x_bin + 1]);
603 auto clusterTimePeak_ClusterEnergy_varBin =
new TH1F(
"clusterTimePeak_ClusterEnergy_varBin",
604 ";ECL cluster energy [GeV];Cluster time fit position [ns]", E_binEdges.size() - 1, &(E_binEdges[0]));
605 auto clusterTimePeakWidth_ClusterEnergy_varBin =
new TH1F(
"clusterTimePeakWidth_ClusterEnergy_varBin",
606 ";ECL cluster energy [GeV];Cluster time fit width [ns]", E_binEdges.size() - 1, &(E_binEdges[0]));
608 int Ebin_counter = 1 ;
611 for (
long unsigned int x_bin = 0; x_bin < binProjectionLeft.size(); x_bin++) {
612 double clusterTime_mean = 0;
613 double clusterTime_mean_unc = 0;
614 double clusterTime_sigma = 0;
616 B2INFO(
"x_bin = " << x_bin);
620 TH1D* h_time = h2->ProjectionY((
"h_time_E_slice_" + std::to_string(x_bin)).c_str(), binProjectionLeft[x_bin],
621 binProjectionRight[x_bin]) ;
624 TH1D* h_E_t_slice = h2->ProjectionX(
"h_E_t_slice", 1, 1) ;
625 double lowE = h_E_t_slice->GetXaxis()->GetBinLowEdge(binProjectionLeft[x_bin]) ;
626 double highE = h_E_t_slice->GetXaxis()->GetBinUpEdge(binProjectionRight[x_bin]) ;
627 double meanE = (lowE + highE) / 2.0 ;
629 B2INFO(
"bin " << Ebin_counter <<
": low E = " << lowE <<
", high E = " << highE <<
" GeV");
631 TH1D* h_timeMask = (TH1D*)h_time->Clone();
632 TH1D* h_timeMasked = (TH1D*)h_time->Clone((std::string(
"h_time_E_slice_masked__") + std::to_string(meanE)).c_str());
633 TH1D* h_timeRebin = (TH1D*)h_time->Clone();
640 h_timeMask->Scale(0.0);
642 time_fit_min = hist_tmax;
643 time_fit_max = hist_tmin;
646 double histRebin_max = h_timeRebin->GetMaximum();
648 bool maskedOutNonZeroBin =
false;
650 for (
int bin = 1; bin <= h_timeRebin->GetNbinsX(); bin++) {
653 if (nonRebinnedBinNumber < h_time->GetNbinsX()) {
655 h_timeMask->SetBinContent(nonRebinnedBinNumber, 1);
658 double x_lower = h_timeRebin->GetXaxis()->GetBinLowEdge(bin);
659 double x_upper = h_timeRebin->GetXaxis()->GetBinUpEdge(bin);
660 if (x_lower < time_fit_min) {
661 time_fit_min = x_lower;
663 if (x_upper > time_fit_max) {
664 time_fit_max = x_upper;
668 if (h_time->GetBinContent(nonRebinnedBinNumber) > 0) {
669 B2DEBUG(22,
"Setting bin " << nonRebinnedBinNumber <<
" from " << h_timeMasked->GetBinContent(nonRebinnedBinNumber) <<
" to 0");
670 maskedOutNonZeroBin =
true;
672 h_timeMasked->SetBinContent(nonRebinnedBinNumber, 0);
677 B2INFO(
"Bins with non-zero values have been masked out: " << maskedOutNonZeroBin);
678 h_timeMasked->ResetStats();
679 h_timeMask->ResetStats();
685 double default_meanMasked = h_timeMasked->GetMean();
687 B2INFO(
"default_meanMasked = " << default_meanMasked);
691 double default_mean = h_time->GetMean();
692 double default_mean_unc = h_time->GetMeanError();
693 double default_sigma = h_time->GetStdDev();
695 B2INFO(
"Fitting crystal between " << time_fit_min <<
" and " << time_fit_max);
698 TF1* gaus =
new TF1(
"func",
"gaus(0)", time_fit_min, time_fit_max);
699 gaus->SetParNames(
"numCrystalHitsNormalization",
"mean",
"sigma");
706 double hist_max = h_time->GetMaximum();
709 double stddev = h_time->GetStdDev();
714 gaus->SetParameter(0, hist_max / 2.);
715 gaus->SetParameter(1, mean);
716 gaus->SetParameter(2, sigma);
723 h_timeMasked->Fit(gaus,
"LQR");
725 double fit_mean = gaus->GetParameter(1);
726 double fit_mean_unc = gaus->GetParError(1);
727 double fit_sigma = gaus->GetParameter(2);
729 double meanDiff = fit_mean - default_mean;
730 double meanUncDiff = fit_mean_unc - default_mean_unc;
731 double sigmaDiff = fit_sigma - default_sigma;
733 bool good_fit =
false;
735 if ((fabs(meanDiff) > 10) ||
736 (fabs(meanUncDiff) > 10) ||
737 (fabs(sigmaDiff) > 10) ||
738 (fit_mean_unc > 3) ||
740 (fit_mean < time_fit_min) ||
741 (fit_mean > time_fit_max)) {
742 B2INFO(
"x_bin = " << x_bin);
743 B2INFO(
"fit mean, default mean = " << fit_mean <<
", " << default_mean);
744 B2INFO(
"fit mean unc, default mean unc = " << fit_mean_unc <<
", " << default_mean_unc);
745 B2INFO(
"fit sigma, default sigma = " << fit_sigma <<
", " << default_sigma);
747 B2INFO(
"crystal fit mean - hist mean = " << meanDiff);
748 B2INFO(
"fit mean unc. - hist mean unc. = " << meanUncDiff);
749 B2INFO(
"fit sigma - hist sigma = " << sigmaDiff);
751 B2INFO(
"fit_mean = " << fit_mean);
752 B2INFO(
"time_fit_min = " << time_fit_min);
753 B2INFO(
"time_fit_max = " << time_fit_max);
755 if (fabs(meanDiff) > 10) B2INFO(
"fit mean diff too large");
756 if (fabs(meanUncDiff) > 10) B2INFO(
"fit mean unc diff too large");
757 if (fabs(sigmaDiff) > 10) B2INFO(
"fit mean sigma diff too large");
758 if (fit_mean_unc > 3) B2INFO(
"fit mean unc too large");
759 if (fit_sigma < 0.1) B2INFO(
"fit sigma too small");
766 int numEntries = h_time->GetEntries();
770 if ((numEntries >= minNumEntries) && good_fit) {
771 clusterTime_mean = fit_mean;
772 clusterTime_mean_unc = fit_mean_unc;
773 clusterTime_sigma = fit_sigma;
775 clusterTime_mean = default_mean;
776 clusterTime_mean_unc = default_mean_unc;
777 clusterTime_sigma = default_sigma;
780 if (numEntries < minNumEntries) B2INFO(
"Number of entries less than minimum");
781 if (numEntries == 0) B2INFO(
"Number of entries == 0");
783 histfile->WriteTObject(h_time, (std::string(
"h_time_E_slice") + std::to_string(meanE)).c_str());
784 histfile->WriteTObject(h_timeMasked, (std::string(
"h_time_E_slice_masked") + std::to_string(meanE)).c_str());
787 clusterTimePeak_ClusterEnergy_varBin->SetBinContent(Ebin_counter, clusterTime_mean);
788 clusterTimePeak_ClusterEnergy_varBin->SetBinError(Ebin_counter, clusterTime_mean_unc);
790 clusterTimePeakWidth_ClusterEnergy_varBin->SetBinContent(Ebin_counter, clusterTime_sigma);
791 clusterTimePeakWidth_ClusterEnergy_varBin->SetBinError(Ebin_counter, 0);
807 vector< pair<double, int> > fitClusterTime__crystalIDBase0__pairs;
811 fitClusterTime__crystalIDBase0__pairs.push_back(make_pair(0.0, cid));
816 fitClusterTime__crystalIDBase0__pairs[crys_id - 1] = make_pair(fabs(t_offsets[crys_id - 1]), crys_id - 1) ;
820 sort(fitClusterTime__crystalIDBase0__pairs.begin(), fitClusterTime__crystalIDBase0__pairs.end());
824 B2INFO(
"-------- List of the (fitted) peak cluster times sorted by their absolute value ----------");
825 B2INFO(
"------------------------------------------------------------------------------------------");
826 B2INFO(
"------------------------------------------------------------------------------------------");
827 B2INFO(
"Quoted # of clusters is before the cutting off of the distribution tails, cellID=1..ECLElementNumbers::c_NCrystals, crysID=0..8735");
829 bool hasHitThresholdBadTimes = false ;
831 int cid = fitClusterTime__crystalIDBase0__pairs[iSortedTimes].second ;
832 if (!hasHitThresholdBadTimes && fitClusterTime__crystalIDBase0__pairs[iSortedTimes].first > 2) {
833 B2INFO(
"======== |t_fit| > Xns threshold ======");
834 hasHitThresholdBadTimes =
true;
837 B2INFO(
"cid = " << cid <<
", peak clust t = " << t_offsets[cid] <<
" +- " << t_offsets_unc[cid] <<
" ns, # clust = " <<
838 numClusterPerCrys[cid] <<
", good fit = " << crysHasGoodFit[cid] <<
", good fit & stats = " << crysHasGoodFitandStats[cid]);
846 B2INFO(
"######## List of poor (fitted) peak cluster times sorted by their absolute value #########");
847 B2INFO(
"##########################################################################################");
848 B2INFO(
"##########################################################################################");
851 int cid = fitClusterTime__crystalIDBase0__pairs[iSortedTimes].second ;
852 if (fitClusterTime__crystalIDBase0__pairs[iSortedTimes].first > 2 && crysHasGoodFitandStats[cid]) {
853 B2INFO(
"WARNING: cid = " << cid <<
", peak clust t = " << t_offsets[cid] <<
" +- " << t_offsets_unc[cid] <<
" ns, # clust = " <<
854 numClusterPerCrys[cid] <<
", good fit = " << crysHasGoodFit[cid] <<
", good fit & stats = " << crysHasGoodFitandStats[cid]);
860 B2INFO(
"Number of crystals with non-zero number of hits = " << numCrysWithNonZeroEntries);
861 B2INFO(
"Number of crystals with good quality fits = " << numCrysWithGoodFit);
864 clusterTimePeak_ClusterEnergy_varBin->ResetStats();
865 clusterTimePeakWidth_ClusterEnergy_varBin->ResetStats();
867 histfile->WriteTObject(clusterTimePeak_ClusterEnergy_varBin,
"clusterTimePeak_ClusterEnergy_varBin");
868 histfile->WriteTObject(clusterTimePeakWidth_ClusterEnergy_varBin,
"clusterTimePeakWidth_ClusterEnergy_varBin");
874 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.");
876 double tsDiffCustomOld_ns = -999;
878 tsDiffCustomOld_ns = (currentValuesCrys[crys_id - 1] - prevValuesCrys[crys_id - 1]) * TICKS_TO_NS;
879 B2INFO(
"Crystal " << crys_id <<
": ts new merged - 'before 1st iter' merged = (" <<
880 currentValuesCrys[crys_id - 1] <<
" - " << prevValuesCrys[crys_id - 1] <<
881 ") ticks * " << TICKS_TO_NS <<
" ns/tick = " << tsDiffCustomOld_ns <<
" ns");
884 tsNew_MINUS_tsCustomPrev__cid->SetBinContent(crys_id, tsDiffCustomOld_ns);
885 tsNew_MINUS_tsCustomPrev__cid->SetBinError(crys_id, 0);
886 tsNew_MINUS_tsCustomPrev__cid->ResetStats();
888 tsNew_MINUS_tsCustomPrev->Fill(tsDiffCustomOld_ns);
889 tsNew_MINUS_tsCustomPrev->ResetStats();
892 histfile->WriteTObject(tsNew_MINUS_tsCustomPrev__cid,
"tsNew_MINUS_tsCustomPrev__cid");
893 histfile->WriteTObject(tsNew_MINUS_tsCustomPrev,
"tsNew_MINUS_tsCustomPrev");
898 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.