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