Belle II Software  release-05-02-19
eclTimeShiftsAlgorithm.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 
10 #include <ecl/calibration/eclTimeShiftsAlgorithm.h>
11 #include <ecl/dbobjects/ECLCrystalCalib.h>
12 #include <ecl/dbobjects/ECLReferenceCrystalPerCrateCalib.h>
13 #include <ecl/digitization/EclConfiguration.h>
14 #include <ecl/utility/ECLChannelMapper.h>
15 #include "TH1F.h"
16 #include "TString.h"
17 #include "TFile.h"
18 #include "TDirectory.h"
19 #include <TCanvas.h>
20 #include <TGraphErrors.h>
21 #include <TLatex.h>
22 #include <sstream>
23 #include <iomanip>
24 #include <TLatex.h>
25 
26 using namespace std;
27 using namespace Belle2;
28 using namespace ECL;
29 using namespace Calibration;
30 
32 //eclTimeShiftsAlgorithm::eclTimeShiftsAlgorithm(): CalibrationAlgorithm("DummyCollector"),
33 eclTimeShiftsAlgorithm::eclTimeShiftsAlgorithm():
34  CalibrationAlgorithm("eclTimeShiftsPlottingCollector"),
35  debugFilenameBase("ECL_time_offsets"),
36  timeShiftForPlotStyle{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
37  crysCrateShift_min(-20),
38  crysCrateShift_max(20),
39  algorithmReadPayloads(false),
40  m_ECLCrystalTimeOffset("ECLCrystalTimeOffset"),
41  m_ECLCrateTimeOffset("ECLCrateTimeOffset"),
42  m_refCrysIDzeroingCrate("ECLReferenceCrystalPerCrateCalib")//,
43 {
44  setDescription(
45  "Plots the ecl crystal and crate time calibations."
46  );
47 }
48 
49 CalibrationAlgorithm::EResult eclTimeShiftsAlgorithm::calibrate()
50 {
52  gROOT->SetBatch();
53 
54  B2INFO("eclTimeShiftsAlgorithm parameters:");
55  B2INFO("debugFilenameBase = " << debugFilenameBase);
56  B2INFO("algorithmReadPayloads = " << algorithmReadPayloads);
57  B2INFO("timeShiftForPlotStyle = {");
58  for (int crateTest = 0; crateTest < 51; crateTest++) {
59  B2INFO(timeShiftForPlotStyle[crateTest] << ",");
60  }
61  B2INFO(timeShiftForPlotStyle[51] << "}");
62 
63 
64  //------------------------------------------------------------------------
65  /* Conversion coefficient from ADC ticks to nanoseconds
66  1/(4fRF) = 0.4913 ns/clock tick, where fRF is the accelerator RF frequency,
67  fRF=508.889 MHz. Same for all crystals. Proper accurate value */
68  const double TICKS_TO_NS = 1.0 / (4.0 * EclConfiguration::m_rf) * 1e3;
69 
70 
71  //------------------------------------------------------------------------
72  /* Set up variables for storing timing information and cutting on
73  timing quality */
74 
75  vector< vector<double> > allCrates_crate_times ;
76  vector< vector<double> > allCrates_run_nums ; // not an integer for plotting purposes
77  vector< vector<double> > allCrates_time_unc ;
78  vector< vector<double> > allCrates_crystalCrate_times ;
79  vector< vector<double> > allCrates_crystalCrate_times_unc ;
80 
81  vector<int> allRunNums;
82 
83  vector<double> mean_crystalCrate_time_ns(m_numCrates, 0);
84 
85  vector< double > blank_vector = {} ;
86  vector< int > blank_vector_int = {} ;
87  for (int temp_crate_id = 0; temp_crate_id < m_numCrates; temp_crate_id++) {
88  allCrates_crate_times.push_back(blank_vector) ;
89  allCrates_run_nums.push_back(blank_vector) ;
90  allCrates_time_unc.push_back(blank_vector) ;
91  allCrates_crystalCrate_times.push_back(blank_vector) ;
92  allCrates_crystalCrate_times_unc.push_back(blank_vector) ;
93  }
94  // This results in : allCrates_crate_time[index for crate number][index for run number]
95 
96 
97 
98  //------------------------------------------------------------------------
99  /* Extract the crystal and crate calibration constant information from the
100  tree as extracted by the collector. */
101 
102  // Pulling in data from collector output. It now returns shared_ptr<T> so the underlying pointer
103  // will delete itself automatically at the end of this scope unless you do something
104  auto tree_perCrys = getObjectPtr<TTree>("tree_perCrystal");
105  if (!tree_perCrys) {
106  B2ERROR("Tree of calibration constants does not exist.");
107  return c_Failure;
108  }
109  B2INFO("Number of Entries in tree_perCrystal was " << tree_perCrys->GetEntries());
110  B2INFO("Number of Entries in tree_perCrystal / 8736 = " << float(tree_perCrys->GetEntries()) / 8736.0);
111 
112 
113  // Define the variables to be read in from the tree
114  tree_perCrys->SetBranchAddress("run", &m_run_perCrystal);
115  tree_perCrys->SetBranchAddress("exp", &m_exp_perCrystal);
116  tree_perCrys->SetBranchAddress("crystalID", &m_crystalID);
117  tree_perCrys->SetBranchAddress("crateID", &m_crateID);
118  tree_perCrys->SetBranchAddress("crateTimeConst", &m_crateTimeConst);
119  tree_perCrys->SetBranchAddress("crateTimeUnc", &m_crateTimeUnc);
120  tree_perCrys->SetBranchAddress("crystalTimeConst", &m_crystalTimeConst);
121  tree_perCrys->SetBranchAddress("crystalTimeUnc", &m_crystalTimeUnc);
122  tree_perCrys->SetBranchAddress("refCrystalID", &m_refCrystalID);
123 
124 
125  int referenceRunNum = -1;
126  int referenceExpNum = -1;
127  //int numAnalysedRuns = 0 ;
128  int previousRunNumTree = -1 ;
129  vector<double> Crate_time_ns_tree(m_numCrates) ;
130  vector<double> Crate_time_tick_tree(m_numCrates) ;
131  vector<double> Crate_time_unc_ns_tree(m_numCrates) ;
132  vector<double> crystalCrate_time_ns_tree(m_numCrates);
133  vector<double> crystalCrate_time_unc_ns_tree(m_numCrates);
134 
135 
136  Int_t numEntriesCrysTree = (Int_t)tree_perCrys->GetEntries();
137 
138  // Loop through the entire tree
139  for (Int_t tree_crys_i = 0; tree_crys_i < numEntriesCrysTree; tree_crys_i++) {
140  for (Int_t tree_crys_j = 0; tree_crys_j < m_numCrystals; tree_crys_j++) {
141  tree_perCrys->GetEntry(tree_crys_i);
142  //B2INFO("tree_crys_i, tree_crys_j = " << tree_crys_i << ", " << tree_crys_j);
143  if (tree_crys_j != m_numCrystals - 1) {
144  tree_crys_i++;
145  }
146 
147  // Make sure that all the information read in for 8736 crystals are all from one (exp,run).
148  if (tree_crys_j == 0) {
149  referenceExpNum = m_exp_perCrystal;
150  referenceRunNum = m_run_perCrystal;
151  B2INFO("Looking at exp,run " << m_exp_perCrystal << ", " << m_run_perCrystal);
152  }
153  if ((m_exp_perCrystal != referenceExpNum) or
154  (m_run_perCrystal != referenceRunNum) or
155  (m_run_perCrystal == previousRunNumTree)) {
156 
157  B2ERROR("m_exp_perCrystal, referenceExpNum" << m_exp_perCrystal << ", " << referenceExpNum);
158  B2ERROR("m_run_perCrystal, referenceRunNum" << m_run_perCrystal << ", " << referenceRunNum);
159  B2ERROR("m_run_perCrystal, previousRunNumTree" << m_run_perCrystal << ", " << previousRunNumTree);
160  B2ERROR("Exp/run number problem");
161  return c_Failure;
162  }
163 
164 
165  int crateID_temp = m_crateID;
166  Crate_time_ns_tree[crateID_temp - 1] = m_crateTimeConst * TICKS_TO_NS ;
167  Crate_time_tick_tree[crateID_temp - 1] = m_crateTimeConst ;
168  Crate_time_unc_ns_tree[crateID_temp - 1] = m_crateTimeUnc * TICKS_TO_NS ;
169 
170  if (m_crystalID == m_refCrystalID) {
171  B2INFO("m_exp_perCrystal, m_run_perCrystal, cell ID (0..8735), m_crateID, m_crateTimeConst = " << m_exp_perCrystal << ", " <<
172  m_run_perCrystal << ", " << tree_crys_j << ", " << m_crateID << ", " << m_crateTimeConst << " ticks") ;
173  crystalCrate_time_ns_tree[crateID_temp - 1] = (m_crystalTimeConst + m_crateTimeConst) * TICKS_TO_NS;
174 
175  crystalCrate_time_unc_ns_tree[crateID_temp - 1] = TICKS_TO_NS * sqrt(
176  (m_crateTimeUnc * m_crateTimeUnc) +
177  (m_crystalTimeUnc * m_crystalTimeUnc)) ;
178  } else if (tree_crys_j == 0 || tree_crys_j == 8735) {
179  B2INFO("m_exp_perCrystal, m_run_perCrystal, cell ID (0..8735), m_crateID, m_crateTimeConst = " << m_exp_perCrystal << ", " <<
180  m_run_perCrystal << ", " << tree_crys_j << ", " << m_crateID << ", " << m_crateTimeConst << " ns") ;
181  } else {
182  B2DEBUG(22, "m_exp_perCrystal, m_run_perCrystal, cell ID (0..8735), m_crateID, m_crateTimeConst = " << m_exp_perCrystal << ", " <<
183  m_run_perCrystal << ", " << tree_crys_j << ", " << m_crateID << ", " << m_crateTimeConst << " ns") ;
184  }
185 
186  }
187 
188  //------------------------------------------------------------------------
191  bool savedThisRunNum = false;
192  for (int iCrate = 0; iCrate < m_numCrates; iCrate++) {
193  double tcrate = Crate_time_ns_tree[iCrate] ;
194  double tcrate_unc = Crate_time_unc_ns_tree[iCrate];
195  double tcrystalCrate = crystalCrate_time_ns_tree[iCrate];
196  double tcrystalCrate_unc = crystalCrate_time_unc_ns_tree[iCrate];
197 
198  if ((tcrate < m_tcrate_max_cut) &&
199  (tcrate > m_tcrate_min_cut) &&
200  (fabs(tcrate_unc) > m_tcrate_unc_min_cut) &&
201  (fabs(tcrate_unc) < m_tcrate_unc_max_cut)) {
202  allCrates_crate_times[iCrate].push_back(tcrate) ;
203  allCrates_run_nums[iCrate].push_back(m_run_perCrystal) ;
204  allCrates_time_unc[iCrate].push_back(tcrate_unc) ;
205  allCrates_crystalCrate_times[iCrate].push_back(tcrystalCrate) ;
206  allCrates_crystalCrate_times_unc[iCrate].push_back(tcrystalCrate_unc) ;
207 
208  mean_crystalCrate_time_ns[iCrate] += tcrystalCrate ;
209 
210  if (!savedThisRunNum) {
211  allRunNums.push_back(m_run_perCrystal);
212  savedThisRunNum = true;
213  }
214  }
215  }
216 
217  //------------------------------------------------------------------------
219  for (int ic = 0; ic < m_numCrates; ic++) {
220  B2INFO("Crate " << ic + 1 << ", t_crate = " << Crate_time_tick_tree[ic] << " ticks = "
221  << Crate_time_ns_tree[ic] << " +- " << Crate_time_unc_ns_tree[ic]
222  << " ns; t crys+crate (no shifts) = " << crystalCrate_time_ns_tree[ic] << " +- "
223  << crystalCrate_time_unc_ns_tree[ic] << " ns") ;
224  }
225 
226  previousRunNumTree = m_run_perCrystal;
227 
228  }
229 
230 
231  B2INFO("Finished reading tree calibration constants. Now extracting here by stepping through runs.");
232 
233 
234 
235 
236 
237  //------------------------------------------------------------------------
241  bool minRunNumBool = false;
242  bool maxRunNumBool = false;
243  int minRunNum = -1;
244  int maxRunNum = -1;
245  int minExpNum = -1;
246  int maxExpNum = -1;
247  for (auto expRun : getRunList()) {
248  int expNumber = expRun.first;
249  int runNumber = expRun.second;
250  if (!minRunNumBool) {
251  minExpNum = expNumber;
252  minRunNum = runNumber;
253  minRunNumBool = true;
254  }
255  if (!maxRunNumBool) {
256  maxExpNum = expNumber;
257  maxRunNum = runNumber;
258  maxRunNumBool = true;
259  }
260  if (((minRunNum > runNumber) && (minExpNum >= expNumber)) ||
261  (minExpNum > expNumber)) {
262  minExpNum = expNumber;
263  minRunNum = runNumber;
264  }
265  if (((maxRunNum < runNumber) && (maxExpNum <= expNumber)) ||
266  (maxExpNum < expNumber)) {
267  maxExpNum = expNumber;
268  maxRunNum = runNumber;
269  }
270  }
271 
272  B2INFO("minExpNum = " << minExpNum) ;
273  B2INFO("minRunNum = " << minRunNum) ;
274  B2INFO("maxExpNum = " << maxExpNum) ;
275  B2INFO("maxRunNum = " << maxRunNum) ;
276 
277 
278  if (minExpNum != maxExpNum) {
279  B2ERROR("The runs must all come from the same experiment");
280  return c_Failure;
281  }
282 
283  int experiment = minExpNum;
284 
285 
286  //------------------------------------------------------------------------
287  //------------------------------------------------------------------------
288  //------------------------------------------------------------------------
289  //------------------------------------------------------------------------
290  /* Extract out the time offset information from the database directly.
291  This method loops over all run numbers so it can more easiy pick up
292  old payloads. It is not the preferred method to use if the payloads
293  have iov gaps.*/
294 
295  if (algorithmReadPayloads) {
296  //------------------------------------------------------------------------
297  // Get the input run list (should be only 1) for us to use to update the DBObjectPtrs
298  auto runs = getRunList();
299  /* Take the first run. For the crystal cosmic calibrations, because of the crate
300  calibrations, there is not a known correct run to use within the range. */
301  ExpRun chosenRun = runs.front();
302  B2INFO("merging using the ExpRun (" << chosenRun.second << "," << chosenRun.first << ")");
303  // After here your DBObjPtrs are correct
304  updateDBObjPtrs(1, chosenRun.second, chosenRun.first);
305 
306  //------------------------------------------------------------------------
307  // Test the DBObjects we want to exist and fail if not all of them do.
308  bool allObjectsFound = true;
309 
311  // Check that the payloads we want to merge are sufficiently loaded
312  if (!m_ECLCrystalTimeOffset) {
313  allObjectsFound = false;
314  B2ERROR("No valid DBObject found for 'ECLCrystalTimeOffset'");
315  }
316 
317  // Check that the crate payload is loaded (used for transforming cosmic payload)
318  if (!m_ECLCrateTimeOffset) {
319  allObjectsFound = false;
320  B2ERROR("No valid DBObject found for 'ECLCrateTimeOffset'");
321  }
322 
323  if (!m_refCrysIDzeroingCrate) {
324  allObjectsFound = false;
325  B2ERROR("No valid DBObject found for 'refCrysIDzeroingCrate'");
326  }
327 
328 
329  if (allObjectsFound) {
330  B2INFO("Valid objects found for 'ECLCrystalTimeOffset'");
331  B2INFO("Valid object found for 'ECLCrateTimeOffset'");
332  B2INFO("Valid object found for 'refCrysIDzeroingCrate'");
333  } else {
334  B2INFO("eclTimeShiftsAlgorithm: Exiting with failure. Some missing valid objects.");
335  return c_Failure;
336  }
337 
338 
339  //------------------------------------------------------------------------
341  vector<float> crystalCalib = m_ECLCrystalTimeOffset->getCalibVector();
342  vector<float> crystalCalibUnc = m_ECLCrystalTimeOffset->getCalibUncVector();
343  B2INFO("Loaded 'ECLCrystalTimeOffset' calibrations");
344 
345  vector<float> crateCalib = m_ECLCrateTimeOffset->getCalibVector();
346  vector<float> crateCalibUnc = m_ECLCrateTimeOffset->getCalibUncVector();
347 
348  B2INFO("Loaded 'ECLCrateTimeOffset' calibration with default exp/run");
349 
350  B2INFO("eclTimeShiftsAlgorithm:: loaded ECLCrateTimeOffset from the database"
351  << LogVar("IoV", m_ECLCrateTimeOffset.getIoV())
352  << LogVar("Checksum", m_ECLCrateTimeOffset.getChecksum()));
353 
354  for (int cellID = 1; cellID <= m_numCrystals; cellID += 511) {
355  B2INFO("crystalCalib = " << crystalCalib[cellID - 1]);
356  B2INFO("crateCalib = " << crateCalib[cellID - 1]);
357  }
358 
359  vector<short> refCrystals = m_refCrysIDzeroingCrate->getReferenceCrystals();
360  for (int icrate = 0; icrate < m_numCrates; icrate++) {
361  B2INFO("reference crystal for crate " << icrate + 1 << " = " << refCrystals[icrate]);
362  }
363 
364 
365 
366  //------------------------------------------------------------------------
368  for (int run = minRunNum; run <= maxRunNum; run++) {
369  B2INFO("---------") ;
370  B2INFO("Looking at run " << run) ;
371 
372  vector<int>::iterator it = find(allRunNums.begin(), allRunNums.end(), run);
373  if (it != allRunNums.end()) {
374  int pos = it - allRunNums.begin() ;
375  B2INFO("allRunNums[" << pos << "] = " << allRunNums[pos]);
376  B2INFO("Run " << run << " already processed so skipping it.");
377  continue;
378  } else {
379  B2INFO("New run. Starting to extract information");
380  }
381 
382  // Forloading database for a specific run
383  int eventNumberForCrates = 1;
384 
385  StoreObjPtr<EventMetaData> evtPtr;
386  // simulate the initialize() phase where we can register objects in the DataStore
387  DataStore::Instance().setInitializeActive(true);
388  evtPtr.registerInDataStore();
389  DataStore::Instance().setInitializeActive(false);
390  // now construct the event metadata
391  evtPtr.construct(eventNumberForCrates, run, experiment);
392  // and update the database contents
393  DBStore& dbstore = DBStore::Instance();
394  dbstore.update();
395  // this is only needed it the payload might be intra-run dependent,
396  // that is if it might change during one run as well
397  dbstore.updateEvent();
398  updateDBObjPtrs(eventNumberForCrates, run, experiment);
399 
400 
401  //------------------------------------------------------------------------
403  shared_ptr< ECL::ECLChannelMapper > crystalMapper(new ECL::ECLChannelMapper()) ;
404  crystalMapper->initFromDB();
405 
407  B2INFO("eclTimeShiftsAlgorithm:: loaded ECLCrystalTimeOffset from the database"
408  << LogVar("IoV", m_ECLCrystalTimeOffset.getIoV())
409  << LogVar("Checksum", m_ECLCrystalTimeOffset.getChecksum()));
410  B2INFO("eclTimeShiftsAlgorithm:: loaded ECLCrateTimeOffset from the database"
411  << LogVar("IoV", m_ECLCrateTimeOffset.getIoV())
412  << LogVar("Checksum", m_ECLCrateTimeOffset.getChecksum()));
413 
414 
415  //------------------------------------------------------------------------
418  vector<float> crystalTimeOffsetsCalib;
419  vector<float> crystalTimeOffsetsCalibUnc;
420  crystalTimeOffsetsCalib = m_ECLCrystalTimeOffset->getCalibVector();
421  crystalTimeOffsetsCalibUnc = m_ECLCrystalTimeOffset->getCalibUncVector();
422 
423  vector<float> crateTimeOffsetsCalib;
424  vector<float> crateTimeOffsetsCalibUnc;
425  crateTimeOffsetsCalib = m_ECLCrateTimeOffset->getCalibVector();
426  crateTimeOffsetsCalibUnc = m_ECLCrateTimeOffset->getCalibUncVector();
427 
428  //------------------------------------------------------------------------
432  vector<double> Crate_time_ns(m_numCrates) ;
433  vector<double> Crate_time_tick(m_numCrates) ;
434  vector<double> Crate_time_unc_ns(m_numCrates) ;
435  vector<double> crystalCrate_time_ns(m_numCrates);
436  vector<double> crystalCrate_time_unc_ns(m_numCrates);
437 
438  for (int crysID = 1; crysID <= m_numCrystals; crysID++) {
439  int crateID_temp = crystalMapper->getCrateID(crysID) ;
440  Crate_time_ns[crateID_temp - 1] = crateTimeOffsetsCalib[crysID - 1] * TICKS_TO_NS ;
441  Crate_time_tick[crateID_temp - 1] = crateTimeOffsetsCalib[crysID - 1] ;
442  Crate_time_unc_ns[crateID_temp - 1] = crateTimeOffsetsCalibUnc[crysID - 1] * TICKS_TO_NS ;
443 
444  if (crysID == refCrystals[crateID_temp - 1]) {
445  crystalCrate_time_ns[crateID_temp - 1] = (crystalTimeOffsetsCalib[crysID - 1] +
446  crateTimeOffsetsCalib[crysID - 1]) * TICKS_TO_NS;
447 
448  crystalCrate_time_unc_ns[crateID_temp - 1] = TICKS_TO_NS * sqrt(
449  (crateTimeOffsetsCalibUnc[crysID - 1] * crateTimeOffsetsCalibUnc[crysID - 1]) +
450  (crystalTimeOffsetsCalibUnc[crysID - 1] * crystalTimeOffsetsCalibUnc[crysID - 1])) ;
451  }
452  }
453 
454 
455  for (int iCrate = 0; iCrate < m_numCrates; iCrate++) {
456  double tcrate = Crate_time_ns[iCrate] ;
457  double tcrate_unc = Crate_time_unc_ns[iCrate];
458  double tcrystalCrate = crystalCrate_time_ns[iCrate];
459  double tcrystalCrate_unc = crystalCrate_time_unc_ns[iCrate];
460 
461  if ((tcrate < m_tcrate_max_cut) &&
462  (tcrate > m_tcrate_min_cut) &&
463  (fabs(tcrate_unc) > m_tcrate_unc_min_cut) &&
464  (fabs(tcrate_unc) < m_tcrate_unc_max_cut)) {
465  allCrates_crate_times[iCrate].push_back(tcrate) ;
466  allCrates_run_nums[iCrate].push_back(run) ;
467  allCrates_time_unc[iCrate].push_back(tcrate_unc) ;
468  allCrates_crystalCrate_times[iCrate].push_back(tcrystalCrate) ;
469  allCrates_crystalCrate_times_unc[iCrate].push_back(tcrystalCrate_unc) ;
470 
471  mean_crystalCrate_time_ns[iCrate] += tcrystalCrate ;
472  }
473  }
474 
475 
476  //------------------------------------------------------------------------
478  for (int ic = 0; ic < m_numCrates; ic++) {
479  B2INFO("Crate " << ic + 1 << ", t_crate = " << Crate_time_tick[ic] << " ticks = "
480  << Crate_time_ns[ic] << " +- " << Crate_time_unc_ns[ic]
481  << " ns; t crys+crate (no shift) = " << crystalCrate_time_ns[ic] << " +- "
482  << crystalCrate_time_unc_ns[ic] << " ns") ;
483  }
484 
485  /* Shift the run number to the end of the iov so that we can skip runs
486  that have the payload with the same revision number */
487  int IOV_exp_high = m_ECLCrateTimeOffset.getIoV().getExperimentHigh() ;
488  int IOV_run_high = m_ECLCrateTimeOffset.getIoV().getRunHigh() ;
489  B2INFO("IOV_exp_high = " << IOV_exp_high);
490  B2INFO("IOV_run_high = " << IOV_run_high);
491  if (IOV_run_high == -1) {
492  B2INFO("IOV_run_high is -1 so stop looping over all runs");
493  break;
494  } else {
495  B2INFO("Set run number to higher iov run number");
496  run = IOV_run_high;
497  }
498  B2INFO("now set run = " << run);
499  }
500  }
501 
502 
503 
504 
505  //------------------------------------------------------------------------
506  //------------------------------------------------------------------------
507  //------------------------------------------------------------------------
508  //------------------------------------------------------------------------
511  B2INFO("Shift all run crys+crate+off times. Show the results for a subset of crates/runs:");
512  for (int iCrate = 0; iCrate < m_numCrates; iCrate++) {
513  double mean_time = mean_crystalCrate_time_ns[iCrate] / allCrates_crate_times[iCrate].size() ;
514  B2INFO("Mean crys+crate times for all runs used as offset (crate " << iCrate + 1 << ") = " << mean_time);
515 
516  for (long unsigned int jRun = 0; jRun < allCrates_crate_times[iCrate].size(); jRun++) {
517  allCrates_crystalCrate_times[iCrate][jRun] += -mean_time + timeShiftForPlotStyle[iCrate] ;
518  if (jRun < 50 || iCrate == 1 || iCrate == 40 || iCrate == 51) {
519  B2INFO("allCrates_crystalCrate_times(crate " << iCrate + 1 << ", run counter " << jRun + 1 << ", runNum " <<
520  allCrates_run_nums[iCrate][jRun] << " | after shifting mean) = " <<
521  allCrates_crystalCrate_times[iCrate][jRun]);
522  }
523  }
524  }
525 
526 
527 
528  //------------------------------------------------------------------------
529  //------------------------------------------------------------------------
532  TFile* tcratefile = 0;
533 
534  B2INFO("Debug output rootfile: " << debugFilenameBase);
535  string runNumsString = string("_") + to_string(minExpNum) + "_" + to_string(minRunNum) + string("-") +
536  to_string(maxExpNum) + "_" + to_string(maxRunNum);
537  string debugFilename = debugFilenameBase + runNumsString + string(".root");
538  TString fname = debugFilename;
539 
540  tcratefile = new TFile(fname, "recreate");
541  tcratefile->cd();
542  B2INFO("Debugging histograms written to " << fname);
543 
544  for (int i = 0; i < m_numCrates; i++) {
545  B2INFO("Starting to make crate time jump plots for crate " << i + 1);
546  shared_ptr< TCanvas > cSmart(new TCanvas);
547 
548  Double_t* single_crate_crate_times = &allCrates_crate_times[i][0] ;
549  Double_t* single_crate_run_nums = &allCrates_run_nums[i][0] ;
550  Double_t* single_crate_time_unc = &allCrates_time_unc[i][0] ;
551  Double_t* single_crate_crystalCrate_times = &allCrates_crystalCrate_times[i][0] ;
552  Double_t* single_crate_crystalCrate_times_unc = &allCrates_crystalCrate_times_unc[i][0] ;
553  B2INFO("Done setting up the arrays for the crate " << i + 1);
554 
555  ostringstream ss;
556  ss << setw(2) << setfill('0') << i + 1 ;
557  string paddedCrateID(ss.str());
558 
559  // ----- crate time constants vs run number ------
560  shared_ptr< TGraphErrors > g_tcrate_vs_runNum(new TGraphErrors(allCrates_crate_times[i].size(), single_crate_run_nums,
561  single_crate_crate_times, NULL, single_crate_time_unc)) ;
562  // NULL for run number errors = 0 for all
563 
564  string tgraph_title = string("e") + to_string(minExpNum) + string("r") + to_string(minRunNum) +
565  string("-e") + to_string(maxExpNum) + string("r") + to_string(maxRunNum) ;
566 
567  string tgraph_name_short = "crateTimeVSrunNum_" ;
568  tgraph_name_short = tgraph_name_short + runNumsString + "_crate";
569 
570  tgraph_title = tgraph_title + string("_crate") + paddedCrateID ;
571  tgraph_name_short = tgraph_name_short + paddedCrateID ;
572  tgraph_title = tgraph_title + string(" (") + to_string(m_tcrate_min_cut) + string(" < tcrate < ") +
573  to_string(m_tcrate_max_cut) + string(" ns, ") + to_string(m_tcrate_unc_min_cut) +
574  string(" < tcrate unc. < ") + to_string(m_tcrate_unc_max_cut) + string(" ns cuts)") ;
575 
576  g_tcrate_vs_runNum->SetName(tgraph_name_short.c_str()) ;
577  g_tcrate_vs_runNum->SetTitle(tgraph_title.c_str()) ;
578  g_tcrate_vs_runNum->GetXaxis()->SetTitle("Run number") ;
579  g_tcrate_vs_runNum->GetYaxis()->SetTitle("Crate time [ns]") ;
580 
581  g_tcrate_vs_runNum->GetYaxis()->SetRangeUser(m_tcrate_min_cut, m_tcrate_max_cut) ;
582 
583  g_tcrate_vs_runNum->Draw("AP") ;
584  g_tcrate_vs_runNum->SetMarkerSize(0.8) ;
585  g_tcrate_vs_runNum->Draw("AP") ;
586 
587  shared_ptr< TLatex > Leg1(new TLatex);
588  Leg1->SetNDC();
589  Leg1->SetTextAlign(11);
590  Leg1->SetTextFont(42);
591  Leg1->SetTextSize(0.035);
592  Leg1->SetTextColor(1);
593  Leg1->AppendPad();
594 
595  g_tcrate_vs_runNum->Write() ;
596  cSmart->SaveAs((tgraph_name_short + string(".pdf")).c_str()) ;
597 
598  B2INFO("Saved pdf: " << tgraph_name_short << ".pdf");
599 
600 
601  // ----- crystal + crate time constants + offset vs run number ------
602  shared_ptr< TGraphErrors > g_crateCrystalTime_vs_runNum(new TGraphErrors(allCrates_crystalCrate_times[i].size(),
603  single_crate_run_nums,
604  single_crate_crystalCrate_times, NULL, single_crate_crystalCrate_times_unc)) ;
605 
606  tgraph_title = string("e") + to_string(minExpNum) + string("r") + to_string(minRunNum) +
607  string("-e") + to_string(maxExpNum) + string("r") + to_string(maxRunNum) ;
608 
609  tgraph_name_short = "crystalCrateTimeVSrunNum_" ;
610  tgraph_name_short = tgraph_name_short + runNumsString + "_crate";
611 
612  tgraph_title = tgraph_title + string("_crate") + paddedCrateID ;
613  tgraph_name_short = tgraph_name_short + paddedCrateID ;
614  tgraph_title = tgraph_title + string(" (") + to_string(m_tcrate_min_cut) + string(" < tcrate < ") +
615  to_string(m_tcrate_max_cut) + string(" ns, ") + to_string(m_tcrate_unc_min_cut) +
616  string(" < tcrate unc. < ") + to_string(m_tcrate_unc_max_cut) + string(" ns cuts)") ;
617 
618 
619  g_crateCrystalTime_vs_runNum->SetName(tgraph_name_short.c_str()) ;
620  g_crateCrystalTime_vs_runNum->SetTitle(tgraph_title.c_str()) ;
621  g_crateCrystalTime_vs_runNum->GetXaxis()->SetTitle("Run number") ;
622  g_crateCrystalTime_vs_runNum->GetYaxis()->SetTitle("Crate time + Crystal time + centring overall offset [ns]") ;
623 
624  g_crateCrystalTime_vs_runNum->GetYaxis()->SetRangeUser(crysCrateShift_min, crysCrateShift_max) ;
625 
626  g_crateCrystalTime_vs_runNum->Draw("AP") ;
627  g_crateCrystalTime_vs_runNum->SetMarkerSize(0.8) ;
628  g_crateCrystalTime_vs_runNum->Draw("AP") ;
629 
630  g_crateCrystalTime_vs_runNum->Write() ;
631  cSmart->SaveAs((tgraph_name_short + string(".pdf")).c_str()) ;
632 
633  B2INFO("Saved pdf: " << tgraph_name_short << ".pdf");
634 
635  // ----- crystal + crate time constants + offset vs run counter------
636  // This will remove gaps and ignore the actual run number
637 
638  /* Define a vector to store a renumbering of the run numbers, incrementing
639  by +1 so that there are no gaps. The runs are not in order so the
640  run numbers&indices first have to be sorted before the "run counter"
641  numbers can used.*/
642  int numRunsWithCrateTimes = allCrates_crystalCrate_times[i].size();
643  vector<Double_t> counterVec(numRunsWithCrateTimes);
644 
645 
646  // Vector to store element
647  // with respective present index
648  vector<pair<int, double> > runNum_index_pairs;
649 
650  // Inserting element in pair vector
651  // to keep track of previous indexes
652  for (int pairIndex = 0; pairIndex < numRunsWithCrateTimes; pairIndex++) {
653  runNum_index_pairs.push_back(make_pair(allCrates_run_nums[i][pairIndex], pairIndex));
654  }
655 
656  B2INFO("Crate id = " << i + 1);
657  B2INFO("Unsorted run numbers");
658  for (int runCounter = 0; runCounter < numRunsWithCrateTimes; runCounter++) {
659  B2INFO("Run number, run number vector index = " << runNum_index_pairs[runCounter].first << ", " <<
660  runNum_index_pairs[runCounter].second);
661  }
662 
663  // Sorting pair vector
664  sort(runNum_index_pairs.begin(), runNum_index_pairs.end());
665 
666  // Fill the run counter vector
667  for (int runCounter = 0; runCounter < numRunsWithCrateTimes; runCounter++) {
668  counterVec[runNum_index_pairs[runCounter].second] = runCounter + 1;
669  }
670 
671  B2INFO("Run numbers with index and times");
672  for (int runCounter = 0; runCounter < numRunsWithCrateTimes; runCounter++) {
673  int idx = (int) round(counterVec[runCounter]);
674  B2INFO("Vector index, Run number, run number sorting order index, tcrystal+tcrate+shifts = " << runCounter << ", " <<
675  allCrates_run_nums[i][runCounter] << ", " << idx << ", " << single_crate_crystalCrate_times[idx - 1] << " ns");
676  }
677 
678 
679  if (numRunsWithCrateTimes > 0) {
680  shared_ptr< TGraphErrors > g_crateCrystalTime_vs_runCounter(new TGraphErrors(numRunsWithCrateTimes, &counterVec[0],
681  single_crate_crystalCrate_times, NULL, single_crate_crystalCrate_times_unc)) ;
682 
683  tgraph_title = string("e") + to_string(minExpNum) + string("r") + to_string(minRunNum) +
684  string("-e") + to_string(maxExpNum) + string("r") + to_string(maxRunNum) ;
685 
686 
687  tgraph_name_short = "crystalCrateTimeVSrunCounter_" ;
688  tgraph_name_short = tgraph_name_short + runNumsString + "_crate";
689 
690 
691  tgraph_title = tgraph_title + string("_crate") + paddedCrateID ;
692  tgraph_name_short = tgraph_name_short + paddedCrateID ;
693  tgraph_title = tgraph_title + string(" (") + to_string(m_tcrate_min_cut) + string(" < tcrate < ") +
694  to_string(m_tcrate_max_cut) + string(" ns, ") + to_string(m_tcrate_unc_min_cut) +
695  string(" < tcrate unc. < ") + to_string(m_tcrate_unc_max_cut) + string(" ns cuts)") ;
696 
697 
698  g_crateCrystalTime_vs_runCounter->SetName(tgraph_name_short.c_str()) ;
699  g_crateCrystalTime_vs_runCounter->SetTitle(tgraph_title.c_str()) ;
700  g_crateCrystalTime_vs_runCounter->GetXaxis()->SetTitle("Run counter (remove gaps from run numbers)") ;
701  g_crateCrystalTime_vs_runCounter->GetYaxis()->SetTitle("Crate time + Crystal time + centring overall offset [ns]") ;
702 
703  g_crateCrystalTime_vs_runCounter->GetYaxis()->SetRangeUser(crysCrateShift_min, crysCrateShift_max) ;
704  g_crateCrystalTime_vs_runCounter->GetXaxis()->SetRangeUser(0, numRunsWithCrateTimes + 1) ;
705 
706  g_crateCrystalTime_vs_runCounter->Draw("AP") ;
707  g_crateCrystalTime_vs_runCounter->SetMarkerSize(0.8) ;
708  g_crateCrystalTime_vs_runCounter->Draw("AP") ;
709 
710  g_crateCrystalTime_vs_runCounter->Write() ;
711  cSmart->SaveAs((tgraph_name_short + string(".pdf")).c_str()) ;
712  B2INFO("Saved pdf: " << tgraph_name_short << ".pdf");
713 
714  B2INFO("Finished making crate time jump plots for crate " << i + 1);
715  } else {
716  B2INFO("Crate " << i + 1 << " has no entries that pass all the cuts so no crystalCrateTimeVSrunCounter_crate plot will be made.");
717  }
718  }
719 
720 
721 
722 
723  /* Loop over all the runs and crates and let the user know when a crate time jump
724  has occurred. Jumps can be of various sizes so have different thresholds. */
725  double smallThreshold = 1 ; //ns
726  double largeThreshold = 6.5 ; //ns
727 
728  B2INFO("======================= Crate time jumps =========================");
729  B2INFO("======================= Small threshold jumps ====================");
730  B2INFO("Crate ID = 1..52");
731  B2INFO("==================================================================");
732 
733  for (int i = 0; i < m_numCrates; i++) {
734  int numRunsWithCrateTimes = allCrates_crystalCrate_times[i].size();
735  for (int runCounter = 0; runCounter < numRunsWithCrateTimes - 1; runCounter++) {
736  int run_i = allCrates_run_nums[i][runCounter] ;
737  int run_f = allCrates_run_nums[i][runCounter + 1] ;
738  double time_i = allCrates_crystalCrate_times[i][runCounter] ;
739  double time_f = allCrates_crystalCrate_times[i][runCounter + 1] ;
740 
741  if (fabs(time_f - time_i) > smallThreshold) {
742  B2INFO("Crate " << i + 1 << " has crate time jump > " << smallThreshold << " ns: t(run " << run_f << ") = " << time_f <<
743  " ns - t(run " << run_i << ") = " << time_i << " ns = " << time_f - time_i);
744  }
745  }
746  }
747 
748 
749  B2INFO("~~~~~~~~~~~~~~~~~~~~~~~ Large threshold jumps ~~~~~~~~~~~~~~~~~~~~");
750  B2INFO("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
751 
752  for (int i = 0; i < m_numCrates; i++) {
753  int numRunsWithCrateTimes = allCrates_crystalCrate_times[i].size();
754  for (int runCounter = 0; runCounter < numRunsWithCrateTimes - 1; runCounter++) {
755  int run_i = allCrates_run_nums[i][runCounter] ;
756  int run_f = allCrates_run_nums[i][runCounter + 1] ;
757  double time_i = allCrates_crystalCrate_times[i][runCounter] ;
758  double time_f = allCrates_crystalCrate_times[i][runCounter + 1] ;
759 
760  if (fabs(time_f - time_i) > largeThreshold) {
761  B2INFO("WARNING: Crate " << i + 1 << " has crate time jump > " << largeThreshold << " ns: t(run " << run_f << ") = " << time_f <<
762  " ns - t(run " << run_i << ") = " << time_i << " ns = " << time_f - time_i);
763  }
764  }
765  }
766 
767 
768 
769 
770  // Just in case, we remember the current TDirectory so we can return to it
771  TDirectory* executeDir = gDirectory;
772 
773  tcratefile->Write();
774  tcratefile->Close();
775  // Go back to original TDirectory
776  executeDir->cd();
777 
778  return c_OK;
779 }
Belle2::ECL::eclTimeShiftsAlgorithm::m_crystalTimeConst
Double_t m_crystalTimeConst
Crystal time calibration constant.
Definition: eclTimeShiftsAlgorithm.h:91
Belle2::ECL::eclTimeShiftsAlgorithm::m_tcrate_max_cut
double m_tcrate_max_cut
Maximum value cut for the crate time calibration constant for plotting
Definition: eclTimeShiftsAlgorithm.h:103
Belle2::ECL::eclTimeShiftsAlgorithm::m_tcrate_unc_min_cut
double m_tcrate_unc_min_cut
Minimum value cut for the crate time calibration constant uncertainty for plotting.
Definition: eclTimeShiftsAlgorithm.h:104
Belle2::ECL::eclTimeShiftsAlgorithm::crysCrateShift_max
double crysCrateShift_max
Plotting time max for crystal+crate shift plots.
Definition: eclTimeShiftsAlgorithm.h:58
Belle2::ECL::eclTimeShiftsAlgorithm::m_refCrystalID
Int_t m_refCrystalID
Crystal ID number for the reference crystal.
Definition: eclTimeShiftsAlgorithm.h:99
Belle2::ECL::eclTimeShiftsAlgorithm::m_ECLCrateTimeOffset
DBObjPtr< ECLCrystalCalib > m_ECLCrateTimeOffset
ECLCrateTimeOffset payload that we want to read from the DB.
Definition: eclTimeShiftsAlgorithm.h:75
Belle2::DataStore::Instance
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:54
Belle2::ECL::ECLChannelMapper::getCrateID
int getCrateID(int iCOPPERNode, int iFINESSE)
get crate number by given COPPER node number and FINESSE number
Definition: ECLChannelMapper.cc:211
Belle2::DataStore::setInitializeActive
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:94
Belle2::ECL::eclTimeShiftsAlgorithm::m_exp_perCrystal
Int_t m_exp_perCrystal
Experiment number
Definition: eclTimeShiftsAlgorithm.h:88
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::updateDBObjPtrs
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
Definition: CalibrationAlgorithm.cc:398
Belle2::ECL::eclTimeShiftsAlgorithm::crysCrateShift_min
double crysCrateShift_min
Plotting time min for crystal+crate shift plots.
Definition: eclTimeShiftsAlgorithm.h:57
Belle2::ECL::eclTimeShiftsAlgorithm::m_run_perCrystal
Int_t m_run_perCrystal
Run number.
Definition: eclTimeShiftsAlgorithm.h:87
Belle2::ECL::eclTimeShiftsAlgorithm::m_crystalID
Int_t m_crystalID
Crystal ID number.
Definition: eclTimeShiftsAlgorithm.h:89
Belle2::ECL::eclTimeShiftsAlgorithm::algorithmReadPayloads
bool algorithmReadPayloads
Whether or not to have the algorithm code to loop over all the runs and read the payloads itself.
Definition: eclTimeShiftsAlgorithm.h:62
Belle2::ECL::eclTimeShiftsAlgorithm::m_tcrate_unc_max_cut
double m_tcrate_unc_max_cut
Maximum value cut for the crate time calibration constant uncertainty for plotting.
Definition: eclTimeShiftsAlgorithm.h:105
Belle2::ECL::eclTimeShiftsAlgorithm::m_tcrate_min_cut
double m_tcrate_min_cut
Minimum value cut for the crate time calibration constant for plotting.
Definition: eclTimeShiftsAlgorithm.h:102
Belle2::ECL::ECLChannelMapper
This class provides access to ECL channel map that is either a) Loaded from the database (see ecl/dbo...
Definition: ECLChannelMapper.h:36
Belle2::DBStore
Singleton class to cache database objects.
Definition: DBStore.h:42
Belle2::StoreObjPtr::construct
bool construct(Args &&... params)
Construct an object of type T in this StoreObjPtr, using the provided constructor arguments.
Definition: StoreObjPtr.h:128
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECL::eclTimeShiftsAlgorithm::m_ECLCrystalTimeOffset
DBObjPtr< ECLCrystalCalib > m_ECLCrystalTimeOffset
ECLCrystalTimeOffset payload that we want to read from the DB.
Definition: eclTimeShiftsAlgorithm.h:72
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
Belle2::CalibrationAlgorithm::c_Failure
@ c_Failure
Failed =3 in Python.
Definition: CalibrationAlgorithm.h:54
Belle2::DBStore::Instance
static DBStore & Instance()
Instance of a singleton DBStore.
Definition: DBStore.cc:36
Belle2::ECL::EclConfiguration::m_rf
static constexpr double m_rf
accelerating RF, http://ptep.oxfordjournals.org/content/2013/3/03A006.full.pdf
Definition: EclConfiguration.h:45
Belle2::ExpRun
Struct containing exp number and run number.
Definition: Splitter.h:55
Belle2::ECL::eclTimeShiftsAlgorithm::timeShiftForPlotStyle
double timeShiftForPlotStyle[52]
List of time offsets, one per crate, used just to centre the time constants around zero.
Definition: eclTimeShiftsAlgorithm.h:55
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::DBStore::updateEvent
void updateEvent()
Updates all intra-run dependent objects.
Definition: DBStore.cc:150
Belle2::ECL::eclTimeShiftsAlgorithm::m_numCrates
const int m_numCrates
Number of Crates expected.
Definition: eclTimeShiftsAlgorithm.h:84
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition: CalibrationAlgorithm.h:276
Belle2::ECL::ECLChannelMapper::initFromDB
bool initFromDB()
Initialize channel mapper from the conditions database.
Definition: ECLChannelMapper.cc:105
Belle2::ECL::eclTimeShiftsAlgorithm::debugFilenameBase
std::string debugFilenameBase
Name of file with debug output, eclTimeShiftsAlgorithm.root by default.
Definition: eclTimeShiftsAlgorithm.h:51
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::ECL::eclTimeShiftsAlgorithm::calibrate
EResult calibrate() override
..Run algorithm
Definition: eclTimeShiftsAlgorithm.cc:49
Belle2::ECL::eclTimeShiftsAlgorithm::m_crateID
Int_t m_crateID
Crate ID number.
Definition: eclTimeShiftsAlgorithm.h:94
Belle2::ECL::eclTimeShiftsAlgorithm::m_crystalTimeUnc
Double_t m_crystalTimeUnc
Uncertainty on the crystal time calibration constant.
Definition: eclTimeShiftsAlgorithm.h:93
Belle2::DBStore::update
void update()
Updates all objects that are outside their interval of validity.
Definition: DBStore.cc:87
Belle2::ECL::eclTimeShiftsAlgorithm::m_crateTimeConst
Double_t m_crateTimeConst
Crate time calibration constant.
Definition: eclTimeShiftsAlgorithm.h:90
Belle2::ECL::eclTimeShiftsAlgorithm::m_numCrystals
const int m_numCrystals
Number of Crystals expected.
Definition: eclTimeShiftsAlgorithm.h:81
Belle2::ECL::eclTimeShiftsAlgorithm::m_crateTimeUnc
Double_t m_crateTimeUnc
Uncertainty on the crate time calibration constant.
Definition: eclTimeShiftsAlgorithm.h:92
Belle2::ECL::eclTimeShiftsAlgorithm::m_refCrysIDzeroingCrate
DBObjPtr< ECLReferenceCrystalPerCrateCalib > m_refCrysIDzeroingCrate
payload that we want to read from the DB
Definition: eclTimeShiftsAlgorithm.h:78