MCTruth/src/MCTruthManager.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file eventgenerator/HepMC/MCTruth/src/MCTruthManager.cc
0027 /// \brief Implementation of the MCTruthManager class
0028 //
0029 //
0030 //
0031 //
0032 // --------------------------------------------------------------
0033 //      GEANT 4 - MCTruthManager class
0034 // --------------------------------------------------------------
0035 //
0036 // Author: Witold POKORSKI (Witold.Pokorski@cern.ch)
0037 //
0038 // --------------------------------------------------------------
0039 //
0040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0041
0042 #include "MCTruthManager.hh"
0043
0044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0045
0046 static MCTruthManager* instance = 0;
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0049
0050 MCTruthManager::MCTruthManager() : fEvent(0), fConfig(0) {}
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0053
0054 MCTruthManager::~MCTruthManager() {}
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0057
0058 MCTruthManager* MCTruthManager::GetInstance()
0059 {
0060   if (instance == 0) {
0061     instance = new MCTruthManager();
0062   }
0063   return instance;
0064 }
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0067
0068 void MCTruthManager::NewEvent()
0069 {
0070   // first delete the old event
0071   delete fEvent;
0072   // and now instaciate a new one
0073   fEvent = new HepMC::GenEvent();
0074 }
0075
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0077
0078 void MCTruthManager::AddParticle(G4LorentzVector& momentum, G4LorentzVector& prodpos,
0079                                  G4LorentzVector& endpos, G4int pdg_id, G4int partID,
0080                                  G4int motherID, G4bool directParent)
0081 {
0082   // we create a new particle with barcode = partID
0083   HepMC::GenParticle* particle = new HepMC::GenParticle(momentum, pdg_id);
0084   particle->suggest_barcode(partID);
0085   // we initialize the 'segmentations' map
0086   // for the moment particle is not 'segmented'
0087   fSegmentations[partID] = 1;
0088
0089   // we create the GenVertex corresponding to the end point of the track
0090   HepMC::GenVertex* endvertex = new HepMC::GenVertex(endpos);
0091
0092   // barcode of the endvertex = - barcode of the track
0093   endvertex->suggest_barcode(-partID);
0094   endvertex->add_particle_in(particle);
0095   fEvent->add_vertex(endvertex);
0096
0097   if (motherID)  // not a primary
0098   {
0099     // here we could try to improve speed by searching only through particles which
0100     // belong to the given primary tree
0101     HepMC::GenParticle* mother = fEvent->barcode_to_particle(motherID);
0102     //
0103     if (mother) {
0104       // we first check whether the mother's end vertex corresponds to the particle's
0105       // production vertex
0106       HepMC::GenVertex* motherendvtx = mother->end_vertex();
0107       HepMC::FourVector mp0 = motherendvtx->position();
0108       G4LorentzVector motherendpos(mp0.x(), mp0.y(), mp0.z(), mp0.t());
0109
0110       if (motherendpos.x() == prodpos.x() && motherendpos.y() == prodpos.y()
0111           && motherendpos.z() == prodpos.z())  // if yes, we attach the particle
0112       {
0113         motherendvtx->add_particle_out(particle);
0114       }
0115       else  // if not, we check whether the mother is biological or adopted
0116       {
0117         if (!directParent)  // adopted
0118         {
0119           G4bool found = false;
0120
0121           // first check if any of the dummy particles
0122           // has the end vertex at the right place
0123           //
0124           for (HepMC::GenVertex::particles_out_const_iterator it =
0125                  motherendvtx->particles_out_const_begin();
0126                it != motherendvtx->particles_out_const_end(); it++)
0127           {
0128             if ((*it)->pdg_id() == -999999) {
0129               HepMC::FourVector dp0 = (*it)->end_vertex()->position();
0130               G4LorentzVector dummypos(dp0.x(), dp0.y(), dp0.z(), dp0.t());
0131               ;
0132
0133               if (dummypos.x() == prodpos.x() && dummypos.y() == prodpos.y()
0134                   && dummypos.z() == prodpos.z())
0135               {
0136                 (*it)->end_vertex()->add_particle_out(particle);
0137                 found = true;
0138                 break;
0139               }
0140             }
0141           }
0142
0143           // and if not, create a dummy particle connecting
0144           // to the end vertex of the mother
0145           //
0146           if (!found) {
0147             HepMC::GenVertex* childvtx = new HepMC::GenVertex(prodpos);
0148             childvtx->add_particle_out(particle);
0149
0150             // the dummy vertex gets the barcode -500000
0151             // minus the daughter particle barcode
0152             //
0153             childvtx->suggest_barcode(-500000 - partID);
0154             fEvent->add_vertex(childvtx);
0155
0156             HepMC::GenParticle* dummypart = new HepMC::GenParticle(G4LorentzVector(), -999999);
0157
0158             // the dummy particle gets the barcode 500000
0159             // plus the daughter particle barcode
0160             //
0161             dummypart->suggest_barcode(500000 + partID);
0162             childvtx->add_particle_in(dummypart);
0163             motherendvtx->add_particle_out(dummypart);
0164           }
0165         }
0166         else  // biological
0167         {
0168           // in case mother was already 'split' we need to look for
0169           // the right 'segment' to add the new daugther.
0170           // We use Time coordinate to locate the place for the new vertex
0171
0172           G4int number_of_segments = fSegmentations[motherID];
0173           G4int segment = 0;
0174
0175           // we loop through the segments
0176           //
0177           while (!((mother->end_vertex()->position().t() > prodpos.t())
0178                    && (mother->production_vertex()->position().t() < prodpos.t())))
0179           {
0180             segment++;
0181             if (segment == number_of_segments)
0182               G4cerr << "Problem!!!! Time coordinates incompatible!" << G4endl;
0183
0184             mother = fEvent->barcode_to_particle(segment * 10000000 + motherID);
0185           }
0186
0187           // now, we 'split' the appropriate 'segment' of the mother particle
0188           // into two particles and create a new vertex
0189           //
0190           HepMC::GenVertex* childvtx = new HepMC::GenVertex(prodpos);
0191           childvtx->add_particle_out(particle);
0192           fEvent->add_vertex(childvtx);
0193
0194           // we first detach the mother from its original vertex
0195           //
0196           HepMC::GenVertex* orig_mother_end_vtx = mother->end_vertex();
0197           orig_mother_end_vtx->remove_particle(mother);
0198
0199           // and attach it to the new vertex
0200           //
0201           childvtx->add_particle_in(mother);
0202
0203           // now we create a new particle representing the mother after
0204           // interaction the barcode of the new particle is 10000000 + the
0205           // original barcode
0206           //
0207           HepMC::GenParticle* mothertwo = new HepMC::GenParticle(*mother);
0208           mothertwo->suggest_barcode(fSegmentations[motherID] * 10000000 + mother->barcode());
0209
0210           // we also reset the barcodes of the vertices
0211           //
0212           orig_mother_end_vtx->suggest_barcode(-fSegmentations[motherID] * 10000000
0213                                                - mother->barcode());
0214           childvtx->suggest_barcode(-mother->barcode());
0215
0216           // we attach it to the new vertex where interaction took place
0217           //
0218           childvtx->add_particle_out(mothertwo);
0219
0220           // and we attach it to the original endvertex
0221           //
0222           orig_mother_end_vtx->add_particle_in(mothertwo);
0223
0224           // and finally ... the increase the 'segmentation counter'
0225           //
0226           fSegmentations[motherID] = fSegmentations[motherID] + 1;
0227         }
0228       }
0229     }
0230     else
0231     // mother GenParticle is not there for some reason...
0232     // if this happens, we need to revise the philosophy...
0233     // a solution would be to create HepMC particles
0234     // at the begining of each track
0235     {
0236       G4cerr << "barcode " << motherID << " mother not there! " << G4endl;
0237     }
0238   }
0239   else  // primary
0240   {
0241     HepMC::GenVertex* primaryvtx = new HepMC::GenVertex(prodpos);
0242     primaryvtx->add_particle_out(particle);
0243     fEvent->add_vertex(primaryvtx);
0244
0245     // add id to the list of primaries
0246     //
0247     fPrimarybarcodes.push_back(partID);
0248   }
0249 }
0250
0251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0252
0253 void MCTruthManager::PrintEvent()
0254 {
0255   fEvent->print();
0256
0257   // looping over primaries and print the decay tree for each of them
0258   //
0259   for (std::vector<int>::const_iterator primarybar = fPrimarybarcodes.begin();
0260        primarybar != fPrimarybarcodes.end(); primarybar++)
0261   {
0262     PrintTree(fEvent->barcode_to_particle(*primarybar), " | ");
0263   }
0264 }
0265
0266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0267
0268 void MCTruthManager::PrintTree(HepMC::GenParticle* particle, G4String offset)
0269 {
0270   G4cout << offset << "---  barcode: " << particle->barcode() << " pdg: " << particle->pdg_id()
0271          << " energy: " << particle->momentum().e()
0272          << " production vertex: " << particle->production_vertex()->position().x() << ", "
0273          << particle->production_vertex()->position().y() << ", "
0274          << particle->production_vertex()->position().z() << ", "
0275          << particle->production_vertex()->position().t() << G4endl;
0276
0277   for (HepMC::GenVertex::particles_out_const_iterator it =
0278          particle->end_vertex()->particles_out_const_begin();
0279        it != particle->end_vertex()->particles_out_const_end(); it++)
0280   {
0281     G4String deltaoffset = "";
0282
0283     G4int curr = std::fmod(double((*it)->barcode()), 10000000.);
0284     G4int part = std::fmod(double(particle->barcode()), 10000000.);
0285     if (curr != part) {
0286       deltaoffset = " | ";
0287     }
0288
0289     PrintTree((*it), offset + deltaoffset);
0290   }
0291 }
0292
0293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....