Merge pull request #196 from mirguest/master

Support geantino/chargedgeantino in the simulation framework

Generator/src/GtGunTool.cpp

@@ -86,10 +86,22 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
         double charge = 0;
 
         TParticlePDG* particle = db_pdg->GetParticle(particle_name.c_str());
+        info() << "According to the " << particle_name << " get " << particle << endmsg;
         if (particle) {
             pdgcode = particle->PdgCode();
             mass = particle->Mass(); // GeV
             charge = particle->Charge()/3; // in ROOT, it is in units of |e|/3
+
+        } else if (particle_name == "geantino") {
+            info() << "Create geantino ... " << endmsg;
+            pdgcode = 0;
+            mass = 0;
+            charge = 0;
+        } else if (particle_name == "chargedgeantino") {
+            info() << "Create chargedgeantino ... " << endmsg;
+            pdgcode = 0;
+            mass = 0;
+            charge = 1; // according to the definition in G4ChargedGeantino.cc
         } else {
             // guess it is pdg code
             pdgcode = atol(particle_name.c_str());
@@ -99,7 +111,18 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
             }
         }
 
-        double energy = m_energymins.value()[i]==m_energymaxs.value()[i] ? m_energymins.value()[i] : CLHEP::RandFlat::shoot(m_energymins.value()[i], m_energymaxs.value()[i]);
+        if (i>=m_energymins.value().size()) {
+            error() << "Mismatch EnergyMins" << endmsg;
+            return false;
+        }
+        if (i>=m_energymaxs.value().size()) {
+            error() << "Mismatch EnergyMaxs" << endmsg;
+            return false;
+        }
+        double energy_min = m_energymins.value()[i];
+        double energy_max = m_energymaxs.value()[i];
+
+        double energy = energy_min==energy_max ? energy_min : CLHEP::RandFlat::shoot(energy_min, energy_max);
 
         // create the MC particle
         edm4hep::MCParticle mcp = event.m_mc_vec.create();
@@ -126,6 +149,26 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
         
         // direction
         // by default, randomize the direction
+
+
+        if (i>=m_thetamins.value().size()) {
+            error() << "Mismatch ThetaMins" << endmsg;
+            return false;
+        }
+        if (i>=m_thetamaxs.value().size()) {
+            error() << "Mismatch ThetaMaxs" << endmsg;
+            return false;
+        }
+        if (i>=m_phimins.value().size()) {
+            error() << "Mismatch PhiMins" << endmsg;
+            return false;
+        }
+        if (i>=m_phimaxs.value().size()) {
+            error() << "Mismatch PhiMaxs" << endmsg;
+            return false;
+        }
+
+
         double theta = m_thetamins.value()[i]==m_thetamaxs.value()[i] ? m_thetamins.value()[i] : CLHEP::RandFlat::shoot(m_thetamins.value()[i], m_thetamaxs.value()[i]);
         double phi =   m_phimins  .value()[i]==m_phimaxs  .value()[i] ? m_phimins  .value()[i] : CLHEP::RandFlat::shoot(m_phimins  .value()[i], m_phimaxs  .value()[i]);
         double costheta = cos(theta*acos(-1)/180);

Simulation/DetSimCore/src/G4PrimaryCnvTool.cpp

@@ -40,8 +40,18 @@ bool G4PrimaryCnvTool::mutate(G4Event* anEvent) {
         // pdg/particle
         int pdgcode = p.getPDG();
         G4ParticleTable* particletbl = G4ParticleTable::GetParticleTable();
-        G4ParticleDefinition* particle_def = particletbl->FindParticle(pdgcode);
+        G4ParticleDefinition* particle_def = nullptr;
 
+        // handle the several exceptions
+        if (pdgcode == 0 && p.getCharge() == 0) {
+            // this is geantino
+            particle_def = particletbl->FindParticle("geantino");
+        } else if (pdgcode == 0 && p.getCharge() == 1) {
+            // this is chargedgeantino
+            particle_def = particletbl->FindParticle("chargedgeantino");
+        } else {
+            particle_def = particletbl->FindParticle(pdgcode);
+        }
         // momentum
         const edm4hep::Vector3f& momentum = p.getMomentum();
         G4PrimaryParticle* g4prim = new G4PrimaryParticle(particle_def,
@@ -49,6 +59,14 @@ bool G4PrimaryCnvTool::mutate(G4Event* anEvent) {
                                                           momentum.y*CLHEP::GeV,
                                                           momentum.z*CLHEP::GeV);
 
+        // modify the mass of the chargedgeantino
+        if (pdgcode == 0 && p.getCharge() == 1) {
+            info() << "The mass of G4ChargedGeantino is "
+                   << m_chargedgeantino_mass.value() 
+                   << endmsg;
+            g4prim->SetMass(m_chargedgeantino_mass.value());
+        }
+
         g4vtx->SetPrimary(g4prim);
 
         anEvent->AddPrimaryVertex(g4vtx);

Simulation/DetSimCore/src/G4PrimaryCnvTool.h

@@ -18,6 +18,7 @@ public:
 private:
     DataHandle<edm4hep::MCParticleCollection> m_mcParCol{"MCParticle", Gaudi::DataHandle::Reader, this};
 
+    Gaudi::Property<double> m_chargedgeantino_mass{this, "ChargedGeantinoMass"};
 };
 
 #endif