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