| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file eventgenerator/HepMC/MCTruth/src/MCT
27 /// \brief Implementation of the MCTruthManage
28 //
29 //
30 //
31 //
32 // -------------------------------------------
33 // GEANT 4 - MCTruthManager class
34 // -------------------------------------------
35 //
36 // Author: Witold POKORSKI (Witold.Pokorski@ce
37 //
38 // -------------------------------------------
39 //
40 //....oooOO0OOooo........oooOO0OOooo........oo
41
42 #include "MCTruthManager.hh"
43
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 static MCTruthManager* instance = 0;
47
48 //....oooOO0OOooo........oooOO0OOooo........oo
49
50 MCTruthManager::MCTruthManager() : fEvent(0),
51
52 //....oooOO0OOooo........oooOO0OOooo........oo
53
54 MCTruthManager::~MCTruthManager() {}
55
56 //....oooOO0OOooo........oooOO0OOooo........oo
57
58 MCTruthManager* MCTruthManager::GetInstance()
59 {
60 if (instance == 0) {
61 instance = new MCTruthManager();
62 }
63 return instance;
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oo
67
68 void MCTruthManager::NewEvent()
69 {
70 // first delete the old event
71 delete fEvent;
72 // and now instaciate a new one
73 fEvent = new HepMC::GenEvent();
74 }
75
76 //....oooOO0OOooo........oooOO0OOooo........oo
77
78 void MCTruthManager::AddParticle(G4LorentzVect
79 G4LorentzVect
80 G4int motherI
81 {
82 // we create a new particle with barcode = p
83 HepMC::GenParticle* particle = new HepMC::Ge
84 particle->suggest_barcode(partID);
85 // we initialize the 'segmentations' map
86 // for the moment particle is not 'segmented
87 fSegmentations[partID] = 1;
88
89 // we create the GenVertex corresponding to
90 HepMC::GenVertex* endvertex = new HepMC::Gen
91
92 // barcode of the endvertex = - barcode of t
93 endvertex->suggest_barcode(-partID);
94 endvertex->add_particle_in(particle);
95 fEvent->add_vertex(endvertex);
96
97 if (motherID) // not a primary
98 {
99 // here we could try to improve speed by s
100 // belong to the given primary tree
101 HepMC::GenParticle* mother = fEvent->barco
102 //
103 if (mother) {
104 // we first check whether the mother's e
105 // production vertex
106 HepMC::GenVertex* motherendvtx = mother-
107 HepMC::FourVector mp0 = motherendvtx->po
108 G4LorentzVector motherendpos(mp0.x(), mp
109
110 if (motherendpos.x() == prodpos.x() && m
111 && motherendpos.z() == prodpos.z())
112 {
113 motherendvtx->add_particle_out(particl
114 }
115 else // if not, we check whether the mo
116 {
117 if (!directParent) // adopted
118 {
119 G4bool found = false;
120
121 // first check if any of the dummy p
122 // has the end vertex at the right p
123 //
124 for (HepMC::GenVertex::particles_out
125 motherendvtx->particles_out_c
126 it != motherendvtx->particles_o
127 {
128 if ((*it)->pdg_id() == -999999) {
129 HepMC::FourVector dp0 = (*it)->e
130 G4LorentzVector dummypos(dp0.x()
131 ;
132
133 if (dummypos.x() == prodpos.x()
134 && dummypos.z() == prodpos.z
135 {
136 (*it)->end_vertex()->add_parti
137 found = true;
138 break;
139 }
140 }
141 }
142
143 // and if not, create a dummy partic
144 // to the end vertex of the mother
145 //
146 if (!found) {
147 HepMC::GenVertex* childvtx = new H
148 childvtx->add_particle_out(particl
149
150 // the dummy vertex gets the barco
151 // minus the daughter particle bar
152 //
153 childvtx->suggest_barcode(-500000
154 fEvent->add_vertex(childvtx);
155
156 HepMC::GenParticle* dummypart = ne
157
158 // the dummy particle gets the bar
159 // plus the daughter particle barc
160 //
161 dummypart->suggest_barcode(500000
162 childvtx->add_particle_in(dummypar
163 motherendvtx->add_particle_out(dum
164 }
165 }
166 else // biological
167 {
168 // in case mother was already 'split
169 // the right 'segment' to add the ne
170 // We use Time coordinate to locate
171
172 G4int number_of_segments = fSegmenta
173 G4int segment = 0;
174
175 // we loop through the segments
176 //
177 while (!((mother->end_vertex()->posi
178 && (mother->production_vert
179 {
180 segment++;
181 if (segment == number_of_segments)
182 G4cerr << "Problem!!!! Time coor
183
184 mother = fEvent->barcode_to_partic
185 }
186
187 // now, we 'split' the appropriate '
188 // into two particles and create a n
189 //
190 HepMC::GenVertex* childvtx = new Hep
191 childvtx->add_particle_out(particle)
192 fEvent->add_vertex(childvtx);
193
194 // we first detach the mother from i
195 //
196 HepMC::GenVertex* orig_mother_end_vt
197 orig_mother_end_vtx->remove_particle
198
199 // and attach it to the new vertex
200 //
201 childvtx->add_particle_in(mother);
202
203 // now we create a new particle repr
204 // interaction the barcode of the ne
205 // original barcode
206 //
207 HepMC::GenParticle* mothertwo = new
208 mothertwo->suggest_barcode(fSegmenta
209
210 // we also reset the barcodes of the
211 //
212 orig_mother_end_vtx->suggest_barcode
213
214 childvtx->suggest_barcode(-mother->b
215
216 // we attach it to the new vertex wh
217 //
218 childvtx->add_particle_out(mothertwo
219
220 // and we attach it to the original
221 //
222 orig_mother_end_vtx->add_particle_in
223
224 // and finally ... the increase the
225 //
226 fSegmentations[motherID] = fSegmenta
227 }
228 }
229 }
230 else
231 // mother GenParticle is not there for som
232 // if this happens, we need to revise the
233 // a solution would be to create HepMC par
234 // at the begining of each track
235 {
236 G4cerr << "barcode " << motherID << " mo
237 }
238 }
239 else // primary
240 {
241 HepMC::GenVertex* primaryvtx = new HepMC::
242 primaryvtx->add_particle_out(particle);
243 fEvent->add_vertex(primaryvtx);
244
245 // add id to the list of primaries
246 //
247 fPrimarybarcodes.push_back(partID);
248 }
249 }
250
251 //....oooOO0OOooo........oooOO0OOooo........oo
252
253 void MCTruthManager::PrintEvent()
254 {
255 fEvent->print();
256
257 // looping over primaries and print the deca
258 //
259 for (std::vector<int>::const_iterator primar
260 primarybar != fPrimarybarcodes.end(); p
261 {
262 PrintTree(fEvent->barcode_to_particle(*pri
263 }
264 }
265
266 //....oooOO0OOooo........oooOO0OOooo........oo
267
268 void MCTruthManager::PrintTree(HepMC::GenParti
269 {
270 G4cout << offset << "--- barcode: " << part
271 << " energy: " << particle->momentum(
272 << " production vertex: " << particle
273 << particle->production_vertex()->pos
274 << particle->production_vertex()->pos
275 << particle->production_vertex()->pos
276
277 for (HepMC::GenVertex::particles_out_const_i
278 particle->end_vertex()->particles_out
279 it != particle->end_vertex()->particles
280 {
281 G4String deltaoffset = "";
282
283 G4int curr = std::fmod(double((*it)->barco
284 G4int part = std::fmod(double(particle->ba
285 if (curr != part) {
286 deltaoffset = " | ";
287 }
288
289 PrintTree((*it), offset + deltaoffset);
290 }
291 }
292
293 //....oooOO0OOooo........oooOO0OOooo........oo
294