Belle II Software development
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
31using namespace std;
32using namespace ROOT::Math;
33
34namespace 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 getModuleID(const Particle* particle)
447 {
448 const auto* topLikelihood = TOPVariable::getTOPLikelihood(particle);
449 return topLikelihood ? topLikelihood->getModuleID() : std::numeric_limits<double>::quiet_NaN();
450 }
451
452 double getEmissionX(const Particle* particle)
453 {
454 const auto* topLikelihood = TOPVariable::getTOPLikelihood(particle);
455 return topLikelihood ? topLikelihood->getX() : std::numeric_limits<double>::quiet_NaN();
456 }
457
458 double getEmissionZ(const Particle* particle)
459 {
460 const auto* topLikelihood = TOPVariable::getTOPLikelihood(particle);
461 return topLikelihood ? topLikelihood->getZ() : std::numeric_limits<double>::quiet_NaN();
462 }
463
464 double getTOPPhotonCount(const Particle* particle)
465 {
466 const auto* topLikelihood = TOPVariable::getTOPLikelihood(particle);
467 return topLikelihood ? topLikelihood->getNphot() : std::numeric_limits<double>::quiet_NaN();
468 }
469
470 double expectedPhotonCount(const Particle* particle, int pdg)
471 {
472 const auto* topLikelihood = getTOPLikelihood(particle);
473 if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
474
475 pdg = pdg != 0 ? pdg : particle->getPDGCode();
476 const auto& chargedStable = Const::chargedStableSet.find(abs(pdg));
477 if (chargedStable == Const::invalidParticle) return std::numeric_limits<double>::quiet_NaN();
478
479 return topLikelihood->getEstPhot(chargedStable);
480 }
481
482 double getExpectedPhotonCount(const Particle* particle)
483 {
484 return TOPVariable::expectedPhotonCount(particle, 0);
485 }
486
487 double getExpectedPhotonCountExpert(const Particle* particle, const vector<double>& vars)
488 {
489 if (vars.size() != 1) {
490 B2FATAL("Need exactly one parameter (pdg id).");
491 }
492 return TOPVariable::expectedPhotonCount(particle, static_cast<int>(vars[0]));
493 }
494
495 double effectiveSignalYield(const Particle* particle, int pdg)
496 {
497 const auto* topLikelihood = getTOPLikelihood(particle);
498 if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
499
500 pdg = pdg != 0 ? pdg : particle->getPDGCode();
501 const auto& chargedStable = Const::chargedStableSet.find(abs(pdg));
502 if (chargedStable == Const::invalidParticle) return std::numeric_limits<double>::quiet_NaN();
503
504 return topLikelihood->getEffectiveSignalYield(chargedStable);
505 }
506
507 double getEffectiveSignalYield(const Particle* particle)
508 {
509 return TOPVariable::effectiveSignalYield(particle, 0);
510 }
511
512 double getEffectiveSignalYieldExpert(const Particle* particle, const vector<double>& vars)
513 {
514 if (vars.size() != 1) {
515 B2FATAL("Need exactly one parameter (pdg id).");
516 }
517 return TOPVariable::effectiveSignalYield(particle, static_cast<int>(vars[0]));
518 }
519
520 double getEstimatedBkgCount(const Particle* particle)
521 {
522 const auto* topLikelihood = TOPVariable::getTOPLikelihood(particle);
523 return topLikelihood ? topLikelihood->getEstBkg() : std::numeric_limits<double>::quiet_NaN();
524 }
525
526 double getElectronLogL(const Particle* particle)
527 {
528 const auto* topLikelihood = getTOPLikelihood(particle);
529 if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
530 return topLikelihood->getLogL_e();
531 }
532
533 double getMuonLogL(const Particle* particle)
534 {
535 const auto* topLikelihood = getTOPLikelihood(particle);
536 if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
537 return topLikelihood->getLogL_mu();
538 }
539
540 double getPionLogL(const Particle* particle)
541 {
542 const auto* topLikelihood = getTOPLikelihood(particle);
543 if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
544 return topLikelihood->getLogL_pi();
545 }
546
547 double getKaonLogL(const Particle* particle)
548 {
549 const auto* topLikelihood = getTOPLikelihood(particle);
550 if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
551 return topLikelihood->getLogL_K();
552 }
553
554 double getProtonLogL(const Particle* particle)
555 {
556 const auto* topLikelihood = getTOPLikelihood(particle);
557 if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
558 return topLikelihood->getLogL_p();
559 }
560
561 double getDeuteronLogL(const Particle* particle)
562 {
563 const auto* topLikelihood = getTOPLikelihood(particle);
564 if (not topLikelihood) return std::numeric_limits<double>::quiet_NaN();
565 return topLikelihood->getLogL(Const::deuteron);
566 }
567
568 //---------------- TOPLikelihoodScanResult based --------------------
569
570 double getLogLScanMass(const Particle* particle)
571 {
572 auto* scanRes = getTOPLikelihoodScanResult(particle);
573 return scanRes ? scanRes->getMostLikelyMass() : std::numeric_limits<double>::quiet_NaN();
574 }
575
576 double getLogLScanMassUpperInterval(const Particle* particle)
577 {
578 auto* scanRes = getTOPLikelihoodScanResult(particle);
579 return scanRes ? scanRes->getMostLikelyMassIntervalUp() : std::numeric_limits<double>::quiet_NaN();
580 }
581
582 double getLogLScanMassLowerInterval(const Particle* particle)
583 {
584 auto* scanRes = getTOPLikelihoodScanResult(particle);
585 return scanRes ? scanRes->getMostLikelyMassIntervalLow() : std::numeric_limits<double>::quiet_NaN();
586 }
587
588 double getLogLScanThreshold(const Particle* particle)
589 {
590 auto* scanRes = getTOPLikelihoodScanResult(particle);
591 return scanRes ? scanRes->getThreshold() : std::numeric_limits<double>::quiet_NaN();
592 }
593
594 double getLogLScanExpectedSignalPhotons(const Particle* particle)
595 {
596 auto* scanRes = getTOPLikelihoodScanResult(particle);
597 return scanRes ? scanRes->getMostLikelySignalPhotonCount() : std::numeric_limits<double>::quiet_NaN();
598 }
599
600 //---------------- TOPRecBunch based --------------------
601
602 double isTOPRecBunchReconstructed([[maybe_unused]] const Particle* particle)
603 {
604 StoreObjPtr<TOPRecBunch> recBunch;
605 if (not recBunch.isValid()) return 0;
606 return recBunch->isReconstructed();
607 }
608
609 double TOPRecBunchNumber([[maybe_unused]] const Particle* particle)
610 {
611 StoreObjPtr<TOPRecBunch> recBunch;
612 if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
613 if (not recBunch->isReconstructed()) return std::numeric_limits<double>::quiet_NaN();
614 return recBunch->getBunchNo();
615 }
616
617 double TOPRecBucketNumber([[maybe_unused]] const Particle* particle)
618 {
619 StoreObjPtr<TOPRecBunch> recBunch;
620 if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
621 auto bucket = recBunch->getBucketNumber();
622 if (bucket == TOPRecBunch::c_Unknown) return std::numeric_limits<double>::quiet_NaN();
623 return bucket;
624 }
625
626 double isTOPRecBunchFilled([[maybe_unused]] const Particle* particle)
627 {
628 StoreObjPtr<TOPRecBunch> recBunch;
629 if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
630 return recBunch->getBucketFillStatus();
631 }
632
633 double isTOPRecBunchNumberEQsim([[maybe_unused]] const Particle* particle)
634 {
635 StoreObjPtr<TOPRecBunch> recBunch;
636 if (not recBunch.isValid()) return 0;
637 if (not recBunch->isReconstructed()) return 0;
638 if (not recBunch->isSimulated()) return 0;
639 return (recBunch->getBunchNo() == recBunch->getMCBunchNo());
640 }
641
642 double TOPRecBunchCurrentOffset([[maybe_unused]] const Particle* particle)
643 {
644 StoreObjPtr<TOPRecBunch> recBunch;
645 if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
646 if (not recBunch->isReconstructed()) return std::numeric_limits<double>::quiet_NaN();
647 return recBunch->getCurrentOffset();
648 }
649
650 double TOPRecBunchTrackCount([[maybe_unused]] const Particle* particle)
651 {
652 StoreObjPtr<TOPRecBunch> recBunch;
653 if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
654 return recBunch->getNumTracks();
655 }
656
657 double TOPRecBunchUsedTrackCount([[maybe_unused]] const Particle* particle)
658 {
659 StoreObjPtr<TOPRecBunch> recBunch;
660 if (not recBunch.isValid()) return std::numeric_limits<double>::quiet_NaN();
661 return recBunch->getUsedTracks();
662 }
663
664 //-------------- Event based -----------------------------------
665
666 double TOPRawPhotonsInSlot([[maybe_unused]] const Particle* particle, const vector<double>& vars)
667 {
668 if (vars.size() != 1) { B2FATAL("Need exactly one parameter (slot id).");}
669 StoreArray<TOPDigit> topDigits;
670 int slotID = static_cast<int>(vars[0]);
671 size_t count = 0;
672 for (const auto& t : topDigits) {
673 if (t.getModuleID() != slotID) continue;
674 count += 1;
675 }
676 return count;
677 }
678
679 double TOPGoodPhotonsInSlot([[maybe_unused]] const Particle* particle, const vector<double>& vars)
680 {
681 if (vars.size() != 1) { B2FATAL("Need exactly one parameter (slot id).");}
682 StoreArray<TOPDigit> topDigits;
683 int slotID = static_cast<int>(vars[0]);
684 size_t count = 0;
685 for (const auto& t : topDigits) {
686 if (t.getModuleID() != slotID) continue;
687 if (t.getHitQuality() != TOPDigit::c_Good) continue;
688 count += 1;
689 }
690 return count;
691 }
692
693 double TOPTracksInSlot(const Particle* particle)
694 {
695 int slotID = getSlotID(particle);
696 if (slotID == 0) return std::numeric_limits<double>::quiet_NaN();
697
698 StoreArray<Track> tracks;
699 int nTracks = 0;
700 for (const auto& t : tracks) {
701 const auto* tl = t.getRelated<TOPLikelihood>();
702 if (not tl) continue;
703 const auto* te = tl->getRelated<ExtHit>();
704 if (not te) continue;
705 if (te->getCopyID() != slotID) continue;
706 nTracks += 1;
707 }
708 return nTracks;
709 }
710
711 } // TOPVariable
712
713
714 VARIABLE_GROUP("TOP Variables (cdst needed)");
715 REGISTER_VARIABLE("topSlotID", TOPVariable::topSlotID,
716 "slot ID of the particle");
717 REGISTER_VARIABLE("topSlotIDMCMatch", TOPVariable::topSlotIDMCMatch,
718 "slot ID of the matched MC particle");
719
720 REGISTER_VARIABLE("topLocalX", TOPVariable::getTOPLocalX,
721 "x coordinate of the particle entry point to TOP in the local frame");
722 REGISTER_VARIABLE("topLocalY", TOPVariable::getTOPLocalY,
723 "y coordinate of the particle entry point to TOP in the local frame");
724 REGISTER_VARIABLE("topLocalZ", TOPVariable::getTOPLocalZ,
725 "z coordinate of the particle entry point to TOP in the local frame");
726 REGISTER_VARIABLE("topLocalXMCMatch", TOPVariable::getTOPLocalXMCMatch,
727 "x coordinate of the matched MC particle entry point to TOP in the local frame");
728 REGISTER_VARIABLE("topLocalYMCMatch", TOPVariable::getTOPLocalYMCMatch,
729 "y coordinate of the matched MC particle entry point to TOP in the local frame");
730 REGISTER_VARIABLE("topLocalZMCMatch", TOPVariable::getTOPLocalZMCMatch,
731 "z coordinate of the matched MC particle entry point to TOP in the local frame");
732 REGISTER_VARIABLE("topLocalPhi", TOPVariable::getTOPLocalPhi,
733 "momentum azimuthal angle of the particle at entry point to TOP in the local frame");
734 REGISTER_VARIABLE("topLocalTheta", TOPVariable::getTOPLocalTheta,
735 "momentum polar angle of the particle at entry point to the TOP in the local frame");
736 REGISTER_VARIABLE("topLocalPhiMCMatch", TOPVariable::getTOPLocalPhiMCMatch,
737 "momentum azimuthal angle of the matched MC particle at entry point to TOP in the local frame");
738 REGISTER_VARIABLE("topLocalThetaMCMatch", TOPVariable::getTOPLocalThetaMCMatch,
739 "momentum polar angle of the matched MC particle at entry point to TOP in the local frame");
740 REGISTER_VARIABLE("topTOF", TOPVariable::getTOF,
741 "time-of-flight of the particle from the origin to TOP");
742 REGISTER_VARIABLE("topTOFMCMatch", TOPVariable::getTOFMCMatch,
743 "time-of-flight of the matched MC particle from the origin to TOP");
744 REGISTER_VARIABLE("topTOFExpert(pdg)", TOPVariable::getTOFExpert,
745 "time-of-flight from the origin to TOP for a given PDG code");
746
747 REGISTER_VARIABLE("extrapTrackToTOPimpactZ", TOPVariable::extrapTrackToTOPz,
748 "z coordinate of the track extrapolated to R = 120 cm using helix from TrackFitResult");
749 REGISTER_VARIABLE("extrapTrackToTOPimpactTheta", TOPVariable::extrapTrackToTOPtheta,
750 "theta coordinate of the track extrapolated to R = 120 cm using helix from TrackFitResult");
751 REGISTER_VARIABLE("extrapTrackToTOPimpactPhi", TOPVariable::extrapTrackToTOPphi,
752 "phi coordinate of the track extrapolated to R = 120 cm using helix from TrackFitResult");
753
754 REGISTER_VARIABLE("topDigitCount", TOPVariable::topDigitCount,
755 "number of good digits in the same module as particle");
756 REGISTER_VARIABLE("topDigitCountSignal", TOPVariable::topSignalDigitCount,
757 "number of good, background-subtracted digits in interval [0, 50 ns] of the same module as particle");
758 REGISTER_VARIABLE("topDigitCountBkg", TOPVariable::topBackgroundDigitCount,
759 "number of good digits in interval [-50 ns, 0] of the same module as particle");
760 REGISTER_VARIABLE("topDigitCountRaw", TOPVariable::topRawDigitCount,
761 "number of all digits (regardless of hit quality) in the same module as particle");
762 REGISTER_VARIABLE("topDigitCountInterval(tmin, tmax)", TOPVariable::countTOPHitsInInterval,
763 "number of good digits in a given time interval of the same module as particle");
764 REGISTER_VARIABLE("topDigitCountIntervalRaw(tmin, tmax)", TOPVariable::countRawTOPHitsInInterval,
765 "number of all digits (regardless of hit quality) in a given time interval of the same module as particle");
766 REGISTER_VARIABLE("topDigitCountMCMatch", TOPVariable::topDigitCountMCMatch,
767 "number of good digits associated with the matched MC particle and in the same module as particle");
768 REGISTER_VARIABLE("topDigitCountIntervalMCMatch(tmin, tmax)", TOPVariable::countTOPHitsInIntervalMCMatch,
769 "number of good digits associated with the matched MC particle in a given time interval and in the same module as particle");
770
771 REGISTER_VARIABLE("topLogLFlag", TOPVariable::getFlag,
772 "reconstruction flag: 1 if log likelihoods available, 0 otherwise");
773 REGISTER_VARIABLE("topLogLSlotID", TOPVariable::getModuleID,
774 "slot ID of the particle (from TOPLikelihoods)");
775 REGISTER_VARIABLE("topLogLEmiX", TOPVariable::getEmissionX,
776 "photon emission point assumed in PDF construction: coordinate x in local frame");
777 REGISTER_VARIABLE("topLogLEmiZ", TOPVariable::getEmissionZ,
778 "photon emission point assumed in PDF construction: coordinate z in local frame");
779 REGISTER_VARIABLE("topLogLPhotonCount", TOPVariable::getTOPPhotonCount,
780 "number of photons used for log likelihood calculation");
781 REGISTER_VARIABLE("topLogLExpectedPhotonCount", TOPVariable::getExpectedPhotonCount,
782 "expected number of photons (including bkg) for this particle");
783 REGISTER_VARIABLE("topLogLExpectedPhotonCountExpert(pdg)", TOPVariable::getExpectedPhotonCountExpert,
784 "expected number of photons (including bkg) for a given PDG code");
785 REGISTER_VARIABLE("topLogLEstimatedBkgCount", TOPVariable::getEstimatedBkgCount,
786 "estimated number of background photons");
787 REGISTER_VARIABLE("topLogLSignalYield", TOPVariable::getEffectiveSignalYield,
788 "effective number of signal photons for this particle (sum of sPlot weights)");
789 REGISTER_VARIABLE("topLogLSignalYieldExpert(pdg)", TOPVariable::getEffectiveSignalYieldExpert,
790 "effective number of signal photons for a given PDG code (sum of sPlot weights)");
791 REGISTER_VARIABLE("topLogLElectron", TOPVariable::getElectronLogL,
792 "electron log likelihood");
793 REGISTER_VARIABLE("topLogLMuon", TOPVariable::getMuonLogL,
794 "muon log likelihood");
795 REGISTER_VARIABLE("topLogLPion", TOPVariable::getPionLogL,
796 "pion log likelihood");
797 REGISTER_VARIABLE("topLogLKaon", TOPVariable::getKaonLogL,
798 "kaon log likelihood");
799 REGISTER_VARIABLE("topLogLProton", TOPVariable::getProtonLogL,
800 "proton log likelihood");
801 REGISTER_VARIABLE("topLogLDeuteron", TOPVariable::getDeuteronLogL,
802 "deuteron log likelihood");
803
804 REGISTER_VARIABLE("logLScanMass", TOPVariable::getLogLScanMass,
805 "mass at the logL maximum from the LL scan. Requires TOPLLScanner in the processing path.");
806 REGISTER_VARIABLE("logLScanMassUpperInterval", TOPVariable::getLogLScanMassUpperInterval,
807 "Upper edge of the mass interval determined by the LL scan. Requires TOPLLScanner in the processing path.");
808 REGISTER_VARIABLE("logLScanMassLowerInterval", TOPVariable::getLogLScanMassLowerInterval,
809 "Lower edge of the mass interval determined by the LL scan. Requires TOPLLScanner in the processing path.");
810 REGISTER_VARIABLE("logLScanThreshold", TOPVariable::getLogLScanThreshold,
811 "Cherenkov threshold determind by the LL scan. Requires TOPLLScanner in the processing path.");
812 REGISTER_VARIABLE("logLScanExpectedSignalPhotons", TOPVariable::getLogLScanExpectedSignalPhotons,
813 "Expected signal photon yield at the LL maximum. Requires TOPLLScanner in the processing path.");
814
815 REGISTER_VARIABLE("topBunchIsReconstructed", TOPVariable::isTOPRecBunchReconstructed,
816 "reconstruction flag: 1 if reconstructed, 0 otherwise");
817 REGISTER_VARIABLE("topBunchIsFilled", TOPVariable::isTOPRecBunchFilled,
818 "bunch fill status: 0 empty, 1 filled, -1 unknown");
819 REGISTER_VARIABLE("topBunchNumber", TOPVariable::TOPRecBunchNumber,
820 "reconstructed bunch number relative to L1 trigger");
821 REGISTER_VARIABLE("topBucketNumber", TOPVariable::TOPRecBucketNumber,
822 "reconstructed bucket number within the ring");
823 REGISTER_VARIABLE("topBunchMCMatch", TOPVariable::isTOPRecBunchNumberEQsim,
824 "MC matching status: 1 if reconstructed bunch equal to simulated bunch, 0 otherwise");
825 REGISTER_VARIABLE("topBunchOffset", TOPVariable::TOPRecBunchCurrentOffset,
826 "current offset to the reconstructed bunch crossing time");
827 REGISTER_VARIABLE("topBunchTrackCount", TOPVariable::TOPRecBunchTrackCount,
828 "number of tracks selected for the bunch reconstruction");
829 REGISTER_VARIABLE("topBunchUsedTrackCount", TOPVariable::TOPRecBunchUsedTrackCount,
830 "number of tracks actually used in the bunch reconstruction");
831
832 REGISTER_VARIABLE("topRawPhotonsInSlot(id)", TOPVariable::TOPRawPhotonsInSlot,
833 "number of all photons in the given slot id");
834 REGISTER_VARIABLE("topGoodPhotonsInSlot(id)", TOPVariable::TOPGoodPhotonsInSlot,
835 "number of good photons in the given slot id");
836 REGISTER_VARIABLE("topTracksInSlot", TOPVariable::TOPTracksInSlot,
837 "number of tracks in the same slot as particle");
838 } // Variable
840} // Belle2
Belle2::Const::ParticleSet::find
const ParticleType & find(int pdg) const
Returns particle in set with given PDG code, or invalidParticle if not found.
Definition: Const.h:571
Belle2::Const::chargedStableSet
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:618
Belle2::Const::invalidParticle
static const ParticleType invalidParticle
Invalid particle, used internally.
Definition: Const.h:681
Belle2::Const::deuteron
static const ChargedStable deuteron
deuteron particle
Definition: Const.h:664
Belle2::TOPRecBunch::c_Unknown
@ c_Unknown
not known
Definition: TOPRecBunch.h:33
Belle2::TOP::TOPGeometryPar::getGeometry
const TOPGeometry * getGeometry() const
Returns pointer to geometry object using basf2 units.
Definition: TOPGeometryPar.cc:165
Belle2::TOP::TOPGeometryPar::Instance
static TOPGeometryPar * Instance()
Static method to obtain the pointer to its instance.
Definition: TOPGeometryPar.cc:43
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
std
STL namespace.