File indexing completed on 2025-01-30 09:17:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <G4GlobalConfig.hh>
0015 #include <G4ParticleTable.hh>
0016
0017
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
0050 HEPMC3EventReader::HEPMC3EventReader(const string& fileName): Geant4EventReader(fileName) {}
0051
0052
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
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
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
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};
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;
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
0127 p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK;
0128 m_inputAction->setGeneratorStatus(genStatus, status);
0129
0130 if ( p->parents.size() == 0 ) {
0131
0132
0133
0134
0135
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 }