Belle II Software  release-06-02-00
ECLCalibrationVariables.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 include
10 #include <ecl/variables/ECLCalibrationVariables.h>
11 
12 // analysis
13 #include <analysis/dataobjects/Particle.h>
14 #include <analysis/dataobjects/ParticleList.h>
15 #include <analysis/dataobjects/ECLEnergyCloseToTrack.h>
16 #include <analysis/dataobjects/ECLTRGInformation.h>
17 #include <analysis/dataobjects/ECLTriggerCell.h>
18 #include <analysis/utility/ReferenceFrame.h>
19 #include <analysis/ClusterUtility/ClusterUtils.h>
20 #include <analysis/VariableManager/Utility.h>
21 
22 // framework
23 #include <framework/core/Module.h>
24 #include <framework/logging/Logger.h>
25 #include <framework/datastore/StoreArray.h>
26 
27 // dataobjects
28 #include <mdst/dataobjects/ECLCluster.h>
29 #include <mdst/dataobjects/Track.h>
30 
31 using namespace std;
32 
33 namespace Belle2 {
39  namespace Variable {
40 
41  double eclEnergy3FWDBarrel(const Particle* particle)
42  {
43 
44  const Track* track = particle->getTrack();
45  if (track) {
46 
47  auto* eclinfo = track->getRelatedTo<ECLEnergyCloseToTrack>();
48 
49  if (eclinfo) {
50  return eclinfo->getEnergy3FWDBarrel();
51  } else {
52  B2WARNING("Relation to ECLEnergyCloseToTrack not found, did you forget to run ECLTrackCalDigitMatchModule?");
53  return std::numeric_limits<float>::quiet_NaN();
54  }
55  }
56 
57  return std::numeric_limits<float>::quiet_NaN();
58  }
59 
60  double eclEnergy3FWDEndcap(const Particle* particle)
61  {
62 
63  const Track* track = particle->getTrack();
64  if (track) {
65 
66  auto* eclinfo = track->getRelatedTo<ECLEnergyCloseToTrack>();
67 
68  if (eclinfo) {
69  return eclinfo->getEnergy3FWDEndcap();
70  } else {
71  B2WARNING("Relation to ECLEnergyCloseToTrack not found, did you forget to run ECLTrackCalDigitMatchModule?");
72  return std::numeric_limits<float>::quiet_NaN();
73  }
74  }
75 
76  return std::numeric_limits<float>::quiet_NaN();
77  }
78 
79  double eclEnergy3BWDBarrel(const Particle* particle)
80  {
81 
82  const Track* track = particle->getTrack();
83  if (track) {
84 
85  auto* eclinfo = track->getRelatedTo<ECLEnergyCloseToTrack>();
86 
87  if (eclinfo) {
88  return eclinfo->getEnergy3BWDBarrel();
89  } else {
90  B2WARNING("Relation to ECLEnergyCloseToTrack not found, did you forget to run ECLTrackCalDigitMatchModule?");
91  return std::numeric_limits<float>::quiet_NaN();
92  }
93  }
94 
95  return std::numeric_limits<float>::quiet_NaN();
96  }
97 
98  double eclEnergy3BWDEndcap(const Particle* particle)
99  {
100  const Track* track = particle->getTrack();
101  if (track) {
102 
103  auto* eclinfo = track->getRelatedTo<ECLEnergyCloseToTrack>();
104 
105  if (eclinfo) {
106  return eclinfo->getEnergy3BWDEndcap();
107  } else {
108  B2WARNING("Relation to ECLEnergyCloseToTrack not found, did you forget to run ECLTrackCalDigitMatchModule?");
109  return std::numeric_limits<float>::quiet_NaN();
110  }
111  }
112 
113  return std::numeric_limits<float>::quiet_NaN();
114  }
115 
116  double eclNumberOfTCsForCluster(const Particle* particle, const std::vector<double>& vars)
117  {
118  if (vars.size() != 4) {
119  B2FATAL("Need exactly four parameters (minthetaid, maxthetaid, minhitwindow, maxhitwindow).");
120  }
121 
122  // if we did not run the ECLTRGInformation module, return NaN
123  // cppcheck-suppress unassignedVariable
124  StoreArray<ECLTriggerCell> ecltc;
125  if (!ecltc) return std::numeric_limits<double>::quiet_NaN();
126 
127  // if theta range makes no sense, return NaN
128  const int minTheta = int(std::lround(vars[0]));
129  const int maxTheta = int(std::lround(vars[1]));
130  if (maxTheta < minTheta) {
131  B2WARNING("minTheta i (vars[0]) must be equal or less than maxTheta j (vars[1]).");
132  return std::numeric_limits<double>::quiet_NaN();
133  }
134  // if hitwin range makes no sense, return NaN
135  const int minHitWin = int(std::lround(vars[2]));
136  const int maxHitWin = int(std::lround(vars[3]));
137  if (maxHitWin < minHitWin) {
138  B2WARNING("minHitWin k (vars[2]) must be equal or less than maxHitWin l (vars[3]).");
139  return std::numeric_limits<double>::quiet_NaN();
140  }
141 
142  double result = 0.;
143  const ECLCluster* cluster = particle->getECLCluster();
144 
145  // if everything else is fine, but we don't have a cluster, return 0
146  if (cluster) {
147  auto relationsTCs = cluster->getRelationsWith<ECLTriggerCell>();
148  for (unsigned int idxTC = 0; idxTC < relationsTCs.size(); ++idxTC) {
149  const auto tc = relationsTCs.object(idxTC);
150  if (tc->getThetaId() >= minTheta and tc->getThetaId() <= maxTheta
151  and tc->getHitWin() >= minHitWin and tc->getHitWin() <= maxHitWin) result += 1.0;
152  }
153  }
154  return result;
155  }
156 
157  double eclTCFADCForCluster(const Particle* particle, const std::vector<double>& vars)
158  {
159  if (vars.size() != 4) {
160  B2FATAL("Need exactly four parameters (minthetaid, maxthetaid, minhitwindow, maxhitwindow).");
161  }
162 
163  // if we did not run the ECLTRGInformation module, return NaN
164  // cppcheck-suppress unassignedVariable
165  StoreArray<ECLTriggerCell> ecltc;
166  if (!ecltc) return std::numeric_limits<double>::quiet_NaN();
167 
168  // if theta range makes no sense, return NaN
169  const int minTheta = int(std::lround(vars[0]));
170  const int maxTheta = int(std::lround(vars[1]));
171  if (maxTheta < minTheta) {
172  B2WARNING("minTheta i (vars[0]) must be equal or less than maxTheta j (vars[1]).");
173  return std::numeric_limits<double>::quiet_NaN();
174  }
175 
176  // if hitwin range makes no sense, return NaN
177  const int minHitWin = int(std::lround(vars[2]));
178  const int maxHitWin = int(std::lround(vars[3]));
179  if (maxHitWin < minHitWin) {
180  B2WARNING("minHitWin k (vars[2]) must be equal or less than maxHitWin l (vars[3]).");
181  return std::numeric_limits<double>::quiet_NaN();
182  }
183 
184  double result = 0.;
185  const ECLCluster* cluster = particle->getECLCluster();
186 
187  // if everything else is fine, but we don't have a cluster, return 0
188  if (cluster) {
189  auto relationsTCs = cluster->getRelationsTo<ECLTriggerCell>();
190  for (unsigned int idxTC = 0; idxTC < relationsTCs.size(); ++idxTC) {
191  const auto tc = relationsTCs.object(idxTC);
192  if (tc->getThetaId() >= minTheta and tc->getThetaId() <= maxTheta
193  and tc->getHitWin() >= minHitWin and tc->getHitWin() <= maxHitWin) result += tc->getFADC();
194  }
195  }
196  return result;
197  }
198 
199  double eclTCIsMaximumForCluster(const Particle* particle)
200  {
201 
202  // if we did not run the ECLTRGInformation module, return NaN
203  StoreArray<ECLTriggerCell> ecltc;
204  if (!ecltc) return std::numeric_limits<double>::quiet_NaN();
205 
206  double result = 0.;
207  const ECLCluster* cluster = particle->getECLCluster();
208 
209  // if everything else is fine, but we don't have a cluster, return 0
210  if (cluster) {
211  auto relationsTCs = cluster->getRelationsTo<ECLTriggerCell>();
212  for (unsigned int idxTC = 0; idxTC < relationsTCs.size(); ++idxTC) {
213  const auto tc = relationsTCs.object(idxTC);
214  if (tc->isHighestFADC()) result = 1.0;
215  }
216  }
217  return result;
218  }
219 
220  double eclClusterTrigger(const Particle* particle)
221  {
222  const ECLCluster* cluster = particle->getECLCluster();
223  if (cluster) {
224  const bool matcher = cluster->hasTriggerClusterMatching();
225 
226  if (matcher) {
227  return cluster->isTriggerCluster();
228  } else {
229  B2WARNING("Particle has an associated ECLCluster but the ECLTriggerClusterMatcher module has not been run!");
230  return std::numeric_limits<float>::quiet_NaN();
231  }
232  }
233  return std::numeric_limits<float>::quiet_NaN();
234  }
235 
236  double getEnergyTC(const Particle*, const std::vector<double>& vars)
237  {
238  if (vars.size() != 1) {
239  B2FATAL("Need exactly one parameter (tcid).");
240  }
241 
242  StoreObjPtr<ECLTRGInformation> tce;
243  const int tcid = int(std::lround(vars[0]));
244 
245  if (!tce) return std::numeric_limits<double>::quiet_NaN();
246  return tce->getEnergyTC(tcid);
247  }
248 
249  double getEnergyTCECLCalDigit(const Particle*, const std::vector<double>& vars)
250  {
251  if (vars.size() != 1) {
252  B2FATAL("Need exactly one parameter (tcid).");
253  }
254 
255  StoreObjPtr<ECLTRGInformation> tce;
256  const int tcid = int(std::lround(vars[0]));
257 
258  if (!tce) return std::numeric_limits<double>::quiet_NaN();
259  return tce->getEnergyTCECLCalDigit(tcid);
260  }
261 
262  double getTimingTC(const Particle*, const std::vector<double>& vars)
263  {
264  if (vars.size() != 1) {
265  B2FATAL("Need exactly one parameter (tcid).");
266  }
267 
268  StoreObjPtr<ECLTRGInformation> tce;
269  const int tcid = int(std::lround(vars[0]));
270 
271  if (!tce) return std::numeric_limits<double>::quiet_NaN();
272  return tce->getTimingTC(tcid);
273  }
274 
275  double eclHitWindowTC(const Particle*, const std::vector<double>& vars)
276  {
277  if (vars.size() != 1) {
278  B2FATAL("Need exactly one parameter (tcid).");
279  }
280 
281  StoreObjPtr<ECLTRGInformation> tce;
282  const int tcid = int(std::lround(vars[0]));
283 
284  if (!tce) return std::numeric_limits<double>::quiet_NaN();
285  return tce->getHitWinTC(tcid);
286  }
287 
288  double getEvtTimingTC(const Particle*)
289  {
290  StoreObjPtr<ECLTRGInformation> tce;
291  if (!tce) return std::numeric_limits<double>::quiet_NaN();
292  return tce->getEvtTiming();
293  }
294 
295  double getMaximumTCId(const Particle*)
296  {
297  StoreObjPtr<ECLTRGInformation> tce;
298  if (!tce) return std::numeric_limits<double>::quiet_NaN();
299  return tce->getMaximumTCId();
300  }
301 
302  double getTimingTCECLCalDigit(const Particle*, const std::vector<double>& vars)
303  {
304  if (vars.size() != 1) {
305  B2FATAL("Need exactly one parameter (tcid).");
306  }
307 
308  StoreObjPtr<ECLTRGInformation> tce;
309  const int tcid = int(std::lround(vars[0]));
310 
311  if (!tce) return std::numeric_limits<double>::quiet_NaN();
312  return tce->getTimingTCECLCalDigit(tcid);
313  }
314 
315  double getNumberOfTCs(const Particle*, const std::vector<double>& vars)
316  {
317  if (vars.size() != 3) {
318  B2FATAL("Need exactly three parameters (fadccut, minthetaid, maxthetaid).");
319  }
320 
321  StoreArray<ECLTriggerCell> ecltcs;
322  const int fadccut = int(std::lround(vars[0]));
323 
324  if (!ecltcs) return std::numeric_limits<double>::quiet_NaN();
325  if (fadccut == 0) return ecltcs.getEntries();
326  else {
327  int minTheta = int(std::lround(vars[1]));
328  int maxTheta = int(std::lround(vars[2]));
329 
330  unsigned nTCs = 0;
331  for (const auto& tc : ecltcs) {
332  if (tc.getFADC() >= fadccut and
333  tc.getThetaId() >= minTheta and
334  tc.getThetaId() <= maxTheta) nTCs++;
335  }
336  return nTCs;
337  }
338  return 0.0;
339  }
340 
341  double eclEnergySumTC(const Particle*, const std::vector<double>& vars)
342  {
343  if (vars.size() != 3) {
344  B2FATAL("Need exactly three parameters (fadccut, minthetaid, maxthetaid).");
345  }
346 
347  StoreObjPtr<ECLTRGInformation> tce;
348  if (!tce) return std::numeric_limits<double>::quiet_NaN();
349 
350  const int fadccut = int(std::lround(vars[0]));
351  const int minTheta = int(std::lround(vars[1]));
352  const int maxTheta = int(std::lround(vars[2]));
353 
354  if (maxTheta < minTheta) {
355  B2WARNING("minTheta i (vars[1]) must be equal or less than maxTheta j (vars[2]).");
356  return std::numeric_limits<double>::quiet_NaN();
357  }
358 
359  double energySum = 0.;
360  for (unsigned idx = 1; idx <= 576; idx++) {
361  if (tce->getThetaIdTC(idx) >= minTheta and tce->getThetaIdTC(idx) <= maxTheta and tce->getEnergyTC(idx) >= fadccut) {
362  energySum += tce->getEnergyTC(idx);
363  }
364  }
365 
366  return energySum;
367  }
368 
369  double eclEnergySumTCECLCalDigit(const Particle*, const std::vector<double>& vars)
370  {
371  if (vars.size() != 3) {
372  B2FATAL("Need exactly three parameters (minthetaid, maxthetaid, option).");
373  }
374 
375  StoreObjPtr<ECLTRGInformation> tce;
376  if (!tce) return std::numeric_limits<double>::quiet_NaN();
377 
378  int minTheta = int(std::lround(vars[0]));
379  int maxTheta = int(std::lround(vars[1]));
380  int option = int(std::lround(vars[2]));
381  double par = vars[3];
382 
383  if (maxTheta < minTheta) {
384  B2WARNING("minTheta i (vars[0]) must be equal or less than maxTheta j (vars[1]).");
385  return std::numeric_limits<double>::quiet_NaN();
386  }
387  if (option < 0 or option > 2) {
388  B2WARNING("Third parameters k (vars[2]) must be >=0 and <=2.");
389  return std::numeric_limits<double>::quiet_NaN();
390  }
391 
392  double energySum = 0.;
393  for (unsigned idx = 1; idx <= 576; idx++) {
394  if (tce->getThetaIdTC(idx) >= minTheta and
395  tce->getThetaIdTC(idx) <= maxTheta) {
396 
397  if (option == 0) {
398  energySum += tce->getEnergyTCECLCalDigit(idx);
399  } else if (option == 1 and tce->getEnergyTC(idx)) {
400  energySum += tce->getEnergyTCECLCalDigit(idx);
401  } else if (option == 2 and tce->getEnergyTCECLCalDigit(idx) > par) { // TCECLCalDigits > par
402  energySum += tce->getEnergyTCECLCalDigit(idx);
403  }
404  }
405  }
406 
407  return energySum;
408  }
409 
410  double eclEnergySumTCECLCalDigitInECLCluster(const Particle*)
411  {
412  StoreObjPtr<ECLTRGInformation> tce;
413  if (!tce) return std::numeric_limits<double>::quiet_NaN();
414  return tce->getSumEnergyTCECLCalDigitInECLCluster();
415  }
416 
417  double eclEnergySumECLCalDigitInECLCluster(const Particle*)
418  {
419  StoreObjPtr<ECLTRGInformation> tce;
420  if (!tce) return std::numeric_limits<double>::quiet_NaN();
421  return tce->getSumEnergyECLCalDigitInECLCluster();
422  }
423 
424  double eclEnergySumTCECLCalDigitInECLClusterThreshold(const Particle*)
425  {
426  StoreObjPtr<ECLTRGInformation> tce;
427  if (!tce) return std::numeric_limits<double>::quiet_NaN();
428  return tce->getClusterEnergyThreshold();
429  }
430 
431  // These variables require cDST inputs and the eclTrackCalDigitMatch module run first
432  VARIABLE_GROUP("ECL calibration");
433 
434  REGISTER_VARIABLE("eclEnergy3FWDBarrel", eclEnergy3FWDBarrel, R"DOC(
435 [Calibration] Returns energy sum of three crystals in forward barrel.
436 )DOC");
437 
438  REGISTER_VARIABLE("eclEnergy3FWDEndcap", eclEnergy3FWDEndcap, R"DOC(
439 [Calibration] Returns energy sum of three crystals in forward endcap.
440 )DOC");
441 
442  REGISTER_VARIABLE("eclEnergy3BWDBarrel", eclEnergy3BWDBarrel, R"DOC(
443 [Calibration] Returns energy sum of three crystals in backward barrel.
444 )DOC");
445 
446  REGISTER_VARIABLE("eclEnergy3BWDEndcap", eclEnergy3BWDEndcap, R"DOC(
447 [Calibration] Returns energy sum of three crystals in backward endcap.
448 )DOC");
449 
450  // These variables require cDST inputs and the eclTRGInformation module run first
451  VARIABLE_GROUP("ECL trigger calibration");
452 
453  REGISTER_VARIABLE("clusterNumberOfTCs(i, j, k, l)", eclNumberOfTCsForCluster, R"DOC(
454 [Calibration] Returns the number of TCs for this ECL cluster for a given TC theta ID range
455 :math:`(i, j)` and hit window :math:`(k, l)`.
456 )DOC");
457  REGISTER_VARIABLE("clusterTCFADC(i, j, k, l)", eclTCFADCForCluster, R"DOC(
458 [Calibration] Returns the total FADC sum related to this ECL cluster for a given TC theta ID
459 range :math:`(i, j)` and hit window :math:`(k, l)`.
460 )DOC");
461  REGISTER_VARIABLE("clusterTCIsMaximum", eclTCIsMaximumForCluster, R"DOC(
462 [Calibration] Returns True if cluster is related to maximum TC.
463 )DOC");
464  REGISTER_VARIABLE("clusterTrigger", eclClusterTrigger, R"DOC(
465 [Calibration] Returns 1.0 if ECL cluster is matched to a trigger cluster (requires to run eclTriggerClusterMatcher
466 (which requires TRGECLClusters in the input file)) and 0 otherwise. Returns -1 if the matching code was not run.
467 NOT FOR PHASE2 DATA!
468 )DOC");
469  REGISTER_VARIABLE("eclEnergyTC(i)", getEnergyTC, R"DOC(
470 [Eventbased][Calibration] Returns the energy (in FADC counts) for the :math:`i`-th trigger cell (TC), 1 based (1..576).
471 )DOC");
472  REGISTER_VARIABLE("eclEnergyTCECLCalDigit(i)", getEnergyTCECLCalDigit, R"DOC(
473 [Eventbased][Calibration] Returns the energy (in GeV) for the :math:`i`-th trigger cell (TC)
474 based on ECLCalDigits, 1 based (1..576).
475 )DOC");
476  REGISTER_VARIABLE("eclTimingTC(i)", getTimingTC, R"DOC(
477 [Eventbased][Calibration] Returns the time (in ns) for the :math:`i`-th trigger cell (TC), 1 based (1..576).
478 )DOC");
479  REGISTER_VARIABLE("eclHitWindowTC(i)", eclHitWindowTC, R"DOC(
480 [Eventbased][Calibration] Returns the hit window for the :math:`i`-th trigger cell (TC), 1 based (1..576).
481 )DOC");
482  REGISTER_VARIABLE("eclEventTimingTC", getEvtTimingTC, R"DOC(
483 [Eventbased][Calibration] Returns the ECL TC event time (in ns).
484 )DOC");
485  REGISTER_VARIABLE("eclMaximumTCId", getMaximumTCId, R"DOC(
486 [Eventbased][Calibration] Returns the TC ID with maximum FADC value.
487 )DOC");
488 
489  REGISTER_VARIABLE("eclTimingTCECLCalDigit(i)", getTimingTCECLCalDigit, R"DOC(
490 [Eventbased][Calibration] Returns the time (in ns) for the :math:`i`-th trigger cell (TC) based
491 on ECLCalDigits, 1 based (1..576)
492 )DOC");
493 
494  REGISTER_VARIABLE("eclNumberOfTCs(i, j, k)", getNumberOfTCs, R"DOC(
495 [Eventbased][Calibration] Returns the number of TCs above threshold (i=FADC counts) for this event
496 for a given theta range (j-k)
497 )DOC");
498  REGISTER_VARIABLE("eclEnergySumTC(i, j)", eclEnergySumTC, R"DOC(
499 [Eventbased][Calibration] Returns energy sum (in FADC counts) of all TC cells between two
500 theta ids i<=thetaid<=j, 1 based (1..17)
501 )DOC");
502  REGISTER_VARIABLE("eclEnergySumTCECLCalDigit(i, j, k, l)", eclEnergySumTCECLCalDigit, R"DOC(
503 [Eventbased][Calibration] Returns energy sum (in GeV) of all TC cells between two theta ids i<=thetaid<=j,
504 1 based (1..17). k is the sum option: 0 (all), 1 (those with actual TC entries), 2 (sum of ECLCalDigit energy
505 in this TC above threshold). l is the threshold parameter for the option k 2.
506 )DOC");
507  REGISTER_VARIABLE("eclEnergySumTCECLCalDigitInECLCluster", eclEnergySumTCECLCalDigitInECLCluster, R"DOC(
508 [Eventbased][Calibration] Returns energy sum (in GeV) of all ECLCalDigits if TC is above threshold
509 that are part of an ECLCluster above eclEnergySumTCECLCalDigitInECLClusterThreshold within TC thetaid 2-15.
510 )DOC");
511  REGISTER_VARIABLE("eclEnergySumECLCalDigitInECLCluster", eclEnergySumECLCalDigitInECLCluster, R"DOC(
512 [Eventbased][Calibration] Returns energy sum (in GeV) of all ECLCalDigits that are part of an ECL cluster
513 above eclEnergySumTCECLCalDigitInECLClusterThreshold within TC thetaid 2-15.
514 )DOC");
515  REGISTER_VARIABLE("eclEnergySumTCECLCalDigitInECLClusterThreshold", eclEnergySumTCECLCalDigitInECLClusterThreshold, R"DOC(
516 [Eventbased][Calibration] Returns threshold used to calculate eclEnergySumTCECLCalDigitInECLCluster.
517 )DOC");
518 
519  }
521 }
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17