ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCTruthManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MCTruthManager.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 //
30 //
31 //
32 // --------------------------------------------------------------
33 // GEANT 4 - MCTruthManager class
34 // --------------------------------------------------------------
35 //
36 // Author: Witold POKORSKI (Witold.Pokorski@cern.ch)
37 //
38 // --------------------------------------------------------------
39 //
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 
42 #include "MCTruthManager.hh"
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 
50 MCTruthManager::MCTruthManager() : fEvent(0), fConfig(0)
51 {}
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 
56 {}
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59 
61 {
62  if( instance == 0 )
63  {
64  instance = new MCTruthManager();
65  }
66  return instance;
67 }
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70 
72 {
73  // first delete the old event
74  delete fEvent;
75  // and now instaciate a new one
76  fEvent = new HepMC::GenEvent();
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
82  G4LorentzVector& prodpos,
83  G4LorentzVector& endpos,
84  G4int pdg_id, G4int partID, G4int motherID,
85  G4bool directParent)
86 {
87  // we create a new particle with barcode = partID
88  HepMC::GenParticle* particle = new HepMC::GenParticle(momentum, pdg_id);
89  particle->suggest_barcode(partID);
90  // we initialize the 'segmentations' map
91  // for the moment particle is not 'segmented'
92  fSegmentations[partID] = 1;
93 
94  // we create the GenVertex corresponding to the end point of the track
95  HepMC::GenVertex* endvertex = new HepMC::GenVertex(endpos);
96 
97  // barcode of the endvertex = - barcode of the track
98  endvertex->suggest_barcode(-partID);
99  endvertex->add_particle_in(particle);
100  fEvent->add_vertex(endvertex);
101 
102  if(motherID) // not a primary
103  {
104  // here we could try to improve speed by searching only through particles which
105  // belong to the given primary tree
106  HepMC::GenParticle* mother = fEvent->barcode_to_particle(motherID);
107  //
108  if(mother)
109  {
110  // we first check whether the mother's end vertex corresponds to the particle's
111  // production vertex
112  HepMC::GenVertex* motherendvtx = mother->end_vertex();
113  HepMC::FourVector mp0 = motherendvtx->position();
114  G4LorentzVector motherendpos(mp0.x(), mp0.y(), mp0.z(), mp0.t());
115 
116  if( motherendpos.x() == prodpos.x() &&
117  motherendpos.y() == prodpos.y() &&
118  motherendpos.z() == prodpos.z() ) // if yes, we attach the particle
119  {
120  motherendvtx->add_particle_out(particle);
121  }
122  else // if not, we check whether the mother is biological or adopted
123  {
124  if(!directParent) // adopted
125  {
126  G4bool found = false;
127 
128  // first check if any of the dummy particles
129  // has the end vertex at the right place
130  //
131  for(HepMC::GenVertex::particles_out_const_iterator
132  it=motherendvtx->particles_out_const_begin();
133  it!=motherendvtx->particles_out_const_end();it++)
134  {
135  if((*it)->pdg_id()==-999999)
136  {
137  HepMC::FourVector dp0 = (*it)->end_vertex()->position();
138  G4LorentzVector dummypos(dp0.x(), dp0.y(), dp0.z(), dp0.t());;
139 
140  if( dummypos.x() == prodpos.x() &&
141  dummypos.y() == prodpos.y() &&
142  dummypos.z() == prodpos.z() )
143  {
144  (*it)->end_vertex()->add_particle_out(particle);
145  found = true;
146  break;
147  }
148  }
149  }
150 
151  // and if not, create a dummy particle connecting
152  // to the end vertex of the mother
153  //
154  if(!found)
155  {
156  HepMC::GenVertex* childvtx = new HepMC::GenVertex(prodpos);
157  childvtx->add_particle_out(particle);
158 
159  // the dummy vertex gets the barcode -500000
160  // minus the daughter particle barcode
161  //
162  childvtx->suggest_barcode(-500000-partID);
163  fEvent->add_vertex(childvtx);
164 
165  HepMC::GenParticle* dummypart =
166  new HepMC::GenParticle(G4LorentzVector(),-999999);
167 
168  // the dummy particle gets the barcode 500000
169  // plus the daughter particle barcode
170  //
171  dummypart->suggest_barcode(500000+partID);
172  childvtx->add_particle_in(dummypart);
173  motherendvtx->add_particle_out(dummypart);
174  }
175  }
176  else // biological
177  {
178  // in case mother was already 'split' we need to look for
179  // the right 'segment' to add the new daugther.
180  // We use Time coordinate to locate the place for the new vertex
181 
182  G4int number_of_segments = fSegmentations[motherID];
183  G4int segment = 0;
184 
185  // we loop through the segments
186  //
187  while ( !((mother->end_vertex()->position().t()>prodpos.t()) &&
188  (mother->production_vertex()->position().t()<prodpos.t())) )
189  {
190  segment++;
191  if (segment == number_of_segments)
192  G4cerr << "Problem!!!! Time coordinates incompatible!" << G4endl;
193 
194  mother = fEvent->barcode_to_particle(segment*10000000 + motherID);
195  }
196 
197  // now, we 'split' the appropriate 'segment' of the mother particle
198  // into two particles and create a new vertex
199  //
200  HepMC::GenVertex* childvtx = new HepMC::GenVertex(prodpos);
201  childvtx->add_particle_out(particle);
202  fEvent->add_vertex(childvtx);
203 
204  // we first detach the mother from its original vertex
205  //
206  HepMC::GenVertex* orig_mother_end_vtx = mother->end_vertex();
207  orig_mother_end_vtx->remove_particle(mother);
208 
209  // and attach it to the new vertex
210  //
211  childvtx->add_particle_in(mother);
212 
213  // now we create a new particle representing the mother after
214  // interaction the barcode of the new particle is 10000000 + the
215  // original barcode
216  //
217  HepMC::GenParticle* mothertwo = new HepMC::GenParticle(*mother);
218  mothertwo->suggest_barcode(fSegmentations[motherID]*10000000
219  + mother->barcode());
220 
221  // we also reset the barcodes of the vertices
222  //
223  orig_mother_end_vtx->suggest_barcode(-fSegmentations[motherID]
224  *10000000 - mother->barcode());
225  childvtx->suggest_barcode(-mother->barcode());
226 
227  // we attach it to the new vertex where interaction took place
228  //
229  childvtx->add_particle_out(mothertwo);
230 
231  // and we attach it to the original endvertex
232  //
233  orig_mother_end_vtx->add_particle_in(mothertwo);
234 
235  // and finally ... the increase the 'segmentation counter'
236  //
237  fSegmentations[motherID] = fSegmentations[motherID]+1;
238  }
239  }
240  }
241  else
242  // mother GenParticle is not there for some reason...
243  // if this happens, we need to revise the philosophy...
244  // a solution would be to create HepMC particles
245  // at the begining of each track
246  {
247  G4cerr << "barcode " << motherID << " mother not there! "<< G4endl;
248  }
249  }
250  else // primary
251  {
252  HepMC::GenVertex* primaryvtx = new HepMC::GenVertex(prodpos);
253  primaryvtx->add_particle_out(particle);
254  fEvent->add_vertex(primaryvtx);
255 
256  // add id to the list of primaries
257  //
258  fPrimarybarcodes.push_back(partID);
259  }
260 }
261 
262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263 
265 {
266  fEvent->print();
267 
268  // looping over primaries and print the decay tree for each of them
269  //
270  for(std::vector<int>::const_iterator primarybar=fPrimarybarcodes.begin();
271  primarybar!=fPrimarybarcodes.end();primarybar++)
272  {
273  PrintTree(fEvent->barcode_to_particle(*primarybar), " | ");
274  }
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278 
280 {
281  G4cout << offset << "--- barcode: " << particle->barcode() << " pdg: "
282  << particle->pdg_id() << " energy: " << particle->momentum().e()
283  << " production vertex: "
284  << particle->production_vertex()->position().x() << ", "
285  << particle->production_vertex()->position().y() << ", "
286  << particle->production_vertex()->position().z() << ", "
287  << particle->production_vertex()->position().t()
288  << G4endl;
289 
290  for(HepMC::GenVertex::particles_out_const_iterator
291  it=particle->end_vertex()->particles_out_const_begin();
292  it!=particle->end_vertex()->particles_out_const_end();
293  it++)
294  {
295  G4String deltaoffset = "";
296 
297  G4int curr = std::fmod(double((*it)->barcode()),10000000.);
298  G4int part = std::fmod(double(particle->barcode()),10000000.);
299  if( curr != part )
300  {
301  deltaoffset = " | ";
302  }
303 
304  PrintTree((*it), offset + deltaoffset);
305  }
306 }
307 
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....