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