Belle II Software development
MCTruthVariables.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9// Own header.
10#include <analysis/variables/MCTruthVariables.h>
11
12// include VariableManager
13#include <analysis/VariableManager/Manager.h>
14
15#include <analysis/dataobjects/Particle.h>
16#include <analysis/dataobjects/TauPairDecay.h>
17#include <analysis/utility/MCMatching.h>
18#include <analysis/utility/ReferenceFrame.h>
19#include <analysis/utility/ValueIndexPairSorting.h>
20
21#include <mdst/dataobjects/MCParticle.h>
22#include <mdst/dataobjects/ECLCluster.h>
23#include <mdst/dataobjects/Track.h>
24
25
26#include <framework/datastore/StoreArray.h>
27#include <framework/datastore/StoreObjPtr.h>
28#include <framework/dataobjects/EventMetaData.h>
29#include <framework/gearbox/Const.h>
30#include <framework/logging/Logger.h>
31#include <framework/database/DBObjPtr.h>
32#include <framework/dbobjects/BeamParameters.h>
33
34#include <queue>
35
36namespace Belle2 {
41 namespace Variable {
42
43 double isSignal(const Particle* part)
44 {
45 const MCParticle* mcparticle = part->getMCParticle();
46 if (!mcparticle) return Const::doubleNaN;
47
48 int status = MCMatching::getMCErrors(part, mcparticle);
49 return (status == MCMatching::c_Correct);
50 }
51
52 double isSignalAcceptWrongFSPs(const Particle* part)
53 {
54 const MCParticle* mcparticle = part->getMCParticle();
55 if (!mcparticle) return Const::doubleNaN;
56
57 int status = MCMatching::getMCErrors(part, mcparticle);
58 //remove the following bits
59 status &= (~MCMatching::c_MisID);
60 status &= (~MCMatching::c_AddedWrongParticle);
61
62 return (status == MCMatching::c_Correct);
63 }
64
65 double isPrimarySignal(const Particle* part)
66 {
67 return (isSignal(part) > 0.5 and particleMCPrimaryParticle(part) > 0.5);
68 }
69
70 double isMisidentified(const Particle* part)
71 {
72 const MCParticle* mcp = part->getMCParticle();
73 if (!mcp) return Const::doubleNaN;
74 int st = MCMatching::getMCErrors(part, mcp);
75 return ((st & MCMatching::c_MisID) != 0);
76 }
77
78 double isWrongCharge(const Particle* part)
79 {
80 const MCParticle* mcp = part->getMCParticle();
81 if (!mcp) return Const::doubleNaN;
82 return (part->getCharge() != mcp->getCharge());
83 }
84
85 double isCloneTrack(const Particle* particle)
86 {
87 // neutrals and composites don't make sense
88 if (!Const::chargedStableSet.contains(Const::ParticleType(abs(particle->getPDGCode()))))
89 return Const::doubleNaN;
90 // get mcparticle weight (mcmatch weight)
91 const auto mcpww = particle->getRelatedToWithWeight<MCParticle>();
92 if (!mcpww.first) return Const::doubleNaN;
93 return (mcpww.second < 0);
94 }
95
96 double isOrHasCloneTrack(const Particle* particle)
97 {
98 // use std::queue to check daughters-- granddaughters etc recursively
99 std::queue<const Particle*> qq;
100 qq.push(particle);
101 while (!qq.empty()) {
102 const auto d = qq.front(); // get daughter
103 qq.pop(); // remove the daughter from the queue
104 if (isCloneTrack(d) == 1.0) return 1.0;
105 size_t nDau = d->getNDaughters(); // number of daughters of daughters
106 for (size_t iDau = 0; iDau < nDau; ++iDau)
107 qq.push(d->getDaughter(iDau));
108 }
109 return 0.0;
110 }
111
112 double genNthMotherPDG(const Particle* part, const std::vector<double>& args)
113 {
114 const MCParticle* mcparticle = part->getMCParticle();
115 if (!mcparticle) return 0.0;
116
117 unsigned int nLevels = args.empty() ? 0 : args[0];
118
119 const MCParticle* curMCParticle = mcparticle;
120 for (unsigned int i = 0; i <= nLevels; ++i) {
121 const MCParticle* curMCMother = curMCParticle->getMother();
122 if (!curMCMother) return 0.0;
123 curMCParticle = curMCMother;
124 }
125 return curMCParticle->getPDG();
126 }
127
128 double genNthMotherIndex(const Particle* part, const std::vector<double>& args)
129 {
130 const MCParticle* mcparticle = part->getMCParticle();
131 if (!mcparticle) return 0.0;
132
133 unsigned int nLevels = args.empty() ? 0 : args[0];
134
135 const MCParticle* curMCParticle = mcparticle;
136 for (unsigned int i = 0; i <= nLevels; ++i) {
137 const MCParticle* curMCMother = curMCParticle->getMother();
138 if (!curMCMother) return 0.0;
139 curMCParticle = curMCMother;
140 }
141 return curMCParticle->getArrayIndex();
142 }
143
144 double genQ2PmPd(const Particle* part, const std::vector<double>& daughter_indices)
145 {
146 const MCParticle* mcparticle = part->getMCParticle();
147 if (!mcparticle) return Const::doubleNaN;
148
149 auto daughters = mcparticle->getDaughters();
150
151 ROOT::Math::PxPyPzEVector p4Daughters;
152 for (auto& double_daughter : daughter_indices) {
153 unsigned long daughter = std::lround(double_daughter);
154 if (daughter >= daughters.size()) return Const::doubleNaN;
155
156 p4Daughters += daughters[daughter]->get4Vector();
157 }
158 auto p4Mother = mcparticle->get4Vector();
159 return (p4Mother - p4Daughters).mag2();
160 }
161
162 double genMotherPDG(const Particle* part)
163 {
164 return genNthMotherPDG(part, {});
165 }
166
167 double genMotherP(const Particle* part)
168 {
169 const MCParticle* mcparticle = part->getMCParticle();
170 if (!mcparticle) return Const::doubleNaN;
171
172 const MCParticle* mcmother = mcparticle->getMother();
173 if (!mcmother) return Const::doubleNaN;
174
175 return mcmother->getMomentum().R();
176 }
177
178 double genMotherIndex(const Particle* part)
179 {
180 return genNthMotherIndex(part, {});
181 }
182
183 double genParticleIndex(const Particle* part)
184 {
185 const MCParticle* mcparticle = part->getMCParticle();
186 if (!mcparticle) return Const::doubleNaN;
187 return mcparticle->getArrayIndex();
188 }
189
190 double isSignalAcceptMissingNeutrino(const Particle* part)
191 {
192 const MCParticle* mcparticle = part->getMCParticle();
193 if (!mcparticle) return Const::doubleNaN;
194
195 int status = MCMatching::getMCErrors(part, mcparticle);
196 //remove the following bits
197 status &= (~MCMatching::c_MissNeutrino);
198
199 return (status == MCMatching::c_Correct);
200 }
201
202 double isSignalAcceptMissingMassive(const Particle* part)
203 {
204 const MCParticle* mcparticle = part->getMCParticle();
205 if (!mcparticle) return Const::doubleNaN;
206
207 int status = MCMatching::getMCErrors(part, mcparticle);
208 //remove the following bits
209 status &= (~MCMatching::c_MissMassiveParticle);
210 status &= (~MCMatching::c_MissKlong);
211
212 return (status == MCMatching::c_Correct);
213 }
214
215 double isSignalAcceptMissingGamma(const Particle* part)
216 {
217 const MCParticle* mcparticle = part->getMCParticle();
218 if (!mcparticle) return Const::doubleNaN;
219
220 int status = MCMatching::getMCErrors(part, mcparticle);
221 //remove the following bits
222 status &= (~MCMatching::c_MissGamma);
223
224 return (status == MCMatching::c_Correct);
225 }
226
227 double isSignalAcceptMissing(const Particle* part)
228 {
229 const MCParticle* mcparticle = part->getMCParticle();
230 if (!mcparticle) return Const::doubleNaN;
231
232 int status = MCMatching::getMCErrors(part, mcparticle);
233 //remove the following bits
234 status &= (~MCMatching::c_MissGamma);
235 status &= (~MCMatching::c_MissMassiveParticle);
236 status &= (~MCMatching::c_MissKlong);
237 status &= (~MCMatching::c_MissNeutrino);
238
239 return (status == MCMatching::c_Correct);
240 }
241
242 double isSignalAcceptBremsPhotons(const Particle* part)
243 {
244 const MCParticle* mcparticle = part->getMCParticle();
245 if (!mcparticle) return Const::doubleNaN;
246
247 int status = MCMatching::getMCErrors(part, mcparticle);
248 //remove the following bits
249 status &= (~MCMatching::c_AddedRecoBremsPhoton);
250
251 return (status == MCMatching::c_Correct);
252 }
253
254 double particleMCMatchPDGCode(const Particle* part)
255 {
256 const MCParticle* mcparticle = part->getMCParticle();
257 if (!mcparticle) return Const::doubleNaN;
258 return mcparticle->getPDG();
259 }
260
261 double particleMCErrors(const Particle* part)
262 {
263 return MCMatching::getMCErrors(part);
264 }
265
266 double particleNumberOfMCMatch(const Particle* particle)
267 {
268 RelationVector<MCParticle> mcRelations = particle->getRelationsTo<MCParticle>();
269 return (mcRelations.size());
270 }
271
272 double particleMCMatchWeight(const Particle* particle)
273 {
274 auto relWithWeight = particle->getRelatedToWithWeight<MCParticle>();
275 if (!relWithWeight.first) return Const::doubleNaN;
276 return relWithWeight.second;
277 }
278
279 double particleMCMatchDecayTime(const Particle* part)
280 {
281 const MCParticle* mcparticle = part->getMCParticle();
282 if (!mcparticle) return Const::doubleNaN;
283 return mcparticle->getDecayTime();
284 }
285
286 double particleMCMatchLifeTime(const Particle* part)
287 {
288 const MCParticle* mcparticle = part->getMCParticle();
289 if (!mcparticle) return Const::doubleNaN;
290 return mcparticle->getLifetime();
291 }
292
293 double particleMCMatchPX(const Particle* part)
294 {
295 const MCParticle* mcparticle = part->getMCParticle();
296 if (!mcparticle) return Const::doubleNaN;
297
298 const auto& frame = ReferenceFrame::GetCurrent();
299 ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
300 return frame.getMomentum(mcpP4).Px();
301 }
302
303 double particleMCMatchPY(const Particle* part)
304 {
305 const MCParticle* mcparticle = part->getMCParticle();
306 if (!mcparticle) return Const::doubleNaN;
307
308 const auto& frame = ReferenceFrame::GetCurrent();
309 ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
310 return frame.getMomentum(mcpP4).Py();
311 }
312
313 double particleMCMatchPZ(const Particle* part)
314 {
315 const MCParticle* mcparticle = part->getMCParticle();
316 if (!mcparticle) return Const::doubleNaN;
317
318 const auto& frame = ReferenceFrame::GetCurrent();
319 ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
320 return frame.getMomentum(mcpP4).Pz();
321 }
322
323 double particleMCMatchPT(const Particle* part)
324 {
325 const MCParticle* mcparticle = part->getMCParticle();
326 if (!mcparticle) return Const::doubleNaN;
327
328 const auto& frame = ReferenceFrame::GetCurrent();
329 ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
330 return frame.getMomentum(mcpP4).Pt();
331 }
332
333 double particleMCMatchE(const Particle* part)
334 {
335 const MCParticle* mcparticle = part->getMCParticle();
336 if (!mcparticle) return Const::doubleNaN;
337
338 const auto& frame = ReferenceFrame::GetCurrent();
339 ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
340 return frame.getMomentum(mcpP4).E();
341 }
342
343 double particleMCMatchP(const Particle* part)
344 {
345 const MCParticle* mcparticle = part->getMCParticle();
346 if (!mcparticle) return Const::doubleNaN;
347
348 const auto& frame = ReferenceFrame::GetCurrent();
349 ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
350 return frame.getMomentum(mcpP4).P();
351 }
352
353 double particleMCMatchTheta(const Particle* part)
354 {
355 const MCParticle* mcparticle = part->getMCParticle();
356 if (!mcparticle) return Const::doubleNaN;
357
358 const auto& frame = ReferenceFrame::GetCurrent();
359 ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
360 return frame.getMomentum(mcpP4).Theta();
361 }
362
363 double particleMCMatchPhi(const Particle* part)
364 {
365 const MCParticle* mcparticle = part->getMCParticle();
366 if (!mcparticle) return Const::doubleNaN;
367
368 const auto& frame = ReferenceFrame::GetCurrent();
369 ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
370 return frame.getMomentum(mcpP4).Phi();
371 }
372
373 double mcParticleNDaughters(const Particle* part)
374 {
375 const MCParticle* mcparticle = part->getMCParticle();
376
377 if (!mcparticle) return Const::doubleNaN;
378 return mcparticle->getNDaughters();
379 }
380
381 double particleMCRecoilMass(const Particle* part)
382 {
383 StoreArray<MCParticle> mcparticles;
384 if (mcparticles.getEntries() < 1) return Const::doubleNaN;
385
386 ROOT::Math::PxPyPzEVector pInitial = mcparticles[0]->get4Vector();
387 ROOT::Math::PxPyPzEVector pDaughters;
388 const std::vector<Particle*> daughters = part->getDaughters();
389 for (auto daughter : daughters) {
390 const MCParticle* mcD = daughter->getMCParticle();
391 if (!mcD) return Const::doubleNaN;
392
393 pDaughters += mcD->get4Vector();
394 }
395 return (pInitial - pDaughters).M();
396 }
397
398 ROOT::Math::PxPyPzEVector MCInvisibleP4(const MCParticle* mcparticle)
399 {
400 ROOT::Math::PxPyPzEVector ResultP4;
401 int pdg = abs(mcparticle->getPDG());
402 bool isNeutrino = (pdg == 12 or pdg == 14 or pdg == 16);
403
404 if (mcparticle->getNDaughters() > 0) {
405 const std::vector<MCParticle*> daughters = mcparticle->getDaughters();
406 for (auto daughter : daughters)
407 ResultP4 += MCInvisibleP4(daughter);
408 } else if (isNeutrino)
409 ResultP4 += mcparticle->get4Vector();
410
411 return ResultP4;
412 }
413
414 double particleMCCosThetaBetweenParticleAndNominalB(const Particle* part)
415 {
416 int particlePDG = abs(part->getPDGCode());
417 if (particlePDG != 511 and particlePDG != 521)
418 B2FATAL("The variable mcCosThetaBetweenParticleAndNominalB is only meant to be used on B mesons!");
419
420 PCmsLabTransform T;
421 double e_Beam = T.getCMSEnergy() / 2.0; // GeV
422 double m_B = part->getPDGMass();
423
424 // Y(4S) mass according PDG (https://pdg.lbl.gov/2020/listings/rpp2020-list-upsilon-4S.pdf)
425 const double mY4S = 10.5794; // GeV
426
427 // if this is a continuum run, use an approximate Y(4S) CMS energy
428 if (e_Beam * e_Beam - m_B * m_B < 0) {
429 e_Beam = mY4S / 2.0;
430 }
431 double p_B = std::sqrt(e_Beam * e_Beam - m_B * m_B);
432
433 // Calculate cosThetaBY with daughter neutrino momenta subtracted
434 const MCParticle* mcB = part->getMCParticle();
435 if (!mcB) return Const::doubleNaN;
436
437 int mcParticlePDG = abs(mcB->getPDG());
438 if (mcParticlePDG != 511 and mcParticlePDG != 521)
439 return Const::doubleNaN;
440
441 ROOT::Math::PxPyPzEVector p = T.rotateLabToCms() * (mcB->get4Vector() - MCInvisibleP4(mcB));
442 double e_d = p.E();
443 double m_d = p.M();
444 double p_d = p.P();
445
446 double theta_BY = (2 * e_Beam * e_d - m_B * m_B - m_d * m_d)
447 / (2 * p_B * p_d);
448 return theta_BY;
449 }
450
451 double mcParticleSecondaryPhysicsProcess(const Particle* p)
452 {
453 const MCParticle* mcp = p->getMCParticle();
454 if (!mcp) return Const::doubleNaN;
455 return mcp->getSecondaryPhysicsProcess();
456 }
457
458 double mcParticleStatus(const Particle* p)
459 {
460 const MCParticle* mcp = p->getMCParticle();
461 if (!mcp) return Const::doubleNaN;
462 return mcp->getStatus();
463 }
464
465 double particleMCPrimaryParticle(const Particle* p)
466 {
467 const MCParticle* mcp = p->getMCParticle();
468 if (!mcp) return Const::doubleNaN;
469
470 unsigned int bitmask = MCParticle::c_PrimaryParticle;
471 return mcp->hasStatus(bitmask);
472 }
473
474 double particleMCVirtualParticle(const Particle* p)
475 {
476 const MCParticle* mcp = p->getMCParticle();
477 if (!mcp) return Const::doubleNaN;
478
479 unsigned int bitmask = MCParticle::c_IsVirtual;
480 return mcp->hasStatus(bitmask);
481 }
482
483 double particleMCInitialParticle(const Particle* p)
484 {
485 const MCParticle* mcp = p->getMCParticle();
486 if (!mcp) return Const::doubleNaN;
487
488 unsigned int bitmask = MCParticle::c_Initial;
489 return mcp->hasStatus(bitmask);
490 }
491
492 double particleMCISRParticle(const Particle* p)
493 {
494 const MCParticle* mcp = p->getMCParticle();
495 if (!mcp) return Const::doubleNaN;
496
497 unsigned int bitmask = MCParticle::c_IsISRPhoton;
498 return mcp->hasStatus(bitmask);
499 }
500
501 double particleMCFSRParticle(const Particle* p)
502 {
503 const MCParticle* mcp = p->getMCParticle();
504 if (!mcp) return Const::doubleNaN;
505
506 unsigned int bitmask = MCParticle::c_IsFSRPhoton;
507 return mcp->hasStatus(bitmask);
508 }
509
510 double particleMCPhotosParticle(const Particle* p)
511 {
512 const MCParticle* mcp = p->getMCParticle();
513 if (!mcp) return Const::doubleNaN;
514
515 unsigned int bitmask = MCParticle::c_IsPHOTOSPhoton;
516 return mcp->hasStatus(bitmask);
517 }
518
519 double generatorEventWeight(const Particle*)
520 {
521 StoreObjPtr<EventMetaData> evtMetaData;
522 if (!evtMetaData) return Const::doubleNaN;
523 return evtMetaData->getGeneratedWeight();
524 }
525
526 int tauPlusMcMode(const Particle*)
527 {
528 StoreObjPtr<TauPairDecay> tauDecay;
529 if (!tauDecay) {
530 B2WARNING("Cannot find tau decay ID, did you forget to run TauDecayMarkerModule?");
531 return 0;
532 }
533 return tauDecay->getTauPlusIdMode();
534 }
535
536 int tauMinusMcMode(const Particle*)
537 {
538 StoreObjPtr<TauPairDecay> tauDecay;
539 if (!tauDecay) {
540 B2WARNING("Cannot find tau decay ID, did you forget to run TauDecayMarkerModule?");
541 return 0;
542 }
543 return tauDecay->getTauMinusIdMode();
544 }
545
546 int tauPlusMcProng(const Particle*)
547 {
548 StoreObjPtr<TauPairDecay> tauDecay;
549 if (!tauDecay) {
550 B2WARNING("Cannot find tau prong, did you forget to run TauDecayMarkerModule?");
551 return 0;
552 }
553 return tauDecay->getTauPlusMcProng();
554 }
555
556 int tauMinusMcProng(const Particle*)
557 {
558 StoreObjPtr<TauPairDecay> tauDecay;
559 if (!tauDecay) {
560 B2WARNING("Cannot find tau prong, did you forget to run TauDecayMarkerModule?");
561 return 0;
562 }
563 return tauDecay->getTauMinusMcProng();
564 }
565
566 double tauPlusEgstar(const Particle*)
567 {
568 StoreObjPtr<TauPairDecay> tauDecay;
569 if (!tauDecay) {
570 B2WARNING("Cannot find tau prong, did you forget to run TauDecayMarkerModule?");
571 return 0;
572 }
573 return tauDecay->getTauPlusEgstar();
574 }
575
576 double tauMinusEgstar(const Particle*)
577 {
578 StoreObjPtr<TauPairDecay> tauDecay;
579 if (!tauDecay) {
580 B2WARNING("Cannot find tau prong, did you forget to run TauDecayMarkerModule?");
581 return 0;
582 }
583 return tauDecay->getTauMinusEgstar();
584 }
585
586 double isReconstructible(const Particle* p)
587 {
588 if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
589 return Const::doubleNaN;
590 const MCParticle* mcp = p->getMCParticle();
591 if (!mcp) return Const::doubleNaN;
592
593 // If charged: make sure it was seen in the SVD.
594 // If neutral: make sure it was seen in the ECL.
595 return (abs(mcp->getCharge()) > 0) ? seenInSVD(p) : seenInECL(p);
596 }
597
598 double isTrackFound(const Particle* p)
599 {
600 if (p->getParticleSource() != Particle::EParticleSourceObject::c_MCParticle)
601 return Const::doubleNaN;
602 const MCParticle* tmp_mcP = p->getMCParticle();
603 if (!Const::chargedStableSet.contains(Const::ParticleType(abs(tmp_mcP->getPDG()))))
604 return Const::doubleNaN;
605 Track* tmp_track = tmp_mcP->getRelated<Track>();
606 if (tmp_track) {
607 const TrackFitResult* tmp_tfr = tmp_track->getTrackFitResultWithClosestMass(Const::ChargedStable(abs(tmp_mcP->getPDG())));
608 if (tmp_tfr->getChargeSign()*tmp_mcP->getCharge() > 0)
609 return 1;
610 else
611 return -1;
612 }
613 return 0;
614 }
615
616 double seenInPXD(const Particle* p)
617 {
618 if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
619 return Const::doubleNaN;
620 const MCParticle* mcp = p->getMCParticle();
621 if (!mcp) return Const::doubleNaN;
622 return mcp->hasSeenInDetector(Const::PXD);
623 }
624
625 double seenInSVD(const Particle* p)
626 {
627 if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
628 return Const::doubleNaN;
629 const MCParticle* mcp = p->getMCParticle();
630 if (!mcp) return Const::doubleNaN;
631 return mcp->hasSeenInDetector(Const::SVD);
632 }
633
634 double seenInCDC(const Particle* p)
635 {
636 if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
637 return Const::doubleNaN;
638 const MCParticle* mcp = p->getMCParticle();
639 if (!mcp) return Const::doubleNaN;
640 return mcp->hasSeenInDetector(Const::CDC);
641 }
642
643 double seenInTOP(const Particle* p)
644 {
645 if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
646 return Const::doubleNaN;
647 const MCParticle* mcp = p->getMCParticle();
648 if (!mcp) return Const::doubleNaN;
649 return mcp->hasSeenInDetector(Const::TOP);
650 }
651
652 double seenInECL(const Particle* p)
653 {
654 if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
655 return Const::doubleNaN;
656 const MCParticle* mcp = p->getMCParticle();
657 if (!mcp) return Const::doubleNaN;
658 return mcp->hasSeenInDetector(Const::ECL);
659 }
660
661 double seenInARICH(const Particle* p)
662 {
663 if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
664 return Const::doubleNaN;
665 const MCParticle* mcp = p->getMCParticle();
666 if (!mcp) return Const::doubleNaN;
667 return mcp->hasSeenInDetector(Const::ARICH);
668 }
669
670 double seenInKLM(const Particle* p)
671 {
672 if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
673 return Const::doubleNaN;
674 const MCParticle* mcp = p->getMCParticle();
675 if (!mcp) return Const::doubleNaN;
676 return mcp->hasSeenInDetector(Const::KLM);
677 }
678
679 int genNStepsToDaughter(const Particle* p, const std::vector<double>& arguments)
680 {
681 if (arguments.size() != 1)
682 B2FATAL("Wrong number of arguments for genNStepsToDaughter");
683
684 const MCParticle* mcp = p->getMCParticle();
685 if (!mcp) {
686 B2WARNING("No MCParticle is associated to the particle");
687 return 0;
688 }
689
690 int nChildren = p->getNDaughters();
691 if (arguments[0] >= nChildren) {
692 return 0;
693 }
694
695 const Particle* daugP = p->getDaughter(arguments[0]);
696 const MCParticle* daugMCP = daugP->getMCParticle();
697 if (!daugMCP) {
698 // This is a strange case.
699 // The particle, p, has the related MC particle, but i-th daughter does not have the related MC Particle.
700 B2WARNING("No MCParticle is associated to the i-th daughter");
701 return 0;
702 }
703
704 if (nChildren == 1) return 1;
705
706 std::vector<int> genMothers;
707 MCMatching::fillGenMothers(daugMCP, genMothers);
708 auto match = std::find(genMothers.begin(), genMothers.end(), mcp->getIndex());
709 return match - genMothers.begin();
710 }
711
712 int genNMissingDaughter(const Particle* p, const std::vector<double>& arguments)
713 {
714 if (arguments.size() < 1)
715 B2FATAL("Wrong number of arguments for genNMissingDaughter");
716
717 const std::vector<int> PDGcodes(arguments.begin(), arguments.end());
718
719 const MCParticle* mcp = p->getMCParticle();
720 if (!mcp) {
721 B2WARNING("No MCParticle is associated to the particle");
722 return 0;
723 }
724
725 return MCMatching::countMissingParticle(p, mcp, PDGcodes);
726 }
727
728 double getHEREnergy(const Particle*)
729 {
730 static DBObjPtr<BeamParameters> beamParamsDB;
731 return (beamParamsDB->getHER()).E();
732 }
733
734 double getLEREnergy(const Particle*)
735 {
736 static DBObjPtr<BeamParameters> beamParamsDB;
737 return (beamParamsDB->getLER()).E();
738 }
739
740 double getCrossingAngleX(const Particle*)
741 {
742 // get the beam momenta from the DB
743 static DBObjPtr<BeamParameters> beamParamsDB;
744 B2Vector3D herVec = beamParamsDB->getHER().Vect();
745 B2Vector3D lerVec = beamParamsDB->getLER().Vect();
746
747 // only looking at the horizontal (XZ plane) -> set y-coordinates to zero
748 herVec.SetY(0);
749 lerVec.SetY(0);
750
751 //calculate the crossing angle
752 return herVec.Angle(-lerVec);
753 }
754
755 double getCrossingAngleY(const Particle*)
756 {
757 // get the beam momenta from the DB
758 static DBObjPtr<BeamParameters> beamParamsDB;
759 B2Vector3D herVec = beamParamsDB->getHER().Vect();
760 B2Vector3D lerVec = beamParamsDB->getLER().Vect();
761
762 // only looking at the vertical (YZ plane) -> set x-coordinates to zero
763 herVec.SetX(0);
764 lerVec.SetX(0);
765
766 //calculate the crossing angle
767 return herVec.Angle(-lerVec);
768 }
769
770
771 double particleClusterMatchWeight(const Particle* particle)
772 {
773 /* Get the weight of the *cluster* mc match for the mcparticle matched to
774 * this particle.
775 *
776 * Note that for track-based particles this is different from the mc match
777 * of the particle (which it inherits from the mc match of the track)
778 */
779 const MCParticle* matchedToParticle = particle->getMCParticle();
780 if (!matchedToParticle) return Const::doubleNaN;
781 int matchedToIndex = matchedToParticle->getArrayIndex();
782
783 const ECLCluster* cluster = particle->getECLCluster();
784 if (!cluster) return Const::doubleNaN;
785
786 const auto mcps = cluster->getRelationsTo<MCParticle>();
787 for (unsigned int i = 0; i < mcps.size(); ++i)
788 if (mcps[i]->getArrayIndex() == matchedToIndex)
789 return mcps.weight(i);
790
791 return Const::doubleNaN;
792 }
793
794 double particleClusterBestMCMatchWeight(const Particle* particle)
795 {
796 /* Get the weight of the best mc match of the cluster associated to
797 * this particle.
798 *
799 * Note for electrons (or any track-based particle) this may not be
800 * the same thing as the mc match of the particle (which is taken
801 * from the track).
802 *
803 * For photons (or any ECL-based particle) this will be the same as the
804 * mcMatchWeight
805 */
806 const ECLCluster* cluster = particle->getECLCluster();
807 if (!cluster) return Const::doubleNaN;
808
809 /* loop over all mcparticles related to this cluster, find the largest
810 * weight by std::sort-ing the doubles
811 */
812 auto mcps = cluster->getRelationsTo<MCParticle>();
813 if (mcps.size() == 0) return Const::doubleNaN;
814
815 std::vector<double> weights;
816 for (unsigned int i = 0; i < mcps.size(); ++i)
817 weights.emplace_back(mcps.weight(i));
818
819 // sort descending by weight
820 std::sort(weights.begin(), weights.end());
821 std::reverse(weights.begin(), weights.end());
822 return weights[0];
823 }
824
825 double particleClusterBestMCPDGCode(const Particle* particle)
826 {
827 /* Get the PDG code of the best mc match of the cluster associated to this
828 * particle.
829 *
830 * Note for electrons (or any track-based particle) this may not be the
831 * same thing as the mc match of the particle (which is taken from the track).
832 *
833 * For photons (or any ECL-based particle) this will be the same as the mcPDG
834 */
835 const ECLCluster* cluster = particle->getECLCluster();
836 if (!cluster) return Const::doubleNaN;
837
838 auto mcps = cluster->getRelationsTo<MCParticle>();
839 if (mcps.size() == 0) return Const::doubleNaN;
840
841 std::vector<std::pair<double, int>> weightsAndIndices;
842 for (unsigned int i = 0; i < mcps.size(); ++i)
843 weightsAndIndices.emplace_back(mcps.weight(i), i);
844
845 // sort descending by weight
846 std::sort(weightsAndIndices.begin(), weightsAndIndices.end(),
847 ValueIndexPairSorting::higherPair<decltype(weightsAndIndices)::value_type>);
848 // cppcheck-suppress containerOutOfBounds
849 return mcps.object(weightsAndIndices[0].second)->getPDG();
850 }
851
852 double particleClusterTotalMCMatchWeight(const Particle* particle)
853 {
854 const ECLCluster* cluster = particle->getECLCluster();
855 if (!cluster) return Const::doubleNaN;
856
857 auto mcps = cluster->getRelationsTo<MCParticle>();
858
859 // if there are no relations to any MCParticles, we return 0!
860 double weightsum = 0;
861 for (unsigned int i = 0; i < mcps.size(); ++i)
862 weightsum += mcps.weight(i);
863
864 return weightsum;
865 }
866
867 // Helper function for particleClusterTotalMCMatchWeightForKlong
868 void getKlongWeightMap(const Particle* particle, std::map<int, double>& mapMCParticleIndxAndWeight)
869 {
870 const ECLCluster* cluster = particle->getECLCluster();
871 auto mcps = cluster->getRelationsTo<MCParticle>();
872
873 for (unsigned int i = 0; i < mcps.size(); ++i) {
874 double weight = mcps.weight(i);
875 const MCParticle* mcp = mcps[i];
876
877 while (mcp) {
878 if (mcp->getPDG() == 130) {
879 int index = mcp->getArrayIndex();
880 if (mapMCParticleIndxAndWeight.find(index) != mapMCParticleIndxAndWeight.end()) {
881 mapMCParticleIndxAndWeight.at(index) = mapMCParticleIndxAndWeight.at(index) + weight;
882 } else {
883 mapMCParticleIndxAndWeight.insert({index, weight});
884 }
885 break;
886 } else {
887 mcp = mcp->getMother();
888 }
889 }
890 }
891 }
892
893 double particleClusterTotalMCMatchWeightForKlong(const Particle* particle)
894 {
895 const ECLCluster* cluster = particle->getECLCluster();
896 if (!cluster) return Const::doubleNaN;
897
898 auto mcps = cluster->getRelationsTo<MCParticle>();
899 if (mcps.size() == 0) return Const::doubleNaN;
900
901 std::map<int, double> mapMCParticleIndxAndWeight;
902 getKlongWeightMap(particle, mapMCParticleIndxAndWeight);
903
904 double totalWeight = 0;
905 for (const auto& map : mapMCParticleIndxAndWeight) {
906 totalWeight += map.second;
907 }
908
909 return totalWeight;
910 }
911
912 double particleClusterTotalMCMatchWeightForBestKlong(const Particle* particle)
913 {
914 const ECLCluster* cluster = particle->getECLCluster();
915 if (!cluster) return Const::doubleNaN;
916
917 auto mcps = cluster->getRelationsTo<MCParticle>();
918 if (mcps.size() == 0) return Const::doubleNaN;
919
920 std::map<int, double> mapMCParticleIndxAndWeight;
921 getKlongWeightMap(particle, mapMCParticleIndxAndWeight);
922
923 if (mapMCParticleIndxAndWeight.size() == 0)
924 return 0.0;
925
926 auto maxMap = std::max_element(mapMCParticleIndxAndWeight.begin(), mapMCParticleIndxAndWeight.end(),
927 [](const auto & x, const auto & y) { return x.second < y.second; }
928 );
929
930 return maxMap->second;
931 }
932
933 double isBBCrossfeed(const Particle* particle)
934 {
935 if (particle == nullptr)
936 return Const::doubleNaN;
937
938 int pdg = particle->getPDGCode();
939 if (abs(pdg) != 511 && abs(pdg) != 521 && abs(pdg) != 531)
940 return Const::doubleNaN;
941
942 std::vector<const Particle*> daughters = particle->getFinalStateDaughters();
943 int nDaughters = daughters.size();
944 if (nDaughters <= 1)
945 return 0;
946 std::vector<int> mother_ids;
947
948 for (int j = 0; j < nDaughters; ++j) {
949 const MCParticle* curMCParticle = daughters[j]->getMCParticle();
950 while (curMCParticle != nullptr) {
951 pdg = curMCParticle->getPDG();
952 if (abs(pdg) == 511 || abs(pdg) == 521 || abs(pdg) == 531) {
953 mother_ids.emplace_back(curMCParticle->getArrayIndex());
954 break;
955 }
956 const MCParticle* curMCMother = curMCParticle->getMother();
957 curMCParticle = curMCMother;
958 }
959 if (curMCParticle == nullptr) {
960 return Const::doubleNaN;
961 }
962 }
963
964 std::set<int> distinctIDs = std::set(mother_ids.begin(), mother_ids.end());
965 if (distinctIDs.size() == 1)
966 return 0;
967 else
968 return 1;
969 }
970
971 int ancestorBIndex(const Particle* particle)
972 {
973 const MCParticle* mcpart = particle->getMCParticle();
974
975 while (mcpart) {
976 int pdg = std::abs(mcpart->getPDG());
977
978 if ((pdg == 521) || (pdg == 511))
979 return mcpart->getArrayIndex();
980
981 mcpart = mcpart->getMother();
982 }
983
984 return -1;
985 }
986
987 VARIABLE_GROUP("MC matching and MC truth");
988 REGISTER_VARIABLE("isSignal", isSignal,
989 "1.0 if Particle is correctly reconstructed (SIGNAL), 0.0 if not, and NaN if no related MCParticle could be found. \n"
990 "It behaves according to DecayStringGrammar.");
991 REGISTER_VARIABLE("isSignalAcceptWrongFSPs", isSignalAcceptWrongFSPs,
992 "1.0 if Particle is almost correctly reconstructed (SIGNAL), 0.0 if not, and NaN if no related MCParticle could be found.\n"
993 "Misidentification of charged FSP is allowed.");
994 REGISTER_VARIABLE("isPrimarySignal", isPrimarySignal,
995 "1.0 if Particle is correctly reconstructed (SIGNAL) and primary, 0.0 if not, and NaN if no related MCParticle could be found");
996 REGISTER_VARIABLE("isSignalAcceptBremsPhotons", isSignalAcceptBremsPhotons,
997 "1.0 if Particle is correctly reconstructed (SIGNAL), 0.0 if not, and NaN if no related MCParticle could be found.\n"
998 "Particles with gamma daughters attached through the bremsstrahlung recovery modules are allowed.");
999 REGISTER_VARIABLE("genMotherPDG", genMotherPDG,
1000 "Check the PDG code of a particles MC mother particle");
1001 REGISTER_VARIABLE("genMotherPDG(i)", genNthMotherPDG,
1002 "Check the PDG code of a particles n-th MC mother particle by providing an argument. 0 is first mother, 1 is grandmother etc. :noindex:");
1003
1004 REGISTER_VARIABLE("genQ2PmPd(i,j,...)", genQ2PmPd, R"DOC(
1005 Returns the generated four momentum transfer squared :math:`q^2` calculated as :math:`q^2 = (p_m - p_{d_i} - p_{d_j} - ...)^2`.
1006
1007 Here :math:`p_m` is the four momentum of the given (mother) particle,
1008 and :math:`p_{d_{i,j,...}}` are the daughter particles with indices given as arguments .
1009 The ordering of daughters is as defined in the DECAY.DEC file used in the generation, with the numbering starting at :math:`N=0`.
1010
1011 Returns NaN if no related MCParticle could be found.
1012 Returns NaN if any of the given indices is larger than the number of daughters of the given particle.
1013
1014 )DOC", ":math:`[\\text{GeV}/\\text{c}]^2`");
1015
1016 REGISTER_VARIABLE("genMotherID", genMotherIndex,
1017 "Check the array index of a particles generated mother");
1018 REGISTER_VARIABLE("genMotherID(i)", genNthMotherIndex,
1019 "Check the array index of a particle n-th MC mother particle by providing an argument. 0 is first mother, 1 is grandmother etc. :noindex:");
1020 // genMotherPDG and genMotherID are overloaded (each are two C++ functions
1021 // sharing one variable name) so one of the two needs to be made the indexed
1022 // variable in sphinx
1023 REGISTER_VARIABLE("isBBCrossfeed", isBBCrossfeed,
1024 "Returns 1 for crossfeed in reconstruction of given B meson, 0 for no crossfeed and NaN for no true B meson or failed truthmatching.");
1025 REGISTER_VARIABLE("ancestorBIndex", ancestorBIndex,
1026 "Returns array index of B ancestor, or -1 if no B or no MC-matching is found.");
1027 REGISTER_VARIABLE("genMotherP", genMotherP,
1028 "Generated momentum of a particles MC mother particle\n\n", "GeV/c");
1029 REGISTER_VARIABLE("genParticleID", genParticleIndex,
1030 "Check the array index of a particle's related MCParticle");
1031 REGISTER_VARIABLE("isSignalAcceptMissingNeutrino",
1032 isSignalAcceptMissingNeutrino,
1033 "Same as isSignal, but also accept missing neutrino");
1034 REGISTER_VARIABLE("isSignalAcceptMissingMassive",
1035 isSignalAcceptMissingMassive,
1036 "Same as isSignal, but also accept missing massive particle");
1037 REGISTER_VARIABLE("isSignalAcceptMissingGamma",
1038 isSignalAcceptMissingGamma,
1039 "Same as isSignal, but also accept missing gamma, such as B -> K* gamma, pi0 -> gamma gamma");
1040 REGISTER_VARIABLE("isSignalAcceptMissing",
1041 isSignalAcceptMissing,
1042 "Same as isSignal, but also accept missing particle");
1043 REGISTER_VARIABLE("isMisidentified", isMisidentified,
1044 "Return 1 if the particle is misidentified: at least one of the final state particles has the wrong PDG code assignment (including wrong charge), 0 if PDG code is fine, and NaN if no related MCParticle could be found.");
1045 REGISTER_VARIABLE("isWrongCharge", isWrongCharge,
1046 "Return 1 if the charge of the particle is wrongly assigned, 0 if it's the correct charge, and NaN if no related MCParticle could be found.");
1047 REGISTER_VARIABLE("isCloneTrack", isCloneTrack,
1048 "Return 1 if the charged final state particle comes from a cloned track, 0 if not a clone. Returns NAN if neutral, composite, or MCParticle not found (like for data or if not MCMatched)");
1049 REGISTER_VARIABLE("isOrHasCloneTrack", isOrHasCloneTrack,
1050 "Return 1 if the particle is a clone track or has a clone track as a daughter, 0 otherwise.");
1051 REGISTER_VARIABLE("mcPDG", particleMCMatchPDGCode,
1052 "The PDG code of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).");
1053 REGISTER_VARIABLE("mcErrors", particleMCErrors,
1054 "The bit pattern indicating the quality of MC match (see MCMatching::MCErrorFlags)");
1055 REGISTER_VARIABLE("mcMatchWeight", particleMCMatchWeight,
1056 "The weight of the Particle -> MCParticle relation (only for the first Relation = largest weight).");
1057 REGISTER_VARIABLE("nMCMatches", particleNumberOfMCMatch,
1058 "The number of relations of this Particle to MCParticle.");
1059 REGISTER_VARIABLE("mcDecayTime", particleMCMatchDecayTime,
1060 "The decay time of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1061 "ns");
1062 REGISTER_VARIABLE("mcLifeTime", particleMCMatchLifeTime,
1063 "The life time of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1064 "ns");
1065 REGISTER_VARIABLE("mcPX", particleMCMatchPX,
1066 "The px of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1067 "GeV/c");
1068 REGISTER_VARIABLE("mcPY", particleMCMatchPY,
1069 "The py of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1070 "GeV/c");
1071 REGISTER_VARIABLE("mcPZ", particleMCMatchPZ,
1072 "The pz of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1073 "GeV/c");
1074 REGISTER_VARIABLE("mcPT", particleMCMatchPT,
1075 "The pt of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1076 "GeV/c");
1077 REGISTER_VARIABLE("mcE", particleMCMatchE,
1078 "The energy of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1079 "GeV");
1080 REGISTER_VARIABLE("mcP", particleMCMatchP,
1081 "The total momentum of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1082 "GeV/c");
1083 REGISTER_VARIABLE("mcPhi", particleMCMatchPhi,
1084 "The phi of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1085 "rad");
1086 REGISTER_VARIABLE("mcTheta", particleMCMatchTheta,
1087 "The theta of matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).\n\n",
1088 "rad");
1089 REGISTER_VARIABLE("nMCDaughters", mcParticleNDaughters,
1090 "The number of daughters of the matched MCParticle, NaN if no match. Requires running matchMCTruth() on the reconstructed particles, or a particle list filled with generator particles (MCParticle objects).");
1091 REGISTER_VARIABLE("mcRecoilMass", particleMCRecoilMass,
1092 "The mass recoiling against the particles attached as particle's daughters calculated using MC truth values.\n\n",
1093 "GeV/:math:`\\text{c}^2`");
1094 REGISTER_VARIABLE("mcCosThetaBetweenParticleAndNominalB",
1095 particleMCCosThetaBetweenParticleAndNominalB,
1096 "Cosine of the angle in CMS between momentum the particle and a nominal B particle. In this calculation, the momenta of all descendant neutrinos are subtracted from the B momentum.");
1097
1098
1099 REGISTER_VARIABLE("mcSecPhysProc", mcParticleSecondaryPhysicsProcess,
1100 R"DOC(
1101Returns the secondary physics process flag, which is set by Geant4 on secondary particles. It indicates the type of process that produced the particle.
1102
1103Returns NaN if the particle is not matched to a MCParticle.
1104
1105Returns -1 in case of unknown process.
1106
1107Returns 0 if the particle is primary, i.e. produced by the event generator and not Geant4. Particles produced by Geant4 (i.e. secondary particles) include those produced in interaction with detector material, Bremsstrahlung, and the decay products of long-lived particles (e.g. muons, pions, K_S0, K_L0, Lambdas, ...).
1108
1109List of possible values (taken from the Geant4 source of
1110`G4DecayProcessType <https://github.com/Geant4/geant4/blob/v10.6.3/source/processes/decay/include/G4DecayProcessType.hh>`_,
1111`G4HadronicProcessType <https://github.com/Geant4/geant4/blob/v10.6.3/source/processes/hadronic/management/include/G4HadronicProcessType.hh>`_,
1112`G4TransportationProcessType <https://github.com/Geant4/geant4/blob/v10.6.3/source/processes/transportation/include/G4TransportationProcessType.hh>`_ and
1113`G4EmProcessSubType <https://github.com/Geant4/geant4/blob/v10.6.3/source/processes/electromagnetic/utils/include/G4EmProcessSubType.hh>`_)
1114
1115* 1 Coulomb scattering
1116* 2 Ionisation
1117* 3 Bremsstrahlung
1118* 4 Pair production by charged
1119* 5 Annihilation
1120* 6 Annihilation to mu mu
1121* 7 Annihilation to hadrons
1122* 8 Nuclear stopping
1123* 9 Electron general process
1124* 10 Multiple scattering
1125* 11 Rayleigh
1126* 12 Photo-electric effect
1127* 13 Compton scattering
1128* 14 Gamma conversion
1129* 15 Gamma conversion to mu mu
1130* 16 Gamma general process
1131* 21 Cerenkov
1132* 22 Scintillation
1133* 23 Synchrotron radiation
1134* 24 Transition radiation
1135* 91 Transportation
1136* 92 Coupled transportation
1137* 111 Hadron elastic
1138* 121 Hadron inelastic
1139* 131 Capture
1140* 132 Mu atomic capture
1141* 141 Fission
1142* 151 Hadron at rest
1143* 152 Lepton at rest
1144* 161 Charge exchange
1145* 201 Decay
1146* 202 Decay with spin
1147* 203 Decay (pion make spin)
1148* 210 Radioactive decay
1149* 211 Unknown decay
1150* 221 Mu atom decay
1151* 231 External decay
1152
1153.. note:: This is what `modularAnalysis.printMCParticles` shows as ``creation process`` when ``showStatus`` is set to ``True``.
1154)DOC");
1155 REGISTER_VARIABLE("mcParticleStatus", mcParticleStatus,
1156 "Returns status bits of related MCParticle or NaN if MCParticle relation is not set.");
1157 REGISTER_VARIABLE("mcPrimary", particleMCPrimaryParticle,
1158 "Returns 1 if Particle is related to primary MCParticle, 0 if Particle is related to non - primary MCParticle, "
1159 "NaN if Particle is not related to MCParticle.");
1160 REGISTER_VARIABLE("mcVirtual", particleMCVirtualParticle,
1161 "Returns 1 if Particle is related to virtual MCParticle, 0 if Particle is related to non - virtual MCParticle, "
1162 "NaN if Particle is not related to MCParticle.")
1163 REGISTER_VARIABLE("mcInitial", particleMCInitialParticle,
1164 "Returns 1 if Particle is related to initial MCParticle, 0 if Particle is related to non - initial MCParticle, "
1165 "NaN if Particle is not related to MCParticle.")
1166 REGISTER_VARIABLE("mcISR", particleMCISRParticle,
1167 "Returns 1 if Particle is related to ISR MCParticle, 0 if Particle is related to non - ISR MCParticle, "
1168 "NaN if Particle is not related to MCParticle.")
1169 REGISTER_VARIABLE("mcFSR", particleMCFSRParticle,
1170 "Returns 1 if Particle is related to FSR MCParticle, 0 if Particle is related to non - FSR MCParticle ,"
1171 "NaN if Particle is not related to MCParticle.")
1172 REGISTER_VARIABLE("mcPhotos", particleMCPhotosParticle,
1173 "Returns 1 if Particle is related to Photos MCParticle, 0 if Particle is related to non - Photos MCParticle, "
1174 "NaN if Particle is not related to MCParticle.")
1175 REGISTER_VARIABLE("generatorEventWeight", generatorEventWeight,
1176 "[Eventbased] Returns the event weight produced by the event generator")
1177
1178 REGISTER_VARIABLE("genNStepsToDaughter(i)", genNStepsToDaughter,
1179 "Returns number of steps to i-th daughter from the particle at generator level. "
1180 "NaN if no MCParticle is associated to the particle or i-th daughter. "
1181 "NaN if i-th daughter does not exist.");
1182 REGISTER_VARIABLE("genNMissingDaughter(PDG)", genNMissingDaughter,
1183 "Returns the number of missing daughters having assigned PDG codes. "
1184 "NaN if no MCParticle is associated to the particle.");
1185 REGISTER_VARIABLE("Eher", getHEREnergy, R"DOC(
1186[Eventbased] The nominal HER energy used by the generator.
1187
1188.. warning:: This variable does not make sense for data.
1189
1190)DOC","GeV");
1191 REGISTER_VARIABLE("Eler", getLEREnergy, R"DOC(
1192[Eventbased] The nominal LER energy used by the generator.
1193
1194.. warning:: This variable does not make sense for data.
1195
1196)DOC","GeV");
1197 REGISTER_VARIABLE("XAngle", getCrossingAngleX, R"DOC(
1198[Eventbased] The nominal beam crossing angle in the x-z plane from generator level beam kinematics.
1199
1200.. warning:: This variable does not make sense for data.
1201
1202)DOC","rad");
1203 REGISTER_VARIABLE("YAngle", getCrossingAngleY, R"DOC(
1204[Eventbased] The nominal beam crossing angle in the y-z plane from generator level beam kinematics.
1205
1206.. warning:: This variable does not make sense for data.
1207
1208)DOC","rad");
1209
1210 VARIABLE_GROUP("Generated tau decay information");
1211 REGISTER_VARIABLE("tauPlusMCMode", tauPlusMcMode,
1212 "[Eventbased] Decay ID for the positive tau lepton in a tau pair generated event.");
1213 REGISTER_VARIABLE("tauMinusMCMode", tauMinusMcMode,
1214 "[Eventbased] Decay ID for the negative tau lepton in a tau pair generated event.");
1215 REGISTER_VARIABLE("tauPlusMCProng", tauPlusMcProng,
1216 "[Eventbased] Prong for the positive tau lepton in a tau pair generated event.");
1217 REGISTER_VARIABLE("tauMinusMCProng", tauMinusMcProng,
1218 "[Eventbased] Prong for the negative tau lepton in a tau pair generated event.");
1219 REGISTER_VARIABLE("tauPlusEgstar", tauPlusEgstar,
1220 "[Eventbased] Energy of radiated gamma from positive tau lepton in a tau pair generated event.");
1221 REGISTER_VARIABLE("tauMinusEgstar", tauMinusEgstar,
1222 "[Eventbased] Energy of radiated gamma from negative tau lepton in a tau pair generated event.");
1223
1224 VARIABLE_GROUP("MC particle seen in subdetectors");
1225 REGISTER_VARIABLE("isReconstructible", isReconstructible,
1226 "Checks charged particles were seen in the SVD and neutrals in the ECL, returns 1.0 if so, 0.0 if not, NaN for composite particles or if no related MCParticle could be found. Useful for generator studies, not for reconstructed particles.");
1227 REGISTER_VARIABLE("seenInPXD", seenInPXD,
1228 "Returns 1.0 if the MC particle was seen in the PXD, 0.0 if not, NaN for composite particles or if no related MCParticle could be found. Useful for generator studies, not for reconstructed particles.");
1229 REGISTER_VARIABLE("isTrackFound", isTrackFound,
1230 "works on charged stable particle list created from MCParticles, returns NaN if not ; returns 1.0 if there is a reconstructed track related to the charged stable MCParticle with the correct charge, return -1.0 if the reconstructed track has the wrong charge, return 0.0 when no reconstructed track is found.");
1231 REGISTER_VARIABLE("seenInSVD", seenInSVD,
1232 "Returns 1.0 if the MC particle was seen in the SVD, 0.0 if not, NaN for composite particles or if no related MCParticle could be found. Useful for generator studies, not for reconstructed particles.");
1233 REGISTER_VARIABLE("seenInCDC", seenInCDC,
1234 "Returns 1.0 if the MC particle was seen in the CDC, 0.0 if not, NaN for composite particles or if no related MCParticle could be found. Useful for generator studies, not for reconstructed particles.");
1235 REGISTER_VARIABLE("seenInTOP", seenInTOP,
1236 "Returns 1.0 if the MC particle was seen in the TOP, 0.0 if not, NaN for composite particles or if no related MCParticle could be found. Useful for generator studies, not for reconstructed particles.");
1237 REGISTER_VARIABLE("seenInECL", seenInECL,
1238 "Returns 1.0 if the MC particle was seen in the ECL, 0.0 if not, NaN for composite particles or if no related MCParticle could be found. Useful for generator studies, not for reconstructed particles.");
1239 REGISTER_VARIABLE("seenInARICH", seenInARICH,
1240 "Returns 1.0 if the MC particle was seen in the ARICH, 0.0 if not, NaN for composite particles or if no related MCParticle could be found. Useful for generator studies, not for reconstructed particles.");
1241 REGISTER_VARIABLE("seenInKLM", seenInKLM,
1242 "Returns 1.0 if the MC particle was seen in the KLM, 0.0 if not, NaN for composite particles or if no related MCParticle could be found. Useful for generator studies, not for reconstructed particles.");
1243
1244 VARIABLE_GROUP("MC Matching for ECLClusters");
1245 REGISTER_VARIABLE("clusterMCMatchWeight", particleClusterMatchWeight,
1246 "Returns the weight of the ECLCluster -> MCParticle relation for the MCParticle matched to the particle. "
1247 "Returns NaN if: no cluster is related to the particle, the particle is not MC matched, or if there are no mcmatches for the cluster. "
1248 "Returns -1 if the cluster *was* matched to particles, but not the match of the particle provided.");
1249 REGISTER_VARIABLE("clusterBestMCMatchWeight", particleClusterBestMCMatchWeight,
1250 "Returns the weight of the ECLCluster -> MCParticle relation for the relation with the largest weight.");
1251 REGISTER_VARIABLE("clusterBestMCPDG", particleClusterBestMCPDGCode,
1252 "Returns the PDG code of the MCParticle for the ECLCluster -> MCParticle relation with the largest weight.");
1253 REGISTER_VARIABLE("clusterTotalMCMatchWeight", particleClusterTotalMCMatchWeight,
1254 "Returns the sum of all weights of the ECLCluster -> MCParticles relations.");
1255
1256 REGISTER_VARIABLE("clusterTotalMCMatchWeightForKlong", particleClusterTotalMCMatchWeightForKlong,
1257 "Returns the sum of all weights of the ECLCluster -> MCParticles relations when MCParticle is a Klong or daughter of a Klong");
1258 REGISTER_VARIABLE("clusterTotalMCMatchWeightForBestKlong", particleClusterTotalMCMatchWeightForBestKlong,
1259 "Returns the sum of all weights of the ECLCluster -> MCParticles relations when MCParticle is the same Klong or daughter of the Klong. If multiple MC Klongs are related to the ECLCluster, returns the sum of weights for the best matched Klong.");
1260
1261
1262 }
1264}
double R
typedef autogenerated by FFTW
void SetX(DataType x)
set X/1st-coordinate
Definition: B2Vector3.h:457
void SetY(DataType y)
set Y/2nd-coordinate
Definition: B2Vector3.h:459
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:618
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
@ c_IsFSRPhoton
bit 7: Particle is from finial state radiation
Definition: MCParticle.h:61
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition: MCParticle.h:57
@ c_IsPHOTOSPhoton
bit 8: Particle is an radiative photon from PHOTOS
Definition: MCParticle.h:63
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:47
@ c_IsVirtual
bit 4: Particle is virtual and not going to Geant4.
Definition: MCParticle.h:55
@ c_IsISRPhoton
bit 6: Particle is from initial state radiation
Definition: MCParticle.h:59
const MCParticle * getMCParticle() const
Returns the pointer to the MCParticle object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:942
static const ReferenceFrame & GetCurrent()
Get current rest frame.
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Abstract base class for different kinds of events.
static int countMissingParticle(const Belle2::Particle *particle, const Belle2::MCParticle *mcParticle, const std::vector< int > &daughterPDG)
Count the number of missing daughters of the 'particle'.
Definition: MCMatching.cc:561
@ c_Correct
This Particle and all its daughters are perfectly reconstructed.
Definition: MCMatching.h:34
@ c_MisID
One of the charged final state particles is mis-identified, i.e.
Definition: MCMatching.h:42
static int getMCErrors(const Belle2::Particle *particle, const Belle2::MCParticle *mcParticle=nullptr)
Returns quality indicator of the match as a bit pattern where the individual bits indicate the the ty...
Definition: MCMatching.cc:282
static void fillGenMothers(const Belle2::MCParticle *mcP, std::vector< int > &genMCPMothers)
Fills vector with array (1-based) indices of all generator ancestors of given MCParticle.
Definition: MCMatching.cc:61