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