Belle II Software  light-2212-foldex
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 include
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 
20 #include <mdst/dataobjects/MCParticle.h>
21 #include <mdst/dataobjects/ECLCluster.h>
22 
23 #include <framework/datastore/StoreArray.h>
24 #include <framework/datastore/StoreObjPtr.h>
25 #include <framework/dataobjects/EventMetaData.h>
26 #include <framework/gearbox/Const.h>
27 #include <framework/logging/Logger.h>
28 #include <framework/database/DBObjPtr.h>
29 #include <framework/dbobjects/BeamParameters.h>
30 
31 #include <queue>
32 
33 namespace Belle2 {
38  namespace Variable {
39 
40  static const double realNaN = std::numeric_limits<double>::quiet_NaN();
41 
42 
43  double isSignal(const Particle* part)
44  {
45  const MCParticle* mcparticle = part->getMCParticle();
46  if (!mcparticle) return realNaN;
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 realNaN;
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 realNaN;
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 realNaN;
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 realNaN;
90  // get mcparticle weight (mcmatch weight)
91  const auto mcpww = particle->getRelatedToWithWeight<MCParticle>();
92  if (!mcpww.first) return realNaN;
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 genMotherPDG(const Particle* part)
145  {
146  return genNthMotherPDG(part, {});
147  }
148 
149  double genMotherP(const Particle* part)
150  {
151  const MCParticle* mcparticle = part->getMCParticle();
152  if (!mcparticle) return realNaN;
153 
154  const MCParticle* mcmother = mcparticle->getMother();
155  if (!mcmother) return realNaN;
156 
157  return mcmother->getMomentum().R();
158  }
159 
160  double genMotherIndex(const Particle* part)
161  {
162  return genNthMotherIndex(part, {});
163  }
164 
165  double genParticleIndex(const Particle* part)
166  {
167  const MCParticle* mcparticle = part->getMCParticle();
168  if (!mcparticle) return realNaN;
169  return mcparticle->getArrayIndex();
170  }
171 
172  double isSignalAcceptMissingNeutrino(const Particle* part)
173  {
174  const MCParticle* mcparticle = part->getMCParticle();
175  if (!mcparticle) return realNaN;
176 
177  int status = MCMatching::getMCErrors(part, mcparticle);
178  //remove the following bits
179  status &= (~MCMatching::c_MissNeutrino);
180 
181  return (status == MCMatching::c_Correct);
182  }
183 
184  double isSignalAcceptMissingMassive(const Particle* part)
185  {
186  const MCParticle* mcparticle = part->getMCParticle();
187  if (!mcparticle) return realNaN;
188 
189  int status = MCMatching::getMCErrors(part, mcparticle);
190  //remove the following bits
191  status &= (~MCMatching::c_MissMassiveParticle);
192  status &= (~MCMatching::c_MissKlong);
193 
194  return (status == MCMatching::c_Correct);
195  }
196 
197  double isSignalAcceptMissingGamma(const Particle* part)
198  {
199  const MCParticle* mcparticle = part->getMCParticle();
200  if (!mcparticle) return realNaN;
201 
202  int status = MCMatching::getMCErrors(part, mcparticle);
203  //remove the following bits
204  status &= (~MCMatching::c_MissGamma);
205 
206  return (status == MCMatching::c_Correct);
207  }
208 
209  double isSignalAcceptMissing(const Particle* part)
210  {
211  const MCParticle* mcparticle = part->getMCParticle();
212  if (!mcparticle) return realNaN;
213 
214  int status = MCMatching::getMCErrors(part, mcparticle);
215  //remove the following bits
216  status &= (~MCMatching::c_MissGamma);
217  status &= (~MCMatching::c_MissMassiveParticle);
218  status &= (~MCMatching::c_MissKlong);
219  status &= (~MCMatching::c_MissNeutrino);
220 
221  return (status == MCMatching::c_Correct);
222  }
223 
224  double isSignalAcceptBremsPhotons(const Particle* part)
225  {
226  const MCParticle* mcparticle = part->getMCParticle();
227  if (!mcparticle) return realNaN;
228 
229  int status = MCMatching::getMCErrors(part, mcparticle);
230  //remove the following bits
231  status &= (~MCMatching::c_AddedRecoBremsPhoton);
232 
233  return (status == MCMatching::c_Correct);
234  }
235 
236  double particleMCMatchPDGCode(const Particle* part)
237  {
238  const MCParticle* mcparticle = part->getMCParticle();
239  if (!mcparticle) return realNaN;
240  return mcparticle->getPDG();
241  }
242 
243  double particleMCErrors(const Particle* part)
244  {
245  return MCMatching::getMCErrors(part);
246  }
247 
248  double particleNumberOfMCMatch(const Particle* particle)
249  {
250  RelationVector<MCParticle> mcRelations = particle->getRelationsTo<MCParticle>();
251  return (mcRelations.size());
252  }
253 
254  double particleMCMatchWeight(const Particle* particle)
255  {
256  auto relWithWeight = particle->getRelatedToWithWeight<MCParticle>();
257  if (!relWithWeight.first) return realNaN;
258  return relWithWeight.second;
259  }
260 
261  double particleMCMatchDecayTime(const Particle* part)
262  {
263  const MCParticle* mcparticle = part->getMCParticle();
264  if (!mcparticle) return realNaN;
265  return mcparticle->getDecayTime();
266  }
267 
268  double particleMCMatchLifeTime(const Particle* part)
269  {
270  const MCParticle* mcparticle = part->getMCParticle();
271  if (!mcparticle) return realNaN;
272  return mcparticle->getLifetime();
273  }
274 
275  double particleMCMatchPX(const Particle* part)
276  {
277  const MCParticle* mcparticle = part->getMCParticle();
278  if (!mcparticle) return realNaN;
279 
280  const auto& frame = ReferenceFrame::GetCurrent();
281  ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
282  return frame.getMomentum(mcpP4).Px();
283  }
284 
285  double particleMCMatchPY(const Particle* part)
286  {
287  const MCParticle* mcparticle = part->getMCParticle();
288  if (!mcparticle) return realNaN;
289 
290  const auto& frame = ReferenceFrame::GetCurrent();
291  ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
292  return frame.getMomentum(mcpP4).Py();
293  }
294 
295  double particleMCMatchPZ(const Particle* part)
296  {
297  const MCParticle* mcparticle = part->getMCParticle();
298  if (!mcparticle) return realNaN;
299 
300  const auto& frame = ReferenceFrame::GetCurrent();
301  ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
302  return frame.getMomentum(mcpP4).Pz();
303  }
304 
305  double particleMCMatchPT(const Particle* part)
306  {
307  const MCParticle* mcparticle = part->getMCParticle();
308  if (!mcparticle) return realNaN;
309 
310  const auto& frame = ReferenceFrame::GetCurrent();
311  ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
312  return frame.getMomentum(mcpP4).Pt();
313  }
314 
315  double particleMCMatchE(const Particle* part)
316  {
317  const MCParticle* mcparticle = part->getMCParticle();
318  if (!mcparticle) return realNaN;
319 
320  const auto& frame = ReferenceFrame::GetCurrent();
321  ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
322  return frame.getMomentum(mcpP4).E();
323  }
324 
325  double particleMCMatchP(const Particle* part)
326  {
327  const MCParticle* mcparticle = part->getMCParticle();
328  if (!mcparticle) return realNaN;
329 
330  const auto& frame = ReferenceFrame::GetCurrent();
331  ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
332  return frame.getMomentum(mcpP4).P();
333  }
334 
335  double particleMCMatchTheta(const Particle* part)
336  {
337  const MCParticle* mcparticle = part->getMCParticle();
338  if (!mcparticle) return realNaN;
339 
340  const auto& frame = ReferenceFrame::GetCurrent();
341  ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
342  return frame.getMomentum(mcpP4).Theta();
343  }
344 
345  double particleMCMatchPhi(const Particle* part)
346  {
347  const MCParticle* mcparticle = part->getMCParticle();
348  if (!mcparticle) return realNaN;
349 
350  const auto& frame = ReferenceFrame::GetCurrent();
351  ROOT::Math::PxPyPzEVector mcpP4 = mcparticle->get4Vector();
352  return frame.getMomentum(mcpP4).Phi();
353  }
354 
355  double mcParticleNDaughters(const Particle* part)
356  {
357  const MCParticle* mcparticle = part->getMCParticle();
358 
359  if (!mcparticle) return realNaN;
360  return mcparticle->getNDaughters();
361  }
362 
363  double particleMCRecoilMass(const Particle* part)
364  {
365  StoreArray<MCParticle> mcparticles;
366  if (mcparticles.getEntries() < 1) return realNaN;
367 
368  ROOT::Math::PxPyPzEVector pInitial = mcparticles[0]->get4Vector();
369  ROOT::Math::PxPyPzEVector pDaughters;
370  const std::vector<Particle*> daughters = part->getDaughters();
371  for (auto daughter : daughters) {
372  const MCParticle* mcD = daughter->getMCParticle();
373  if (!mcD) return realNaN;
374 
375  pDaughters += mcD->get4Vector();
376  }
377  return (pInitial - pDaughters).M();
378  }
379 
380  ROOT::Math::PxPyPzEVector MCInvisibleP4(const MCParticle* mcparticle)
381  {
382  ROOT::Math::PxPyPzEVector ResultP4;
383  int pdg = abs(mcparticle->getPDG());
384  bool isNeutrino = (pdg == 12 or pdg == 14 or pdg == 16);
385 
386  if (mcparticle->getNDaughters() > 0) {
387  const std::vector<MCParticle*> daughters = mcparticle->getDaughters();
388  for (auto daughter : daughters)
389  ResultP4 += MCInvisibleP4(daughter);
390  } else if (isNeutrino)
391  ResultP4 += mcparticle->get4Vector();
392 
393  return ResultP4;
394  }
395 
396  double particleMCCosThetaBetweenParticleAndNominalB(const Particle* part)
397  {
398  int particlePDG = abs(part->getPDGCode());
399  if (particlePDG != 511 and particlePDG != 521)
400  B2FATAL("The variable mcCosThetaBetweenParticleAndNominalB is only meant to be used on B mesons!");
401 
402  PCmsLabTransform T;
403  double e_Beam = T.getCMSEnergy() / 2.0; // GeV
404  double m_B = part->getPDGMass();
405 
406  // Y(4S) mass according PDG (https://pdg.lbl.gov/2020/listings/rpp2020-list-upsilon-4S.pdf)
407  const double mY4S = 10.5794; // GeV
408 
409  // if this is a continuum run, use an approximate Y(4S) CMS energy
410  if (e_Beam * e_Beam - m_B * m_B < 0) {
411  e_Beam = mY4S / 2.0;
412  }
413  double p_B = std::sqrt(e_Beam * e_Beam - m_B * m_B);
414 
415  // Calculate cosThetaBY with daughter neutrino momenta subtracted
416  const MCParticle* mcB = part->getMCParticle();
417  if (!mcB) return realNaN;
418 
419  int mcParticlePDG = abs(mcB->getPDG());
420  if (mcParticlePDG != 511 and mcParticlePDG != 521)
421  return realNaN;
422 
423  ROOT::Math::PxPyPzEVector p = T.rotateLabToCms() * (mcB->get4Vector() - MCInvisibleP4(mcB));
424  double e_d = p.E();
425  double m_d = p.M();
426  double p_d = p.P();
427 
428  double theta_BY = (2 * e_Beam * e_d - m_B * m_B - m_d * m_d)
429  / (2 * p_B * p_d);
430  return theta_BY;
431  }
432 
433  double mcParticleSecondaryPhysicsProcess(const Particle* p)
434  {
435  const MCParticle* mcp = p->getMCParticle();
436  if (!mcp) return realNaN;
437  return mcp->getSecondaryPhysicsProcess();
438  }
439 
440  double mcParticleStatus(const Particle* p)
441  {
442  const MCParticle* mcp = p->getMCParticle();
443  if (!mcp) return realNaN;
444  return mcp->getStatus();
445  }
446 
447  double particleMCPrimaryParticle(const Particle* p)
448  {
449  const MCParticle* mcp = p->getMCParticle();
450  if (!mcp) return realNaN;
451 
452  unsigned int bitmask = MCParticle::c_PrimaryParticle;
453  return mcp->hasStatus(bitmask);
454  }
455 
456  double particleMCVirtualParticle(const Particle* p)
457  {
458  const MCParticle* mcp = p->getMCParticle();
459  if (!mcp) return realNaN;
460 
461  unsigned int bitmask = MCParticle::c_IsVirtual;
462  return mcp->hasStatus(bitmask);
463  }
464 
465  double particleMCInitialParticle(const Particle* p)
466  {
467  const MCParticle* mcp = p->getMCParticle();
468  if (!mcp) return realNaN;
469 
470  unsigned int bitmask = MCParticle::c_Initial;
471  return mcp->hasStatus(bitmask);
472  }
473 
474  double particleMCISRParticle(const Particle* p)
475  {
476  const MCParticle* mcp = p->getMCParticle();
477  if (!mcp) return realNaN;
478 
479  unsigned int bitmask = MCParticle::c_IsISRPhoton;
480  return mcp->hasStatus(bitmask);
481  }
482 
483  double particleMCFSRParticle(const Particle* p)
484  {
485  const MCParticle* mcp = p->getMCParticle();
486  if (!mcp) return realNaN;
487 
488  unsigned int bitmask = MCParticle::c_IsFSRPhoton;
489  return mcp->hasStatus(bitmask);
490  }
491 
492  double particleMCPhotosParticle(const Particle* p)
493  {
494  const MCParticle* mcp = p->getMCParticle();
495  if (!mcp) return realNaN;
496 
497  unsigned int bitmask = MCParticle::c_IsPHOTOSPhoton;
498  return mcp->hasStatus(bitmask);
499  }
500 
501  double generatorEventWeight(const Particle*)
502  {
503  StoreObjPtr<EventMetaData> evtMetaData;
504  if (!evtMetaData) return realNaN;
505  return evtMetaData->getGeneratedWeight();
506  }
507 
508  int tauPlusMcMode(const Particle*)
509  {
510  StoreObjPtr<TauPairDecay> tauDecay;
511  if (!tauDecay) {
512  B2WARNING("Cannot find tau decay ID, did you forget to run TauDecayMarkerModule?");
513  return 0;
514  }
515  return tauDecay->getTauPlusIdMode();
516  }
517 
518  int tauMinusMcMode(const Particle*)
519  {
520  StoreObjPtr<TauPairDecay> tauDecay;
521  if (!tauDecay) {
522  B2WARNING("Cannot find tau decay ID, did you forget to run TauDecayMarkerModule?");
523  return 0;
524  }
525  return tauDecay->getTauMinusIdMode();
526  }
527 
528  int tauPlusMcProng(const Particle*)
529  {
530  StoreObjPtr<TauPairDecay> tauDecay;
531  if (!tauDecay) {
532  B2WARNING("Cannot find tau prong, did you forget to run TauDecayMarkerModule?");
533  return 0;
534  }
535  return tauDecay->getTauPlusMcProng();
536  }
537 
538  int tauMinusMcProng(const Particle*)
539  {
540  StoreObjPtr<TauPairDecay> tauDecay;
541  if (!tauDecay) {
542  B2WARNING("Cannot find tau prong, did you forget to run TauDecayMarkerModule?");
543  return 0;
544  }
545  return tauDecay->getTauMinusMcProng();
546  }
547 
548  double isReconstructible(const Particle* p)
549  {
550  if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
551  return realNaN;
552  const MCParticle* mcp = p->getMCParticle();
553  if (!mcp) return realNaN;
554 
555  // If charged: make sure it was seen in the SVD.
556  // If neutral: make sure it was seen in the ECL.
557  return (abs(mcp->getCharge()) > 0) ? seenInSVD(p) : seenInECL(p);
558  }
559 
560  double seenInPXD(const Particle* p)
561  {
562  if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
563  return realNaN;
564  const MCParticle* mcp = p->getMCParticle();
565  if (!mcp) return realNaN;
566  return mcp->hasSeenInDetector(Const::PXD);
567  }
568 
569  double seenInSVD(const Particle* p)
570  {
571  if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
572  return realNaN;
573  const MCParticle* mcp = p->getMCParticle();
574  if (!mcp) return realNaN;
575  return mcp->hasSeenInDetector(Const::SVD);
576  }
577 
578  double seenInCDC(const Particle* p)
579  {
580  if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
581  return realNaN;
582  const MCParticle* mcp = p->getMCParticle();
583  if (!mcp) return realNaN;
584  return mcp->hasSeenInDetector(Const::CDC);
585  }
586 
587  double seenInTOP(const Particle* p)
588  {
589  if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
590  return realNaN;
591  const MCParticle* mcp = p->getMCParticle();
592  if (!mcp) return realNaN;
593  return mcp->hasSeenInDetector(Const::TOP);
594  }
595 
596  double seenInECL(const Particle* p)
597  {
598  if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
599  return realNaN;
600  const MCParticle* mcp = p->getMCParticle();
601  if (!mcp) return realNaN;
602  return mcp->hasSeenInDetector(Const::ECL);
603  }
604 
605  double seenInARICH(const Particle* p)
606  {
607  if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
608  return realNaN;
609  const MCParticle* mcp = p->getMCParticle();
610  if (!mcp) return realNaN;
611  return mcp->hasSeenInDetector(Const::ARICH);
612  }
613 
614  double seenInKLM(const Particle* p)
615  {
616  if (p->getParticleSource() == Particle::EParticleSourceObject::c_Composite)
617  return realNaN;
618  const MCParticle* mcp = p->getMCParticle();
619  if (!mcp) return realNaN;
620  return mcp->hasSeenInDetector(Const::KLM);
621  }
622 
623  int genNStepsToDaughter(const Particle* p, const std::vector<double>& arguments)
624  {
625  if (arguments.size() != 1)
626  B2FATAL("Wrong number of arguments for genNStepsToDaughter");
627 
628  const MCParticle* mcp = p->getMCParticle();
629  if (!mcp) {
630  B2WARNING("No MCParticle is associated to the particle");
631  return 0;
632  }
633 
634  int nChildren = p->getNDaughters();
635  if (arguments[0] >= nChildren) {
636  return 0;
637  }
638 
639  const Particle* daugP = p->getDaughter(arguments[0]);
640  const MCParticle* daugMCP = daugP->getMCParticle();
641  if (!daugMCP) {
642  // This is a strange case.
643  // The particle, p, has the related MC particle, but i-th daughter does not have the related MC Particle.
644  B2WARNING("No MCParticle is associated to the i-th daughter");
645  return 0;
646  }
647 
648  if (nChildren == 1) return 1;
649 
650  std::vector<int> genMothers;
651  MCMatching::fillGenMothers(daugMCP, genMothers);
652  auto match = std::find(genMothers.begin(), genMothers.end(), mcp->getIndex());
653  return match - genMothers.begin();
654  }
655 
656  int genNMissingDaughter(const Particle* p, const std::vector<double>& arguments)
657  {
658  if (arguments.size() < 1)
659  B2FATAL("Wrong number of arguments for genNMissingDaughter");
660 
661  const std::vector<int> PDGcodes(arguments.begin(), arguments.end());
662 
663  const MCParticle* mcp = p->getMCParticle();
664  if (!mcp) {
665  B2WARNING("No MCParticle is associated to the particle");
666  return 0;
667  }
668 
669  return MCMatching::countMissingParticle(p, mcp, PDGcodes);
670  }
671 
672  double getHEREnergy(const Particle*)
673  {
674  static DBObjPtr<BeamParameters> beamParamsDB;
675  return (beamParamsDB->getHER()).E();
676  }
677 
678  double getLEREnergy(const Particle*)
679  {
680  static DBObjPtr<BeamParameters> beamParamsDB;
681  return (beamParamsDB->getLER()).E();
682  }
683 
684  double getCrossingAngleX(const Particle*)
685  {
686  // get the beam momenta from the DB
687  static DBObjPtr<BeamParameters> beamParamsDB;
688  B2Vector3D herVec = beamParamsDB->getHER().Vect();
689  B2Vector3D lerVec = beamParamsDB->getLER().Vect();
690 
691  // only looking at the horizontal (XZ plane) -> set y-coordinates to zero
692  herVec.SetY(0);
693  lerVec.SetY(0);
694 
695  //calculate the crossing angle
696  return herVec.Angle(-lerVec);
697  }
698 
699  double getCrossingAngleY(const Particle*)
700  {
701  // get the beam momenta from the DB
702  static DBObjPtr<BeamParameters> beamParamsDB;
703  B2Vector3D herVec = beamParamsDB->getHER().Vect();
704  B2Vector3D lerVec = beamParamsDB->getLER().Vect();
705 
706  // only looking at the vertical (YZ plane) -> set x-coordinates to zero
707  herVec.SetX(0);
708  lerVec.SetX(0);
709 
710  //calculate the crossing angle
711  return herVec.Angle(-lerVec);
712  }
713 
714 
715  double particleClusterMatchWeight(const Particle* particle)
716  {
717  /* Get the weight of the *cluster* mc match for the mcparticle matched to
718  * this particle.
719  *
720  * Note that for track-based particles this is different from the mc match
721  * of the particle (which it inherits from the mc match of the track)
722  */
723  const MCParticle* matchedToParticle = particle->getMCParticle();
724  if (!matchedToParticle) return realNaN;
725  int matchedToIndex = matchedToParticle->getArrayIndex();
726 
727  const ECLCluster* cluster = particle->getECLCluster();
728  if (!cluster) return realNaN;
729 
730  const auto mcps = cluster->getRelationsTo<MCParticle>();
731  for (unsigned int i = 0; i < mcps.size(); ++i)
732  if (mcps[i]->getArrayIndex() == matchedToIndex)
733  return mcps.weight(i);
734 
735  return realNaN;
736  }
737 
738  double particleClusterBestMCMatchWeight(const Particle* particle)
739  {
740  /* Get the weight of the best mc match of the cluster associated to
741  * this particle.
742  *
743  * Note for electrons (or any track-based particle) this may not be
744  * the same thing as the mc match of the particle (which is taken
745  * from the track).
746  *
747  * For photons (or any ECL-based particle) this will be the same as the
748  * mcMatchWeight
749  */
750  const ECLCluster* cluster = particle->getECLCluster();
751  if (!cluster) return realNaN;
752 
753  /* loop over all mcparticles related to this cluster, find the largest
754  * weight by std::sort-ing the doubles
755  */
756  auto mcps = cluster->getRelationsTo<MCParticle>();
757  if (mcps.size() == 0) return realNaN;
758 
759  std::vector<double> weights;
760  for (unsigned int i = 0; i < mcps.size(); ++i)
761  weights.emplace_back(mcps.weight(i));
762 
763  // sort descending by weight
764  std::sort(weights.begin(), weights.end());
765  std::reverse(weights.begin(), weights.end());
766  return weights[0];
767  }
768 
769  double particleClusterBestMCPDGCode(const Particle* particle)
770  {
771  /* Get the PDG code of the best mc match of the cluster associated to this
772  * particle.
773  *
774  * Note for electrons (or any track-based particle) this may not be the
775  * same thing as the mc match of the particle (which is taken from the track).
776  *
777  * For photons (or any ECL-based particle) this will be the same as the mcPDG
778  */
779  const ECLCluster* cluster = particle->getECLCluster();
780  if (!cluster) return realNaN;
781 
782  auto mcps = cluster->getRelationsTo<MCParticle>();
783  if (mcps.size() == 0) return realNaN;
784 
785  std::vector<std::pair<double, int>> weightsAndIndices;
786  for (unsigned int i = 0; i < mcps.size(); ++i)
787  weightsAndIndices.emplace_back(mcps.weight(i), i);
788 
789  // sort descending by weight
790  std::sort(
791  weightsAndIndices.begin(), weightsAndIndices.end(),
792  [](const std::pair<double, int>& l, const std::pair<double, int>& r) {
793  return l.first > r.first;
794  });
795  // cppcheck-suppress containerOutOfBounds
796  return mcps.object(weightsAndIndices[0].second)->getPDG();
797  }
798 
799  double particleClusterTotalMCMatchWeight(const Particle* particle)
800  {
801  const ECLCluster* cluster = particle->getECLCluster();
802  if (!cluster) return realNaN;
803 
804  auto mcps = cluster->getRelationsTo<MCParticle>();
805 
806  // if there are no relations to any MCParticles, we return 0!
807  double weightsum = 0;
808  for (unsigned int i = 0; i < mcps.size(); ++i)
809  weightsum += mcps.weight(i);
810 
811  return weightsum;
812  }
813 
814  double isBBCrossfeed(const Particle* particle)
815  {
816  if (particle == nullptr)
817  return std::numeric_limits<double>::quiet_NaN();
818 
819  int pdg = particle->getPDGCode();
820  if (abs(pdg) != 511 && abs(pdg) != 521 && abs(pdg) != 531)
821  return std::numeric_limits<double>::quiet_NaN();
822 
823  std::vector<const Particle*> daughters = particle->getFinalStateDaughters();
824  int nDaughters = daughters.size();
825  if (nDaughters <= 1)
826  return 0;
827  std::vector<int> mother_ids;
828 
829  for (int j = 0; j < nDaughters; ++j) {
830  const MCParticle* curMCParticle = daughters[j]->getMCParticle();
831  while (curMCParticle != nullptr) {
832  pdg = curMCParticle->getPDG();
833  if (abs(pdg) == 511 || abs(pdg) == 521 || abs(pdg) == 531) {
834  mother_ids.emplace_back(curMCParticle->getArrayIndex());
835  break;
836  }
837  const MCParticle* curMCMother = curMCParticle->getMother();
838  curMCParticle = curMCMother;
839  }
840  if (curMCParticle == nullptr) {
841  return std::numeric_limits<double>::quiet_NaN();
842  }
843  }
844 
845  std::set<int> distinctIDs = std::set(mother_ids.begin(), mother_ids.end());
846  if (distinctIDs.size() == 1)
847  return 0;
848  else
849  return 1;
850  }
851 
852  VARIABLE_GROUP("MC matching and MC truth");
853  REGISTER_VARIABLE("isSignal", isSignal,
854  "1.0 if Particle is correctly reconstructed (SIGNAL), 0.0 if not, and NaN if no related MCParticle could be found. \n"
855  "It behaves according to DecayStringGrammar.");
856  REGISTER_VARIABLE("isSignalAcceptWrongFSPs", isSignalAcceptWrongFSPs,
857  "1.0 if Particle is almost correctly reconstructed (SIGNAL), 0.0 if not, and NaN if no related MCParticle could be found.\n"
858  "Misidentification of charged FSP is allowed.");
859  REGISTER_VARIABLE("isPrimarySignal", isPrimarySignal,
860  "1.0 if Particle is correctly reconstructed (SIGNAL) and primary, 0.0 if not, and NaN if no related MCParticle could be found");
861  REGISTER_VARIABLE("isSignalAcceptBremsPhotons", isSignalAcceptBremsPhotons,
862  "1.0 if Particle is correctly reconstructed (SIGNAL), 0.0 if not, and NaN if no related MCParticle could be found.\n"
863  "Particles with gamma daughters attached through the bremsstrahlung recovery modules are allowed.");
864  REGISTER_VARIABLE("genMotherPDG", genMotherPDG,
865  "Check the PDG code of a particles MC mother particle");
866  REGISTER_VARIABLE("genMotherPDG(i)", genNthMotherPDG,
867  "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:");
868 
869  REGISTER_VARIABLE("genMotherID", genMotherIndex,
870  "Check the array index of a particles generated mother");
871  REGISTER_VARIABLE("genMotherID(i)", genNthMotherIndex,
872  "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:");
873  // genMotherPDG and genMotherID are overloaded (each are two C++ functions
874  // sharing one variable name) so one of the two needs to be made the indexed
875  // variable in sphinx
876  REGISTER_VARIABLE("isBBCrossfeed", isBBCrossfeed,
877  "Returns 1 for crossfeed in reconstruction of given B meson, 0 for no crossfeed and NaN for no true B meson or failed truthmatching.");
878  REGISTER_VARIABLE("genMotherP", genMotherP,
879  "Generated momentum of a particles MC mother particle", "GeV/c");
880  REGISTER_VARIABLE("genParticleID", genParticleIndex,
881  "Check the array index of a particle's related MCParticle");
882  REGISTER_VARIABLE("isSignalAcceptMissingNeutrino",
883  isSignalAcceptMissingNeutrino,
884  "Same as isSignal, but also accept missing neutrino");
885  REGISTER_VARIABLE("isSignalAcceptMissingMassive",
886  isSignalAcceptMissingMassive,
887  "Same as isSignal, but also accept missing massive particle");
888  REGISTER_VARIABLE("isSignalAcceptMissingGamma",
889  isSignalAcceptMissingGamma,
890  "Same as isSignal, but also accept missing gamma, such as B -> K* gamma, pi0 -> gamma gamma");
891  REGISTER_VARIABLE("isSignalAcceptMissing",
892  isSignalAcceptMissing,
893  "Same as isSignal, but also accept missing particle");
894  REGISTER_VARIABLE("isMisidentified", isMisidentified,
895  "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.");
896  REGISTER_VARIABLE("isWrongCharge", isWrongCharge,
897  "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.");
898  REGISTER_VARIABLE("isCloneTrack", isCloneTrack,
899  "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)");
900  REGISTER_VARIABLE("isOrHasCloneTrack", isOrHasCloneTrack,
901  "Return 1 if the particle is a clone track or has a clone track as a daughter, 0 otherwise.");
902  REGISTER_VARIABLE("mcPDG", particleMCMatchPDGCode,
903  "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).");
904  REGISTER_VARIABLE("mcErrors", particleMCErrors,
905  "The bit pattern indicating the quality of MC match (see MCMatching::MCErrorFlags)");
906  REGISTER_VARIABLE("mcMatchWeight", particleMCMatchWeight,
907  "The weight of the Particle -> MCParticle relation (only for the first Relation = largest weight).");
908  REGISTER_VARIABLE("nMCMatches", particleNumberOfMCMatch,
909  "The number of relations of this Particle to MCParticle.");
910  REGISTER_VARIABLE("mcDecayTime", particleMCMatchDecayTime,
911  "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).",
912  "ns");
913  REGISTER_VARIABLE("mcLifeTime", particleMCMatchLifeTime,
914  "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).",
915  "ns");
916  REGISTER_VARIABLE("mcPX", particleMCMatchPX,
917  "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).",
918  "GeV/c");
919  REGISTER_VARIABLE("mcPY", particleMCMatchPY,
920  "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).",
921  "GeV/c");
922  REGISTER_VARIABLE("mcPZ", particleMCMatchPZ,
923  "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).",
924  "GeV/c");
925  REGISTER_VARIABLE("mcPT", particleMCMatchPT,
926  "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).",
927  "GeV/c");
928  REGISTER_VARIABLE("mcE", particleMCMatchE,
929  "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).",
930  "GeV");
931  REGISTER_VARIABLE("mcP", particleMCMatchP,
932  "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).",
933  "GeV/c");
934  REGISTER_VARIABLE("mcPhi", particleMCMatchPhi,
935  "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).",
936  "rad");
937  REGISTER_VARIABLE("mcTheta", particleMCMatchTheta,
938  "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).",
939  "rad");
940  REGISTER_VARIABLE("nMCDaughters", mcParticleNDaughters,
941  "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).");
942  REGISTER_VARIABLE("mcRecoilMass", particleMCRecoilMass,
943  "The mass recoiling against the particles attached as particle's daughters calculated using MC truth values.",
944  "GeV/:math:`\\text{c}^2`");
945  REGISTER_VARIABLE("mcCosThetaBetweenParticleAndNominalB",
946  particleMCCosThetaBetweenParticleAndNominalB,
947  "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.");
948 
949 
950  REGISTER_VARIABLE("mcSecPhysProc", mcParticleSecondaryPhysicsProcess,
951  "Returns the secondary physics process flag.");
952  REGISTER_VARIABLE("mcParticleStatus", mcParticleStatus,
953  "Returns status bits of related MCParticle or NaN if MCParticle relation is not set.");
954  REGISTER_VARIABLE("mcPrimary", particleMCPrimaryParticle,
955  "Returns 1 if Particle is related to primary MCParticle, 0 if Particle is related to non - primary MCParticle, "
956  "NaN if Particle is not related to MCParticle.");
957  REGISTER_VARIABLE("mcVirtual", particleMCVirtualParticle,
958  "Returns 1 if Particle is related to virtual MCParticle, 0 if Particle is related to non - virtual MCParticle, "
959  "NaN if Particle is not related to MCParticle.")
960  REGISTER_VARIABLE("mcInitial", particleMCInitialParticle,
961  "Returns 1 if Particle is related to initial MCParticle, 0 if Particle is related to non - initial MCParticle, "
962  "NaN if Particle is not related to MCParticle.")
963  REGISTER_VARIABLE("mcISR", particleMCISRParticle,
964  "Returns 1 if Particle is related to ISR MCParticle, 0 if Particle is related to non - ISR MCParticle, "
965  "NaN if Particle is not related to MCParticle.")
966  REGISTER_VARIABLE("mcFSR", particleMCFSRParticle,
967  "Returns 1 if Particle is related to FSR MCParticle, 0 if Particle is related to non - FSR MCParticle ,"
968  "NaN if Particle is not related to MCParticle.")
969  REGISTER_VARIABLE("mcPhotos", particleMCPhotosParticle,
970  "Returns 1 if Particle is related to Photos MCParticle, 0 if Particle is related to non - Photos MCParticle, "
971  "NaN if Particle is not related to MCParticle.")
972  REGISTER_VARIABLE("generatorEventWeight", generatorEventWeight,
973  "[Eventbased] Returns the event weight produced by the event generator")
974 
975  REGISTER_VARIABLE("genNStepsToDaughter(i)", genNStepsToDaughter,
976  "Returns number of steps to i-th daughter from the particle at generator level. "
977  "NaN if no MCParticle is associated to the particle or i-th daughter. "
978  "NaN if i-th daughter does not exist.");
979  REGISTER_VARIABLE("genNMissingDaughter(PDG)", genNMissingDaughter,
980  "Returns the number of missing daughters having assigned PDG codes. "
981  "NaN if no MCParticle is associated to the particle.");
982  REGISTER_VARIABLE("Eher", getHEREnergy, R"DOC(
983 [Eventbased] The nominal HER energy used by the generator.
984 
985 .. warning:: This variable does not make sense for data.
986 )DOC","GeV");
987  REGISTER_VARIABLE("Eler", getLEREnergy, R"DOC(
988 [Eventbased] The nominal LER energy used by the generator.
989 
990 .. warning:: This variable does not make sense for data.
991 )DOC","GeV");
992  REGISTER_VARIABLE("XAngle", getCrossingAngleX, R"DOC(
993 [Eventbased] The nominal beam crossing angle in the x-z plane from generator level beam kinematics.
994 
995 .. warning:: This variable does not make sense for data.
996 )DOC","rad");
997  REGISTER_VARIABLE("YAngle", getCrossingAngleY, R"DOC(
998 [Eventbased] The nominal beam crossing angle in the y-z plane from generator level beam kinematics.
999 
1000 .. warning:: This variable does not make sense for data.
1001 )DOC","rad");
1002 
1003  VARIABLE_GROUP("Generated tau decay information");
1004  REGISTER_VARIABLE("tauPlusMCMode", tauPlusMcMode,
1005  "[Eventbased] Decay ID for the positive tau lepton in a tau pair generated event.");
1006  REGISTER_VARIABLE("tauMinusMCMode", tauMinusMcMode,
1007  "[Eventbased] Decay ID for the negative tau lepton in a tau pair generated event.");
1008  REGISTER_VARIABLE("tauPlusMCProng", tauPlusMcProng,
1009  "[Eventbased] Prong for the positive tau lepton in a tau pair generated event.");
1010  REGISTER_VARIABLE("tauMinusMCProng", tauMinusMcProng,
1011  "[Eventbased] Prong for the negative tau lepton in a tau pair generated event.");
1012 
1013  VARIABLE_GROUP("MC particle seen in subdetectors");
1014  REGISTER_VARIABLE("isReconstructible", isReconstructible,
1015  "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.");
1016  REGISTER_VARIABLE("seenInPXD", seenInPXD,
1017  "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.");
1018  REGISTER_VARIABLE("seenInSVD", seenInSVD,
1019  "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.");
1020  REGISTER_VARIABLE("seenInCDC", seenInCDC,
1021  "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.");
1022  REGISTER_VARIABLE("seenInTOP", seenInTOP,
1023  "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.");
1024  REGISTER_VARIABLE("seenInECL", seenInECL,
1025  "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.");
1026  REGISTER_VARIABLE("seenInARICH", seenInARICH,
1027  "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.");
1028  REGISTER_VARIABLE("seenInKLM", seenInKLM,
1029  "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.");
1030 
1031  VARIABLE_GROUP("MC Matching for ECLClusters");
1032  REGISTER_VARIABLE("clusterMCMatchWeight", particleClusterMatchWeight,
1033  "Returns the weight of the ECLCluster -> MCParticle relation for the MCParticle matched to the particle. "
1034  "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. "
1035  "Returns -1 if the cluster *was* matched to particles, but not the match of the particle provided.");
1036  REGISTER_VARIABLE("clusterBestMCMatchWeight", particleClusterBestMCMatchWeight,
1037  "Returns the weight of the ECLCluster -> MCParticle relation for the relation with the largest weight.");
1038  REGISTER_VARIABLE("clusterBestMCPDG", particleClusterBestMCPDGCode,
1039  "Returns the PDG code of the MCParticle for the ECLCluster -> MCParticle relation with the largest weight.");
1040  REGISTER_VARIABLE("clusterTotalMCMatchWeight", particleClusterTotalMCMatchWeight,
1041  "Returns the sum of all weights of the ECLCluster -> MCParticles relations.");
1042 
1043  }
1045 }
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:609
@ 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:972
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.
Definition: ClusterUtils.h:23
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:558
@ 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.
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:279
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