ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHHepMCParticleSelectorDecayProductChain.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHHepMCParticleSelectorDecayProductChain.cc
2 
3 #include "PHHepMCGenEvent.h"
4 #include "PHHepMCGenEventMap.h"
5 
7 #include <fun4all/SubsysReco.h> // for SubsysReco
8 
9 #include <phool/getClass.h>
10 #include <phool/phool.h> // for PHWHERE
11 
12 #include <HepMC/GenEvent.h> // for GenEvent::particle_const_ite...
13 #include <HepMC/GenParticle.h> // for GenParticle
14 #include <HepMC/GenVertex.h> // for GenVertex, GenVertex::partic...
15 #include <HepMC/IteratorRange.h> // for ancestors, children, descend...
16 #include <HepMC/SimpleVector.h> // for FourVector
17 
18 #include <cstdlib> // for abs
19 #include <iostream> // for operator<<, basic_ostream::o...
20 
21 class PHCompositeNode;
22 
23 using namespace std;
24 
26  : SubsysReco(name)
27  , _embedding_id(0)
28 {
29  _theParticle = 11;
30  return;
31 }
32 
34 {
36 }
37 
38 HepMC::GenParticle* PHHepMCParticleSelectorDecayProductChain::GetParent(HepMC::GenParticle* p, HepMC::GenEvent* /*event*/)
39 {
40  HepMC::GenParticle* parent = nullptr;
41  if (!p->production_vertex()) return parent;
42 
43  for (HepMC::GenVertex::particle_iterator mother = p->production_vertex()->particles_begin(HepMC::ancestors);
44  mother != p->production_vertex()->particles_end(HepMC::ancestors);
45  ++mother)
46  {
47  for (unsigned int i = 0; i < _theAncestors.size(); i++)
48  {
49  if (abs((*mother)->pdg_id()) == _theAncestors[i])
50  {
51  parent = *mother;
52  break;
53  }
54  }
55  if (parent != nullptr) break;
56  }
57 
58  return parent;
59 }
60 
62 {
63  if (_theParticle == 0 && _theDaughters.size() == 0)
64  {
65  cout << PHWHERE << "Doing nothing." << endl;
67  }
68 
69  PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
70  if (!geneventmap)
71  {
72  cerr << "ERROR: PHHepMCGenEventMap node not found!" << endl;
74  }
75  PHHepMCGenEvent* inEvent = geneventmap->get(_embedding_id);
76  if (!inEvent)
77  {
78  cerr << "PHHepMCParticleSelectorDecayProductChain::process_event - WARNING: PHHepMCGenEvent with embedding ID " << _embedding_id << " is not found! Move on." << endl;
80  }
81 
82  HepMC::GenEvent* event = inEvent->getEvent();
83  int npart = event->particles_size();
84  int nvert = event->vertices_size();
85  if (Verbosity() > 0) cout << "=========== Event " << event->event_number() << " contains " << npart << " particles and " << nvert << " vertices." << endl;
86 
87  // list of vertices to keep
88  vector<HepMC::GenVertex> vkeep;
89 
90  if (_theParticle != 0) // keep _theParticle and all its daughters (if any)
91  {
92  for (HepMC::GenEvent::particle_const_iterator p = event->particles_begin(); p != event->particles_end(); ++p)
93  {
94  // find _thePartcle
95  if (abs((*p)->pdg_id()) == _theParticle)
96  {
97  // do we need to check for ancestors?
98  if (_theAncestors.size() > 0)
99  {
100  HepMC::GenParticle* parent = GetParent(*p, event);
101  if (parent)
102  {
103  vkeep.push_back(*(*p)->production_vertex());
104  }
105  }
106  else
107  {
108  vkeep.push_back(*(*p)->production_vertex());
109  }
110 
111  // do we need to keep the daughters?
112  if (_theDaughters.size() > 0)
113  {
114  if ((*p)->end_vertex())
115  {
116  for (HepMC::GenVertex::particle_iterator des = (*p)->end_vertex()->particles_begin(HepMC::descendants);
117  des != (*p)->end_vertex()->particles_end(HepMC::descendants);
118  ++des)
119  {
120  for (unsigned int i = 0; i < _theDaughters.size(); i++)
121  {
122  if (abs((*des)->pdg_id()) == _theDaughters[i])
123  {
124  vkeep.push_back(*(*p)->end_vertex());
125  break;
126  }
127  }
128  }
129  }
130  } // there are daughters
131 
132  } // this is _theParticle
133  } // end loop over particles
134  }
135  else // save only particles in _theDaughters list no matter where they came from
136  {
137  for (HepMC::GenEvent::particle_const_iterator p = event->particles_begin(); p != event->particles_end(); ++p)
138  {
139  for (unsigned int ip = 0; ip < _theDaughters.size(); ip++)
140  {
141  if (abs((*p)->pdg_id()) == _theDaughters[ip])
142  {
143  vkeep.push_back(*(*p)->production_vertex());
144  }
145  }
146  }
147  }
148 
149  // loop over vertices and keep only selected ones.
150  for (HepMC::GenEvent::vertex_const_iterator v = event->vertices_begin(); v != event->vertices_end(); ++v)
151  {
152  bool goodvertex = false;
153  for (unsigned int i = 0; i < vkeep.size(); i++)
154  {
155  HepMC::GenVertex tmp1 = (*(*v));
156  HepMC::GenVertex tmp2 = vkeep[i];
157  if (tmp1 == tmp2)
158  {
159  goodvertex = true;
160  }
161  }
162  if (!goodvertex)
163  {
164  bool tmp = event->remove_vertex((*v));
165  if (Verbosity() > 10 && tmp)
166  {
167  cout << PHWHERE << " Erasing empty vertex." << endl;
168  }
169  }
170  }
171 
172  // clean up the vertices
173  for (HepMC::GenEvent::vertex_const_iterator v = event->vertices_begin(); v != event->vertices_end(); ++v)
174  {
175  std::vector<HepMC::GenParticle*> removep;
176 
177  for (HepMC::GenVertex::particle_iterator itpart = (*v)->particles_begin(HepMC::children);
178  itpart != (*v)->particles_end(HepMC::children);
179  ++itpart)
180  {
181  bool keepparticle = false;
182  if (abs((*itpart)->pdg_id()) == _theParticle)
183  {
184  keepparticle = true;
185  }
186  for (unsigned int j = 0; j < _theDaughters.size(); j++)
187  {
188  if (abs((*itpart)->pdg_id()) == _theDaughters[j] && (*itpart)->status() == 1)
189  {
190  keepparticle = true;
191  }
192  }
193  if (!keepparticle)
194  {
195  removep.push_back((*itpart));
196  }
197  } // end loop over particles in this vertex
198 
199  for (HepMC::GenVertex::particle_iterator itpart = (*v)->particles_begin(HepMC::parents);
200  itpart != (*v)->particles_end(HepMC::parents);
201  ++itpart)
202  {
203  bool keepparticle = false;
204  if (abs((*itpart)->pdg_id()) == _theParticle)
205  {
206  keepparticle = true;
207  }
208  for (unsigned int j = 0; j < _theDaughters.size(); j++)
209  {
210  if (abs((*itpart)->pdg_id()) == _theDaughters[j] && (*itpart)->status() == 1)
211  {
212  keepparticle = true;
213  }
214  }
215  if (!keepparticle)
216  {
217  removep.push_back((*itpart));
218  }
219  } // end loop over particles in this vertex
220 
221  for (unsigned int k = 0; k < removep.size(); k++)
222  {
223  HepMC::GenParticle* tmp = (*v)->remove_particle(removep[k]);
224  if (tmp->end_vertex())
225  {
226  delete tmp->end_vertex();
227  }
228  }
229  }
230 
231  int partcount = 0;
232  if (Verbosity() > 0)
233  {
234  cout << "FINAL Event " << event->event_number() << " contains " << event->particles_size() << " particles and " << event->vertices_size() << " vertices." << endl;
235  cout << "FINAL LIST OF PARTICLES:" << endl;
236  }
237  for (HepMC::GenEvent::particle_const_iterator p = event->particles_begin(); p != event->particles_end(); ++p)
238  {
239  int pid = (*p)->pdg_id();
240  int status = (*p)->status();
241  double pz = ((*p)->momentum()).pz();
242  double pt = ((*p)->momentum()).perp();
243  double eta = ((*p)->momentum()).eta();
244  double mass = ((*p)->momentum()).m();
245  if (Verbosity() > 0)
246  {
247  cout << pid << " " << mass << " " << status << " " << pt << " " << pz << " " << eta << " " << (*p)->production_vertex() << " " << (*p)->end_vertex() << endl;
248  }
249  partcount++;
250  }
251 
252  // if there is nothing to write out the code crashes
253  if (partcount == 0)
254  {
255  if (Verbosity() > 0)
256  {
257  cout << "EVENT ABORTED: No particles to write out." << endl;
258  }
260  }
261  else
262  {
264  }
265 }
266 
268 {
269  _theParticle = pid;
270  return;
271 }
272 
274 {
275  _theAncestors.push_back(pid);
276  return;
277 }
278 
280 {
281  _theDaughters.push_back(pid);
282  return;
283 }