source: trunk/examples/extended/eventgenerator/HepMC/MCTruth/src/MCTruthManager.cc@ 1036

Last change on this file since 1036 was 807, checked in by garnier, 17 years ago

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