g4simulation/g4decayer/G4EvtGenDecayer.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 // $Id: G4EvtGenDecayer.cc,v 1.2 2014/10/07 03:06:54 mccumber Exp $
0027 //
0028 /// \file eventgenerator/pythia/decayer6/src/G4EvtGenDecayer.cc
0029 /// \brief Implementation of the G4EvtGenDecayer class
0030
0031 // ----------------------------------------------------------------------------
0032 // According to TPythia6Decayer class in Root:
0033 // http://root.cern.ch/
0034 // see http://root.cern.ch/root/License.html
0035 // ----------------------------------------------------------------------------
0036
0037 #include "G4EvtGenDecayer.hh"
0038 #include "PHEvtGenRandomEngine.hh"  // for PHEvtGenRandomEngine
0039
0040 #include <EvtGen/EvtGen.hh>                      // for EvtGen
0041 #include <EvtGenBase/EvtAbsRadCorr.hh>           // for EvtAbsRadCorr
0042 #include <EvtGenBase/EvtHepMCEvent.hh>           // for EvtHepMCEvent
0043 #include <EvtGenBase/EvtId.hh>                   // for EvtId
0044 #include <EvtGenBase/EvtPDL.hh>                  // for EvtPDL
0045 #include <EvtGenBase/EvtParticle.hh>             // for EvtParticle
0046 #include <EvtGenBase/EvtParticleFactory.hh>      // for EvtParticleFactory
0047 #include <EvtGenBase/EvtRandom.hh>               // for EvtRandom
0048 #include <EvtGenBase/EvtVector4R.hh>             // for EvtVector4R
0049 #include <EvtGenExternal/EvtExternalGenList.hh>  // for EvtExternalGenList
0050
0051 #include <HepMC3/Attribute.h>        // for string
0052 #include <HepMC3/FourVector.h>       // for FourVector
0053 #include <HepMC3/GenEvent.h>         // for GenEvent
0054 #include <HepMC3/GenParticle.h>      // for GenParticle
0055 #include <HepMC3/GenParticle_fwd.h>  // for GenParticlePtr
0056 #include <HepMC3/GenVertex.h>        // for GenVertex
0057 #include <HepMC3/GenVertex_fwd.h>    // for GenVertexPtr
0058 #include <HepMC3/Print.h>            // for Print
0059
0060 #include <Geant4/G4DecayProducts.hh>
0061 #include <Geant4/G4DynamicParticle.hh>
0062 #include <Geant4/G4ParticleDefinition.hh>  // for G4ParticleDefinition
0063 #include <Geant4/G4ParticleTable.hh>
0064 #include <Geant4/G4String.hh>  // for G4String
0065 #include <Geant4/G4SystemOfUnits.hh>
0066 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0067 #include <Geant4/G4Track.hh>
0068 #include <Geant4/G4Types.hh>        // for G4int, G4double
0069 #include <Geant4/G4VExtDecayer.hh>  // for G4VExtDecayer
0070
0071 #include <algorithm>  // for copy, max
0072 #include <cstdlib>    // for abs
0073 #include <iostream>   // for operator<<, basic_ostream:...
0074 #include <memory>     // for __shared_ptr_access
0075 #include <stack>      // for operator<<
0076 #include <string>     // for operator<<
0077 #include <vector>     // for vector
0078
0079 class EvtDecayBase;
0080
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082
0083 G4EvtGenDecayer::G4EvtGenDecayer()
0084   : G4VExtDecayer("G4EvtGenDecayer")
0085   , fVerboseLevel(0)
0086 {
0087   mEvtGenRandomEngine = new PHEvtGenRandomEngine();
0088
0089   EvtRandom::setRandomEngine(static_cast<EvtRandomEngine*>(mEvtGenRandomEngine));
0090
0091   radCorrEngine = genList.getPhotosModel();
0092
0093   extraModels = genList.getListOfModels();
0094
0095   // the hardcoded paths are temporary
0096   const char* offline_main = getenv("OFFLINE_MAIN");
0097
0098   if (!offline_main)
0099   {
0100     exit(1);
0101   }
0102
0103   string decay = string(offline_main) + "/share/EvtGen/DECAY.DEC";  // Using PDG 2019 reference as the input for now
0104   string evt = string(offline_main) + "/share/EvtGen/evt.pdl";
0105
0106   mEvtGen = new EvtGen(decay, evt, static_cast<EvtRandomEngine*>(mEvtGenRandomEngine), radCorrEngine, &extraModels);
0107   extraModels.clear();
0108 }
0109
0110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0111
0112 G4EvtGenDecayer::~G4EvtGenDecayer()
0113 {
0114   delete mEvtGen;
0115
0116   delete radCorrEngine;
0117 }
0118
0119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0120
0121 G4ParticleDefinition* G4EvtGenDecayer::
0122     GetParticleDefinition(int ParPDGID) const
0123 {
0124   /// Return G4 particle definition for given TParticle
0125
0126   // get particle definition from G4ParticleTable
0127   G4int pdgEncoding = ParPDGID;
0128   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0129   G4ParticleDefinition* particleDefinition = nullptr;
0130   if (pdgEncoding != 0)
0131   {
0132     particleDefinition = particleTable->FindParticle(pdgEncoding);
0133   }
0134
0135   return particleDefinition;
0136 }
0137
0138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0139
0140 G4DecayProducts* G4EvtGenDecayer::ImportDecayProducts(const G4Track& track)
0141 {
0142   /// Import decay products
0143
0144   G4ThreeVector momentum = track.GetMomentum();
0145   G4double etot = track.GetDynamicParticle()->GetTotalEnergy();
0146
0147   G4ParticleDefinition* particleDef = track.GetDefinition();
0148   G4int pdgEncoding = particleDef->GetPDGEncoding();
0149
0150   EvtVector4R p_init(etot / GeV, momentum.x() / GeV, momentum.y() / GeV, momentum.z() / GeV);
0151
0152   EvtId parentID = EvtPDL::evtIdFromLundKC(pdgEncoding);
0153
0154   auto mParticle = EvtParticleFactory::particleFactory(parentID, p_init);
0155   // Decay the particle
0156   mEvtGen->generateDecay(mParticle);
0157
0158   EvtHepMCEvent theEvent;
0159
0160   theEvent.constructEvent(mParticle);
0161
0162   HepMC3::GenEvent* evt = theEvent.getEvent();  // Directly use HepMC3 records
0163
0164   G4DecayProducts* decayProducts = new G4DecayProducts(*(track.GetDynamicParticle()));
0165
0166   std::stack<HepMC3::GenParticlePtr> stack;
0167
0168   auto part = evt->particles()[0];
0169
0170   if (part->pdg_id() != pdgEncoding)
0171   {
0172     std::cout << "Issue Found: The first particle in the decay chain is NOT the incoming particle decayed by EvtGen" << std::endl;
0173   }
0174
0175   stack.push(part);
0176
0177   while (!stack.empty())
0178   {
0179     auto particle = stack.top();
0180     stack.pop();
0181
0182     for (const auto& children : particle->children())
0183     {
0184       float EvtWidth = 9999;
0185       int pdg = children->pdg_id();
0186
0187       G4ParticleDefinition* particleDefinition = GetParticleDefinition(pdg);
0188       if (EvtPDL::evtIdFromLundKC(pdg).getId() != -1)
0189       {
0190         EvtWidth = EvtPDL::getWidth(EvtPDL::evtIdFromLundKC(pdg)) * 1000000;  // Decay Width Unit: GeV -> keV to define stable particles in EvtGen -> G4
0191       }
0192
0193       if (particleDefinition && EvtWidth < WidthThreshold)
0194       {
0195         bool SameVtxPos = true;
0196
0197         if (children->production_vertex()->position().x() != particle->end_vertex()->position().x() || children->production_vertex()->position().y() != particle->end_vertex()->position().y() || children->production_vertex()->position().z() != particle->end_vertex()->position().z() || children->production_vertex()->position().t() != particle->end_vertex()->position().t())
0198         {
0199           SameVtxPos = false;
0200         }
0201         if (!SameVtxPos)
0202         {
0203           std::cout << "Issue Found: in this vertex, particles pushed to GEANT have different production positions, need to check and understand why!!!" << std::endl;
0204           std::cout << "This is due to the particle PDGID: " << children->pdg_id() << "  from the decay vertex of the parent particle:" << particle->pdg_id() << std::endl;
0205           std::cout << "Here is the Event Content" << std::endl;
0206           const HepMC3::GenEvent& CheckEvent = *evt;
0207           HepMC3::Print::content(CheckEvent);
0208         }
0209
0210         G4ThreeVector G4Mom = G4ThreeVector(children->momentum().px() * GeV, children->momentum().py() * GeV, children->momentum().pz() * GeV);
0211         G4DynamicParticle* dynamicParticle = new G4DynamicParticle(particleDefinition, G4Mom);
0212
0213         if (dynamicParticle)
0214         {
0215           if (fVerboseLevel > 0)
0216           {
0217             std::cout << "  G4 particle name: "
0218                       << dynamicParticle->GetDefinition()->GetParticleName()
0219                       << std::endl;
0220           }
0221
0222           decayProducts->PushProducts(dynamicParticle);
0223         }
0224         else
0225         {
0226           std::cout << "Dynamical particle to be pushed from EvtGen to GEANT 4 is not properly defined: need to check why this happens!" << std::endl;
0227           exit(1);
0228         }
0229       }
0230       else
0231       {
0232         stack.push(children);
0233       }
0234     }
0235   }
0236
0237   evt->clear();
0238   mParticle->deleteTree();
0239   return decayProducts;
0240 }
0241
0242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0243
0244 void G4EvtGenDecayer::SetDecayTable(const string& decayTable, bool useXml)
0245 {
0246   mEvtGen->readUDecay(decayTable, useXml);
0247 }