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