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