Geant4  10.02.p01
MCTruthManager.cc
Go to the documentation of this file.
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 // $Id: MCTruthManager.cc 73446 2013-08-27 11:32:59Z gcosmo $
31 //
32 //
33 // --------------------------------------------------------------
34 // GEANT 4 - MCTruthManager class
35 // --------------------------------------------------------------
36 //
37 // Author: Witold POKORSKI (Witold.Pokorski@cern.ch)
38 //
39 // --------------------------------------------------------------
40 
41 #include "MCTruthManager.hh"
42 
44 
46 {}
47 
49 {}
50 
52 {
53  if( instance == 0 )
54  {
55  instance = new MCTruthManager();
56  }
57  return instance;
58 }
59 
61 {
62  // first delete the old event
63  delete event;
64  // and now instaciate a new one
65  event = new HepMC::GenEvent();
66 }
67 
69  G4LorentzVector& prodpos,
70  G4LorentzVector& endpos,
71  G4int pdg_id, G4int partID, G4int motherID,
72  G4bool directParent)
73 {
74  // we create a new particle with barcode = partID
75  HepMC::GenParticle* particle = new HepMC::GenParticle(momentum, pdg_id);
76  particle->suggest_barcode(partID);
77  // we initialize the 'segmentations' map
78  // for the moment particle is not 'segmented'
79  segmentations[partID] = 1;
80 
81  // we create the GenVertex corresponding to the end point of the track
82  HepMC::GenVertex* endvertex = new HepMC::GenVertex(endpos);
83 
84  // barcode of the endvertex = - barcode of the track
85  endvertex->suggest_barcode(-partID);
86  endvertex->add_particle_in(particle);
87  event->add_vertex(endvertex);
88 
89  if(motherID) // not a primary
90  {
91  // here we could try to improve speed by searching only through particles which
92  // belong to the given primary tree
93  HepMC::GenParticle* mother = event->barcode_to_particle(motherID);
94  //
95  if(mother)
96  {
97  // we first check whether the mother's end vertex corresponds to the particle's
98  // production vertex
99  HepMC::GenVertex* motherendvtx = mother->end_vertex();
100  HepMC::FourVector mp0 = motherendvtx->position();
101  G4LorentzVector motherendpos(mp0.x(), mp0.y(), mp0.z(), mp0.t());
102 
103  if( motherendpos.x() == prodpos.x() &&
104  motherendpos.y() == prodpos.y() &&
105  motherendpos.z() == prodpos.z() ) // if yes, we attach the particle
106  {
107  motherendvtx->add_particle_out(particle);
108  }
109  else // if not, we check whether the mother is biological or adopted
110  {
111  if(!directParent) // adopted
112  {
113  G4bool found = false;
114 
115  // first check if any of the dummy particles
116  // has the end vertex at the right place
117  //
118  for(HepMC::GenVertex::particles_out_const_iterator
119  it=motherendvtx->particles_out_const_begin();
120  it!=motherendvtx->particles_out_const_end();it++)
121  {
122  if((*it)->pdg_id()==-999999)
123  {
124  HepMC::FourVector dp0 = (*it)->end_vertex()->position();
125  G4LorentzVector dummypos(dp0.x(), dp0.y(), dp0.z(), dp0.t());;
126 
127  if( dummypos.x() == prodpos.x() &&
128  dummypos.y() == prodpos.y() &&
129  dummypos.z() == prodpos.z() )
130  {
131  (*it)->end_vertex()->add_particle_out(particle);
132  found = true;
133  break;
134  }
135  }
136  }
137 
138  // and if not, create a dummy particle connecting
139  // to the end vertex of the mother
140  //
141  if(!found)
142  {
143  HepMC::GenVertex* childvtx = new HepMC::GenVertex(prodpos);
144  childvtx->add_particle_out(particle);
145 
146  // the dummy vertex gets the barcode -500000
147  // minus the daughter particle barcode
148  //
149  childvtx->suggest_barcode(-500000-partID);
150  event->add_vertex(childvtx);
151 
152  HepMC::GenParticle* dummypart =
153  new HepMC::GenParticle(G4LorentzVector(),-999999);
154 
155  // the dummy particle gets the barcode 500000
156  // plus the daughter particle barcode
157  //
158  dummypart->suggest_barcode(500000+partID);
159  childvtx->add_particle_in(dummypart);
160  motherendvtx->add_particle_out(dummypart);
161  }
162  }
163  else // biological
164  {
165  // in case mother was already 'split' we need to look for
166  // the right 'segment' to add the new daugther.
167  // We use Time coordinate to locate the place for the new vertex
168 
169  G4int number_of_segments = segmentations[motherID];
170  G4int segment = 0;
171 
172  // we loop through the segments
173  //
174  while ( !((mother->end_vertex()->position().t()>prodpos.t()) &&
175  (mother->production_vertex()->position().t()<prodpos.t())) )
176  {
177  segment++;
178  if (segment == number_of_segments)
179  G4cerr << "Problem!!!! Time coordinates incompatible!" << G4endl;
180 
181  mother = event->barcode_to_particle(segment*10000000 + motherID);
182  }
183 
184  // now, we 'split' the appropriate 'segment' of the mother particle
185  // into two particles and create a new vertex
186  //
187  HepMC::GenVertex* childvtx = new HepMC::GenVertex(prodpos);
188  childvtx->add_particle_out(particle);
189  event->add_vertex(childvtx);
190 
191  // we first detach the mother from its original vertex
192  //
193  HepMC::GenVertex* orig_mother_end_vtx = mother->end_vertex();
194  orig_mother_end_vtx->remove_particle(mother);
195 
196  // and attach it to the new vertex
197  //
198  childvtx->add_particle_in(mother);
199 
200  // now we create a new particle representing the mother after
201  // interaction the barcode of the new particle is 10000000 + the
202  // original barcode
203  //
204  HepMC::GenParticle* mothertwo = new HepMC::GenParticle(*mother);
205  mothertwo->suggest_barcode(segmentations[motherID]*10000000
206  + mother->barcode());
207 
208  // we also reset the barcodes of the vertices
209  //
210  orig_mother_end_vtx->suggest_barcode(-segmentations[motherID]
211  *10000000 - mother->barcode());
212  childvtx->suggest_barcode(-mother->barcode());
213 
214  // we attach it to the new vertex where interaction took place
215  //
216  childvtx->add_particle_out(mothertwo);
217 
218  // and we attach it to the original endvertex
219  //
220  orig_mother_end_vtx->add_particle_in(mothertwo);
221 
222  // and finally ... the increase the 'segmentation counter'
223  //
224  segmentations[motherID] = segmentations[motherID]+1;
225  }
226  }
227  }
228  else
229  // mother GenParticle is not there for some reason...
230  // if this happens, we need to revise the philosophy...
231  // a solution would be to create HepMC particles
232  // at the begining of each track
233  {
234  G4cerr << "barcode " << motherID << " mother not there! "<< G4endl;
235  }
236  }
237  else // primary
238  {
239  HepMC::GenVertex* primaryvtx = new HepMC::GenVertex(prodpos);
240  primaryvtx->add_particle_out(particle);
241  event->add_vertex(primaryvtx);
242 
243  // add id to the list of primaries
244  //
245  primarybarcodes.push_back(partID);
246  }
247 }
248 
250 {
251  event->print();
252 
253  // looping over primaries and print the decay tree for each of them
254  //
255  for(std::vector<int>::const_iterator primarybar=primarybarcodes.begin();
256  primarybar!=primarybarcodes.end();primarybar++)
257  {
258  printTree(event->barcode_to_particle(*primarybar), " | ");
259  }
260 }
261 
262 void MCTruthManager::printTree(HepMC::GenParticle* particle, G4String offset)
263 {
264  G4cout << offset << "--- barcode: " << particle->barcode() << " pdg: "
265  << particle->pdg_id() << " energy: " << particle->momentum().e()
266  << " production vertex: "
267  << particle->production_vertex()->position().x() << ", "
268  << particle->production_vertex()->position().y() << ", "
269  << particle->production_vertex()->position().z() << ", "
270  << particle->production_vertex()->position().t()
271  << G4endl;
272 
273  for(HepMC::GenVertex::particles_out_const_iterator
274  it=particle->end_vertex()->particles_out_const_begin();
275  it!=particle->end_vertex()->particles_out_const_end();
276  it++)
277  {
278  G4String deltaoffset = "";
279 
280  G4int curr = std::fmod(double((*it)->barcode()),10000000.);
281  G4int part = std::fmod(double(particle->barcode()),10000000.);
282  if( curr != part )
283  {
284  deltaoffset = " | ";
285  }
286 
287  printTree((*it), offset + deltaoffset);
288  }
289 }
static MCTruthManager * GetInstance()
void AddParticle(G4LorentzVector &, G4LorentzVector &, G4LorentzVector &, G4int, G4int, G4int, G4bool)
virtual ~MCTruthManager()
Definition of the MCTruthManager class.
int G4int
Definition: G4Types.hh:78
std::map< G4int, G4int > segmentations
G4GLOB_DLL std::ostream G4cout
HepMC::GenEvent * event
bool G4bool
Definition: G4Types.hh:79
std::vector< G4int > primarybarcodes
void printTree(HepMC::GenParticle *, G4String)
#define G4endl
Definition: G4ios.hh:61
static MCTruthManager * instance
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector