feat: enable geantino particle generation (#1888)

Allows shooting geantino particles into the detector with our particle gun Example gun: ``` addParticleGun( s, ParticleConfig(1, acts.PdgParticle.eInvalid, randomizeCharge=False, charge=1, mass=0), rnd=rnd, ) ```
acts-project · Feb 23, 2023 · 386611f · 386611f
1 parent fffb998
commit 386611f
Show file tree

Hide file tree

Showing 10 changed files with 45 additions and 30 deletions.
diff --git a/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/SimParticleTranslation.hpp b/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/SimParticleTranslation.hpp
@@ -36,10 +36,11 @@ class SimParticleTranslation final : public G4VUserPrimaryGeneratorAction {
     /// The input particle collection
     std::string inputParticles = "";
 
-    /// Force pdgCode & mass (this is needed for Geantino simulation)
-    bool forceParticle = false;
-    G4int forcedPdgCode = 998;
-    G4double forcedMass = 0.;
+    /// Force pdgCode & mass & charge in G4 units (this is needed for Geantino
+    /// simulation)
+    std::optional<G4int> forcedPdgCode;
+    std::optional<G4double> forcedCharge;  // e.g. 1 for charged geantino
+    std::optional<G4double> forcedMass;    // e.g. 0 for geantino
 
     /// The number of hits per particle to be expected
     /// @note best to include secondaries for that

diff --git a/Examples/Algorithms/Geant4/src/ParticleTrackingAction.cpp b/Examples/Algorithms/Geant4/src/ParticleTrackingAction.cpp
@@ -45,9 +45,9 @@ ActsExamples::SimParticle ActsExamples::ParticleTrackingAction::convert(
 
   // Get all the information from the Track
   const G4ParticleDefinition* particleDef = aTrack.GetParticleDefinition();
-  G4double mass = particleDef->GetPDGMass();
-  G4double charge = particleDef->GetPDGCharge();
   G4int pdg = particleDef->GetPDGEncoding();
+  G4double charge = particleDef->GetPDGCharge();
+  G4double mass = particleDef->GetPDGMass();
   G4int id = aTrack.GetTrackID();
   G4int parentId = aTrack.GetParentID();
   G4ThreeVector pPosition = convertLength * aTrack.GetPosition();

diff --git a/Examples/Algorithms/Geant4/src/SimParticleTranslation.cpp b/Examples/Algorithms/Geant4/src/SimParticleTranslation.cpp
@@ -95,24 +95,24 @@ void ActsExamples::SimParticleTranslation::GeneratePrimaries(G4Event* anEvent) {
     Acts::Vector4 mom4 = part.fourMomentum() * convertEnergy;
 
     // Particle properties, may be forced to specific value
-    G4int particlePdgCode =
-        m_cfg.forceParticle ? m_cfg.forcedPdgCode : part.pdg();
+    G4int particlePdgCode = m_cfg.forcedPdgCode.value_or(part.pdg());
+    G4double particleCharge = m_cfg.forcedCharge.value_or(part.charge());
     G4double particleMass =
-        m_cfg.forceParticle ? part.mass() * convertEnergy : 0.;
+        m_cfg.forcedMass.value_or(part.mass() * convertEnergy);
 
     // Check if it is a Geantino / ChargedGeantino
     G4ParticleDefinition* particleDefinition =
         particleTable->FindParticle(particlePdgCode);
     if (particleDefinition == nullptr) {
-      switch (particlePdgCode) {
-        case 999: {
-          particleDefinition = G4Geantino::Definition();
-        } break;
-        case 998: {
-          particleDefinition = G4ChargedGeantino::Definition();
-        } break;
-        default:
-          break;
+      if (particlePdgCode == 0 && particleMass == 0 && particleCharge == 0) {
+        particleDefinition = G4Geantino::Definition();
+      }
+      if (particlePdgCode == 0 && particleMass == 0 && particleCharge != 0) {
+        if (particleCharge != 1) {
+          ACTS_ERROR("invalid charged geantino charge " << particleCharge
+                                                        << ". should be 1");
+        }
+        particleDefinition = G4ChargedGeantino::Definition();
       }
     }
 
@@ -127,14 +127,14 @@ void ActsExamples::SimParticleTranslation::GeneratePrimaries(G4Event* anEvent) {
                  << particleDefinition->GetParticleName()
                  << "' and properties:");
     ACTS_VERBOSE(" -> mass: " << particleMass);
+    ACTS_VERBOSE(" -> charge: " << particleCharge);
     ACTS_VERBOSE(" -> momentum: " << mom4.transpose());
-    ACTS_VERBOSE(" -> charge: " << part.charge());
 
     G4PrimaryParticle* particle = new G4PrimaryParticle(particleDefinition);
 
     particle->SetMass(particleMass);
+    particle->SetCharge(particleCharge);
     particle->Set4Momentum(mom4[0], mom4[1], mom4[2], mom4[3]);
-    particle->SetCharge(part.charge());
     particle->SetTrackID(trackId++);
 
     // Add the primary to the vertex

diff --git a/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp b/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp
@@ -16,8 +16,8 @@
 ActsExamples::ParametricParticleGenerator::ParametricParticleGenerator(
     const Config& cfg)
     : m_cfg(cfg),
-      m_charge(ActsFatras::findCharge(m_cfg.pdg)),
-      m_mass(ActsFatras::findMass(m_cfg.pdg)),
+      m_charge(cfg.charge.value_or(ActsFatras::findCharge(m_cfg.pdg))),
+      m_mass(cfg.mass.value_or(ActsFatras::findMass(m_cfg.pdg))),
       // since we want to draw the direction uniform on the unit sphere, we must
       // draw from cos(theta) instead of theta. see e.g.
       // https://mathworld.wolfram.com/SpherePointPicking.html

diff --git a/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp b/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp
@@ -55,6 +55,11 @@ class ParametricParticleGenerator : public EventGenerator::ParticlesGenerator {
     bool randomizeCharge = false;
     /// Number of particles.
     size_t numParticles = 1;
+
+    /// Overrides particle charge.
+    std::optional<double> charge;
+    /// Overrides particle mass.
+    std::optional<double> mass;
   };
 
   ParametricParticleGenerator(const Config& cfg);

diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py
@@ -26,8 +26,8 @@
 PhiConfig = namedtuple("PhiConfig", ["min", "max"], defaults=[None, None])
 ParticleConfig = namedtuple(
     "ParticleConfig",
-    ["num", "pdg", "randomizeCharge"],
-    defaults=[None, None, None],
+    ["num", "pdg", "randomizeCharge", "charge", "mass"],
+    defaults=[None, None, None, None, None],
 )
 ParticleSelectorConfig = namedtuple(
     "ParticleSelectorConfig",
@@ -82,7 +82,7 @@ def addParticleGun(
         pseudorapidity configuration: eta min, eta max, uniform
     phiConfig : PhiConfig(min, max)
         azimuthal angle configuration: phi min, phi max
-    particleConfig : ParticleConfig(num, pdg, randomizeCharge)
+    particleConfig : ParticleConfig(num, pdg, randomizeCharge, charge, mass)
         particle configuration: number of particles, particle type, charge flip
     multiplicity : int, 1
         number of generated vertices
@@ -107,7 +107,8 @@ def addParticleGun(
                 multiplicity=FixedMultiplicityGenerator(n=multiplicity),
                 vertex=vtxGen
                 or acts.examples.GaussianVertexGenerator(
-                    stddev=acts.Vector4(0, 0, 0, 0), mean=acts.Vector4(0, 0, 0, 0)
+                    mean=acts.Vector4(0, 0, 0, 0),
+                    stddev=acts.Vector4(0, 0, 0, 0),
                 ),
                 particles=acts.examples.ParametricParticleGenerator(
                     **acts.examples.defaultKWArgs(
@@ -119,6 +120,8 @@ def addParticleGun(
                         numParticles=particleConfig.num,
                         pdg=particleConfig.pdg,
                         randomizeCharge=particleConfig.randomizeCharge,
+                        charge=particleConfig.charge,
+                        mass=particleConfig.mass,
                     )
                 ),
             )

diff --git a/Examples/Python/src/Geant4Component.cpp b/Examples/Python/src/Geant4Component.cpp
@@ -116,9 +116,9 @@ PYBIND11_MODULE(ActsPythonBindingsGeant4, mod) {
         // Read the particle from the generator
         SimParticleTranslation::Config g4PrCfg;
         g4PrCfg.inputParticles = inputParticles;
-        g4PrCfg.forceParticle = true;
+        g4PrCfg.forcedPdgCode = 0;
+        g4PrCfg.forcedCharge = 0.;
         g4PrCfg.forcedMass = 0.;
-        g4PrCfg.forcedPdgCode = 999;
 
         auto g4Cfg = makeGeant4Config(level, std::move(randomNumbers), detector,
                                       physicsList, g4PrCfg);

diff --git a/Examples/Python/src/Generators.cpp b/Examples/Python/src/Generators.cpp
@@ -131,6 +131,8 @@ void addGenerators(Context& ctx) {
         .def_readwrite("pdg", &Config::pdg)
         .def_readwrite("randomizeCharge", &Config::randomizeCharge)
         .def_readwrite("numParticles", &Config::numParticles)
+        .def_readwrite("mass", &Config::mass)
+        .def_readwrite("charge", &Config::charge)
         .def_property(
             "p",
             [](Config& cfg) {

diff --git a/Examples/Run/Geant4/Common/src/Geant4Common.cpp b/Examples/Run/Geant4/Common/src/Geant4Common.cpp
@@ -75,9 +75,9 @@ void setupGeant4Simulation(
   g4PrCfg.inputParticles = materialRecording ? Simulation::kParticlesInitial
                                              : Simulation::kParticlesSelection;
   if (materialRecording) {
-    g4PrCfg.forceParticle = true;
+    g4PrCfg.forcedPdgCode = 0;
+    g4PrCfg.forcedCharge = 0.;
     g4PrCfg.forcedMass = 0.;
-    g4PrCfg.forcedPdgCode = 999;
     // Set the material tracks at output
     g4Cfg.outputMaterialTracks = Simulation::kMaterialTracks;
   }

diff --git a/Examples/Scripts/Python/material_recording.py b/Examples/Scripts/Python/material_recording.py
@@ -43,6 +43,10 @@ def runMaterialRecording(g4geo, outputDir, tracksPerEvent=10000, s=None):
                     mean=acts.Vector4(0, 0, 0, 0),
                 ),
                 particles=ParametricParticleGenerator(
+                    pdg=acts.PdgParticle.eInvalid,
+                    charge=0,
+                    randomizeCharge=False,
+                    mass=0,
                     p=(1 * u.GeV, 10 * u.GeV),
                     eta=(-4, 4),
                     numParticles=tracksPerEvent,