Belle II Software light-2406-ragdoll
MCParticleGraph.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/MCParticleGraph.h>
10
11#include <framework/datastore/StoreArray.h>
12#include <framework/logging/Logger.h>
13
14#include <boost/graph/graph_traits.hpp>
15#include <boost/graph/adjacency_list.hpp>
16#include <boost/graph/depth_first_search.hpp>
17
18#include <limits>
19#include <vector>
20#include <queue>
21
22
23using namespace std;
24using namespace boost;
25using namespace Belle2;
26
28struct cycle_detector : public dfs_visitor<> {
36 template <class Edge, class Graph> void back_edge(Edge, Graph&) { throw MCParticleGraph::CyclicReferenceError(); }
37};
38
39
50
51public:
52
60 ParticleSorter(MemoryPool<MCParticleGraph::GraphParticle>& particles, TClonesArray* plist, bool setVertex, bool setTime):
61 m_index(0), m_particles(particles), m_plist(plist), m_setVertex(setVertex), m_setTime(setTime) {}
62
68 void setStartIndex(int index) { m_index = index; }
69
74 template <class Graph> void sort(Graph& g)
75 {
76 //Set seen flag for all vertices to false
77 m_seen.clear();
78 m_seen.resize(num_vertices(g), false);
79
80 //Create a dummy GraphParticle, needed only to find all primary particles
82
83 //Add all direct children of the 0th vertex to the queue.
84 //This are the primary particles
85 find_daughters(0, g, dummy);
86
87 //Go through the queue and write out each particle.
88 //Daughters of particles will be added to the queue by find_daughters
89 while (!m_vqueue.empty()) {
90 unsigned int cur = m_vqueue.front();
91 m_vqueue.pop();
92 finish_vertex(cur, g);
93 }
94 }
95
96
102 template <class Vertex, class Graph> void finish_vertex(Vertex v, Graph& g)
103 {
105
106 //Reset daughter information, will be filled by find_daughters
107 p.setFirstDaughter(0);
108 p.setLastDaughter(0);
109 //Find all direct daughters
110 find_daughters(v, g, p);
111
112 //If stable particle, set decaytime to infinity
113 if (out_degree(v, g) == 0 && m_setTime) {
114 p.setDecayTime(numeric_limits<double>::infinity());
115 }
116 //If given a pointer to a TClonesArray, create MCParticle at the appropriate index position
117 if (m_plist) {
118 new (m_plist->AddrAt(p.getIndex() - 1)) MCParticle(m_plist, p);
119 }
120 }
121
122
129 template <class Vertex, class Graph> void find_daughters(Vertex v, Graph& g, MCParticleGraph::GraphParticle& mother)
130 {
131 //References to the daughter information of the mother for easier access
132 int& d1 = mother.m_firstDaughter;
133 int& d2 = mother.m_lastDaughter;
134
135 typename graph_traits<Graph>::out_edge_iterator j, j_end;
136 for (tie(j, j_end) = out_edges(v, g); j != j_end; ++j) {
137 //Get daughter particle from list
138 Vertex nv = target(*j, g);
139 MCParticleGraph::GraphParticle& daughter = *m_particles[nv - 1];
140
141 if (daughter.m_ignore) {
142 //daughter ignored, search its children and treat them as direct children of mother
143 //if we haven't seen this particle yet, set its index to that of its last unignored parent
144 if (!m_seen[nv]) daughter.setIndex(mother.getIndex());
145 find_daughters(nv, g, mother);
146 } else {
147 //If we didn't see this particle already, set its index and add it to the queue for writing out
148 if (!m_seen[nv]) {
149 daughter.setIndex(++m_index);
150 m_vqueue.push(nv);
151 }
152 //Set daughter information of mother. If 0, no daughters yet so just take current daughter as only
153 //daughter. Otherwise allow extension of daughter information in both directions.
154 if (d1 == 0) {
155 d1 = daughter.getIndex();
156 d2 = d1;
157 } else if ((d2 + 1) == daughter.getIndex()) {
158 ++d2;
159 } else if ((d1 - 1) == daughter.getIndex()) {
160 --d1;
161 } else {
162 //Daughter indices are not continous, cannot continue
163 throw MCParticleGraph::NonContinousDaughtersError();
164 }
165 //Set Vertex and time information if requested
166 setVertexTime(mother, daughter);
167 daughter.m_mother = mother.getIndex();
168 }
169 //Mark particle as seen
170 m_seen[nv] = true;
171 }
172 }
173
174
181 {
182 //Only set vertex information if both particles have a valid vertex set
183 m.setValidVertex(m.hasValidVertex() && d.hasValidVertex());
184 if (m.hasValidVertex() && d.getProductionTime() >= m.getDecayTime()) {
185 if (m_setVertex) {
186 m.setDecayVertex(d.getProductionVertex());
187 }
188 if (m_setTime) {
189 m.setDecayTime(d.getProductionTime());
190 }
191 }
192 }
193
194
195protected:
196
199 TClonesArray*
203 vector<bool>
205 std::queue<unsigned int> m_vqueue;
206};
207
208
209void MCParticleGraph::generateList(const string& name, int options)
210{
211 StoreArray<MCParticle> MCParticles(name);
212
213 //Make Graph and connect all primary vertices (particles without mother)
214 //to an artificial 0ths vertex to be able to find them easily
215 typedef adjacency_list<vecS, vecS, directedS> Graph;
216 int num_particles(0);
217 //Determine number of not ignored particles and add an edge from 0ths vertex to any primary
218 //particle
219 for (unsigned int i = 0; i < m_particles.size(); ++i) {
220 if (!m_particles[i]->m_ignore) ++num_particles;
221 if (m_particles[i]->m_primary) m_decays.insert(DecayLine(0, i + 1));
222 }
223 Graph g(m_decays.begin(), m_decays.end(), m_particles.size() + 1);
224
225 //Check for cyclic dependency
226 if (options & c_checkCyclic) {
227 cycle_detector vis;
228 depth_first_search(g, visitor(vis));
229 }
230
231 //Fill TClonesArray in correct order
232 if (options & c_clearParticles) MCParticles.getPtr()->Clear();
233 MCParticles.getPtr()->Expand(num_particles + MCParticles.getEntries());
234 MCParticleGraph::ParticleSorter psorter(m_particles, MCParticles.getPtr(), options & c_setDecayVertex, options & c_setDecayTime);
235 psorter.setStartIndex(MCParticles.getEntries());
236 psorter.sort(g);
237}
238
239void MCParticleGraph::loadList(const string& name)
240{
241 StoreArray<MCParticle> MCParticles(name);
242 if (!MCParticles) {
243 B2ERROR("MCParticle Collection is not valid, cannot load into Graph");
244 return;
245 }
246
247 unsigned numParticles = MCParticles.getEntries();
248 unsigned particleOffset = size();
249 //Here we assume that the MCParticle collection is somehow ordered: All
250 //particles which are product of a decay come in the list after the mother,
251 //thus having a higher index. This is true for all lists generated by the
252 //MCParticleGraph and also for the standard lists produced by Evtgen and
253 //similiar generators.
254 for (unsigned i = 0; i < numParticles; ++i) {
255 GraphParticle& newParticle = addParticle();
256 const MCParticle& oldParticle = *MCParticles[i];
257 //Copy all values
258 newParticle = oldParticle;
259 //If this particle has a mother we just add the decay to this mother
260 const MCParticle* oldMother = oldParticle.getMother();
261 if (oldMother != nullptr) {
262 unsigned motherIndex = oldMother->getArrayIndex() + particleOffset;
263 if (motherIndex >= size())
264 B2FATAL("MCParticle collection \"" << name << "\" not sorted correctly: mother index larger than daughter. Cannot load into Graph");
265 newParticle.comesFrom((*this)[oldMother->getArrayIndex() + particleOffset]);
266 }
267 }
268}
Class to represent Particle data in graph.
void comesFrom(GraphParticle &mother)
Tells the graph that this particle is a decay product of mother.
Class to go over all the particles in the Graph an sort them in a sensible way.
void setVertexTime(MCParticleGraph::GraphParticle &m, const MCParticleGraph::GraphParticle &d)
Set the vertex and time information of the mother particle.
ParticleSorter(MemoryPool< MCParticleGraph::GraphParticle > &particles, TClonesArray *plist, bool setVertex, bool setTime)
ParticleSorter constructor.
MemoryPool< MCParticleGraph::GraphParticle > & m_particles
Reference to the list of particles which should be sorted.
void find_daughters(Vertex v, Graph &g, MCParticleGraph::GraphParticle &mother)
Find the daughters of the given vertex.
void finish_vertex(Vertex v, Graph &g)
Go through the daughters of the vertex.
bool m_setTime
True if the production time information should be saved.
std::queue< unsigned int > m_vqueue
The list of the vertices that will be visited.
vector< bool > m_seen
Vector of the particles that were already seen while sorting the graph.
void setStartIndex(int index)
Set the starting index for the particle graph.
int m_index
The latest index given to a particle.
bool m_setVertex
True if the vertex information should be saved.
TClonesArray * m_plist
The final array of sorted particles which is stored in the DataStore.
void sort(Graph &g)
Sort the particles and generate MCParticle list.
@ c_setDecayVertex
Set the decay vertex to the production vertex of the last daughter (ordered by production time)
@ c_setDecayTime
Set decay time to the largest production time of the daughters.
@ c_checkCyclic
Check for cyclic dependencies.
@ c_clearParticles
Clear the particle list before adding the graph.
size_t size() const
Return the number of particles in the graph.
std::pair< unsigned int, unsigned int > DecayLine
Type representing a decay in the graph.
void loadList(const std::string &name="")
Load the MCParticle list given by name into the Graph.
std::set< DecayLine > m_decays
internal set of decay lines
MemoryPool< GraphParticle > m_particles
internal list of particles
void generateList(const std::string &name="", int options=c_setNothing)
Generates the MCParticle list and stores it in the StoreArray with the given name.
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 getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:230
int m_firstDaughter
1-based index of first daughter particle in collection, 0 if no daughters
Definition: MCParticle.h:561
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
bool hasValidVertex() const
Indication whether vertex and time information is useful or just default.
Definition: MCParticle.h:153
float getProductionTime() const
Return production time in ns.
Definition: MCParticle.h:159
Class to provide a constant access time memory pool for one kind of objects.
Definition: MemoryPool.h:33
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
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
MCParticle * getMother() const
Returns a pointer to the mother particle.
Definition: MCParticle.h:600
GraphParticle & addParticle()
Add new particle to the graph.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
Simple struct to check boost graph for cyclic references.
void back_edge(Edge, Graph &)
This method is invoked on back edges in the graph.