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