Belle II Software development
eclTValidationAlgorithm.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9/* Own header. */
10#include <ecl/calibration/eclTValidationAlgorithm.h>
11
12/* ECL headers. */
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>
18
19/* Basf2 headers. */
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>
25
26/* ROOT headers. */
27#include <TF1.h>
28#include <TFile.h>
29#include <TGraphAsymmErrors.h>
30#include <TH2F.h>
31#include <TROOT.h>
32#include <TString.h>
33
34using namespace std;
35using namespace Belle2;
36using namespace ECL;
37
38
39/* By default assume the timing validation collector used a hadronic event selection
40 but a bhabha event selection also exists and uses the same algorithm to
41 analyse the results:
42 The collector name should be set in the steering script.*/
44 // Parameters
45 CalibrationAlgorithm("eclHadronTimeCalibrationValidationCollector"),
46 cellIDLo(1),
47 cellIDHi(ECLElementNumbers::c_NCrystals),
48 readPrevCrysPayload(false),
49 meanCleanRebinFactor(1),
50 meanCleanCutMinFactor(0),
51 clusterTimesFractionWindow_maxtime(8),
52 debugFilenameBase("eclTValidationAlgorithm")
53{
54 setDescription("Fit gaussian function to the cluster times to validate results.");
55}
56
57
58
59
60/* By default assume the timing validation collector used a hadronic event selection
61 but a bhabha event selection also exists and uses the same algorithm to
62 analyse the results:
63 The collector name should be set in the steering script.*/
64eclTValidationAlgorithm::eclTValidationAlgorithm(string physicsProcessCollectorName):
65 // Parameters
66 CalibrationAlgorithm(physicsProcessCollectorName.c_str()),
67 cellIDLo(1),
68 cellIDHi(ECLElementNumbers::c_NCrystals),
69 readPrevCrysPayload(false),
70 meanCleanRebinFactor(1),
71 meanCleanCutMinFactor(0),
72 clusterTimesFractionWindow_maxtime(8),
73 debugFilenameBase("eclTValidationAlgorithm")
74{
75 setDescription("Fit gaussian function to the cluster times to validate results.");
76}
77
78
79
80
82{
84 gROOT->SetBatch();
85
86
88 B2INFO("eclTValidationAlgorithm parameters:");
89 B2INFO("cellIDLo = " << cellIDLo);
90 B2INFO("cellIDHi = " << cellIDHi);
91 B2INFO("readPrevCrysPayload = " << readPrevCrysPayload);
92 B2INFO("meanCleanRebinFactor = " << meanCleanRebinFactor);
93 B2INFO("meanCleanCutMinFactor = " << meanCleanCutMinFactor);
94 B2INFO("clusterTimesFractionWindow_maxtime = " << clusterTimesFractionWindow_maxtime);
95
96
97 /* Histogram with the data collected by eclTimeCalibrationValidationCollector*/
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");
106
107 // Collect other plots just for reference - combines all the runs for these plots.
108 auto cutflow = getObjectPtr<TH1F>("cutflow");
109
110 vector <int> binProjectionLeft_Time_vs_E_runDep ;
111 vector <int> binProjectionRight_Time_vs_E_runDep ;
112
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);
116 }
117
118 if (!clusterTime_cid) return c_Failure;
119
122 TFile* histfile = 0;
123
124 // Vector of time offsets to track how far from nominal the cluster times are.
125 vector<float> t_offsets(ECLElementNumbers::c_NCrystals, 0.0);
126 vector<float> t_offsets_unc(ECLElementNumbers::c_NCrystals, 0.0);
127 vector<long> numClusterPerCrys(ECLElementNumbers::c_NCrystals, 0);
128 vector<bool> crysHasGoodFitandStats(ECLElementNumbers::c_NCrystals, false);
129 vector<bool> crysHasGoodFit(ECLElementNumbers::c_NCrystals, false);
130 int numCrysWithNonZeroEntries = 0 ;
131 int numCrysWithGoodFit = 0 ;
132
133 int minNumEntries = 40;
134
135 double mean;
136 double sigma;
137
138
139 bool minRunNumBool = false;
140 bool maxRunNumBool = false;
141 int minRunNum = -1;
142 int maxRunNum = -1;
143 int minExpNum = -1;
144 int maxExpNum = -1;
145 for (auto expRun : getRunList()) {
146 int expNumber = expRun.first;
147 int runNumber = expRun.second;
148 if (!minRunNumBool) {
149 minExpNum = expNumber;
150 minRunNum = runNumber;
151 minRunNumBool = true;
152 }
153 if (!maxRunNumBool) {
154 maxExpNum = expNumber;
155 maxRunNum = runNumber;
156 maxRunNumBool = true;
157 }
158 if (((minRunNum > runNumber) && (minExpNum >= expNumber)) ||
159 (minExpNum > expNumber)) {
160 minExpNum = expNumber;
161 minRunNum = runNumber;
162 }
163 if (((maxRunNum < runNumber) && (maxExpNum <= expNumber)) ||
164 (maxExpNum < expNumber))
165
166 {
167 maxExpNum = expNumber;
168 maxRunNum = runNumber;
169 }
170 }
171
172 B2INFO("debugFilenameBase = " << debugFilenameBase);
173 string runNumsString = string("_") + to_string(minExpNum) + "_" + to_string(minRunNum) + string("-") +
174 to_string(maxExpNum) + "_" + to_string(maxRunNum);
175 string debugFilename = debugFilenameBase + runNumsString + string(".root");
176
177
178 // Need to load information about the event/run/experiment to get the right database information
179 // Will be used for:
180 // * ECLChannelMapper (to map crystal to crates)
181 // * crystal payload updating for iterating crystal and crate fits
182 int eventNumberForCrates = 1;
183
184
185 //-------------------------------------------------------------------
186 /* Uploading older payloads for the current set of runs */
187
189 // simulate the initialize() phase where we can register objects in the DataStore
191 evtPtr.registerInDataStore();
193 // now construct the event metadata
194 evtPtr.construct(eventNumberForCrates, minRunNum, minExpNum);
195 // and update the database contents
196 DBStore& dbstore = DBStore::Instance();
197 dbstore.update();
198 // this is only needed it the payload might be intra-run dependent,
199 // that is if it might change during one run as well
200 dbstore.updateEvent();
201
202
203 B2INFO("Uploading payload for exp " << minExpNum << ", run " << minRunNum << ", event " << eventNumberForCrates);
204 updateDBObjPtrs(eventNumberForCrates, minRunNum, minExpNum);
205 unique_ptr<ECLChannelMapper> crystalMapper(new ECL::ECLChannelMapper());
206 crystalMapper->initFromDB();
207
208 /* 1/(4fRF) = 0.4913 ns/clock tick, where fRF is the accelerator RF frequency.
209 Same for all crystals. */
210 const double TICKS_TO_NS = 1.0 / (4.0 * EclConfiguration::getRF()) * 1e3;
211
212 //------------------------------------------------------------------------
213 //..Read payloads from database
214 DBObjPtr<Belle2::ECLCrystalCalib> crystalTimeObject("ECLCrystalTimeOffset");
215 B2INFO("Dumping payload");
216
217 //..Get vectors of values from the payloads
218 std::vector<float> currentValuesCrys = crystalTimeObject->getCalibVector();
219 std::vector<float> currentUncCrys = crystalTimeObject->getCalibUncVector();
220
221 //..Print out a few values for quality control
222 B2INFO("Values read from database. Write out for their values for comparison against those from tcol");
223 for (int ic = 0; ic < ECLElementNumbers::c_NCrystals; ic += 500) {
224 B2INFO("ts: cellID " << ic + 1 << " " << currentValuesCrys[ic] << " +/- " << currentUncCrys[ic]);
225 }
226
227
228 //..Read in the previous crystal payload values
229 DBObjPtr<Belle2::ECLCrystalCalib> customPrevCrystalTimeObject("ECLCrystalTimeOffsetPreviousValues");
230 vector<float> prevValuesCrys(ECLElementNumbers::c_NCrystals);
232 //..Get vectors of values from the payloads
233 prevValuesCrys = customPrevCrystalTimeObject->getCalibVector();
234
235 //..Print out a few values for quality control
236 B2INFO("Previous values read from database. Write out for their values for comparison against those from tcol");
237 for (int ic = 0; ic < ECLElementNumbers::c_NCrystals; ic += 500) {
238 B2INFO("ts custom previous payload: cellID " << ic + 1 << " " << prevValuesCrys[ic]);
239 }
240 }
241
242
243 //------------------------------------------------------------------------
244 //..Start looking at timing information
245
246 B2INFO("Debug output rootfile: " << debugFilename);
247 histfile = new TFile(debugFilename.c_str(), "recreate");
248
249
250 clusterTime ->Write();
251 clusterTime_cid ->Write();
252 clusterTime_run ->Write();
253 clusterTimeClusterE ->Write();
254 dt99_clusterE ->Write();
255 eventT0 ->Write();
256 eventT0Detector ->Write();
257 clusterTimeE0E1diff ->Write();
258
259 cutflow->Write();
260
261
262 double hist_tmin = clusterTime->GetXaxis()->GetXmin();
263 double hist_tmax = clusterTime->GetXaxis()->GetXmax();
264 int hist_nTbins = clusterTime->GetNbinsX();
265
266 B2INFO("hist_tmin = " << hist_tmin);
267 B2INFO("hist_tmax = " << hist_tmax);
268 B2INFO("hist_nTbins = " << hist_nTbins);
269
270 double time_fit_min = hist_tmax; // Set min value to largest possible value so that it gets reset
271 double time_fit_max = hist_tmin; // Set max value to smallest possible value so that it gets reset
272
273
274 // define histogram for keeping track of the peak of the cluster times per crystal
275 auto peakClusterTime_cid = new TH1F("peakClusterTime_cid", ";cell id;Peak cluster time [ns]", ECLElementNumbers::c_NCrystals, 1,
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);
283
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]",
287
288
289 // define histograms to keep track of the difference in the new crystal times vs the old ones
290 auto tsNew_MINUS_tsCustomPrev__cid = new TH1F("TsNew_MINUS_TsCustomPrev__cid",
291 ";cell id; ts(new|merged) - ts(old = 'pre-calib'|merged) [ns]",
293
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);
297
298
299
300 // Histogram to keep track of the fraction of cluster times within a window.
301 double timeWindow_maxTime = clusterTimesFractionWindow_maxtime;
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.");
316
317 TString fracWindowHistTitle = "Fraction of cluster times in window [" + s_lowTime + ", " + s_highTime +
318 "] ns;Fraction of cluster times in window;Number of crystals";
319
320 auto clusterTimeNumberInWindow__cid = new TH1F("clusterTimeNumberInWindow__cid", fracWindowTitle, ECLElementNumbers::c_NCrystals, 1,
322 auto clusterTimeNumberInWindowInGoodECLRings__cid = new TH1F("clusterTimeNumberInWindowInGoodECLRings__cid", fracWindowTitle,
325 auto clusterTimeNumber__cid = new TH1F("clusterTimeNumber_cid", fracWindowTitle, ECLElementNumbers::c_NCrystals, 1,
327 auto clusterTimeFractionInWindow = new TH1F("clusterTimeFractionInWindow", fracWindowHistTitle, 110, 0.0, 1.1);
328
329 clusterTimeNumberInWindow__cid->Sumw2();
330 clusterTimeNumberInWindowInGoodECLRings__cid->Sumw2();
331 clusterTimeNumber__cid->Sumw2();
332
333
334
335 /* CRYSTAL BY CRYSTAL VALIDATION */
336
338
339 // Loop over all the crystals for doing the crystal calibation
340 for (int crys_id = cellIDLo; crys_id <= cellIDHi; crys_id++) {
341 double clusterTime_mean = 0;
342 double clusterTime_mean_unc = 0;
343
344 B2INFO("Crystal cell id = " << crys_id);
345
346 eclgeo->Mapping(crys_id - 1);
347 int thetaID = eclgeo->GetThetaID();
348
349
350 /* Determining which bins to mask out for mean calculation
351 */
352
353 TH1D* h_time = clusterTime_cid->ProjectionY((std::string("h_time_psi__") + std::to_string(crys_id)).c_str(),
354 crys_id, crys_id);
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();
358
359 // Do rebinning and cleaning of some bins but only if the user selection values call for it since it slows the code down
361
362 h_timeRebin->Rebin(meanCleanRebinFactor);
363
364 h_timeMask->Scale(0.0); // set all bins to being masked off
365
366 time_fit_min = hist_tmax; // Set min value to largest possible value so that it gets reset
367 time_fit_max = hist_tmin; // Set max value to smallest possible value so that it gets reset
368
369 // Find value of bin with max value
370 double histRebin_max = h_timeRebin->GetMaximum();
371
372 bool maskedOutNonZeroBin = false;
373 // Loop over all bins to find those with content less than a certain threshold. Mask the non-rebinned histogram for the corresponding bins
374 for (int bin = 1; bin <= h_timeRebin->GetNbinsX(); bin++) {
375 for (int rebinCounter = 1; rebinCounter <= meanCleanRebinFactor; rebinCounter++) {
376 int nonRebinnedBinNumber = (bin - 1) * meanCleanRebinFactor + rebinCounter;
377 if (nonRebinnedBinNumber < h_time->GetNbinsX()) {
378 if (h_timeRebin->GetBinContent(bin) >= histRebin_max * meanCleanCutMinFactor) {
379 h_timeMask->SetBinContent(nonRebinnedBinNumber, 1);
380
381 // Save the lower and upper edges of the rebin histogram time range for fitting purposes
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;
386 }
387 if (x_upper > time_fit_max) {
388 time_fit_max = x_upper;
389 }
390
391 } else {
392 if (h_time->GetBinContent(nonRebinnedBinNumber) > 0) {
393 B2DEBUG(22, "Setting bin " << nonRebinnedBinNumber << " from " << h_timeMasked->GetBinContent(nonRebinnedBinNumber) << " to 0");
394 maskedOutNonZeroBin = true;
395 }
396 h_timeMasked->SetBinContent(nonRebinnedBinNumber, 0);
397 }
398 }
399 }
400 }
401 B2INFO("Bins with non-zero values have been masked out: " << maskedOutNonZeroBin);
402 h_timeMasked->ResetStats();
403 h_timeMask->ResetStats();
404
405 }
406
407 // Calculate mean from masked histogram
408 double default_meanMasked = h_timeMasked->GetMean();
409 //double default_meanMasked_unc = h_timeMasked->GetMeanError();
410 B2INFO("default_meanMasked = " << default_meanMasked);
411
412
413 // Get the overall mean and standard deviation of the distribution within the plot. This doesn't require a fit.
414 double default_mean = h_time->GetMean();
415 double default_mean_unc = h_time->GetMeanError();
416 double default_sigma = h_time->GetStdDev();
417
418 B2INFO("Fitting crystal between " << time_fit_min << " and " << time_fit_max);
419
420 // gaus(0) is a substitute for [0]*exp(-0.5*((x-[1])/[2])**2)
421 TF1* gaus = new TF1("func", "gaus(0)", time_fit_min, time_fit_max);
422 gaus->SetParNames("numCrystalHitsNormalization", "mean", "sigma");
423 /*
424 gaus->ReleaseParameter(0); // number of crystals
425 gaus->ReleaseParameter(1); // mean
426 gaus->ReleaseParameter(2); // standard deviation
427 */
428
429 double hist_max = h_time->GetMaximum();
430
431 //=== Estimate initial value of sigma as std dev.
432 double stddev = h_time->GetStdDev();
433 sigma = stddev;
434 mean = default_mean;
435
436 //=== Setting parameters for initial iteration
437 gaus->SetParameter(0, hist_max / 2.);
438 gaus->SetParameter(1, mean);
439 gaus->SetParameter(2, sigma);
440 // L -- Use log likelihood method
441 // I -- Use integral of function in bin instead of value at bin center // not using
442 // R -- Use the range specified in the function range
443 // B -- Fix one or more parameters with predefined function // not using
444 // Q -- Quiet mode
445
446 h_timeMasked->Fit(gaus, "LQR"); // L for likelihood, R for x-range, Q for fit quiet mode
447
448 double fit_mean = gaus->GetParameter(1);
449 double fit_mean_unc = gaus->GetParError(1);
450 double fit_sigma = gaus->GetParameter(2);
451
452 double meanDiff = fit_mean - default_mean;
453 double meanUncDiff = fit_mean_unc - default_mean_unc;
454 double sigmaDiff = fit_sigma - default_sigma;
455
456 bool good_fit = false;
457
458 if ((fabs(meanDiff) > 10) ||
459 (fabs(meanUncDiff) > 10) ||
460 (fabs(sigmaDiff) > 10) ||
461 (fit_mean_unc > 3) ||
462 (fit_sigma < 0.1) ||
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);
469
470 B2INFO("crystal fit mean - hist mean = " << meanDiff);
471 B2INFO("fit mean unc. - hist mean unc. = " << meanUncDiff);
472 B2INFO("fit sigma - hist sigma = " << sigmaDiff);
473
474 B2INFO("fit_mean = " << fit_mean);
475 B2INFO("time_fit_min = " << time_fit_min);
476 B2INFO("time_fit_max = " << time_fit_max);
477
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");
483
484 } else {
485 good_fit = true;
486 numCrysWithGoodFit++;
487 crysHasGoodFit[crys_id - 1] = true ;
488 }
489
490
491 int numEntries = h_time->GetEntries();
492 // If number of entries in histogram is greater than X then use the statistical information from the data otherwise leave crystal uncalibrated. Histograms are still shown though.
493 // ALSO require the that fits are good.
494 if ((numEntries >= minNumEntries) && good_fit) {
495 clusterTime_mean = fit_mean;
496 clusterTime_mean_unc = fit_mean_unc;
497 crysHasGoodFitandStats[crys_id - 1] = true ;
498 } else {
499 clusterTime_mean = default_mean;
500 clusterTime_mean_unc = default_mean_unc;
501 }
502
503 if (numEntries < minNumEntries) B2INFO("Number of entries less than minimum");
504 if (numEntries == 0) B2INFO("Number of entries == 0");
505
506
507 t_offsets[crys_id - 1] = clusterTime_mean ;
508 t_offsets_unc[crys_id - 1] = clusterTime_mean_unc ;
509 numClusterPerCrys[crys_id - 1] = numEntries ;
510
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());
513
514 // Set this for each crystal even if there are zero entries
515 peakClusterTime_cid->SetBinContent(crys_id, t_offsets[crys_id - 1]);
516 peakClusterTime_cid->SetBinError(crys_id, t_offsets_unc[crys_id - 1]);
517
518 /* Store mean cluster time info in a separate histogram but only if there is at
519 least one entry for that crystal. */
520 if (numEntries > 0) {
521 peakClusterTimes->Fill(t_offsets[crys_id - 1]);
522 numCrysWithNonZeroEntries++ ;
523 }
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]);
528 }
529
530
531 // Find the fraction of cluster times within +-X ns and fill histograms
532 double numClusterTimesWithinWindowFraction = h_time->Integral(binyLeft, binyRight) ;
533 double clusterTimesWithinWindowFraction = numClusterTimesWithinWindowFraction;
534 if (numEntries > 0) {
535 clusterTimesWithinWindowFraction /= numEntries;
536 } else {
537 clusterTimesWithinWindowFraction = -0.1;
538 }
539
540 B2INFO("Crystal cell id = " << crys_id << ", theta id = " <<
541 thetaID << ", clusterTimesWithinWindowFraction = " <<
542 numClusterTimesWithinWindowFraction << " / " << numEntries << " = " <<
543 clusterTimesWithinWindowFraction);
544
545 clusterTimeFractionInWindow->Fill(clusterTimesWithinWindowFraction);
546 clusterTimeNumberInWindow__cid->SetBinContent(crys_id, numClusterTimesWithinWindowFraction);
547 clusterTimeNumber__cid->SetBinContent(crys_id, numEntries);
548
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);
554 }
555
556
557 delete gaus;
558 }
559
560 // Find the fraction of cluster times within +-X ns and fill histogram
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);
566
567
568 peakClusterTime_cid->ResetStats();
569 peakClusterTimesGoodFit__cid->ResetStats();
570
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");
578
579
580
581 /* -----------------------------------------------------------
582 Fit the time histograms for different energy slices */
583
584 vector <int> binProjectionLeft = binProjectionLeft_Time_vs_E_runDep;
585 vector <int> binProjectionRight = binProjectionRight_Time_vs_E_runDep;
586
587 auto h2 = clusterTimeClusterE;
588
589
590 double max_E = h2->GetXaxis()->GetXmax();
591
592 // Determine the energy bins. Save the left edge for histogram purposes
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]);
601 }
602 }
603
604
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]));
609
610 int Ebin_counter = 1 ;
611
612 // Loop over all the energy bins
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;
617
618 B2INFO("x_bin = " << x_bin);
619
620 /* Determining which bins to mask out for mean calculation
621 */
622 TH1D* h_time = h2->ProjectionY(("h_time_E_slice_" + std::to_string(x_bin)).c_str(), binProjectionLeft[x_bin],
623 binProjectionRight[x_bin]) ;
624
625
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 ;
630
631 B2INFO("bin " << Ebin_counter << ": low E = " << lowE << ", high E = " << highE << " GeV");
632
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();
636
637
639
640 h_timeRebin->Rebin(meanCleanRebinFactor);
641
642 h_timeMask->Scale(0.0); // set all bins to being masked off
643
644 time_fit_min = hist_tmax; // Set min value to largest possible value so that it gets reset
645 time_fit_max = hist_tmin; // Set max value to smallest possible value so that it gets reset
646
647 // Find value of bin with max value
648 double histRebin_max = h_timeRebin->GetMaximum();
649
650 bool maskedOutNonZeroBin = false;
651 // Loop over all bins to find those with content less than a certain threshold. Mask the non-rebinned histogram for the corresponding bins
652 for (int bin = 1; bin <= h_timeRebin->GetNbinsX(); bin++) {
653 for (int rebinCounter = 1; rebinCounter <= meanCleanRebinFactor; rebinCounter++) {
654 int nonRebinnedBinNumber = (bin - 1) * meanCleanRebinFactor + rebinCounter;
655 if (nonRebinnedBinNumber < h_time->GetNbinsX()) {
656 if (h_timeRebin->GetBinContent(bin) >= histRebin_max * meanCleanCutMinFactor) {
657 h_timeMask->SetBinContent(nonRebinnedBinNumber, 1);
658
659 // Save the lower and upper edges of the rebin histogram time range for fitting purposes
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;
664 }
665 if (x_upper > time_fit_max) {
666 time_fit_max = x_upper;
667 }
668
669 } else {
670 if (h_time->GetBinContent(nonRebinnedBinNumber) > 0) {
671 B2DEBUG(22, "Setting bin " << nonRebinnedBinNumber << " from " << h_timeMasked->GetBinContent(nonRebinnedBinNumber) << " to 0");
672 maskedOutNonZeroBin = true;
673 }
674 h_timeMasked->SetBinContent(nonRebinnedBinNumber, 0);
675 }
676 }
677 }
678 }
679 B2INFO("Bins with non-zero values have been masked out: " << maskedOutNonZeroBin);
680 h_timeMasked->ResetStats();
681 h_timeMask->ResetStats();
682
683 }
684
685
686 // Calculate mean from masked histogram
687 double default_meanMasked = h_timeMasked->GetMean();
688 //double default_meanMasked_unc = h_timeMasked->GetMeanError();
689 B2INFO("default_meanMasked = " << default_meanMasked);
690
691
692 // Get the overall mean and standard deviation of the distribution within the plot. This doesn't require a fit.
693 double default_mean = h_time->GetMean();
694 double default_mean_unc = h_time->GetMeanError();
695 double default_sigma = h_time->GetStdDev();
696
697 B2INFO("Fitting crystal between " << time_fit_min << " and " << time_fit_max);
698
699 // gaus(0) is a substitute for [0]*exp(-0.5*((x-[1])/[2])**2)
700 TF1* gaus = new TF1("func", "gaus(0)", time_fit_min, time_fit_max);
701 gaus->SetParNames("numCrystalHitsNormalization", "mean", "sigma");
702 /*
703 gaus->ReleaseParameter(0); // number of crystals
704 gaus->ReleaseParameter(1); // mean
705 gaus->ReleaseParameter(2); // standard deviation
706 */
707
708 double hist_max = h_time->GetMaximum();
709
710 //=== Estimate initial value of sigma as std dev.
711 double stddev = h_time->GetStdDev();
712 sigma = stddev;
713 mean = default_mean;
714
715 //=== Setting parameters for initial iteration
716 gaus->SetParameter(0, hist_max / 2.);
717 gaus->SetParameter(1, mean);
718 gaus->SetParameter(2, sigma);
719 // L -- Use log likelihood method
720 // I -- Use integral of function in bin instead of value at bin center // not using
721 // R -- Use the range specified in the function range
722 // B -- Fix one or more parameters with predefined function // not using
723 // Q -- Quiet mode
724
725 h_timeMasked->Fit(gaus, "LQR"); // L for likelihood, R for x-range, Q for fit quiet mode
726
727 double fit_mean = gaus->GetParameter(1);
728 double fit_mean_unc = gaus->GetParError(1);
729 double fit_sigma = gaus->GetParameter(2);
730
731 double meanDiff = fit_mean - default_mean;
732 double meanUncDiff = fit_mean_unc - default_mean_unc;
733 double sigmaDiff = fit_sigma - default_sigma;
734
735 bool good_fit = false;
736
737 if ((fabs(meanDiff) > 10) ||
738 (fabs(meanUncDiff) > 10) ||
739 (fabs(sigmaDiff) > 10) ||
740 (fit_mean_unc > 3) ||
741 (fit_sigma < 0.1) ||
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);
748
749 B2INFO("crystal fit mean - hist mean = " << meanDiff);
750 B2INFO("fit mean unc. - hist mean unc. = " << meanUncDiff);
751 B2INFO("fit sigma - hist sigma = " << sigmaDiff);
752
753 B2INFO("fit_mean = " << fit_mean);
754 B2INFO("time_fit_min = " << time_fit_min);
755 B2INFO("time_fit_max = " << time_fit_max);
756
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");
762
763 } else {
764 good_fit = true;
765 }
766
767
768 int numEntries = h_time->GetEntries();
769 /* If number of entries in histogram is greater than X then use the statistical information
770 from the data otherwise leave crystal uncalibrated. Histograms are still shown though.
771 ALSO require the that fits are good. */
772 if ((numEntries >= minNumEntries) && good_fit) {
773 clusterTime_mean = fit_mean;
774 clusterTime_mean_unc = fit_mean_unc;
775 clusterTime_sigma = fit_sigma;
776 } else {
777 clusterTime_mean = default_mean;
778 clusterTime_mean_unc = default_mean_unc;
779 clusterTime_sigma = default_sigma;
780 }
781
782 if (numEntries < minNumEntries) B2INFO("Number of entries less than minimum");
783 if (numEntries == 0) B2INFO("Number of entries == 0");
784
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());
787
788 // store mean cluster time info in a separate histogram
789 clusterTimePeak_ClusterEnergy_varBin->SetBinContent(Ebin_counter, clusterTime_mean);
790 clusterTimePeak_ClusterEnergy_varBin->SetBinError(Ebin_counter, clusterTime_mean_unc);
791
792 clusterTimePeakWidth_ClusterEnergy_varBin->SetBinContent(Ebin_counter, clusterTime_sigma);
793 clusterTimePeakWidth_ClusterEnergy_varBin->SetBinError(Ebin_counter, 0);
794
795 Ebin_counter++;
796
797 delete gaus;
798 }
799
800
801
802 /***************************************************************************
803 For the user, print out some information about the peak cluster times.
804 It is sorted by the absolute value of the peak cluster time so that the
805 worst times are at the end.
806 ***************************************************************************/
807
808 // Vector to store element with respective present index
809 vector< pair<double, int> > fitClusterTime__crystalIDBase0__pairs;
810
811 // Prepare a vector of pairs containing the fitted cluster time and cell ID (base 0)
812 for (int cid = 0; cid < ECLElementNumbers::c_NCrystals; cid++) {
813 fitClusterTime__crystalIDBase0__pairs.push_back(make_pair(0.0, cid));
814 }
815
816 // Inserting element in pair vector to keep track of crystal id.
817 for (int crys_id = cellIDLo; crys_id <= cellIDHi; crys_id++) {
818 fitClusterTime__crystalIDBase0__pairs[crys_id - 1] = make_pair(fabs(t_offsets[crys_id - 1]), crys_id - 1) ;
819 }
820
821 // Sorting by the absolute value of the fitted cluster time for the crystal
822 sort(fitClusterTime__crystalIDBase0__pairs.begin(), fitClusterTime__crystalIDBase0__pairs.end());
823
824
825 // Print out the fitted peak cluster time values sorted by their absolute value
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");
830
831 bool hasHitThresholdBadTimes = false ;
832 for (int iSortedTimes = 0; iSortedTimes < ECLElementNumbers::c_NCrystals; iSortedTimes++) {
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;
837 }
838 //B2INFO("crystal ID = " << cid << ", peak clust t = " << t_offsets[cid] << " +- " << t_offsets_unc[cid] << ", # clusters = " << numClusterPerCrys[cid] << ", fabs(t) = " << fitClusterTime__crystalIDBase0__pairs[iSortedTimes].first );
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]);
841 }
842
843
844
845
846
847 // Print out just a subset that definitely don't look good even though they have good stats.
848 B2INFO("######## List of poor (fitted) peak cluster times sorted by their absolute value #########");
849 B2INFO("##########################################################################################");
850 B2INFO("##########################################################################################");
851
852 for (int iSortedTimes = 0; iSortedTimes < ECLElementNumbers::c_NCrystals; iSortedTimes++) {
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]);
857 }
858 }
859
860
861 B2INFO("~~~~~~~~");
862 B2INFO("Number of crystals with non-zero number of hits = " << numCrysWithNonZeroEntries);
863 B2INFO("Number of crystals with good quality fits = " << numCrysWithGoodFit);
864
865
866 clusterTimePeak_ClusterEnergy_varBin->ResetStats();
867 clusterTimePeakWidth_ClusterEnergy_varBin->ResetStats();
868
869 histfile->WriteTObject(clusterTimePeak_ClusterEnergy_varBin, "clusterTimePeak_ClusterEnergy_varBin");
870 histfile->WriteTObject(clusterTimePeakWidth_ClusterEnergy_varBin, "clusterTimePeakWidth_ClusterEnergy_varBin");
871
872
873
874 /* Fill histograms with the difference in the ts values from this iteration
875 and the previous values read in from the payload. */
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.");
877 for (int crys_id = 1; crys_id <= ECLElementNumbers::c_NCrystals; crys_id++) {
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");
884
885 }
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();
889
890 tsNew_MINUS_tsCustomPrev->Fill(tsDiffCustomOld_ns);
891 tsNew_MINUS_tsCustomPrev->ResetStats();
892 }
893
894 histfile->WriteTObject(tsNew_MINUS_tsCustomPrev__cid, "tsNew_MINUS_tsCustomPrev__cid");
895 histfile->WriteTObject(tsNew_MINUS_tsCustomPrev, "tsNew_MINUS_tsCustomPrev");
896
897
898 histfile->Close();
899
900 B2INFO("Finished validations algorithm");
901 return c_OK;
902}
903
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.
Definition: DBObjPtr.h:21
Singleton class to cache database objects.
Definition: DBStore.h:31
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:53
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:93
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.
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.
Definition: StoreObjPtr.h:95
bool construct(Args &&... params)
Construct an object of type T in this StoreObjPtr, using the provided constructor arguments.
Definition: StoreObjPtr.h:118
static DBStore & Instance()
Instance of a singleton DBStore.
Definition: DBStore.cc:26
void updateEvent()
Updates all intra-run dependent objects.
Definition: DBStore.cc:140
void update()
Updates all objects that are outside their interval of validity.
Definition: DBStore.cc:77
const int c_NCrystals
Number of crystals.
Abstract base class for different kinds of events.
STL namespace.