13 #include "Geant4/G4Event.hh" 14 #include "Geant4/G4IonTable.hh" 15 #include "Geant4/G4ParticleDefinition.hh" 16 #include "Geant4/G4ParticleTable.hh" 17 #include "Geant4/G4PrimaryParticle.hh" 18 #include "Geant4/G4PrimaryVertex.hh" 19 #include "Geant4/G4Types.hh" 21 #include "CLHEP/Vector/LorentzVector.h" 23 #include "TLorentzVector.h" 32 constexpr
bool genie_specific(
int const pdg)
34 constexpr
int genieLoPdg = 2000000001;
35 constexpr
int genieHiPdg = 2000000300;
36 return pdg >= genieLoPdg && pdg <= genieHiPdg;
43 : PrimaryGeneratorActionBase(p.
get<string>(
"name",
"MCTruthEventActionService"))
50 std::ostringstream badtxt;
52 badtxt <<
"\n Unknown PDG code = " << pdg <<
", seen " <<
n <<
" times.";
53 if (genie_specific(pdg)) { badtxt <<
" (GENIE specific)"; }
56 mf::LogWarning(
"MCTruthEventAction") <<
"The following unknown PDG codes were present in the " 57 <<
"simb::MCTruth input.\n" 58 <<
"They were not processed by Geant4." << badtxt.str();
63 std::ostringstream non1txt;
65 non1txt <<
"\n PDG code = " << pdg <<
", seen " <<
n <<
" times.";
66 if (genie_specific(pdg)) { non1txt <<
" (GENIE specific)"; }
69 mf::LogDebug(
"MCTruthEventAction") <<
"The following PDG codes were present in the " 70 <<
"simb::MCTruth input with Status != 1.\n" 71 <<
"They were not processed by Geant4." << non1txt.str();
76 std::ostringstream goodtxt;
78 goodtxt <<
"\n PDG code = " << pdg <<
", seen " <<
n <<
" times.";
80 mf::LogDebug(
"MCTruthEventAction") <<
"The following PDG codes were present in the " 81 <<
"simb::MCTruth input with Status = 1 " 82 <<
"and were processed by Geant4" << goodtxt.str();
92 auto const& mclistHandles = *
fMCLists;
94 size_t mclSize = mclistHandles.size();
95 mf::LogDebug(
"generatePrimaries") <<
"MCTruth Handles Size: " << mclSize;
99 std::map<CLHEP::HepLorentzVector, G4PrimaryVertex*> vertexMap;
100 for (
size_t mcl = 0; mcl < mclSize; ++mcl) {
102 <<
"MCTruth Handle Number: " << (mcl + 1) <<
" of " << mclSize;
106 for (
size_t i = 0; i < mclistHandle->size(); ++i) {
110 <<
"art::Ptr number " << (i + 1) <<
" of " << mclistHandle->size() <<
", Ptr: " << mclist;
112 MF_LOG_INFO(
"generatePrimaries") <<
"Generating " << nPart <<
" particles";
115 for (
int m = 0; m != nPart; ++m) {
117 G4int
const pdgCode = particle.
PdgCode();
121 <<
"Status code != 1, skipping particle number with MCTruth index = " 122 <<
"(" << mcl <<
"," << i <<
")" 123 <<
" and particle index = " << m;
128 if (((m + 1) % nPart) < 2)
130 mf::LogDebug(
"generatePrimaries") <<
"Particle Number: " << (m + 1) <<
" of " << nPart;
132 mf::LogDebug(
"generatePrimaries") <<
"index: " << index;
135 G4double
x = particle.
Vx() * CLHEP::cm;
136 G4double
y = particle.
Vy() * CLHEP::cm;
137 G4double
z = particle.
Vz() * CLHEP::cm;
138 G4double t = particle.
T() * CLHEP::ns;
141 CLHEP::HepLorentzVector fourpos(x, y, z, t);
144 G4PrimaryVertex*
vertex =
nullptr;
145 auto result = vertexMap.find(fourpos);
146 if (result == vertexMap.end()) {
148 vertex =
new G4PrimaryVertex(x, y, z, t);
149 vertexMap[fourpos] = vertex;
152 anEvent->AddPrimaryVertex(vertex);
156 vertex = result->second;
160 TLorentzVector
const& momentum = particle.
Momentum();
169 G4ParticleDefinition* particleDefinition =
nullptr;
171 if (pdgCode == 0) { particleDefinition =
fParticleTable->FindParticle(
"opticalphoton"); }
176 if (pdgCode >= 2
'000'000
'000) { // If the particle is generator-specific 177 mf::LogDebug("ConvertPrimaryToGeant4") 178 << ": %%% Will skip particle with generator-specific PDG code = " << pdgCode; 179 fUnknownPDG[pdgCode] += 1; 182 if (pdgCode > 1'000
'000'000) {
184 <<
": %%% Nuclear PDG code = " << pdgCode <<
" (x,y,z,t)=(" << x <<
"," << y <<
"," << z
186 <<
" P=" << momentum.P() <<
", E=" << momentum.E();
190 if (!particleDefinition) {
191 int const Z = (pdgCode % 10
'000'000) / 10
'000; // atomic number 192 int const A = (pdgCode % 10'000) / 10;
193 particleDefinition =
fParticleTable->GetIonTable()->GetIon(Z, A, 0.);
199 if (particleDefinition ==
nullptr) {
200 mf::LogDebug(
"ConvertPrimaryToGeant4") <<
": %%% Code not found = " << pdgCode;
208 auto* g4particle =
new G4PrimaryParticle(particleDefinition,
209 momentum.Px() * CLHEP::GeV,
210 momentum.Py() * CLHEP::GeV,
211 momentum.Pz() * CLHEP::GeV);
214 g4particle->SetCharge(particleDefinition->GetPDGCharge());
215 g4particle->SetPolarization(polarization.x(), polarization.y(), polarization.z());
218 vertex->SetPrimary(g4particle);
230 g4particle->SetUserInformation(primaryParticleInfo);
const TVector3 & Polarization() const
~MCTruthEventActionService()
static G4ParticleTable * fParticleTable
Geant4's table of particle definitions.
std::vector< art::Handle< std::vector< simb::MCTruth > > > const * fMCLists
MCTruthCollection input lists.
void generatePrimaries(G4Event *anEvent) override
double T(const int i=0) const
std::map< G4int, G4int > fProcessedPDG
PDG codes processed.
#define MF_LOG_INFO(category)
std::map< G4int, G4int > fNon1StatusPDG
PDG codes skipped because not status 1.
const simb::MCParticle & GetParticle(int i) const
double Vx(const int i=0) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
std::map< G4int, G4int > fUnknownPDG
map of unknown PDG codes to instances
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
const TLorentzVector & Momentum(const int i=0) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double Vz(const int i=0) const
void SetMCTruth(const simb::MCTruth *m, size_t idx=0, GeneratedParticleIndex_t indexInTruth=simb::NoGeneratedParticleIndex)
MCTruthEventActionService(fhicl::ParameterSet const &)
double Vy(const int i=0) const