PrimGenInterface/src/eASTHepMC3Interface.cc

0001 // ********************************************************************
0002 //
0003 // eASTHepMC3Interface.cc
0004 //   HepMC3 interface
0005 //
0006 // History
0007 //  June 23rd, 2021 : first empty class definition - Makoto Asai (SLAC)
0008 //
0009 // ********************************************************************
0010
0011 #include "eASTHepMC3Interface.hh"
0012
0013 #include "G4Event.hh"
0014 #include "G4PrimaryParticle.hh"
0015 #include "G4PrimaryVertex.hh"
0016 #include "G4GenericMessenger.hh"
0017 #include "G4SystemOfUnits.hh"
0018
0019 // The newer ReaderFactory is header-only and can be used for older versions
0020 // This file is copied verbatim from https://gitlab.cern.ch/hepmc/HepMC3
0021 // Copyright (C) 2014-2020 The HepMC collaboration, licensed under GPL3
0022 #include <HepMC3/Version.h>
0023 #if HEPMC3_VERSION_CODE < 3002004
0024 #include <HepMC_3_2_4_ReaderFactory.h>
0025 #else
0026 #include <HepMC3/ReaderFactory.h>
0027 #endif
0028
0029 #include <HepMC3/GenEvent.h>
0030 #include <HepMC3/GenVertex.h>
0031 #include <HepMC3/GenParticle.h>
0032
0033
0034 #include "G4AutoLock.hh"
0035 namespace { G4Mutex myHepMC3Mutex = G4MUTEX_INITIALIZER; }
0036
0037 eASTHepMC3Interface* eASTHepMC3Interface::instance = nullptr;
0038 G4bool eASTHepMC3Interface::instantiated = false;
0039 eASTHepMC3Interface* eASTHepMC3Interface::GetInstance()
0040 {
0041   G4AutoLock mlock(&myHepMC3Mutex);
0042   if(instance==nullptr && !instantiated)
0043   { instance = new eASTHepMC3Interface(); }
0044   return instance;
0045 }
0046
0047 eASTHepMC3Interface::eASTHepMC3Interface()
0048 {
0049   messenger = new G4GenericMessenger(this,"/eAST/generator/HepMC3/",
0050                   "HepMC3 interface commands");
0051
0052   auto& fileCmd = messenger->DeclareMethod("openFile",
0053           &eASTHepMC3Interface::OpenFile, "Open the HepMC3 input file");
0054   fileCmd.SetParameterName("fileName",false);
0055   fileCmd.SetToBeBroadcasted(false);
0056
0057   auto& vtxCmd = messenger->DeclareMethodWithUnit("vertex","mm",
0058           &eASTHepMC3Interface::SetVertexPosition,
0059           "Set vertex position");
0060   vtxCmd.SetParameterName("position",false);
0061   vtxCmd.SetToBeBroadcasted(false);
0062
0063   auto& timeCmd = messenger->DeclareMethodWithUnit("time0","s",
0064           &eASTHepMC3Interface::SetVertexTime,
0065           "Set t0 (start time at the primary vertex)");
0066   timeCmd.SetParameterName("t0",false);
0067   timeCmd.SetToBeBroadcasted(false);
0068
0069   auto& verCmd = messenger->DeclareMethod("verbose",
0070          &eASTHepMC3Interface::SetVerbose,"Verbose level");
0071   verCmd.SetParameterName("verboseLevel",true);
0072   verCmd.SetDefaultValue("0");
0073   verCmd.SetToBeBroadcasted(false);
0074
0075   instantiated = true;
0076 }
0077
0078 eASTHepMC3Interface::~eASTHepMC3Interface()
0079 {
0080   G4AutoLock mlock(&myHepMC3Mutex);
0081   delete messenger;
0082   instance = nullptr;
0083 }
0084
0085 G4bool eASTHepMC3Interface::OpenFile(G4String fName)
0086 {
0087   fileName = fName;
0088   if(verboseLevel>0){
0089     G4cout << "eASTHepMC3Interface - opening " << fileName << G4endl;
0090   }
0091
0092   // Make a new reader for every file.
0093   // Note: This supports a variety of formats (versions 1, 2, 3, HepEct, Root)
0094   //       and sources (files, streams, urls, ...)
0095   HepMC3Reader = HepMC3::deduce_reader(fileName);
0096
0097   if ( !HepMC3Reader || HepMC3Reader->failed() ) {
0098     G4Exception("eASTHepMC3Interface::GeneratePrimaryVertex","Event0201",
0099         FatalException, "eASTHepMC3Interface:: cannot open input.");
0100   }
0101   // if ( !HepMC3Reader || HepMC3Reader->failed() ) return false;
0102
0103   if(verboseLevel>0){
0104     G4cout << "eASTHepMC3Interface - " << fileName << " is open." << G4endl;
0105   }
0106
0107   return true;
0108 }
0109
0110 void eASTHepMC3Interface::GeneratePrimaryVertex(G4Event* g4event)
0111 {
0112   G4AutoLock mlock(&myHepMC3Mutex);
0113
0114   //read the event
0115   HepMC3::GenEvent hepevt(HepMC3::Units::GEV,HepMC3::Units::MM);
0116
0117   if (HepMC3Reader->failed()){
0118     // Todo: Don't know if there's a correct exceptionCode
0119     G4Exception("eASTHepMC3Interface::GeneratePrimaryVertex","Event0201",
0120         FatalException, "eASTHepMC3Interface:: cannot open input.");
0121   }
0122   if ( !HepMC3Reader->read_event(hepevt) ) return; // reached eof
0123
0124   // The root vertex is the default primary vertex
0125   // There can be multiple, unconnected graphs in the event.
0126   // In all cases I've seen, that is done for technical reasons,
0127   // the record still only has one event, i.e. they all share the same primary vertex
0128   auto pos = hepevt.event_pos();
0129   auto* g4vtx  = new G4PrimaryVertex(pos.x()*mm, pos.y()*mm, pos.z()*mm, 0);
0130
0131   // loop over particles
0132   // allowed statuses: 1 - final, 2 - decayed hadron/lepton
0133   // decay daughters will be enrolled by their mothers.
0134   // topological ordering in HepMC guarantees that mothers are
0135   // seen first, so only need to keep track of already created ones
0136   // to either ignore or reuse
0137   // using the std::map created_daughters;
0138   created_daughters.clear();
0139
0140   // better safe than sorry, detect faulty double-counting
0141   bool safetycheck=true;
0142   std::set<int> used;
0143
0144   for(auto hep_p : hepevt.particles()) {
0145
0146     auto status = hep_p->status();
0147     if ( status != 1 && status != 2  ) continue;
0148
0149     // already created?
0150     auto id = hep_p->id();
0151     auto finditer=created_daughters.find( id );
0152     G4PrimaryParticle* g4prim;
0153     if ( finditer != created_daughters.end() ){
0154       g4prim = finditer->second;
0155     } else {
0156       // doesn't exist already -> need to create
0157       // that also means it belongs to the primary vertex
0158       g4prim = MakeParticle ( hep_p, safetycheck, used );
0159       g4vtx->SetPrimary(g4prim);
0160     }
0161
0162     if ( hep_p->status() == 1 ){ // final
0163       // Nothing else to do, but should look for sanity
0164       if ( hep_p->children().size() !=0 ){
0165     // This should never happen
0166     G4Exception("eASTHepMC3Interface::GeneratePrimaryVertex","HepmcStatus",
0167             FatalException, "eASTHepMC3Interface:: status 1 but children.");
0168       }
0169       continue;
0170     }
0171
0172
0173     if ( hep_p->status() == 2 ){ // decayed, has daugthers
0174       if ( hep_p->children().size() ==0 ){
0175     // This should never happen
0176     G4Exception("eASTHepMC3Interface::GeneratePrimaryVertex","HepmcStatus",
0177             FatalException, "eASTHepMC3Interface:: status 2 but no children.");
0178       }
0179       // Now register daughters
0180       for ( auto hep_d : hep_p->children()) {
0181     auto id_d = hep_d->id();
0182     auto finditer_d=created_daughters.find( id_d );
0183     if ( finditer_d != created_daughters.end() ){
0184       // This should never happen
0185       G4cout << " PROBLEM: This primary with id = " << hep_d->id()
0186          << ", mother id = " << hep_p->id()
0187          << ", pid = " << hep_d->pid()
0188          << ", status = " << hep_d->status()
0189          << " and 4-momentum = " << hep_d->momentum().px() << " " << hep_d->momentum().py() << " , " << hep_d->momentum().pz() << " , " << hep_d->momentum().e()
0190        << G4endl;
0191       G4cout << " PROBLEM: already there with: " << G4endl;
0192       finditer_d->second->Print();
0193
0194       G4Exception("eASTHepMC3Interface::GeneratePrimaryVertex","Daughter",
0195               FatalException, "eASTHepMC3Interface:: found daughter out of order.");
0196     }
0197     if ( hep_p->status() == 2 || hep_p->status() == 1 ){
0198       auto g4daughter = MakeParticle ( hep_d, safetycheck, used );
0199       g4prim->SetDaughter(g4daughter);
0200       // and register
0201       created_daughters[id_d] = g4daughter;
0202     }
0203       }
0204     } // end of status == 1
0205   }//particle loop
0206
0207   // clean the HepMC objects
0208   hepevt.clear();
0209
0210   // set the primary vertex to G4Event
0211   g4event->AddPrimaryVertex(g4vtx);
0212 }
0213
0214 G4PrimaryParticle* eASTHepMC3Interface::MakeParticle ( const HepMC3::ConstGenParticlePtr hep_p,
0215                                const bool safetycheck, std::set<int>& used){
0216   // Shouldn't see a particle more than once
0217   if ( safetycheck ){
0218     auto id = hep_p->id();
0219     if ( used.count ( id ) ){ // contains() is c++20...
0220       G4Exception("eASTHepMC3Interface::MakeParticle","Particle",
0221           FatalException, "eASTHepMC3Interface:: Double-counting particles.");
0222     } else {
0223       used.insert( id );
0224     }
0225   }
0226   auto p = hep_p->momentum();
0227   auto* g4prim = new G4PrimaryParticle(hep_p->pid(), p.x()*GeV, p.y()*GeV, p.z()*GeV);
0228   g4prim->SetPolarization(0, 0, 0);
0229
0230   if ( verboseLevel > 1) {
0231     G4cout << " Made primary with pid = " << hep_p->pid()
0232        << ", status = " << hep_p->status()
0233        << " and 4-momentum = " << p.px() << " " << p.py() << " , " << p.pz() << " , " << p.e()
0234        << G4endl;
0235   }
0236
0237   return g4prim;
0238 }
0239