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