Belle II Software development
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
23using namespace std;
24
25namespace 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
449range :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.
457NOT 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)
464based 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
481on 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
486for 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
490theta 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,
4941 based (1..17). k is the sum option: 0 (all), 1 (those with actual TC entries), 2 (sum of ECLCalDigit energy
495in 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
499that 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
503above 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
std
STL namespace.