Belle II Software  release-08-01-10
TOPDigitVariables.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 #include <top/variables/TOPDigitVariables.h>
9 
10 // // framework
11 #include <framework/gearbox/Const.h>
12 #include <framework/datastore/StoreObjPtr.h>
13 #include <framework/datastore/StoreArray.h>
14 
15 // // dataobjects
16 #include <tracking/dataobjects/ExtHit.h>
17 #include <mdst/dataobjects/Track.h>
18 #include <analysis/dataobjects/Particle.h>
19 #include <top/dataobjects/TOPDigit.h>
20 #include <top/geometry/TOPGeometryPar.h>
21 #include <top/dataobjects/TOPLikelihood.h>
22 #include <top/dataobjects/TOPRecBunch.h>
23 #include <top/dataobjects/TOPLikelihoodScanResult.h>
24 
25 #include <top/dataobjects/TOPBarHit.h>
26 #include <mdst/dataobjects/MCParticle.h>
27 #include <framework/dataobjects/MCInitialParticles.h>
28 
29 #include <algorithm> // for sort
30 
31 using namespace std;
32 using namespace ROOT::Math;
33 
34 namespace Belle2 {
39  namespace Variable {
40  namespace TOPVariable {
41 
42  //---------------- helpers --------------------
43 
44  const TOPLikelihood* getTOPLikelihood(const Particle* particle)
45  {
46  if (not particle) return nullptr;
47  const auto* track = particle->getTrack();
48  return track ? track->getRelated<TOPLikelihood>() : nullptr;
49  }
50 
51  const ExtHit* getExtHit(const Particle* particle)
52  {
53  const auto* topLikelihood = getTOPLikelihood(particle);
54  return topLikelihood ? topLikelihood->getRelated<ExtHit>() : nullptr;
55  }
56 
57  const TOPBarHit* getBarHit(const Particle* particle)
58  {
59  const auto* topLikelihood = getTOPLikelihood(particle);
60  return topLikelihood ? topLikelihood->getRelated<TOPBarHit>() : nullptr;
61  }
62 
63  const TOPLikelihoodScanResult* getTOPLikelihoodScanResult(const Particle* particle)
64  {
65  if (not particle) return nullptr;
66  const auto* track = particle->getTrack();
67  if (not track) return nullptr;
68  auto scanRes = track->getRelated<TOPLikelihoodScanResult>();
69  if (not scanRes) {
70  B2WARNING("No TOPLikelihoodScanResult object found. Are you sure you added TOPLLScanner to the path?");
71  return nullptr;
72  }
73  return scanRes;
74  }
75 
76  int getSlotID(const Particle* particle)
77  {
78  const auto* extHit = getExtHit(particle);
79  return extHit ? extHit->getCopyID() : 0;
80  }
81 
82 
83  //---------------- ExtHit (or TOPBarHit) based --------------------
84 
85  double topSlotID(const Particle* particle)
86  {
87  const auto* extHit = getExtHit(particle);
88  return extHit ? extHit->getCopyID() : std::numeric_limits<double>::quiet_NaN();
89  }
90 
91  double topSlotIDMCMatch(const Particle* particle)
92  {
93  const auto* barHit = getBarHit(particle);
94  return barHit ? barHit->getModuleID() : std::numeric_limits<double>::quiet_NaN();
95  }
96 
97  bool getLocalPosition(const Particle* particle, XYZPoint& result)
98  {
99  const auto* extHit = getExtHit(particle);
100  if (not extHit) return false;
101  int slotID = extHit->getCopyID();
102  const auto& position = static_cast<XYZPoint>(extHit->getPosition());
103  const auto* geo = TOP::TOPGeometryPar::Instance()->getGeometry();
104  if (not geo or not geo->isModuleIDValid(slotID)) return false;
105  const auto& module = geo->getModule(slotID);
106  result = module.pointToLocal(position);
107  return true;
108  }
109 
110  double getTOPLocalX(const Particle* particle)
111  {
112  XYZPoint position;
113  bool ok = TOPVariable::getLocalPosition(particle, position);
114  if (not ok) return std::numeric_limits<double>::quiet_NaN();
115  return position.X();
116  }
117 
118  double getTOPLocalY(const Particle* particle)
119  {
120  XYZPoint position;
121  bool ok = TOPVariable::getLocalPosition(particle, position);
122  if (not ok) return std::numeric_limits<double>::quiet_NaN();
123  return position.Y();
124  }
125 
126  double getTOPLocalZ(const Particle* particle)
127  {
128  XYZPoint position;
129  bool ok = TOPVariable::getLocalPosition(particle, position);
130  if (not ok) return std::numeric_limits<double>::quiet_NaN();
131  return position.Z();
132  }
133 
134  bool getLocalPositionMCMatch(const Particle* particle, XYZPoint& result)
135  {
136  const auto* barHit = getBarHit(particle);
137  if (not barHit) return false;
138  int slotID = barHit->getModuleID();
139  const auto& position = barHit->getPosition();
140  const auto* geo = TOP::TOPGeometryPar::Instance()->getGeometry();
141  if (not geo or not geo->isModuleIDValid(slotID)) return false;
142  const auto& module = geo->getModule(slotID);
143  result = module.pointToLocal(position);
144  return true;
145  }
146 
147  double getTOPLocalXMCMatch(const Particle* particle)
148  {
149  XYZPoint position;
150  bool ok = TOPVariable::getLocalPositionMCMatch(particle, position);
151  if (not ok) return std::numeric_limits<double>::quiet_NaN();
152  return position.X();
153  }
154 
155  double getTOPLocalYMCMatch(const Particle* particle)
156  {
157  XYZPoint position;
158  bool ok = TOPVariable::getLocalPositionMCMatch(particle, position);
159  if (not ok) return std::numeric_limits<double>::quiet_NaN();
160  return position.Y();
161  }
162 
163  double getTOPLocalZMCMatch(const Particle* particle)
164  {
165  XYZPoint position;
166  bool ok = TOPVariable::getLocalPositionMCMatch(particle, position);
167  if (not ok) return std::numeric_limits<double>::quiet_NaN();
168  return position.Z();
169  }
170 
171  bool getLocalMomentum(const Particle* particle, XYZVector& result)
172  {
173  const auto* extHit = getExtHit(particle);
174  if (not extHit) return false;
175  int slotID = extHit->getCopyID();
176  const auto& momentum = extHit->getMomentum();
177  const auto* geo = TOP::TOPGeometryPar::Instance()->getGeometry();
178  if ((not geo) or (not geo->isModuleIDValid(slotID))) return false;
179  const auto& module = geo->getModule(slotID);
180  result = module.momentumToLocal(momentum);
181  return true;
182  }
183 
184  double getTOPLocalPhi(const Particle* particle)
185  {
186  XYZVector momentum;
187  bool ok = TOPVariable::getLocalMomentum(particle, momentum);
188  if (not ok) return std::numeric_limits<double>::quiet_NaN();
189  return momentum.Phi();
190  }
191 
192  double getTOPLocalTheta(const Particle* particle)
193  {
194  XYZVector momentum;
195  bool ok = TOPVariable::getLocalMomentum(particle, momentum);
196  if (not ok) return std::numeric_limits<double>::quiet_NaN();
197  return momentum.Theta();
198  }
199 
200  bool getLocalMomentumMCMatch(const Particle* particle, XYZVector& result)
201  {
202  const auto* barHit = getBarHit(particle);
203  if (not barHit) return false;
204  int slotID = barHit->getModuleID();
205  const auto& momentum = barHit->getMomentum();
206  const auto* geo = TOP::TOPGeometryPar::Instance()->getGeometry();
207  if ((not geo) or (not geo->isModuleIDValid(slotID))) return false;
208  const auto& module = geo->getModule(slotID);
209  result = module.momentumToLocal(momentum);
210  return true;
211  }
212 
213  double getTOPLocalPhiMCMatch(const Particle* particle)
214  {
215  XYZVector momentum;
216  bool ok = TOPVariable::getLocalMomentumMCMatch(particle, momentum);
217  if (not ok) return std::numeric_limits<double>::quiet_NaN();
218  return momentum.Phi();
219  }
220 
221  double getTOPLocalThetaMCMatch(const Particle* particle)
222  {
223  XYZVector momentum;
224  bool ok = TOPVariable::getLocalMomentumMCMatch(particle, momentum);
225  if (not ok) return std::numeric_limits<double>::quiet_NaN();
226  return momentum.Theta();
227  }
228 
229  double computeTOF(const Particle* particle, int pdg)
230  {
231  const auto* extHit = getExtHit(particle);
232  if (not extHit) return std::numeric_limits<double>::quiet_NaN();
233  auto extPDGCode = abs(extHit->getPdgCode());
234  double pmom = particle->getMomentumMagnitude();
235  double massExtHit = Const::ChargedStable(extPDGCode).getMass();
236  double betaExtHit = pmom / sqrt(pmom * pmom + massExtHit * massExtHit);
237  double mass = pdg == 0 ? particle->getMass() : Const::ChargedStable(abs(pdg)).getMass();
238  double beta = pmom / sqrt(pmom * pmom + mass * mass);
239  return extHit->getTOF() * betaExtHit / beta;
240  }
241 
242  double getTOF(const Particle* particle)
243  {
244  return computeTOF(particle, 0);
245  }
246 
247  double getTOFMCMatch(const Particle* particle)
248  {
249  const auto* barHit = getBarHit(particle);
250  if (not barHit) return std::numeric_limits<double>::quiet_NaN();
251  StoreObjPtr<MCInitialParticles> mcInitialParticles;
252  double trueEventT0 = 0;
253  if (mcInitialParticles.isValid()) trueEventT0 = mcInitialParticles->getTime();
254  return barHit->getTime() - trueEventT0;
255  }
256 
257  double getTOFExpert(const Particle* particle, const vector<double>& vars)
258  {
259  if (vars.size() != 1) {
260  B2FATAL("topTOFExpert(pdg): Need exactly one parameter (PDG code).");
261  }
262  int pdg = static_cast<int>(vars[0]);
263  return computeTOF(particle, pdg);
264  }
265 
266  //---------------- Helix extrapolation --------------------
267 
268  double extrapTrackToTOPz(const Particle* particle)
269  {
270  auto trk = particle->getTrack();
271  if (not trk) return std::numeric_limits<double>::quiet_NaN();
272  auto trkfit = trk->getTrackFitResultWithClosestMass(Const::ChargedStable(std::abs(particle->getPDGCode())));
273  if (not trkfit) return std::numeric_limits<double>::quiet_NaN();
274  auto helix = trkfit->getHelix();
275  double arcLength = helix.getArcLength2DAtCylindricalR(120);
276  const auto& result = helix.getPositionAtArcLength2D(arcLength);
277  return result.Z();
278  }
279 
280  double extrapTrackToTOPtheta(const Particle* particle)
281  {
282  auto trk = particle->getTrack();
283  if (not trk) return std::numeric_limits<double>::quiet_NaN();
284  auto trkfit = trk->getTrackFitResultWithClosestMass(Const::ChargedStable(std::abs(particle->getPDGCode())));
285  if (not trkfit) return std::numeric_limits<double>::quiet_NaN();
286  auto helix = trkfit->getHelix();
287  double arcLength = helix.getArcLength2DAtCylindricalR(120);
288  const auto& result = helix.getPositionAtArcLength2D(arcLength);
289  return result.Theta();
290  }
291 
292  double extrapTrackToTOPphi(const Particle* particle)
293  {
294  auto trk = particle->getTrack();
295  if (not trk) return std::numeric_limits<double>::quiet_NaN();
296  auto trkfit = trk->getTrackFitResultWithClosestMass(Const::ChargedStable(std::abs(particle->getPDGCode())));
297  if (not trkfit) return std::numeric_limits<double>::quiet_NaN();
298  auto helix = trkfit->getHelix();
299  double arcLength = helix.getArcLength2DAtCylindricalR(120);
300  const auto& result = helix.getPositionAtArcLength2D(arcLength);
301  return result.Phi();
302  }
303 
304  //---------------- TOPDigit based --------------------
305 
306  double countHits(const Particle* particle, double tmin, double tmax, bool clean)
307  {
308  int slotID = getSlotID(particle);
309  if (slotID == 0) return std::numeric_limits<double>::quiet_NaN();
310 
311  StoreArray<TOPDigit> digits;
312  int count = 0;
313  for (const auto& digit : digits) {
314  if (digit.getModuleID() != slotID) continue;
315  // skip bad digits only when we want to clean
316  if (clean and digit.getHitQuality() != TOPDigit::c_Good) continue;
317  if (digit.getTime() < tmin or digit.getTime() > tmax) continue;
318  count++;
319  }
320  return count;
321  }
322 
323  double topDigitCount(const Particle* particle)
324  {
325  int slotID = getSlotID(particle);
326  if (slotID == 0) return std::numeric_limits<double>::quiet_NaN();
327 
328  StoreArray<TOPDigit> topDigits;
329  int count = 0;
330  for (const auto& t : topDigits) {
331  if (t.getModuleID() != slotID) continue;
332  if (t.getHitQuality() != TOPDigit::c_Good) continue;
333  count++;
334  }
335  return count;
336  }
337 
338  double topDigitCountMCMatch(const Particle* particle)
339  {
340  int slotID = getSlotID(particle);
341  if (slotID == 0) return std::numeric_limits<double>::quiet_NaN();
342 
343  auto* trk = particle->getTrack();
344  if (not trk) return std::numeric_limits<double>::quiet_NaN();
345 
346  auto* mcParticle = trk->getRelated<MCParticle>();
347  if (not mcParticle) return std::numeric_limits<double>::quiet_NaN();
348 
349  auto digits = mcParticle->getRelationsWith<TOPDigit>();
350  int count = 0;
351  for (const auto& digit : digits) {
352  if (digit.getModuleID() != slotID) continue;
353  if (digit.getHitQuality() != TOPDigit::c_Good) continue;
354  count++;
355  }
356  return count;
357  }
358 
359  double topSignalDigitCount(const Particle* particle)
360  {
361  int slotID = getSlotID(particle);
362  if (slotID == 0) return std::numeric_limits<double>::quiet_NaN();
363 
364  StoreArray<TOPDigit> digits;
365  int count = 0;
366  for (const auto& digit : digits) {
367  if (digit.getModuleID() != slotID) continue;
368  if (digit.getHitQuality() != TOPDigit::c_Good) continue;
369  if (abs(digit.getTime()) > 50) continue;
370  if (digit.getTime() > 0) count++;
371  else count--;
372  }
373  return count;
374  }
375 
376  double topBackgroundDigitCount(const Particle* particle)
377  {
378  return TOPVariable::countHits(particle, -50, 0, true);
379  }
380 
381  double topRawDigitCount(const Particle* particle)
382  {
383  int slotID = getSlotID(particle);
384  if (slotID == 0) return std::numeric_limits<double>::quiet_NaN();
385 
386  StoreArray<TOPDigit> topDigits;
387  int count = 0;
388  for (const auto& t : topDigits) {
389  if (t.getModuleID() != slotID) continue;
390  count++;
391  }
392  return count;
393  }
394 
395  double countTOPHitsInInterval(const Particle* particle, const vector<double>& vars)
396  {
397  if (vars.size() != 2) {
398  B2FATAL("countTOPHitsInInterval(tmin, tmax): Need exactly two parameters (tmin, tmax)");
399  }
400  return TOPVariable::countHits(particle, vars[0], vars[1], true);
401  }
402 
403  double countRawTOPHitsInInterval(const Particle* particle, const vector<double>& vars)
404  {
405  if (vars.size() != 2) {
406  B2FATAL("countRawTOPHitsInInterval(tmin, tmax): Need exactly two parameters (tmin, tmax)");
407  }
408  return TOPVariable::countHits(particle, vars[0], vars[1], false);
409  }
410 
411  double countTOPHitsInIntervalMCMatch(const Particle* particle, const std::vector<double>& vars)
412  {
413  if (vars.size() != 2) {
414  B2FATAL("countTOPHitsInIntervalMCMatch(tmin, tmax): Need exactly two parameters (tmin, tmax)");
415  }
416  int slotID = getSlotID(particle);
417  if (slotID == 0) return std::numeric_limits<double>::quiet_NaN();
418 
419  auto* trk = particle->getTrack();
420  if (not trk) return std::numeric_limits<double>::quiet_NaN();
421 
422  auto* mcParticle = trk->getRelated<MCParticle>();
423  if (not mcParticle) return std::numeric_limits<double>::quiet_NaN();
424 
425  auto digits = mcParticle->getRelationsWith<TOPDigit>();
426  int count = 0;
427  for (const auto& digit : digits) {
428  if (digit.getModuleID() != slotID) continue;
429  if (digit.getHitQuality() != TOPDigit::c_Good) continue;
430  if (digit.getTime() < vars[0] or digit.getTime() > vars[1]) continue;
431  count++;
432  }
433  return count;
434  }
435 
436 
437  //---------------- TOPLikelihood based --------------------
438 
439  double getFlag(const Particle* particle)
440  {
441  const auto* topLikelihood = getTOPLikelihood(particle);
442  if (not topLikelihood) return 0;
443  return (topLikelihood->getFlag() == 1);
444  }
445 
446  double getTOPPhotonCount(const Particle* particle)
447  {
448  const auto* topLikelihood = TOPVariable::getTOPLikelihood(particle);
449  return topLikelihood ? topLikelihood->getNphot() : std::numeric_limits<double>::quiet_NaN();
450  }
451 
452  double expectedPhotonCount(const Particle* particle, int pdg)
453  {
454  const auto* topLikelihood = getTOPLikelihood(particle);
455  if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
456 
457  pdg = pdg != 0 ? pdg : particle->getPDGCode();
458  const auto& chargedStable = Const::chargedStableSet.find(abs(pdg));
459  if (chargedStable == Const::invalidParticle) return std::numeric_limits<double>::quiet_NaN();
460 
461  return topLikelihood->getEstPhot(chargedStable);
462  }
463 
464  double getExpectedPhotonCount(const Particle* particle)
465  {
466  return TOPVariable::expectedPhotonCount(particle, 0);
467  }
468 
469  double getExpectedPhotonCountExpert(const Particle* particle, const vector<double>& vars)
470  {
471  if (vars.size() != 1) {
472  B2FATAL("Need exactly one parameter (pdg id).");
473  }
474  return TOPVariable::expectedPhotonCount(particle, static_cast<int>(vars[0]));
475  }
476 
477  double getEstimatedBkgCount(const Particle* particle)
478  {
479  const auto* topLikelihood = TOPVariable::getTOPLikelihood(particle);
480  return topLikelihood ? topLikelihood->getEstBkg() : std::numeric_limits<double>::quiet_NaN();
481  }
482 
483  double getElectronLogL(const Particle* particle)
484  {
485  const auto* topLikelihood = getTOPLikelihood(particle);
486  if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
487  return topLikelihood->getLogL_e();
488  }
489 
490  double getMuonLogL(const Particle* particle)
491  {
492  const auto* topLikelihood = getTOPLikelihood(particle);
493  if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
494  return topLikelihood->getLogL_mu();
495  }
496 
497  double getPionLogL(const Particle* particle)
498  {
499  const auto* topLikelihood = getTOPLikelihood(particle);
500  if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
501  return topLikelihood->getLogL_pi();
502  }
503 
504  double getKaonLogL(const Particle* particle)
505  {
506  const auto* topLikelihood = getTOPLikelihood(particle);
507  if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
508  return topLikelihood->getLogL_K();
509  }
510 
511  double getProtonLogL(const Particle* particle)
512  {
513  const auto* topLikelihood = getTOPLikelihood(particle);
514  if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
515  return topLikelihood->getLogL_p();
516  }
517 
518  double getDeuteronLogL(const Particle* particle)
519  {
520  const auto* topLikelihood = getTOPLikelihood(particle);
521  if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
522  return topLikelihood->getLogL(Const::deuteron);
523  }
524 
525  //---------------- TOPLikelihoodScanResult based --------------------
526 
527  double getLogLScanMass(const Particle* particle)
528  {
529  auto* scanRes = getTOPLikelihoodScanResult(particle);
530  return scanRes ? scanRes->getMostLikelyMass() : std::numeric_limits<double>::quiet_NaN();
531  }
532 
533  double getLogLScanMassUpperInterval(const Particle* particle)
534  {
535  auto* scanRes = getTOPLikelihoodScanResult(particle);
536  return scanRes ? scanRes->getMostLikelyMassIntervalUp() : std::numeric_limits<double>::quiet_NaN();
537  }
538 
539  double getLogLScanMassLowerInterval(const Particle* particle)
540  {
541  auto* scanRes = getTOPLikelihoodScanResult(particle);
542  return scanRes ? scanRes->getMostLikelyMassIntervalLow() : std::numeric_limits<double>::quiet_NaN();
543  }
544 
545  double getLogLScanThreshold(const Particle* particle)
546  {
547  auto* scanRes = getTOPLikelihoodScanResult(particle);
548  return scanRes ? scanRes->getThreshold() : std::numeric_limits<double>::quiet_NaN();
549  }
550 
551  double getLogLScanExpectedSignalPhotons(const Particle* particle)
552  {
553  auto* scanRes = getTOPLikelihoodScanResult(particle);
554  return scanRes ? scanRes->getMostLikelySignalPhotonCount() : std::numeric_limits<double>::quiet_NaN();
555  }
556 
557  //---------------- TOPRecBunch based --------------------
558 
559  double isTOPRecBunchReconstructed([[maybe_unused]] const Particle* particle)
560  {
561  StoreObjPtr<TOPRecBunch> recBunch;
562  if (not recBunch.isValid()) return 0;
563  return recBunch->isReconstructed();
564  }
565 
566  double TOPRecBunchNumber([[maybe_unused]] const Particle* particle)
567  {
568  StoreObjPtr<TOPRecBunch> recBunch;
569  if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
570  if (not recBunch->isReconstructed()) return std::numeric_limits<double>::quiet_NaN();
571  return recBunch->getBunchNo();
572  }
573 
574  double TOPRecBucketNumber([[maybe_unused]] const Particle* particle)
575  {
576  StoreObjPtr<TOPRecBunch> recBunch;
577  if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
578  auto bucket = recBunch->getBucketNumber();
579  if (bucket == TOPRecBunch::c_Unknown) return std::numeric_limits<double>::quiet_NaN();
580  return bucket;
581  }
582 
583  double isTOPRecBunchFilled([[maybe_unused]] const Particle* particle)
584  {
585  StoreObjPtr<TOPRecBunch> recBunch;
586  if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
587  return recBunch->getBucketFillStatus();
588  }
589 
590  double isTOPRecBunchNumberEQsim([[maybe_unused]] const Particle* particle)
591  {
592  StoreObjPtr<TOPRecBunch> recBunch;
593  if (not recBunch.isValid()) return 0;
594  if (not recBunch->isReconstructed()) return 0;
595  if (not recBunch->isSimulated()) return 0;
596  return (recBunch->getBunchNo() == recBunch->getMCBunchNo());
597  }
598 
599  double TOPRecBunchCurrentOffset([[maybe_unused]] const Particle* particle)
600  {
601  StoreObjPtr<TOPRecBunch> recBunch;
602  if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
603  if (not recBunch->isReconstructed()) return std::numeric_limits<double>::quiet_NaN();
604  return recBunch->getCurrentOffset();
605  }
606 
607  double TOPRecBunchTrackCount([[maybe_unused]] const Particle* particle)
608  {
609  StoreObjPtr<TOPRecBunch> recBunch;
610  if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
611  return recBunch->getNumTracks();
612  }
613 
614  double TOPRecBunchUsedTrackCount([[maybe_unused]] const Particle* particle)
615  {
616  StoreObjPtr<TOPRecBunch> recBunch;
617  if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
618  return recBunch->getUsedTracks();
619  }
620 
621  //-------------- Event based -----------------------------------
622 
623  double TOPRawPhotonsInSlot([[maybe_unused]] const Particle* particle, const vector<double>& vars)
624  {
625  if (vars.size() != 1) { B2FATAL("Need exactly one parameter (slot id).");}
626  StoreArray<TOPDigit> topDigits;
627  int slotID = static_cast<int>(vars[0]);
628  size_t count = 0;
629  for (const auto& t : topDigits) {
630  if (t.getModuleID() != slotID) continue;
631  count += 1;
632  }
633  return count;
634  }
635 
636  double TOPGoodPhotonsInSlot([[maybe_unused]] const Particle* particle, const vector<double>& vars)
637  {
638  if (vars.size() != 1) { B2FATAL("Need exactly one parameter (slot id).");}
639  StoreArray<TOPDigit> topDigits;
640  int slotID = static_cast<int>(vars[0]);
641  size_t count = 0;
642  for (const auto& t : topDigits) {
643  if (t.getModuleID() != slotID) continue;
644  if (t.getHitQuality() != TOPDigit::c_Good) continue;
645  count += 1;
646  }
647  return count;
648  }
649 
650  double TOPTracksInSlot(const Particle* particle)
651  {
652  int slotID = getSlotID(particle);
653  if (slotID == 0) return std::numeric_limits<double>::quiet_NaN();
654 
655  StoreArray<Track> tracks;
656  int nTracks = 0;
657  for (const auto& t : tracks) {
658  const auto* tl = t.getRelated<TOPLikelihood>();
659  if (not tl) continue;
660  const auto* te = tl->getRelated<ExtHit>();
661  if (not te) continue;
662  if (te->getCopyID() != slotID) continue;
663  nTracks += 1;
664  }
665  return nTracks;
666  }
667 
668  } // TOPVariable
669 
670 
671  VARIABLE_GROUP("TOP Variables (cdst needed)");
672  REGISTER_VARIABLE("topSlotID", TOPVariable::topSlotID,
673  "slot ID of the particle");
674  REGISTER_VARIABLE("topSlotIDMCMatch", TOPVariable::topSlotIDMCMatch,
675  "slot ID of the matched MC particle");
676 
677  REGISTER_VARIABLE("topLocalX", TOPVariable::getTOPLocalX,
678  "x coordinate of the particle entry point to TOP in the local frame");
679  REGISTER_VARIABLE("topLocalY", TOPVariable::getTOPLocalY,
680  "y coordinate of the particle entry point to TOP in the local frame");
681  REGISTER_VARIABLE("topLocalZ", TOPVariable::getTOPLocalZ,
682  "z coordinate of the particle entry point to TOP in the local frame");
683  REGISTER_VARIABLE("topLocalXMCMatch", TOPVariable::getTOPLocalXMCMatch,
684  "x coordinate of the matched MC particle entry point to TOP in the local frame");
685  REGISTER_VARIABLE("topLocalYMCMatch", TOPVariable::getTOPLocalYMCMatch,
686  "y coordinate of the matched MC particle entry point to TOP in the local frame");
687  REGISTER_VARIABLE("topLocalZMCMatch", TOPVariable::getTOPLocalZMCMatch,
688  "z coordinate of the matched MC particle entry point to TOP in the local frame");
689  REGISTER_VARIABLE("topLocalPhi", TOPVariable::getTOPLocalPhi,
690  "momentum azimuthal angle of the particle at entry point to TOP in the local frame");
691  REGISTER_VARIABLE("topLocalTheta", TOPVariable::getTOPLocalTheta,
692  "momentum polar angle of the particle at entry point to the TOP in the local frame");
693  REGISTER_VARIABLE("topLocalPhiMCMatch", TOPVariable::getTOPLocalPhiMCMatch,
694  "momentum azimuthal angle of the matched MC particle at entry point to TOP in the local frame");
695  REGISTER_VARIABLE("topLocalThetaMCMatch", TOPVariable::getTOPLocalThetaMCMatch,
696  "momentum polar angle of the matched MC particle at entry point to TOP in the local frame");
697  REGISTER_VARIABLE("topTOF", TOPVariable::getTOF,
698  "time-of-flight of the particle from the origin to TOP");
699  REGISTER_VARIABLE("topTOFMCMatch", TOPVariable::getTOFMCMatch,
700  "time-of-flight of the matched MC particle from the origin to TOP");
701  REGISTER_VARIABLE("topTOFExpert(pdg)", TOPVariable::getTOFExpert,
702  "time-of-flight from the origin to TOP for a given PDG code");
703 
704  REGISTER_VARIABLE("extrapTrackToTOPimpactZ", TOPVariable::extrapTrackToTOPz,
705  "z coordinate of the track extrapolated to R = 120 cm using helix from TrackFitResult");
706  REGISTER_VARIABLE("extrapTrackToTOPimpactTheta", TOPVariable::extrapTrackToTOPtheta,
707  "theta coordinate of the track extrapolated to R = 120 cm using helix from TrackFitResult");
708  REGISTER_VARIABLE("extrapTrackToTOPimpactPhi", TOPVariable::extrapTrackToTOPphi,
709  "phi coordinate of the track extrapolated to R = 120 cm using helix from TrackFitResult");
710 
711  REGISTER_VARIABLE("topDigitCount", TOPVariable::topDigitCount,
712  "number of good digits in the same module as particle");
713  REGISTER_VARIABLE("topDigitCountSignal", TOPVariable::topSignalDigitCount,
714  "number of good, background-subtracted digits in interval [0, 50 ns] of the same module as particle");
715  REGISTER_VARIABLE("topDigitCountBkg", TOPVariable::topBackgroundDigitCount,
716  "number of good digits in interval [-50 ns, 0] of the same module as particle");
717  REGISTER_VARIABLE("topDigitCountRaw", TOPVariable::topRawDigitCount,
718  "number of all digits (regardless of hit quality) in the same module as particle");
719  REGISTER_VARIABLE("topDigitCountInterval(tmin, tmax)", TOPVariable::countTOPHitsInInterval,
720  "number of good digits in a given time interval of the same module as particle");
721  REGISTER_VARIABLE("topDigitCountIntervalRaw(tmin, tmax)", TOPVariable::countRawTOPHitsInInterval,
722  "number of all digits (regardless of hit quality) in a given time interval of the same module as particle");
723  REGISTER_VARIABLE("topDigitCountMCMatch", TOPVariable::topDigitCountMCMatch,
724  "number of good digits associated with the matched MC particle and in the same module as particle");
725  REGISTER_VARIABLE("topDigitCountIntervalMCMatch(tmin, tmax)", TOPVariable::countTOPHitsInIntervalMCMatch,
726  "number of good digits associated with the matched MC particle in a given time interval and in the same module as particle");
727 
728  REGISTER_VARIABLE("topLogLFlag", TOPVariable::getFlag,
729  "reconstruction flag: 1 if log likelihoods available, 0 otherwise");
730  REGISTER_VARIABLE("topLogLPhotonCount", TOPVariable::getTOPPhotonCount,
731  "number of photons used for log likelihood calculation");
732  REGISTER_VARIABLE("topLogLExpectedPhotonCount", TOPVariable::getExpectedPhotonCount,
733  "expected number of photons (including bkg) for this particle");
734  REGISTER_VARIABLE("topLogLExpectedPhotonCountExpert(pdg)", TOPVariable::getExpectedPhotonCountExpert,
735  "expected number of photons (including bkg) for a given PDG code");
736  REGISTER_VARIABLE("topLogLEstimatedBkgCount", TOPVariable::getEstimatedBkgCount,
737  "estimated number of background photons");
738  REGISTER_VARIABLE("topLogLElectron", TOPVariable::getElectronLogL,
739  "electron log likelihood");
740  REGISTER_VARIABLE("topLogLMuon", TOPVariable::getMuonLogL,
741  "muon log likelihood");
742  REGISTER_VARIABLE("topLogLPion", TOPVariable::getPionLogL,
743  "pion log likelihood");
744  REGISTER_VARIABLE("topLogLKaon", TOPVariable::getKaonLogL,
745  "kaon log likelihood");
746  REGISTER_VARIABLE("topLogLProton", TOPVariable::getProtonLogL,
747  "proton log likelihood");
748  REGISTER_VARIABLE("topLogLDeuteron", TOPVariable::getDeuteronLogL,
749  "deuteron log likelihood");
750 
751  REGISTER_VARIABLE("logLScanMass", TOPVariable::getLogLScanMass,
752  "mass at the logL maximum from the LL scan. Requires TOPLLScanner in the processing path.");
753  REGISTER_VARIABLE("logLScanMassUpperInterval", TOPVariable::getLogLScanMassUpperInterval,
754  "Upper edge of the mass interval determined by the LL scan. Requires TOPLLScanner in the processing path.");
755  REGISTER_VARIABLE("logLScanMassLowerInterval", TOPVariable::getLogLScanMassLowerInterval,
756  "Lower edge of the mass interval determined by the LL scan. Requires TOPLLScanner in the processing path.");
757  REGISTER_VARIABLE("logLScanThreshold", TOPVariable::getLogLScanThreshold,
758  "Cherenkov threshold determind by the LL scan. Requires TOPLLScanner in the processing path.");
759  REGISTER_VARIABLE("logLScanExpectedSignalPhotons", TOPVariable::getLogLScanExpectedSignalPhotons,
760  "Expected signal photon yield at the LL maximum. Requires TOPLLScanner in the processing path.");
761 
762  REGISTER_VARIABLE("topBunchIsReconstructed", TOPVariable::isTOPRecBunchReconstructed,
763  "reconstruction flag: 1 if reconstructed, 0 otherwise");
764  REGISTER_VARIABLE("topBunchIsFilled", TOPVariable::isTOPRecBunchFilled,
765  "bunch fill status: 0 empty, 1 filled, -1 unknown");
766  REGISTER_VARIABLE("topBunchNumber", TOPVariable::TOPRecBunchNumber,
767  "reconstructed bunch number relative to L1 trigger");
768  REGISTER_VARIABLE("topBucketNumber", TOPVariable::TOPRecBucketNumber,
769  "reconstructed bucket number within the ring");
770  REGISTER_VARIABLE("topBunchMCMatch", TOPVariable::isTOPRecBunchNumberEQsim,
771  "MC matching status: 1 if reconstructed bunch equal to simulated bunch, 0 otherwise");
772  REGISTER_VARIABLE("topBunchOffset", TOPVariable::TOPRecBunchCurrentOffset,
773  "current offset to the reconstructed bunch crossing time");
774  REGISTER_VARIABLE("topBunchTrackCount", TOPVariable::TOPRecBunchTrackCount,
775  "number of tracks selected for the bunch reconstruction");
776  REGISTER_VARIABLE("topBunchUsedTrackCount", TOPVariable::TOPRecBunchUsedTrackCount,
777  "number of tracks actually used in the bunch reconstruction");
778 
779  REGISTER_VARIABLE("topRawPhotonsInSlot(id)", TOPVariable::TOPRawPhotonsInSlot,
780  "number of all photons in the given slot id");
781  REGISTER_VARIABLE("topGoodPhotonsInSlot(id)", TOPVariable::TOPGoodPhotonsInSlot,
782  "number of good photons in the given slot id");
783  REGISTER_VARIABLE("topTracksInSlot", TOPVariable::TOPTracksInSlot,
784  "number of tracks in the same slot as particle");
785  } // Variable
787 } // Belle2
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17