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