Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-30 09:17:19

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 //
0011 //====================================================================
0012 
0013 // Geant4 include files
0014 #include <G4GlobalConfig.hh>
0015 #include <G4ParticleTable.hh>
0016 
0017 // Framework include files
0018 #include "HepMC3EventReader.h"
0019 #include <DD4hep/Printout.h>
0020 #include <DDG4/Geant4Primary.h>
0021 #include <DDG4/Geant4Context.h>
0022 #include <DDG4/Factories.h>
0023 
0024 #include <G4ParticleTable.hh>
0025 
0026 #include <HepMC3/FourVector.h>
0027 #include <HepMC3/GenEvent.h>
0028 #include <HepMC3/GenParticle.h>
0029 #include <HepMC3/GenVertex.h>
0030 #include <HepMC3/Units.h>
0031 
0032 #include <G4Event.hh>
0033 
0034 
0035 using dd4hep::sim::HEPMC3EventReader;
0036 using SGenParticle = std::shared_ptr<HepMC3::GenParticle>;
0037 using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
0038 
0039 namespace {
0040   template<class T> inline int GET_ENTRY(const std::map<T,int>& mcparts, T part) {
0041     auto ip = mcparts.find(part);
0042     if (ip == mcparts.end())  {
0043       throw std::runtime_error("Unknown particle identifier look-up!");
0044     }
0045     return (*ip).second;
0046   }
0047 }
0048 
0049 /// Initializing constructor
0050 HEPMC3EventReader::HEPMC3EventReader(const string& fileName): Geant4EventReader(fileName) {}
0051 
0052 /// Read an event and fill a vector of MCParticles.
0053 HEPMC3EventReader::EventReaderStatus
0054 HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles& particles) {
0055   HepMC3::GenEvent genEvent;
0056   std::map<SGenParticle, int> mcparts;
0057   std::vector<SGenParticle>  mcpcoll;
0058   EventReaderStatus ret = readGenEvent(event_number, genEvent);
0059 
0060   if ( ret != EVENT_READER_OK ) return ret;
0061 
0062   int NHEP = genEvent.particles().size();
0063   // check if there is at least one particle
0064   if ( NHEP == 0 ) return EVENT_READER_NO_PRIMARIES;
0065 
0066   mcpcoll.resize(NHEP,0);
0067   for(int i=0; i<NHEP; ++i ) {
0068     auto p = genEvent.particles().at(i);
0069     mcparts[p] = i;
0070     mcpcoll[i] = std::move(p);
0071   }
0072 
0073   //treat event attributes, flow[12]
0074   std::vector<std::map<int, int>> colorFlow(2, std::map<int, int>());
0075   std::map<std::string, std::string> eventAttributes {};
0076   for(auto const& attr: genEvent.attributes()){
0077     for(auto const& inAttr: attr.second){
0078       if(attr.first == m_flow1){
0079         colorFlow[0][inAttr.first] = std::atoi(inAttr.second->unparsed_string().c_str());
0080       } else if(attr.first == m_flow2){
0081         colorFlow[1][inAttr.first] = std::atoi(inAttr.second->unparsed_string().c_str());
0082       }
0083     }
0084   }
0085 
0086   double mom_unit = (genEvent.momentum_unit() == HepMC3::Units::MomentumUnit::GEV) ? CLHEP::GeV : CLHEP::MeV;
0087   double len_unit = (genEvent.length_unit() == HepMC3::Units::LengthUnit::MM) ? CLHEP::mm : CLHEP::cm;
0088   // build collection of MCParticles
0089   for(int i=0; i<NHEP; ++i )   {
0090     auto const& mcp = mcpcoll[i];
0091     Geant4ParticleHandle p(new Particle(i));
0092     auto const& mom = mcp->momentum();
0093     auto const& vsx = mcp->production_vertex()->position();
0094     auto const& vex = (mcp->end_vertex()) ?  mcp->end_vertex()->position() : HepMC3::FourVector();
0095     const float spin[]  = {0.0, 0.0, 0.0}; //mcp->getSpin(); // FIXME
0096     const int   color[] = {colorFlow[0][mcp->id()], colorFlow[1][mcp->id()]};
0097     const int   pdg     = mcp->pid();
0098     p->pdgID        = pdg;
0099     p->charge       = 0; // int(mcp->getCharge()*3.0); // FIXME
0100     p->psx          = mom.get_component(0) * mom_unit;
0101     p->psy          = mom.get_component(1) * mom_unit;
0102     p->psz          = mom.get_component(2) * mom_unit;
0103     p->time         = vsx.get_component(3) * len_unit / CLHEP::c_light;
0104     p->properTime   = vsx.get_component(3) * len_unit / CLHEP::c_light;
0105     p->vsx          = vsx.get_component(0) * len_unit;
0106     p->vsy          = vsx.get_component(1) * len_unit;
0107     p->vsz          = vsx.get_component(2) * len_unit;
0108     p->vex          = vex.get_component(0) * len_unit;
0109     p->vey          = vex.get_component(1) * len_unit;
0110     p->vez          = vex.get_component(2) * len_unit;
0111     p->process      = 0;
0112     p->spin[0]      = spin[0];
0113     p->spin[1]      = spin[1];
0114     p->spin[2]      = spin[2];
0115     p->colorFlow[0] = color[0];
0116     p->colorFlow[1] = color[1];
0117     p->mass         = mcp->generated_mass() * mom_unit;
0118     auto const &par = mcp->parents(), &dau=mcp->children();
0119     for(int num=dau.size(), k=0; k<num; ++k)
0120       p->daughters.insert(GET_ENTRY(mcparts,dau[k]));
0121     for(int num=par.size(), k=0; k<num; ++k)
0122       p->parents.insert(GET_ENTRY(mcparts,par[k]));
0123 
0124     PropertyMask status(p->status);
0125     int genStatus = mcp->status();
0126     // Copy raw generator status
0127     p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK;
0128     m_inputAction->setGeneratorStatus(genStatus, status);
0129 
0130     if ( p->parents.size() == 0 )  {
0131       // A particle without a parent in HepMC3 can only be (something like) a beam particle, and it is attached to the
0132       // root vertex, by default (0,0,0,0) and equal for all parent-less particles.  Therefore we can take the end vertex
0133       // of the parentless particle as the start vertex for outgoing particles.  Note that for a particle without end
0134       // vertex (such as in a particle gun), it defaults to (0,0,0,0). This cannot be fixed, the information simply isn't
0135       // in the HepMC file. Having a parent enforces a vertex, having no parent forbids a vertex.
0136 
0137       Geant4Vertex* vtx = new Geant4Vertex ;
0138       vertices.emplace_back( vtx );
0139 
0140       vtx->x    = vex.get_component(0);
0141       vtx->y    = vex.get_component(1);
0142       vtx->z    = vex.get_component(2);
0143       vtx->time = vex.get_component(3) * len_unit / CLHEP::c_light;
0144 
0145       vtx->out.insert(p->id) ;
0146     }
0147 
0148     particles.emplace_back(p);
0149   }
0150   return EVENT_READER_OK;
0151 }