Belle II Software light-2406-ragdoll
MCParticle.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#include <mdst/dataobjects/MCParticle.h>
10#include <framework/logging/Logger.h>
11#include <framework/datastore/DataStore.h>
12#include <framework/datastore/StoreArray.h>
13#include <framework/utilities/HTML.h>
14#include <framework/utilities/Conversion.h>
15
16#include <TDatabasePDG.h>
17
18#include <iostream>
19#include <sstream>
20#include <boost/algorithm/string.hpp>
21
22using namespace std;
23using namespace Belle2;
24
25const double MCParticle::c_epsilon = 10e-7;
26
27
29{
30 if (TDatabasePDG::Instance()->GetParticle(m_pdg) == nullptr)
31 throw(ParticlePDGNotKnownError() << m_pdg);
32 m_mass = TDatabasePDG::Instance()->GetParticle(m_pdg)->Mass();
33}
34
35
37{
38 // Geant4 "optical photon" (m_pdg == 0) is not known to TDatabasePDG::Instance().
39 if (m_pdg == 0) {
40 return 0.0;
41 }
42
43 if (TDatabasePDG::Instance()->GetParticle(m_pdg) == nullptr) {
44 B2ERROR("PDG=" << m_pdg << " ***code unknown to TDatabasePDG");
45 return 0.0;
46 }
47
48 return TDatabasePDG::Instance()->GetParticle(m_pdg)->Charge() / 3.0;
49}
50
51
52vector<MCParticle*> MCParticle::getDaughters() const
53{
54 vector<MCParticle*> result;
55 if (m_firstDaughter > 0) {
57 if (m_lastDaughter > m_plist->GetEntriesFast()) throw LastChildIndexOutOfRangError();
58 TClonesArray& plist = *m_plist;
59 result.reserve(m_lastDaughter - m_firstDaughter + 1);
60 for (int i = m_firstDaughter - 1; i < m_lastDaughter; i++) {
61 result.push_back(static_cast<MCParticle*>(plist[i]));
62 }
63 }
64 return result;
65}
66
68{
69 if (i >= getNDaughters()) {
70 return nullptr;
71 }
72 return getDaughters().at(i);
73}
74
76{
77 if (getFirstDaughter() == 0) //no daughters
78 return 0;
79 return getLastDaughter() - getFirstDaughter() + 1;
80}
81
83{
84 if (m_plist != 0) return;
85
86 TClonesArray* plist(0);
87
88 //Search default location
89 //TODO: this could be replaced with RelationsObject::getArrayIndex()/getArrayName()
90 StoreArray<MCParticle> MCParticles;
91 if (MCParticles && MCParticles.getPtr()->IndexOf(this) >= 0) {
92 plist = MCParticles.getPtr();
93 } else {
94 //Search all StoreArrays which happen to store MCParticles
96 for (DataStore::StoreEntryConstIter iter = map.begin(); iter != map.end(); ++iter) {
97 TClonesArray* value = dynamic_cast<TClonesArray*>(iter->second.ptr);
98 if (value && value->GetClass() == Class() && value->IndexOf(this) >= 0) {
99 plist = value;
100 break;
101 }
102 }
103 }
104 //Could not find any collection, raise exception
105 if (!plist) {
106 B2ERROR("Could not determine StoreArray the MCParticle belongs to !");
107 throw NoParticleListSetError();
108 }
109
110 //Set plist pointer and index for whole array
111 for (int i = 0; i < plist->GetEntriesFast(); i++) {
112 MCParticle& mc = *(static_cast<MCParticle*>(plist->At(i)));
113 mc.m_plist = plist;
114 mc.m_index = i + 1;
115 }
116}
117std::string MCParticle::getName() const
118{
119 const TParticlePDG* p = TDatabasePDG::Instance()->GetParticle(m_pdg);
120 if (p)
121 return p->GetName();
122 else //handle unknown PDG codes
123 return std::to_string(m_pdg);
124}
125std::string MCParticle::getInfoHTML() const
126{
127 std::stringstream out;
128 out << "<b>Charge</b>=" << (int)getCharge();
129 out << ", <b>PDG</b>=" << getPDG();
130 out << " (" << getName() << ")";
131 out << "<br>";
132 out << "<b>isPrimaryParticle</b>=" << isPrimaryParticle();
133 out << ",<b>isInitial</b>=" << isInitial();
134 out << ",<b>isVirtual</b>=" << isVirtual();
135 out << "<br>";
136
137 out << "<b>pT</b>=" << getMomentum().Rho();
138 out << ", <b>pZ</b>=" << m_momentum_z;
139 out << "<br>";
140 std::string unitType = HTML::chooseUnitOfLength(getProductionVertex());
141 int precision = 3;
142 out << "<b>V</b>=" << HTML::getStringConvertToUnit(getProductionVertex(), precision, unitType);
143
144 const MCParticle* mom = getMother();
145 if (mom) {
146 out << "<br>";
147 out << "<b>Mother</b>: " << mom->getArrayName() << "[" << mom->getArrayIndex() << "] (" << mom->getName() << ")";
148 }
149 return out.str();
150}
151
152const MCParticle* MCParticle::getParticleFromGeneralizedIndexString(const std::string& generalizedIndex) const
153{
154 // Split the generalizedIndex string in a vector of strings.
155 std::vector<std::string> generalizedIndexes;
156 boost::split(generalizedIndexes, generalizedIndex, boost::is_any_of(":"));
157
158 if (generalizedIndexes.empty()) {
159 B2WARNING("Generalized index of MC daughter particle is empty. Skipping.");
160 return nullptr;
161 }
162
163 // To explore a decay tree of unknown depth, we need a place to store
164 // both the root particle and the daughter particle at each iteration
165 const MCParticle* dauPart =
166 nullptr; // This will be eventually returned
167 const MCParticle* currentPart = this; // This is the root particle of the next iteration
168
169 // Loop over the generalizedIndexes until you get to the particle you want
170 for (auto& indexString : generalizedIndexes) {
171 // indexString is a string. First try to convert it into an int
172 int dauIndex = 0;
173 try {
174 dauIndex = Belle2::convertString<int>(indexString);
175 } catch (std::invalid_argument&) {
176 B2WARNING("Found the string " << indexString << "instead of a daughter index.");
177 return nullptr;
178 }
179
180 // Check that the daughter index is smaller than the number of daughters of the current root particle
181 if (dauIndex >= int(currentPart->getNDaughters()) or dauIndex < 0) {
182 B2WARNING("Daughter index out of range" << LogVar("daughter index", dauIndex));
183 B2WARNING("Trying to access non-existing particle.");
184 return nullptr;
185 } else {
186 dauPart = currentPart->getDaughter(dauIndex); // Pick the particle indicated by the generalizedIndex
187 currentPart = dauPart;
188 }
189 }
190 return dauPart;
191}
StoreEntryMap::const_iterator StoreEntryConstIter
const_iterator for a StoreEntry map.
Definition: DataStore.h:89
StoreEntryMap & getStoreEntryMap(EDurability durability)
Get a reference to the object/array map.
Definition: DataStore.h:325
@ c_Event
Different object in each event, all objects/arrays are invalidated after event() function has been ca...
Definition: DataStore.h:59
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:54
std::map< std::string, StoreEntry > StoreEntryMap
Map for StoreEntries.
Definition: DataStore.h:87
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
int m_lastDaughter
1-based index of last daughter particle in collection, 0 if no daughters
Definition: MCParticle.h:562
int m_firstDaughter
1-based index of first daughter particle in collection, 0 if no daughters
Definition: MCParticle.h:561
float m_mass
mass of the particle
Definition: MCParticle.h:542
virtual std::string getName() const override
Return name of this particle.
Definition: MCParticle.cc:117
const MCParticle * getParticleFromGeneralizedIndexString(const std::string &generalizedIndex) const
Explores the decay tree of the MC particle and returns the (grand^n)daughter identified by a generali...
Definition: MCParticle.cc:152
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition: MCParticle.cc:52
static const double c_epsilon
limit of precision for two doubles to be the same.
Definition: MCParticle.h:563
int getArrayIndex() const
Get 0-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:244
ROOT::Math::XYZVector getProductionVertex() const
Return production vertex position.
Definition: MCParticle.h:189
virtual std::string getInfoHTML() const override
Return a short summary of this object's contents in HTML format.
Definition: MCParticle.cc:125
int m_pdg
PDG-Code of the particle.
Definition: MCParticle.h:541
const MCParticle * getDaughter(int i) const
Return i-th daughter.
Definition: MCParticle.cc:67
int getNDaughters() const
Return number of daughter MCParticles.
Definition: MCParticle.cc:75
int getLastDaughter() const
Get 1-based index of last daughter, 0 if no daughters.
Definition: MCParticle.h:259
float getCharge() const
Return the particle charge defined in TDatabasePDG.
Definition: MCParticle.cc:36
float m_momentum_z
momentum of particle, z component
Definition: MCParticle.h:546
void fixParticleList() const
Search the DataStore for the corresponding MCParticle array.
Definition: MCParticle.cc:82
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:112
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition: MCParticle.h:198
TClonesArray * m_plist
Internal pointer to DataStore Array containing particles belonging to this collection.
Definition: MCParticle.h:532
int getFirstDaughter() const
Get 1-based index of first daughter, 0 if no daughters.
Definition: MCParticle.h:251
void setMassFromPDG()
Sets the mass for the particle from the particle's PDG code.
Definition: MCParticle.cc:28
std::string getArrayName() const
Get name of array this object is stored in, or "" if not found.
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
TClonesArray * getPtr() const
Raw access to the underlying TClonesArray.
Definition: StoreArray.h:311
Class to store variables with their name which were sent to the logging service.
bool isInitial() const
Check if particle is an initial particle such as ISR.
Definition: MCParticle.h:590
MCParticle * getMother() const
Returns a pointer to the mother particle.
Definition: MCParticle.h:600
bool isPrimaryParticle() const
Check if particle is a primary particle which was created by the generator (and not,...
Definition: MCParticle.h:595
bool isVirtual() const
Check if particle is virtual.
Definition: MCParticle.h:575
std::string getStringConvertToUnit(const ROOT::Math::XYZVector &vec, int precision=2, const std::string &unitType="cm")
get a string with vector coordinates: (x, y, z).
Definition: HTML.cc:85
std::string chooseUnitOfLength(const ROOT::Math::XYZVector &vec)
get a string with a unit type to convert a vector, so that it is easily readable.
Definition: HTML.cc:102
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.