Belle II Software  light-2212-foldex
ECLVariables.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 <analysis/variables/ECLVariables.h>
11 
12 //framework
13 #include <framework/logging/Logger.h>
14 
15 //analysis
16 #include <analysis/dataobjects/Particle.h>
17 #include <analysis/dataobjects/ParticleList.h>
18 #include <analysis/dataobjects/ECLEnergyCloseToTrack.h>
19 #include <analysis/utility/ReferenceFrame.h>
20 #include <analysis/ClusterUtility/ClusterUtils.h>
21 #include <analysis/VariableManager/Utility.h>
22 
23 //MDST
24 #include <mdst/dataobjects/KlId.h>
25 #include <mdst/dataobjects/ECLCluster.h>
26 #include <mdst/dataobjects/Track.h>
27 #include <mdst/dataobjects/EventLevelClusteringInfo.h>
28 
29 #include <Math/Vector4D.h>
30 
31 #include <cmath>
32 #include <stack>
33 
34 
35 
36 namespace Belle2 {
41  namespace Variable {
42  double beamBackgroundSuppression(const Particle* particle)
43  {
44  if (particle->hasExtraInfo("beamBackgroundSuppression")) {
45  return particle->getExtraInfo("beamBackgroundSuppression");
46  } else {
47  B2WARNING("The extraInfo beamBackgroundSuppression is not registered! \n"
48  "This variable is only available for photons, and you either have to run the function getBeamBackgroundProbability or turn the argument loadPhotonBeamBackgroundMVA to True when using fillParticleList.");
49  return std::numeric_limits<float>::quiet_NaN();
50  }
51  }
52 
53  double hadronicSplitOffSuppression(const Particle* particle)
54  {
55  if (particle->hasExtraInfo("hadronicSplitOffSuppression")) {
56  return particle->getExtraInfo("hadronicSplitOffSuppression");
57  } else {
58  B2WARNING("The extraInfo hadronicSplitOffSuppression is not registered! \n"
59  "This variable is only available for photons, and you either have to run the function getHadronicSplitOffProbability or turn the argument loadPhotonHadronicSplitOffMVA to True when using fillParticleList.");
60  return std::numeric_limits<float>::quiet_NaN();
61  }
62  }
63 
64  double eclClusterKlId(const Particle* particle)
65  {
66  const ECLCluster* cluster = particle->getECLCluster();
67  if (!cluster) {
68  return std::numeric_limits<double>::quiet_NaN();
69  }
70  const KlId* klid = cluster->getRelatedTo<KlId>();
71  if (!klid) {
72  return std::numeric_limits<double>::quiet_NaN();
73  }
74  return klid->getKlId();
75  }
76 
77 
78  double eclPulseShapeDiscriminationMVA(const Particle* particle)
79  {
80  const ECLCluster* cluster = particle->getECLCluster();
81  if (cluster) {
82  if (eclClusterHasPulseShapeDiscrimination(particle)) {
83  return cluster->getPulseShapeDiscriminationMVA();
84  } else {
85  return std::numeric_limits<float>::quiet_NaN();
86  }
87  }
88  return std::numeric_limits<float>::quiet_NaN();
89  }
90 
91  double eclClusterNumberOfHadronDigits(const Particle* particle)
92  {
93 
94  const ECLCluster* cluster = particle->getECLCluster();
95  if (cluster) {
96  if (eclClusterHasPulseShapeDiscrimination(particle)) {
97  return cluster->getNumberOfHadronDigits();
98  } else
99  return std::numeric_limits<float>::quiet_NaN();
100  }
101  return std::numeric_limits<float>::quiet_NaN();
102  }
103 
104  double eclClusterDetectionRegion(const Particle* particle)
105  {
106 
107  const ECLCluster* cluster = particle->getECLCluster();
108  if (cluster)
109  return cluster->getDetectorRegion();
110 
111  return std::numeric_limits<float>::quiet_NaN();
112  }
113 
114  double eclClusterIsolation(const Particle* particle)
115  {
116 
117  const ECLCluster* cluster = particle->getECLCluster();
118  if (cluster)
119  return cluster->getMinTrkDistance();
120 
121  return std::numeric_limits<float>::quiet_NaN();
122  }
123 
124  double eclClusterIsolationID(const Particle* particle)
125  {
126 
127  const ECLCluster* cluster = particle->getECLCluster();
128  if (cluster)
129  return cluster->getMinTrkDistanceID();
130 
131  return std::numeric_limits<float>::quiet_NaN();
132  }
133 
134  Manager::FunctionPtr eclClusterIsolationVar(const std::vector<std::string>& arguments)
135  {
136  if (arguments.size() > 2 or arguments.size() == 0)
137  B2FATAL("Wrong number of arguments (2 required) for meta variable minC2TDistVar");
138  std::string listName = "pi-:all";
139  std::string variableName = arguments[0];
140  if (arguments.size() == 2)
141  listName = arguments[1];
142 
143 
144  auto func = [listName, variableName](const Particle * particle) -> double {
145  StoreObjPtr<ParticleList> particleList(listName);
146  if (!(particleList.isValid()))
147  {
148  B2FATAL("Invalid Listname " << listName << " given to minC2TDistVar!");
149  }
150  const Variable::Manager::Var* var = Manager::Instance().getVariable(variableName);
151  const ECLCluster* cluster = particle->getECLCluster();
152  if (!cluster)
153  return std::numeric_limits<float>::quiet_NaN();
154  auto trackID = cluster->getMinTrkDistanceID();
155  double result = std::numeric_limits<float>::quiet_NaN();
156  // Find particle with that track ID:
157  for (unsigned int i = 0; i < particleList->getListSize(); i++)
158  {
159  const Particle* listParticle = particleList->getParticle(i);
160  if (listParticle and listParticle->getTrack() and listParticle->getTrack()->getArrayIndex() == trackID) {
161  result = std::get<double>(var->function(listParticle));
162  break;
163  }
164  }
165  return result;
166  };
167  return func;
168  }
169 
170  double eclClusterConnectedRegionID(const Particle* particle)
171  {
172 
173  const ECLCluster* cluster = particle->getECLCluster();
174  if (cluster)
175  return cluster->getConnectedRegionId();
176 
177  return std::numeric_limits<float>::quiet_NaN();
178  }
179 
180  double eclClusterDeltaL(const Particle* particle)
181  {
182 
183  const ECLCluster* cluster = particle->getECLCluster();
184  if (cluster)
185  return cluster->getDeltaL();
186 
187  return std::numeric_limits<float>::quiet_NaN();
188  }
189 
190  double eclClusterErrorE(const Particle* particle)
191  {
192 
193  const ECLCluster* cluster = particle->getECLCluster();
194  if (cluster) {
195  return cluster->getUncertaintyEnergy();
196  }
197  return std::numeric_limits<float>::quiet_NaN();
198  }
199 
200  double eclClusterUncorrectedE(const Particle* particle)
201  {
202 
203  const ECLCluster* cluster = particle->getECLCluster();
204  if (cluster) {
205  return cluster->getEnergyRaw();
206  }
207  return std::numeric_limits<float>::quiet_NaN();
208  }
209 
210  double eclClusterE(const Particle* particle)
211  {
212  const auto& frame = ReferenceFrame::GetCurrent();
213 
214  const ECLCluster* cluster = particle->getECLCluster();
215  if (cluster) {
216  ClusterUtils clutls;
217  ROOT::Math::PxPyPzEVector p4Cluster = clutls.GetCluster4MomentumFromCluster(cluster, particle->getECLClusterEHypothesisBit());
218 
219  return frame.getMomentum(p4Cluster).E();
220  }
221  return std::numeric_limits<float>::quiet_NaN();
222  }
223 
224  double eclClusterHighestE(const Particle* particle)
225  {
226 
227  const ECLCluster* cluster = particle->getECLCluster();
228  if (cluster) {
229  return cluster->getEnergyHighestCrystal();
230  }
231  return std::numeric_limits<float>::quiet_NaN();
232  }
233 
234  double eclClusterCellId(const Particle* particle)
235  {
236 
237  const ECLCluster* cluster = particle->getECLCluster();
238  if (cluster) {
239  return cluster->getMaxECellId();
240  }
241  return std::numeric_limits<float>::quiet_NaN();
242  }
243 
244  // An array with each number representing the last number of the cellID per thetaID. There are 69 thetaIDs in total.
245  const std::array<int, 69> lastCellIDperThetaID{48, 96, 160, 224, 288, 384, 480, 576, 672, 768, 864,
246  1008, 1152, 1296, 1440, 1584, 1728, 1872, 2016, 2160, 2304, 2448,
247  2592, 2736, 2880, 3024, 3168, 3312, 3456, 3600, 3744, 3888, 4032,
248  4176, 4320, 4464, 4608, 4752, 4896, 5040, 5184, 5328, 5472, 5616,
249  5760, 5904, 6048, 6192, 6336, 6480, 6624, 6768, 6912, 7056, 7200,
250  7344, 7488, 7632, 7776, 7920, 8064, 8160, 8256, 8352, 8448, 8544,
251  8608, 8672, 8736};
252 
253  double eclClusterThetaId(const Particle* particle)
254  {
255 
256  const ECLCluster* cluster = particle->getECLCluster();
257  if (cluster) {
258  int cellID = cluster->getMaxECellId();
259  return std::distance(lastCellIDperThetaID.begin(), std::lower_bound(lastCellIDperThetaID.begin(), lastCellIDperThetaID.end(),
260  cellID));
261  }
262  return std::numeric_limits<float>::quiet_NaN();
263  }
264 
265  double eclClusterPhiId(const Particle* particle)
266  {
267 
268  const ECLCluster* cluster = particle->getECLCluster();
269  if (cluster) {
270  int cellID = cluster->getMaxECellId();
271  if (cellID <= 48) {
272  return cellID - 1;
273  } else {
274  int closestinlist = lastCellIDperThetaID[std::distance(lastCellIDperThetaID.begin(), std::lower_bound(lastCellIDperThetaID.begin(),
275  lastCellIDperThetaID.end(), cellID)) - 1];
276  return cellID - closestinlist - 1;
277  }
278  }
279  return std::numeric_limits<float>::quiet_NaN();
280  }
281 
282  double eclClusterTiming(const Particle* particle)
283  {
284 
285  const ECLCluster* cluster = particle->getECLCluster();
286  if (cluster) {
287  return cluster->getTime();
288  }
289  return std::numeric_limits<float>::quiet_NaN();
290  }
291 
292  double eclClusterHasFailedTiming(const Particle* particle)
293  {
294  const ECLCluster* cluster = particle->getECLCluster();
295  if (cluster) {
296  return cluster->hasFailedFitTime();
297  }
298  return std::numeric_limits<float>::quiet_NaN();
299  }
300 
301  double eclClusterErrorTiming(const Particle* particle)
302  {
303 
304  const ECLCluster* cluster = particle->getECLCluster();
305  if (cluster) {
306  return cluster->getDeltaTime99();
307  }
308  return std::numeric_limits<float>::quiet_NaN();
309  }
310 
311  double eclClusterHasFailedErrorTiming(const Particle* particle)
312  {
313  const ECLCluster* cluster = particle->getECLCluster();
314  if (cluster) {
315  return cluster->hasFailedTimeResolution();
316  }
317  return std::numeric_limits<float>::quiet_NaN();
318  }
319 
320  double eclClusterTheta(const Particle* particle)
321  {
322  const auto& frame = ReferenceFrame::GetCurrent();
323 
324  const ECLCluster* cluster = particle->getECLCluster();
325  if (cluster) {
326  ClusterUtils clutls;
327  ROOT::Math::PxPyPzEVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, particle->getECLClusterEHypothesisBit());
328 
329  return frame.getMomentum(p4Cluster).Theta();
330  }
331  return std::numeric_limits<float>::quiet_NaN();
332  }
333 
334  double eclClusterErrorTheta(const Particle* particle)
335  {
336 
337  const ECLCluster* cluster = particle->getECLCluster();
338  if (cluster) {
339  return cluster->getUncertaintyTheta();
340  }
341  return std::numeric_limits<float>::quiet_NaN();
342  }
343 
344  double eclClusterErrorPhi(const Particle* particle)
345  {
346 
347  const ECLCluster* cluster = particle->getECLCluster();
348  if (cluster) {
349  return cluster->getUncertaintyPhi();
350  }
351  return std::numeric_limits<float>::quiet_NaN();
352  }
353 
354  double eclClusterPhi(const Particle* particle)
355  {
356  const auto& frame = ReferenceFrame::GetCurrent();
357 
358  const ECLCluster* cluster = particle->getECLCluster();
359  if (cluster) {
360  ClusterUtils clutls;
361  ROOT::Math::PxPyPzEVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, particle->getECLClusterEHypothesisBit());
362 
363  return frame.getMomentum(p4Cluster).Phi();
364  }
365  return std::numeric_limits<float>::quiet_NaN();
366  }
367 
368  double eclClusterR(const Particle* particle)
369  {
370 
371  const ECLCluster* cluster = particle->getECLCluster();
372  if (cluster) {
373  return cluster->getR();
374  }
375  return std::numeric_limits<float>::quiet_NaN();
376  }
377 
378  double eclClusterE1E9(const Particle* particle)
379  {
380 
381  const ECLCluster* cluster = particle->getECLCluster();
382  if (cluster) {
383  return cluster->getE1oE9();
384  }
385  return std::numeric_limits<float>::quiet_NaN();
386  }
387 
388  double eclClusterE9E21(const Particle* particle)
389  {
390 
391  const ECLCluster* cluster = particle->getECLCluster();
392  if (cluster) {
393  return cluster->getE9oE21();
394  }
395  return std::numeric_limits<float>::quiet_NaN();
396  }
397 
398  double eclClusterAbsZernikeMoment40(const Particle* particle)
399  {
400 
401  const ECLCluster* cluster = particle->getECLCluster();
402  if (cluster) {
403  return cluster->getAbsZernike40();
404  }
405  return std::numeric_limits<float>::quiet_NaN();
406  }
407 
408  double eclClusterAbsZernikeMoment51(const Particle* particle)
409  {
410 
411  const ECLCluster* cluster = particle->getECLCluster();
412  if (cluster) {
413  return cluster->getAbsZernike51();
414  }
415  return std::numeric_limits<float>::quiet_NaN();
416  }
417 
418  double eclClusterZernikeMVA(const Particle* particle)
419  {
420 
421  const ECLCluster* cluster = particle->getECLCluster();
422  if (cluster) {
423  return cluster->getZernikeMVA();
424  }
425  return std::numeric_limits<float>::quiet_NaN();
426  }
427 
428  double eclClusterSecondMoment(const Particle* particle)
429  {
430 
431  const ECLCluster* cluster = particle->getECLCluster();
432  if (cluster) {
433  return cluster->getSecondMoment();
434  }
435  return std::numeric_limits<float>::quiet_NaN();
436  }
437 
438  double eclClusterLAT(const Particle* particle)
439  {
440 
441  const ECLCluster* cluster = particle->getECLCluster();
442  if (cluster) {
443  return cluster->getLAT();
444  }
445  return std::numeric_limits<float>::quiet_NaN();
446  }
447 
448  double eclClusterNHits(const Particle* particle)
449  {
450 
451  const ECLCluster* cluster = particle->getECLCluster();
452  if (cluster) {
453  return cluster->getNumberOfCrystals();
454  }
455  return std::numeric_limits<float>::quiet_NaN();
456  }
457 
458  double eclClusterTrackMatched(const Particle* particle)
459  {
460 
461  const ECLCluster* cluster = particle->getECLCluster();
462  if (cluster) {
463  const Track* track = cluster->getRelatedFrom<Track>();
464 
465  if (track)
466  return 1.0;
467  else
468  return 0.0;
469  }
470  return std::numeric_limits<float>::quiet_NaN();
471  }
472 
473  double nECLClusterTrackMatches(const Particle* particle)
474  {
475  // if no ECL cluster then nan
476  const ECLCluster* cluster = particle->getECLCluster();
477  if (!cluster)
478  return std::numeric_limits<double>::quiet_NaN();
479 
480  // one or more tracks may be matched to charged particles
481  size_t out = cluster->getRelationsFrom<Track>().size();
482  return double(out);
483  }
484 
485  double eclClusterConnectedRegionId(const Particle* particle)
486  {
487  const ECLCluster* cluster = particle->getECLCluster();
488  if (cluster) {
489  return cluster->getConnectedRegionId();
490  }
491  return std::numeric_limits<float>::quiet_NaN();
492  }
493 
494  double eclClusterId(const Particle* particle)
495  {
496  const ECLCluster* cluster = particle->getECLCluster();
497  if (cluster) {
498  return cluster->getClusterId();
499  }
500  return std::numeric_limits<float>::quiet_NaN();
501  }
502 
503  double eclClusterHasNPhotonsHypothesis(const Particle* particle)
504  {
505  const ECLCluster* cluster = particle->getECLCluster();
506  if (cluster) {
507  return cluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
508  }
509  return std::numeric_limits<float>::quiet_NaN();
510  }
511 
512  double eclClusterHasNeutralHadronHypothesis(const Particle* particle)
513  {
514  const ECLCluster* cluster = particle->getECLCluster();
515  if (cluster) {
516  return cluster->hasHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron);
517  }
518  return std::numeric_limits<float>::quiet_NaN();
519  }
520 
521  double eclClusterHasPulseShapeDiscrimination(const Particle* particle)
522  {
523  const ECLCluster* cluster = particle->getECLCluster();
524  if (cluster) {
525  return cluster->hasPulseShapeDiscrimination();
526  }
527  return std::numeric_limits<float>::quiet_NaN();
528  }
529 
530  double eclExtTheta(const Particle* particle)
531  {
532 
533  const Track* track = particle->getTrack();
534  if (track) {
535 
536  auto* eclinfo = track->getRelatedTo<ECLEnergyCloseToTrack>();
537 
538  if (eclinfo) {
539  return eclinfo->getExtTheta();
540  } else {
541  B2WARNING("Relation to ECLEnergyCloseToTrack not found, did you forget to run ECLTrackCalDigitMatchModule?");
542  return std::numeric_limits<float>::quiet_NaN();
543  }
544  }
545 
546  return std::numeric_limits<float>::quiet_NaN();
547  }
548 
549  double eclExtPhi(const Particle* particle)
550  {
551 
552  const Track* track = particle->getTrack();
553  if (track) {
554 
555  auto* eclinfo = track->getRelatedTo<ECLEnergyCloseToTrack>();
556 
557  if (eclinfo) {
558  return eclinfo->getExtPhi();
559  } else {
560  B2WARNING("Relation to ECLEnergyCloseToTrack not found, did you forget to run ECLTrackCalDigitMatchModule?");
561  return std::numeric_limits<float>::quiet_NaN();
562  }
563  }
564 
565  return std::numeric_limits<float>::quiet_NaN();
566  }
567 
568  double eclExtPhiId(const Particle* particle)
569  {
570  const Track* track = particle->getTrack();
571  if (track) {
572 
573  auto* eclinfo = track->getRelatedTo<ECLEnergyCloseToTrack>();
574 
575  if (eclinfo) {
576  return eclinfo->getExtPhiId();
577  } else {
578  B2WARNING("Relation to ECLEnergyCloseToTrack not found, did you forget to run ECLTrackCalDigitMatchModule?");
579  return std::numeric_limits<float>::quiet_NaN();
580  }
581  }
582 
583  return std::numeric_limits<float>::quiet_NaN();
584  }
585 
586  double weightedAverageECLTime(const Particle* particle)
587  {
588  int nDaughters = particle->getNDaughters();
589  if (nDaughters < 1) {
590  B2WARNING("The provided particle has no daughters!");
591  return std::numeric_limits<float>::quiet_NaN();
592  }
593 
594  double numer = 0, denom = 0;
595  int numberOfClusterDaughters = 0;
596 
597  auto weightedECLTimeAverage = [&numer, &denom, &numberOfClusterDaughters](const Particle * p) {
598  const ECLCluster* cluster = p->getECLCluster();
599  if (cluster and not cluster->hasFailedFitTime()) {
600  numberOfClusterDaughters ++;
601 
602  double time = cluster->getTime();
603  B2DEBUG(10, "time[" << numberOfClusterDaughters << "] = " << time);
604  double deltatime = cluster->getDeltaTime99();
605  B2DEBUG(10, "deltatime[" << numberOfClusterDaughters << "] = " << deltatime);
606  numer += time / pow(deltatime, 2);
607  B2DEBUG(11, "numer[" << numberOfClusterDaughters << "] = " << numer);
608  denom += 1 / pow(deltatime, 2);
609  B2DEBUG(11, "denom[" << numberOfClusterDaughters << "] = " << denom);
610  }
611  return false;
612  };
613 
614  particle->forEachDaughter(weightedECLTimeAverage, true, true);
615 
616  if (numberOfClusterDaughters < 1) {
617  B2WARNING("There are no clusters or cluster matches amongst the daughters of the provided particle!");
618  return std::numeric_limits<float>::quiet_NaN();
619  }
620 
621  if (denom == 0) {
622  B2WARNING("The denominator of the weighted mean is zero!");
623  return std::numeric_limits<float>::quiet_NaN();
624  } else {
625  B2DEBUG(10, "numer/denom = " << numer / denom);
626  return numer / denom;
627  }
628  }
629 
630  double maxWeightedDistanceFromAverageECLTime(const Particle* particle)
631  {
632  int nDaughters = particle->getNDaughters();
633  if (nDaughters < 1) {
634  B2WARNING("The provided particle has no daughters!");
635  return std::numeric_limits<float>::quiet_NaN();
636  }
637 
638  double maxTimeDiff = -DBL_MAX;
639  int numberOfClusterDaughters = 0;
640 
641  double averageECLTime = weightedAverageECLTime(particle);
642 
643  auto maxTimeDifference = [&maxTimeDiff, &numberOfClusterDaughters, &averageECLTime](const Particle * p) {
644 
645  const ECLCluster* cluster = p->getECLCluster();
646  if (cluster) {
647  numberOfClusterDaughters ++;
648 
649  double time = cluster->getTime();
650  B2DEBUG(10, "time[" << numberOfClusterDaughters << "] = " << time);
651  double deltatime = cluster->getDeltaTime99();
652  B2DEBUG(10, "deltatime[" << numberOfClusterDaughters << "] = " << deltatime);
653  double maxTimeDiff_temp = fabs((time - averageECLTime) / deltatime);
654  B2DEBUG(11, "maxTimeDiff_temp[" << numberOfClusterDaughters << "] = " << maxTimeDiff_temp);
655  if (maxTimeDiff_temp > maxTimeDiff)
656  maxTimeDiff = maxTimeDiff_temp;
657  B2DEBUG(11, "maxTimeDiff[" << numberOfClusterDaughters << "] = " << maxTimeDiff);
658  }
659  return false;
660  };
661 
662  particle->forEachDaughter(maxTimeDifference, true, true);
663 
664  if (numberOfClusterDaughters < 1) {
665  B2WARNING("There are no clusters or cluster matches amongst the daughters of the provided particle!");
666  return std::numeric_limits<float>::quiet_NaN();
667  }
668 
669  if (maxTimeDiff < 0) {
670  B2WARNING("The max time difference is negative!");
671  return std::numeric_limits<float>::quiet_NaN();
672  } else {
673  B2DEBUG(10, "maxTimeDiff = " << maxTimeDiff);
674  return maxTimeDiff;
675  }
676  }
677 
678  double eclClusterMdstIndex(const Particle* particle)
679  {
680  const ECLCluster* cluster = particle->getECLCluster();
681  if (cluster) {
682  return cluster->getArrayIndex();
683  } else return std::numeric_limits<double>::quiet_NaN();
684 
685  return std::numeric_limits<double>::quiet_NaN();
686  }
687 
688 
689  /*************************************************************
690  * Event-based ECL clustering information
691  */
692  double nECLOutOfTimeCrystalsFWDEndcap(const Particle*)
693  {
694  StoreObjPtr<EventLevelClusteringInfo> elci;
695  if (!elci) return std::numeric_limits<double>::quiet_NaN();
696  return (double)elci->getNECLCalDigitsOutOfTimeFWD();
697  }
698 
699  double nECLOutOfTimeCrystalsBarrel(const Particle*)
700  {
701  StoreObjPtr<EventLevelClusteringInfo> elci;
702  if (!elci) return std::numeric_limits<double>::quiet_NaN();
703  return (double)elci->getNECLCalDigitsOutOfTimeBarrel();
704  }
705 
706  double nECLOutOfTimeCrystalsBWDEndcap(const Particle*)
707  {
708  StoreObjPtr<EventLevelClusteringInfo> elci;
709  if (!elci) return std::numeric_limits<double>::quiet_NaN();
710  return (double)elci->getNECLCalDigitsOutOfTimeBWD();
711  }
712 
713  double nECLOutOfTimeCrystals(const Particle*)
714  {
715  StoreObjPtr<EventLevelClusteringInfo> elci;
716  if (!elci) return std::numeric_limits<double>::quiet_NaN();
717  return (double)elci->getNECLCalDigitsOutOfTime();
718  }
719 
720  double nRejectedECLShowersFWDEndcap(const Particle*)
721  {
722  StoreObjPtr<EventLevelClusteringInfo> elci;
723  if (!elci) return std::numeric_limits<double>::quiet_NaN();
724  return (double)elci->getNECLShowersRejectedFWD();
725  }
726 
727  double nRejectedECLShowersBarrel(const Particle*)
728  {
729  StoreObjPtr<EventLevelClusteringInfo> elci;
730  if (!elci) return std::numeric_limits<double>::quiet_NaN();
731  return (double)elci->getNECLShowersRejectedBarrel();
732  }
733 
734  double nRejectedECLShowersBWDEndcap(const Particle*)
735  {
736  StoreObjPtr<EventLevelClusteringInfo> elci;
737  if (!elci) return std::numeric_limits<double>::quiet_NaN();
738  return (double)elci->getNECLShowersRejectedBWD();
739  }
740 
741  double nRejectedECLShowers(const Particle*)
742  {
743  StoreObjPtr<EventLevelClusteringInfo> elci;
744  if (!elci) return std::numeric_limits<double>::quiet_NaN();
745  return (double) elci->getNECLShowersRejected();
746  }
747 
748  double eclClusterEoP(const Particle* part)
749  {
750  double E = eclClusterE(part);
751  if (part->hasExtraInfo("bremsCorrectedPhotonEnergy")) {
752  E += part->getExtraInfo("bremsCorrectedPhotonEnergy");
753  }
754  const double p = part->getMomentumMagnitude();
755  if (0 == p) { return std::numeric_limits<float>::quiet_NaN();}
756  return E / p;
757  }
758 
759  double eclClusterOnlyInvariantMass(const Particle* part)
760  {
761  int nDaughters = part->getNDaughters();
762  ROOT::Math::PxPyPzEVector sum;
763 
764  if (nDaughters < 1) {
765  return part->getMass();
766  } else {
767  int nClusterDaughters = 0;
768  std::stack<const Particle*> stacked;
769  stacked.push(part);
770  while (!stacked.empty()) {
771  const Particle* current = stacked.top();
772  stacked.pop();
773 
774  const ECLCluster* cluster = current->getECLCluster();
775  if (cluster) {
776  const ECLCluster::EHypothesisBit clusterBit = current->getECLClusterEHypothesisBit();
777  nClusterDaughters ++;
778  ClusterUtils clutls;
779  ROOT::Math::PxPyPzEVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, clusterBit);
780  sum += p4Cluster;
781  } else {
782  const std::vector<Particle*> daughters = current->getDaughters();
783  nDaughters = current->getNDaughters();
784  for (int iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
785  stacked.push(daughters[iDaughter]);
786  }
787  }
788  }
789 
790  if (nClusterDaughters < 1) {
791  B2WARNING("There are no clusters amongst the daughters of the provided particle!");
792  return std::numeric_limits<double>::quiet_NaN();
793  }
794  B2DEBUG(10, "Number of daughters with cluster associated = " << nClusterDaughters);
795  return sum.M();
796  }
797  }
798 
799  Manager::FunctionPtr photonHasOverlap(const std::vector<std::string>& arguments)
800  {
801  std::string cutString = "";
802  if (arguments.size() > 0) {
803  cutString = arguments[0];
804  }
805  std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
806 
807  std::string photonlistname = "gamma:all";
808  if (arguments.size() > 1) {
809  photonlistname = arguments[1];
810  }
811 
812  std::string tracklistname = "e-:all";
813  if (arguments.size() > 2) {
814  tracklistname = arguments[2];
815  }
816 
817  auto func = [cut, photonlistname, tracklistname](const Particle * particle) -> double {
818 
819  if (particle->getPDGCode() != Const::photon.getPDGCode())
820  {
821  B2WARNING("The variable photonHasOverlap is supposed to be calculated for photons. Returning NaN.");
822  return std::numeric_limits<double>::quiet_NaN();
823  }
824 
825  StoreObjPtr<ParticleList> photonlist(photonlistname);
826  if (!(photonlist.isValid()))
827  {
828  B2WARNING("The provided particle list " << photonlistname << " does not exist."
829  " Therefore, the variable photonHasOverlap can not be calculated. Returning NaN.");
830  return std::numeric_limits<double>::quiet_NaN();
831  }
832  if (photonlist->getPDGCode() != Const::photon.getPDGCode())
833  {
834  B2WARNING("The list " << photonlistname << " does not contain photons."
835  " Therefore, the variable photonHasOverlap can not be calculated reliably. Returning NaN.");
836  return std::numeric_limits<double>::quiet_NaN();
837  }
838 
839  StoreObjPtr<ParticleList> tracklist(tracklistname);
840  if (!(tracklist.isValid()))
841  {
842  B2WARNING("The provided particle list " << tracklistname << " does not exist."
843  " Therefore, the variable photonHasOverlap can not be calculated. Returning NaN.");
844  return std::numeric_limits<double>::quiet_NaN();
845  }
846  if (!Const::chargedStableSet.contains(Const::ParticleType(abs(tracklist->getPDGCode()))))
847  {
848  B2WARNING("The list " << tracklistname << " does not contain charged final state particles."
849  " Therefore, the variable photonHasOverlap can not be calculated reliably. Returning NaN.");
850  return std::numeric_limits<double>::quiet_NaN();
851  }
852 
853  double connectedRegionID = eclClusterConnectedRegionID(particle);
854  unsigned mdstArrayIndex = particle->getMdstArrayIndex();
855 
856  for (unsigned int i = 0; i < photonlist->getListSize(); i++)
857  {
858  const Particle* part = photonlist->getParticle(i);
859 
860  // skip the particle itself
861  if (part->getMdstArrayIndex() == mdstArrayIndex) {
862  continue;
863  }
864 
865  // skip photons that do not fulfill the provided criteria
866  if (!cut->check(part)) {
867  continue;
868  }
869 
870  if (connectedRegionID == eclClusterConnectedRegionID(part)) {
871  return 1;
872  }
873  }
874 
875  for (unsigned int i = 0; i < tracklist->getListSize(); i++)
876  {
877  const Particle* part = tracklist->getParticle(i);
878 
879  // skip tracks that do not fulfill the provided criteria
880  if (!cut->check(part)) {
881  continue;
882  }
883 
884  if (connectedRegionID == eclClusterConnectedRegionID(part)) {
885  return 1;
886  }
887  }
888 
889  return 0;
890  };
891  return func;
892  }
893 
894  VARIABLE_GROUP("ECL Cluster related");
895  REGISTER_VARIABLE("clusterEoP", eclClusterEoP, R"DOC(
896 Returns ratio of uncorrelated energy E over momentum p, a convenience
897 alias for (clusterE / p).
898 )DOC");
899  REGISTER_VARIABLE("clusterReg", eclClusterDetectionRegion, R"DOC(
900 Returns an integer code for the ECL region of a cluster.
901 
902  - 1: forward, 2: barrel, 3: backward,
903  - 11: between FWD and barrel, 13: between BWD and barrel,
904  - 0: otherwise
905 )DOC");
906  REGISTER_VARIABLE("clusterDeltaLTemp", eclClusterDeltaL, R"DOC(
907 | Returns DeltaL for the shower shape.
908 | A cluster comprises the energy depositions of several crystals. All these crystals have slightly
909  different orientations in space. A shower direction can be constructed by calculating the weighted
910  average of these orientations using the corresponding energy depositions as weights. The intersection
911  (more precisely the point of closest approach) of the vector with this direction originating from the
912  cluster center and an extrapolated track can be used as reference for the calculation of the shower
913  depth. It is defined as the distance between this intersection and the cluster center.
914 
915 .. warning::
916  This distance is calculated on the reconstructed level and is temporarily
917  included to the ECL cluster MDST data format for studying purposes. If it is found
918  that it is not crucial for physics analysis then this variable will be removed
919  in future releases.
920  Therefore, keep in mind that this variable might be removed in the future!
921 
922 .. note::
923  | Please read `this <importantNoteECL>` first.
924  | Lower limit: :math:`-250.0`
925  | Upper limit: :math:`250.0`
926  | Precision: :math:`10` bit
927 ..
928 )DOC","cm");
929 
930 
931  REGISTER_VARIABLE("minC2TDist", eclClusterIsolation, R"DOC(
932 Returns distance between ECL cluster and nearest track hitting the ECL.
933 
934 A cluster comprises the energy depositions of several crystals. All these crystals have slightly
935 different orientations in space. A shower direction can be constructed by calculating the weighted
936 average of these orientations using the corresponding energy depositions as weights. The intersection
937 (more precisely the point of closest approach) of the vector with this direction originating from the
938 cluster center and an extrapolated track can be used as reference for the calculation of the track depth.
939 It is defined as the distance between this intersection and the track hit position on the front face of the ECL.
940 
941 .. note::
942  This distance is calculated on the reconstructed level.
943 
944 .. note::
945  | Please read `this <importantNoteECL>` first.
946  | Lower limit: :math:`0.0`
947  | Upper limit: :math:`250.0`
948  | Precision: :math:`10` bit
949 ..
950 )DOC","cm");
951  REGISTER_VARIABLE("minC2TDistID", eclClusterIsolationID, "Nearest track array index");
952  REGISTER_METAVARIABLE("minC2TDistVar(variable,particleList=pi-:all)", eclClusterIsolationVar, R"DOC(
953 Returns variable value for the nearest track to the given ECL cluster. First argument is a variable name, e.g. nCDCHits.
954 The second argument is the particle list name which will be used to pick up the nearest track, default is pi-:all.
955 )DOC", Manager::VariableDataType::c_double);
956  REGISTER_VARIABLE("clusterE", eclClusterE, R"DOC(
957 Returns ECL cluster's energy corrected for leakage and background.
958 
959 The raw photon energy is given by the weighted sum of all ECL crystal energies within the ECL cluster.
960 The weights per crystals are :math:`\leq 1` after cluster energy splitting in the case of overlapping
961 clusters. The number of crystals that are included in the sum depends on a initial energy estimation
962 and local beam background levels at the highest energy crystal position. It is optimized to minimize
963 the core width (resolution) of true photons. Photon energy distributions always show a low energy tail
964 due to unavoidable longitudinal and transverse leakage that can be further modified by the clustering
965 algorithm and beam backgrounds.The peak position of the photon energy distributions are corrected to
966 match the true photon energy in MC:
967 
968  - Leakage correction: Using large MC samples of mono-energetic single photons, a correction factor
969  :math:`f` as function of reconstructed detector position, reconstructed photon energy and beam backgrounds
970  is determined via :math:`f = \frac{\text{peak_reconstructed}}{\text{energy_true}}`.
971 
972  - Cluster energy calibration (data only): To reach the target precision of :math:`< 1.8\%` energy
973  resolution for high energetic photons, the remaining difference between MC and data must be calibrated
974  using kinematically fit muon pairs. This calibration is only applied to data and not to MC and will
975  take time to develop.
976 
977  - Energy Bias Correction module, sub-percent correction, is NOT applied on clusterE, but on photon energy
978  and momentum. Only applied to data.
979 
980 It is important to note that after perfect leakage correction and cluster energy calibration,
981 the :math:`\pi^{0}` mass peak will be shifted slightly to smaller values than the PDG average
982 due to the low energy tails of photons. The :math:`\pi^{0}` mass peak must not be corrected
983 to the PDG value by adjusting the reconstructed photon energies. Selection criteria based on
984 the mass for :math:`\pi^{0}` candidates must be based on the biased value. Most analysis
985 will used mass constrained :math:`\pi^{0}` s anyhow.
986 
987 .. warning::
988  We only store clusters with :math:`E > 20\,` MeV.
989 
990 .. note::
991  | Please read `this <importantNoteECL>` first.
992  | Lower limit: :math:`-5` (:math:`e^{-5} = 0.00674\,` GeV)
993  | Upper limit: :math:`3.0` (:math:`e^3 = 20.08553\,` GeV)
994  | Precision: :math:`18` bit
995  | This value can be changed to a different reference frame with :b2:var:`useCMSFrame`.
996 ..
997 )DOC","GeV");
998  REGISTER_VARIABLE("clusterErrorE", eclClusterErrorE, R"DOC(
999 Returns ECL cluster's uncertainty on energy
1000 (from background level and energy dependent tabulation).
1001 )DOC","GeV");
1002  REGISTER_VARIABLE("clusterErrorPhi", eclClusterErrorPhi, R"DOC(
1003 Returns ECL cluster's uncertainty on :math:`\phi`
1004 (from background level and energy dependent tabulation).
1005 )DOC","rad");
1006  REGISTER_VARIABLE("clusterErrorTheta", eclClusterErrorTheta, R"DOC(
1007 Returns ECL cluster's uncertainty on :math:`\theta`
1008 (from background level and energy dependent tabulation).
1009 )DOC","rad");
1010 
1011  REGISTER_VARIABLE("clusterR", eclClusterR, R"DOC(
1012 Returns ECL cluster's centroid distance from :math:`(0,0,0)`.
1013 )DOC","cm");
1014  REGISTER_VARIABLE("clusterPhi", eclClusterPhi, R"DOC(
1015 Returns ECL cluster's azimuthal angle :math:`\phi`
1016 (this is not generally equal to a photon azimuthal angle).
1017 
1018 | The direction of a cluster is given by the connecting line of :math:`\,(0,0,0)\,` and
1019  cluster centroid position in the ECL.
1020 | The cluster centroid position is calculated using up to 21 crystals (5x5 excluding corners)
1021  after cluster energy splitting in the case of overlapping clusters.
1022 | The centroid position is the logarithmically weighted average of all crystals evaluated at
1023  the crystal centers. Cluster centroids are generally biased towards the centers of the
1024  highest energetic crystal. This effect is larger for low energetic photons.
1025 | Beam backgrounds slightly decrease the position resolution, mainly for low energetic photons.
1026 
1027 .. note::
1028  Radius of a cluster is almost constant in the barrel and should not be used directly in any selection.
1029 
1030 Unlike for charged tracks, the uncertainty (covariance) of the photon directions is not determined
1031 based on individual cluster properties but taken from on MC-based parametrizations of the resolution
1032 as function of true photon energy, true photon direction and beam background level.
1033 
1034 .. warning::
1035  Users must use the actual particle direction (done automatically in the modularAnalysis using the average
1036  IP position (can be changed if needed)) and not the ECL Cluster direction (position in the ECL measured
1037  from :math:`(0,0,0)`) for particle kinematics.
1038 
1039 .. note::
1040  | Please read `this <importantNoteECL>` first.
1041  | Lower limit: :math:`-\pi`
1042  | Upper limit: :math:`\pi`
1043  | Precision: :math:`16` bit
1044 ..
1045 )DOC","rad");
1046  REGISTER_VARIABLE("clusterConnectedRegionID", eclClusterConnectedRegionID, R"DOC(
1047 Returns ECL cluster's connected region ID.
1048 )DOC");
1049  REGISTER_VARIABLE("clusterTheta", eclClusterTheta, R"DOC(
1050 Returns ECL cluster's polar angle :math:`\theta`
1051 (this is not generally equal to a photon polar angle).
1052 
1053 | The direction of a cluster is given by the connecting line of :math:`\,(0,0,0)\,` and
1054  cluster centroid position in the ECL.
1055 | The cluster centroid position is calculated using up to 21 crystals (5x5 excluding corners)
1056  after cluster energy splitting in the case of overlapping clusters.
1057 | The centroid position is the logarithmically weighted average of all crystals evaluated at
1058  the crystal centers. Cluster centroids are generally biased towards the centers of the
1059  highest energetic crystal. This effect is larger for low energetic photons.
1060 | Beam backgrounds slightly decrease the position resolution, mainly for low energetic photons.
1061 
1062 .. note::
1063  Radius of a cluster is almost constant in the barrel and should not be used directly in any selection.
1064 
1065 Unlike for charged tracks, the uncertainty (covariance) of the photon directions is not determined
1066 based on individual cluster properties but taken from on MC-based parametrizations of the resolution
1067 as function of true photon energy, true photon direction and beam background level.
1068 
1069 .. warning::
1070  Users must use the actual particle direction (done automatically in the modularAnalysis using the average
1071  IP position (can be changed if needed)) and not the ECL Cluster direction (position in the ECL measured
1072  from :math:`(0,0,0)`) for particle kinematics.
1073 
1074 .. note::
1075  | Please read `this <importantNoteECL>` first.
1076  | Lower limit: :math:`0.0`
1077  | Upper limit: :math:`\pi`
1078  | Precision: :math:`16` bit
1079 ..
1080 )DOC","rad");
1081  REGISTER_VARIABLE("clusterTiming", eclClusterTiming, R"DOC(
1082 Returns the time of the ECL cluster. It is calculated as the Photon timing minus the Event t0.
1083 Photon timing is given by the fitted time of the recorded waveform of the highest energy crystal in the
1084 cluster. After all calibrations and corrections (including Time-Of-Flight), photons from the interaction
1085 point (IP) should have a Photon timing that corresponds to the Event t0, :math:`t_{0}`. The Event t0 is the
1086 time of the event and may be measured by a different sub-detector (see Event t0 documentation). For an ECL
1087 cluster produced at the interaction point in time with the event, the cluster time should be consistent with zero
1088 within the uncertainties. Special values are returned if the fit for the Photon timing fails (see
1089 documentation for `clusterHasFailedTiming`). (For MC, the calibrations and corrections are not fully simulated).
1090 
1091 .. note::
1092  | Please read `this <importantNoteECL>` first.
1093  | Lower limit: :math:`-1000.0`
1094  | Upper limit: :math:`1000.0`
1095  | Precision: :math:`12` bit
1096 ..
1097 )DOC","ns");
1098  REGISTER_VARIABLE("clusterHasFailedTiming", eclClusterHasFailedTiming, R"DOC(
1099 Status bit for if the ECL cluster's timing fit failed. Photon timing is given by the fitted time
1100 of the recorded waveform of the highest energetic crystal in a cluster; however, that fit can fail and so
1101 this variable tells the user if that has happened.
1102 )DOC");
1103  REGISTER_VARIABLE("clusterErrorTiming", eclClusterErrorTiming, R"DOC(
1104 Returns ECL cluster's timing uncertainty that contains :math:`99\%` of true photons (dt99).
1105 
1106 The photon timing uncertainty is currently determined using MC. The resulting parametrization depends on
1107 the true energy deposition in the highest energetic crystal and the local beam background level in that crystal.
1108 The resulting timing distribution is non-Gaussian and for each photon the value dt99 is stored,
1109 where :math:`|\text{timing}| / \text{dt99} < 1` is designed to give a :math:`99\%`
1110 timing efficiency for true photons from the IP.
1111 The resulting efficiency is approximately flat in energy and independent of beam background levels.
1112 
1113 Very large values of dt99 are an indication of failed waveform fits in the ECL.
1114 We remove such clusters in most physics photon lists.
1115 
1116 .. note::
1117  | Please read `this <importantNoteECL>` first.
1118  | Lower limit: :math:`0.0`
1119  | Upper limit: :math:`1000.0`
1120  | Precision: :math:`12` bit
1121 
1122 .. warning::
1123  In real data there will be a sizeable number of high energetic Bhabha events
1124  (from previous or later bunch collisions) that can easily be rejected by timing cuts.
1125  However, these events create large ECL clusters that can overlap with other ECL clusters
1126  and it is not clear that a simple rejection is the correction strategy.
1127 ..
1128 )DOC","ns");
1129  REGISTER_VARIABLE("clusterHasFailedErrorTiming", eclClusterHasFailedErrorTiming, R"DOC(
1130 Status bit for if the ECL cluster's timing uncertainty calculation failed. Photon timing is given by the fitted time
1131 of the recorded waveform of the highest energetic crystal in a cluster; however, that fit can fail and so
1132 this variable tells the user if that has happened.
1133 )DOC");
1134  REGISTER_VARIABLE("clusterHighestE", eclClusterHighestE, R"DOC(
1135 Returns energy of the highest energetic crystal in the ECL cluster after reweighting.
1136 
1137 .. warning::
1138  This variable must be used carefully since it can bias shower selection
1139  towards photons that hit crystals in the center and hence have a large energy
1140  deposition in the highest energy crystal.
1141 
1142 .. note::
1143  | Please read `this <importantNoteECL>` first.
1144  | Lower limit: :math:`-5` (:math:`e^{-5} = 0.00674\,` GeV)
1145  | Upper limit: :math:`3.0` (:math:`e^3 = 20.08553\,` GeV)
1146  | Precision: :math:`18` bit
1147 ..
1148 )DOC","GeV");
1149  REGISTER_VARIABLE("clusterCellID", eclClusterCellId,
1150  "Returns cellId of the crystal with highest energy in the ECLCluster.");
1151  REGISTER_VARIABLE("clusterThetaID", eclClusterThetaId,
1152  "Returns thetaId of the crystal with highest energy in the ECLCluster.");
1153  REGISTER_VARIABLE("clusterPhiID", eclClusterPhiId,
1154  "Returns phiId of the crystal with highest energy in the ECLCluster.");
1155  REGISTER_VARIABLE("clusterE1E9", eclClusterE1E9, R"DOC(
1156 Returns ratio of energies of the central crystal, E1, and 3x3 crystals, E9, around the central crystal.
1157 Since :math:`E1 \leq E9`, this ratio is :math:`\leq 1` and tends towards larger values for photons
1158 and smaller values for hadrons.
1159 
1160 .. note::
1161  | Please read `this <importantNoteECL>` first.
1162  | Lower limit: :math:`0.0`
1163  | Upper limit: :math:`1.0`
1164  | Precision: :math:`10` bit
1165 )DOC");
1166  REGISTER_VARIABLE("clusterE9E25", eclClusterE9E25, R"DOC(
1167 Deprecated - kept for backwards compatibility - returns clusterE9E21.
1168 )DOC");
1169  REGISTER_VARIABLE("clusterE9E21", eclClusterE9E21, R"DOC(
1170 Returns ratio of energies in inner 3x3 crystals, E9, and 5x5 crystals around the central crystal without corners.
1171 Since :math:`E9 \leq E21`, this ratio is :math:`\leq 1` and tends towards larger values for photons
1172 and smaller values for hadrons.
1173 
1174 .. note::
1175  | Please read `this <importantNoteECL>` first.
1176  | Lower limit: :math:`0.0`
1177  | Upper limit: :math:`1.0`
1178  | Precision: :math:`10` bit
1179 )DOC");
1180  REGISTER_VARIABLE("clusterAbsZernikeMoment40", eclClusterAbsZernikeMoment40, R"DOC(
1181 Returns absolute value of Zernike moment 40 (:math:`|Z_{40}|`). (shower shape variable).
1182 
1183 .. note::
1184  | Please read `this <importantNoteECL>` first.
1185  | Lower limit: :math:`0.0`
1186  | Upper limit: :math:`1.7`
1187  | Precision: :math:`10` bit
1188 )DOC");
1189  REGISTER_VARIABLE("clusterAbsZernikeMoment51", eclClusterAbsZernikeMoment51, R"DOC(
1190 Returns absolute value of Zernike moment 51 (:math:`|Z_{51}|`). (shower shape variable).
1191 
1192 .. note::
1193  | Please read `this <importantNoteECL>` first.
1194  | Lower limit: :math:`0.0`
1195  | Upper limit: :math:`1.2`
1196  | Precision: :math:`10` bit
1197 )DOC");
1198  REGISTER_VARIABLE("clusterZernikeMVA", eclClusterZernikeMVA, R"DOC(
1199 Returns output of a MVA using eleven Zernike moments of the cluster. Zernike moments are calculated per
1200 shower in a plane perpendicular to the shower direction via
1201 
1202 .. math::
1203  |Z_{nm}| = \frac{n+1}{\pi} \frac{1}{\sum_{i} w_{i} E_{i}} \left|\sum_{i} R_{nm}(\rho_{i}) e^{-im\alpha_{i}} w_{i} E_{i} \right|
1204 
1205 where n, m are the integers, :math:`i` runs over the crystals in the shower,
1206 :math:`E_{i}` is the energy of the i-th crystal in the shower,
1207 :math:`R_{nm}` is a polynomial of degree :math:`n`,
1208 :math:`\rho_{i}` is the radial distance of the :math:`i`-th crystal in the perpendicular plane,
1209 and :math:`\alpha_{i}` is the polar angle of the :math:`i`-th crystal in the perpendicular plane.
1210 As a crystal can be related to more than one shower, :math:`w_{i}` is the fraction of the
1211 energy of the :math:`i`-th crystal associated with the shower.
1212 
1213 More details about the implementation can be found in `BELLE2-NOTE-TE-2017-001 <https://docs.belle2.org/record/454?ln=en>`_ .
1214 
1215 More details about Zernike polynomials can be found in `Wikipedia <https://en.wikipedia.org/wiki/Zernike_polynomials>`_ .
1216 
1217 | For cluster with hypothesisId==N1: raw MVA output.
1218 | For cluster with hypothesisId==N2: 1 - prod{clusterZernikeMVA}, where the product is on all N1 showers
1219  belonging to the same connected region (shower shape variable).
1220 
1221 .. note::
1222  | Please read `this <importantNoteECL>` first.
1223  | Lower limit: :math:`0.0`
1224  | Upper limit: :math:`1.0`
1225  | Precision: :math:`10` bit
1226 )DOC");
1227  REGISTER_VARIABLE("clusterSecondMoment", eclClusterSecondMoment, R"DOC(
1228 Returns second moment :math:`S`. It is defined as:
1229 
1230 .. math::
1231  S = \frac{\sum_{i=0}^{n} w_{i} E_{i} r^2_{i}}{\sum_{i=0}^{n} w_{i} E_{i}}
1232 
1233 where :math:`E_{i} = (E_0, E_1, ...)` are the single crystal energies sorted by energy, :math:`w_{i}` is
1234 the crystal weight, and :math:`r_{i}` is the distance of the :math:`i`-th digit to the shower center projected
1235 to a plane perpendicular to the shower axis.
1236 
1237 .. note::
1238  | Please read `this <importantNoteECL>` first.
1239  | Lower limit: :math:`0.0`
1240  | Upper limit: :math:`40.0`
1241  | Precision: :math:`10` bit
1242 ..
1243 )DOC",":math:`\\text{cm}^2`");
1244  REGISTER_VARIABLE("clusterLAT", eclClusterLAT, R"DOC(
1245 Returns lateral energy distribution (shower variable). It is defined as following:
1246 
1247 .. math::
1248  S = \frac{\sum_{i=2}^{n} w_{i} E_{i} r^2_{i}}{(w_{0} E_{0} + w_{1} E_{1}) r^2_{0} + \sum_{i=2}^{n} w_{i} E_{i} r^2_{i}}
1249 
1250 where :math:`E_{i} = (E_{0}, E_{1}, ...)` are the single crystal energies sorted by energy
1251 (:math:`E_{0}` is the highest energy and :math:`E_{1}` the second highest), :math:`w_{i}`
1252 is the crystal weight, :math:`r_{i}` is the distance of the :math:`i`-th digit to the
1253 shower center projected to a plane perpendicular to the shower axis,
1254 and :math:`r_{0} \approx 5\,cm` is the distance between two crystals.
1255 
1256 clusterLAT peaks around 0.3 for radially symmetrical electromagnetic showers and is larger
1257 for hadronic events, and electrons with a close-by radiative or Bremsstrahlung photon.
1258 
1259 .. note::
1260  | Please read `this <importantNoteECL>` first.
1261  | Lower limit: :math:`0.0`
1262  | Upper limit: :math:`1.0`
1263  | Precision: :math:`10` bit
1264 )DOC");
1265  REGISTER_VARIABLE("clusterNHits", eclClusterNHits, R"DOC(
1266 Returns sum of weights :math:`w_{i}` (:math:`w_{i} \leq 1`) of all crystals in an ECL cluster.
1267 For non-overlapping clusters this is equal to the number of crystals in the cluster.
1268 In case of energy splitting among nearby clusters, this can be a non-integer value.
1269 
1270 .. note::
1271  | Please read `this <importantNoteECL>` first.
1272  | Lower limit: :math:`0.0`
1273  | Upper limit: :math:`200.0`
1274  | Precision: :math:`10` bit
1275  | If fractional weights are not of interest, this value should be cast to the nearest integer.
1276 )DOC");
1277  REGISTER_VARIABLE("clusterTrackMatch", eclClusterTrackMatched, R"DOC(
1278 Returns 1.0 if at least one reconstructed charged track is matched to the ECL cluster.
1279 
1280 Every reconstructed charged track is extrapolated into the ECL.
1281 Every ECL crystal that is crossed by the track extrapolation is marked.
1282 Each ECL cluster that contains any marked crystal is matched to the track.
1283 Multiple tracks can be matched to one cluster and multiple clusters can be matched to one track.
1284 It is conceptually correct to have two tracks matched to the same cluster.
1285 )DOC");
1286  REGISTER_VARIABLE("nECLClusterTrackMatches", nECLClusterTrackMatches, R"DOC(
1287 Returns number of charged tracks matched to this cluster.
1288 
1289 .. note::
1290  Sometimes (perfectly correctly) two tracks are extrapolated into the same cluster.
1291 
1292  - For charged particles, this should return at least 1 (but sometimes 2 or more).
1293  - For neutrals, this should always return 0.
1294  - Returns NaN if there is no cluster.
1295 )DOC");
1296  REGISTER_VARIABLE("clusterHasPulseShapeDiscrimination", eclClusterHasPulseShapeDiscrimination, R"DOC(
1297 Status bit to indicate if cluster has digits with waveforms that passed energy and :math:`\chi^2`
1298 thresholds for computing PSD variables.
1299 )DOC");
1300  REGISTER_VARIABLE("beamBackgroundSuppression", beamBackgroundSuppression, R"DOC(
1301 Returns the output of an MVA classifier that uses shower-related variables to distinguish true photon clusters from beam background clusters.
1302 The classes are:
1303 
1304  - 1 for true photon clusters
1305  - 0 for beam background clusters
1306 
1307 The MVA has been trained using samples of signal photons and beam background photons coming from MC. The features used are (in decreasing order of significance):
1308 
1309  - `clusterTiming`
1310  - `clusterPulseShapeDiscriminationMVA`
1311  - `clusterE`
1312  - `clusterTheta`
1313  - `clusterZernikeMVA`
1314  - `clusterE1E9`
1315  - `clusterLAT`
1316  - `clusterSecondMoment`
1317 )DOC");
1318  REGISTER_VARIABLE("hadronicSplitOffSuppression", hadronicSplitOffSuppression, R"DOC(
1319 Returns the output of an MVA classifier that uses shower-related variables to distinguish true photon clusters from hadronic splitoff clusters.
1320 The classes are:
1321 
1322  - 1 for true photon clusters
1323  - 0 for hadronic splitoff clusters
1324 
1325 The MVA has been trained using samples of signal photons and hadronic splitoff photons coming from MC. The features used are (in decreasing order of significance):
1326 
1327  - `clusterPulseShapeDiscriminationMVA`
1328  - `minC2TDist`
1329  - `clusterZernikeMVA`
1330  - `clusterE`
1331  - `clusterLAT`
1332  - `clusterE1E9`
1333  - `clusterSecondMoment`
1334 )DOC");
1335  REGISTER_VARIABLE("clusterKlId", eclClusterKlId, R"DOC(
1336 Returns MVA classifier that uses ECL clusters variables to discriminate Klong clusters from em background.
1337 
1338  - 1 for Kl
1339  - 0 for background
1340 )DOC");
1341  REGISTER_VARIABLE("clusterPulseShapeDiscriminationMVA", eclPulseShapeDiscriminationMVA, R"DOC(
1342 Returns MVA classifier that uses pulse shape discrimination to identify electromagnetic vs hadronic showers.
1343 
1344  - 1 for electromagnetic showers
1345  - 0 for hadronic showers
1346 )DOC");
1347  REGISTER_VARIABLE("clusterNumberOfHadronDigits", eclClusterNumberOfHadronDigits, R"DOC(
1348 Returns ECL cluster's number of hadron digits in cluster (pulse shape discrimination variable).
1349 Weighted sum of digits in cluster with significant scintillation emission (:math:`> 3\,` MeV)
1350 in the hadronic scintillation component.
1351 Computed only using cluster digits with energy :math:`> 50\,` MeV and good offline waveform fit :math:`\chi^2`.
1352 
1353 .. note::
1354  | Please read `this <importantNoteECL>` first.
1355  | Lower limit: :math:`0.0`
1356  | Upper limit: :math:`255.0`
1357  | Precision: :math:`18` bit
1358 )DOC");
1359  REGISTER_VARIABLE("clusterClusterID", eclClusterId, R"DOC(
1360 Returns ECL cluster ID of this ECL cluster within the connected region (CR) to which it belongs to.
1361 )DOC");
1362  REGISTER_VARIABLE("clusterHasNPhotons", eclClusterHasNPhotonsHypothesis, R"DOC(
1363 Returns 1.0 if cluster has the 'N photons' hypothesis (historically called 'N1'),
1364 0.0 if not, and NaN if no cluster is associated to the particle.
1365 )DOC");
1366  REGISTER_VARIABLE("clusterHasNeutralHadron", eclClusterHasNeutralHadronHypothesis, R"DOC(
1367 Returns 1.0 if the cluster has the 'neutral hadrons' hypothesis (historically called 'N2'),
1368 0.0 if not, and NaN if no cluster is associated to the particle.
1369 )DOC");
1370  REGISTER_VARIABLE("eclExtTheta", eclExtTheta, R"DOC(
1371 Returns extrapolated :math:`\theta` of particle track associated to the cluster (if any). Requires module ECLTrackCalDigitMatch to be executed.
1372 )DOC","rad");
1373  REGISTER_VARIABLE("eclExtPhi", eclExtPhi, R"DOC(
1374 Returns extrapolated :math:`\phi` of particle track associated to the cluster (if any). Requires module ECLTrackCalDigitMatch to be executed..
1375 )DOC","rad");
1376  REGISTER_VARIABLE("eclExtPhiId", eclExtPhiId, R"DOC(
1377 Returns extrapolated :math:`\phi` ID of particle track associated to the cluster (if any). Requires module ECLTrackCalDigitMatch to be executed..
1378 )DOC");
1379  REGISTER_VARIABLE("weightedAverageECLTime", weightedAverageECLTime, R"DOC(
1380 Returns ECL weighted average time of all clusters (neutrals) and matched clusters (charged) of daughters
1381 (of any generation) of the provided particle.
1382 )DOC", "ns");
1383  REGISTER_VARIABLE("maxWeightedDistanceFromAverageECLTime", maxWeightedDistanceFromAverageECLTime, R"DOC(
1384 Returns maximum weighted distance between time of the cluster of a photon and the ECL average time, amongst
1385 the clusters (neutrals) and matched clusters (charged) of daughters (of all generations) of the provided particle.
1386 )DOC", "ns");
1387  REGISTER_VARIABLE("clusterMdstIndex", eclClusterMdstIndex, R"DOC(
1388 StoreArray index(0 - based) of the MDST ECLCluster (useful for track-based particles matched to a cluster).
1389 )DOC");
1390 
1391  REGISTER_VARIABLE("nECLOutOfTimeCrystals", nECLOutOfTimeCrystals, R"DOC(
1392 [Eventbased] Returns the number of crystals (ECLCalDigits) that are out of time.
1393 )DOC");
1394 
1395  REGISTER_VARIABLE("nECLOutOfTimeCrystalsFWDEndcap", nECLOutOfTimeCrystalsFWDEndcap, R"DOC(
1396 [Eventbased] Returns the number of crystals (ECLCalDigits) that are out of time in the forward endcap.
1397 )DOC");
1398 
1399  REGISTER_VARIABLE("nECLOutOfTimeCrystalsBarrel", nECLOutOfTimeCrystalsBarrel, R"DOC(
1400 [Eventbased] Returns the number of crystals (ECLCalDigits) that are out of time in the barrel.
1401 )DOC");
1402 
1403  REGISTER_VARIABLE("nECLOutOfTimeCrystalsBWDEndcap", nECLOutOfTimeCrystalsBWDEndcap, R"DOC(
1404 [Eventbased] Returns the number of crystals (ECLCalDigits) that are out of time in the backward endcap.
1405 )DOC");
1406 
1407  REGISTER_VARIABLE("nRejectedECLShowers", nRejectedECLShowers, R"DOC(
1408 [Eventbased] Returns the number of showers in the ECL that do not become clusters.
1409 )DOC");
1410 
1411  REGISTER_VARIABLE("nRejectedECLShowersFWDEndcap", nRejectedECLShowersFWDEndcap, R"DOC(
1412 [Eventbased] Returns the number of showers in the ECL that do not become clusters, from the forward endcap.
1413 If the number exceeds 255 (uint8_t maximum value) the variable is set to 255.
1414 )DOC");
1415 
1416  REGISTER_VARIABLE("nRejectedECLShowersBarrel", nRejectedECLShowersBarrel, R"DOC(
1417 [Eventbased] Returns the number of showers in the ECL that do not become clusters, from the barrel.
1418 If the number exceeds 255 (uint8_t maximum value) the variable is set to 255.
1419 )DOC");
1420 
1421  REGISTER_VARIABLE("nRejectedECLShowersBWDEndcap", nRejectedECLShowersBWDEndcap, R"DOC(
1422 [Eventbased] Returns the number of showers in the ECL that do not become clusters, from the backward endcap.
1423 If the number exceeds 255 (uint8_t maximum value) the variable is set to 255.
1424 )DOC");
1425 
1426  REGISTER_VARIABLE("eclClusterOnlyInvariantMass", eclClusterOnlyInvariantMass, R"DOC(
1427 [Expert] The invariant mass calculated from all ECLCluster daughters (i.e. photons) and
1428 cluster-matched tracks using the cluster 4-momenta.
1429 
1430 Used for ECL-based dark sector physics and debugging track-cluster matching.
1431 )DOC","GeV/:math:`\\text{c}^2`");
1432 
1433  REGISTER_METAVARIABLE("photonHasOverlap(cutString, photonlistname, tracklistname)", photonHasOverlap, R"DOC(
1434  Returns true if the connected ECL region of the particle's cluster is shared by another particle's cluster.
1435  Neutral and charged cluster are considered.
1436  A cut string can be provided to ignore cluster that do not satisfy the given criteria.
1437  By default, the ParticleList ``gamma:all`` is used for the check of neutral ECL cluster,
1438  and the ParticleList ``e-:all`` for the check of charged ECL cluster.
1439  However, one can customize the name of the ParticleLists via additional arguments.
1440  If no argument or only a cut string is provided and ``gamma:all`` or ``e-:all`` does not exist
1441  or if the variable is requested for a particle that is not a photon, NaN is returned.
1442  )DOC", Manager::VariableDataType::c_double);
1443 
1444  REGISTER_VARIABLE("clusterUncorrE", eclClusterUncorrectedE, R"DOC(
1445 [Expert] [Calibration] Returns ECL cluster's uncorrected energy. That is, before leakage corrections.
1446 This variable should only be used for study of the ECL. Please see :b2:var:`clusterE`.
1447 )DOC","GeV");
1448 
1449  }
1451 }
int getPDGCode() const
PDG code.
Definition: Const.h:464
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:609
static const ParticleType photon
photon particle
Definition: Const.h:664
EHypothesisBit
The hypothesis bits for this ECLCluster (Connected region (CR) is split using this hypothesis.
Definition: ECLCluster.h:31
@ c_nPhotons
CR is split into n photons (N1)
@ c_neutralHadron
CR is reconstructed as a neutral hadron (N2)
static std::unique_ptr< GeneralCut > compile(const std::string &cut)
Creates an instance of a cut and returns a unique_ptr to it, if you need a copy-able object instead y...
Definition: GeneralCut.h:84
static const ReferenceFrame & GetCurrent()
Get current rest frame.
std::function< VarVariant(const Particle *)> FunctionPtr
functions stored take a const Particle* and return VarVariant.
Definition: Manager.h:113
const Var * getVariable(std::string name)
Get the variable belonging to the given key.
Definition: Manager.cc:57
static Manager & Instance()
get singleton instance.
Definition: Manager.cc:25
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:23