Make monopoles mass configurable via parameter

Author: jackal1-66

Common/SimConfig/include/SimConfig/G4Params.h

@@ -14,6 +14,7 @@
 
 #include "CommonUtils/ConfigurableParam.h"
 #include "CommonUtils/ConfigurableParamHelper.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
 
 namespace o2
 {
@@ -67,6 +68,9 @@ struct G4Params : public o2::conf::ConfigurableParamHelper<G4Params> {
   // the only volume with an electric field, and its magnitude is the TPC's own
   // TPCGEMParam.ElectricField[0]
   float monopoleMagneticCharge = 1.f;
+  // Monopole mass in GeV, preliminarly applied to all the monopoles species
+  // To-do: allow different monopoles masses in the same run
+  float monopoleMass = static_cast<float>(o2::sim::MonopoleMassDefaultGeV);
 
   // Fast simulation. Empty fastSimModels (the default) disables the feature
   // entirely; see Detectors/gconfig/include/SimSetup/G4FastSimulation.h.

DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h

@@ -29,6 +29,9 @@ constexpr int MonopolePdgSymm = 4110000;
 /// monopole carrying opposite electric and magnetic charge
 constexpr int MonopolePdgAsymm = 4120000;
 
+/// Default monopole mass in GeV.
+constexpr double MonopoleMassDefaultGeV = 100.;
+
 /// true for the monopole and antimonopole of both species
 constexpr bool isMonopole(int pdg) noexcept
 {

DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h

@@ -19,6 +19,7 @@
 #include <string>
 #include "TDatabasePDG.h"
 #include "TParticlePDG.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
 
 namespace o2
 {
@@ -44,7 +45,13 @@ class O2DatabasePDG
   }
 
   // adds ALICE particles to a given TDatabasePDG instance
-  static void addALICEParticles(TDatabasePDG* db = TDatabasePDG::Instance());
+  // monopoleMass is the mass of all monopole species in GeV; it has to match
+  // G4Params.monopoleMass or the generator and the transport will disagree
+  // To fix: all quantities in the MCTrack dependent from GetMass() will get the default 100 GeV.
+  //         A possible fix would be to stop deriving the masses from the compiled-in table,
+  //         or a better solution should be found in the future
+  static void addALICEParticles(TDatabasePDG* db = TDatabasePDG::Instance(),
+                                double monopoleMass = o2::sim::MonopoleMassDefaultGeV);
   static void addParticlesFromExternalFile(TDatabasePDG* db);
 
   // get particle's (if any) mass
@@ -84,7 +91,7 @@ class O2DatabasePDG
 
 // by keeping this inline, we can use it in other parts of the code, for instance Framework or Analysis,
 // without needing to link against this library
-inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db)
+inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMass)
 {
   //
   // Add ALICE particles to the ROOT PDG data base
@@ -191,11 +198,11 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db)
 
   // BSM targeted inclusions
   // Monopoles with same electric and magnetic charge
-  db->AddParticle("Monopole_symm", "Monopole_symm", 100., kTRUE, 0.0, 0, "BSM", 4110000);
-  db->AddParticle("AntiMonopole_symm", "AntiMonopole_symm", 100., kTRUE, 0.0, 0, "BSM", -4110000);
+  db->AddParticle("Monopole_symm", "Monopole_symm", monopoleMass, kTRUE, 0.0, 0, "BSM", 4110000);
+  db->AddParticle("AntiMonopole_symm", "AntiMonopole_symm", monopoleMass, kTRUE, 0.0, 0, "BSM", -4110000);
   // Monopoles with opposite electric and magnetic charge
-  db->AddParticle("Monopole_asymm", "Monopole_asymm", 100., kTRUE, 0.0, 0, "BSM", 4120000);
-  db->AddParticle("AntiMonopole_asymm", "AntiMonopole_asymm", 100., kTRUE, 0.0, 0, "BSM", -4120000);
+  db->AddParticle("Monopole_asymm", "Monopole_asymm", monopoleMass, kTRUE, 0.0, 0, "BSM", 4120000);
+  db->AddParticle("AntiMonopole_asymm", "AntiMonopole_asymm", monopoleMass, kTRUE, 0.0, 0, "BSM", -4120000);
 
   // IONS
   //

Detectors/TPC/simulation/src/Detector.cxx

@@ -18,6 +18,7 @@
 #include "DetectorsBase/Stack.h"
 #include "SimulationDataFormat/TrackReference.h"
 #include "SimulationDataFormat/MonopoleParticles.h"
+#include "SimConfig/G4Params.h"
 
 #include "FairVolume.h" // for FairVolume
 
@@ -127,15 +128,16 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
   // To-do: add dyons case
   bool isMonopole = false;
   if (static_cast<int>(trackCharge) == 0) {
-    isMonopole = o2::sim::isMonopole(fMC->TrackPid());
+    static const bool sMonopoleIonisation = o2::conf::G4Params::Instance().monopole;
+    isMonopole = sMonopoleIonisation && o2::sim::isMonopole(fMC->TrackPid());
     if (!isMonopole) {
       // set a very large step size for neutral particles
       fMC->SetMaxStep(1.e10);
       return kFALSE; // take only charged particles
     }
     if (fMC->Edep() <= 0.) {
-      // The monopole deposits nothing when no ionisation process is attached to
-      // it (G4.monopole=0), so no hit to make.
+      // An ionising monopole depositing nothing on this step: there
+      // is no hit to create, and the step limit is deliberately left alone
       return kFALSE;
     }
   }

Detectors/gconfig/src/O2MonopolePhysics.cxx

@@ -420,11 +420,11 @@ class O2MonopolePhysics : public G4VUserPhysicsList
   {
     auto* table = G4ParticleTable::GetParticleTable();
 
-    // Built once and shared by all monopoles species: O2MonopoleEquation reads the
-    // sign of the magnetic charge off the track, so one chord finder serves
-    // monopoles and anti-monopoles
-    const MonopoleFieldSetup fieldSetup =
-      buildMonopoleFieldSetup(mMagneticCharge / CLHEP::eplus, tpcDriftFieldMagnitude());
+    // Built on the first monopole species actually found, then shared by the rest:
+    // O2MonopoleEquation reads the sign of the magnetic charge off the track, so
+    // one chord finder serves monopoles and anti-monopoles
+    // No monopoles in the particles table == no monopole ionisation attached
+    MonopoleFieldSetup fieldSetup;
 
     int nAttached = 0;
     for (int pdg : gMonopolePDGs) {
@@ -477,6 +477,9 @@ class O2MonopolePhysics : public G4VUserPhysicsList
       // Deflect the monopole in the field as well; without this only the energy
       // loss above would act and the monopole would fly straight through, since
       // its electric charge (and hence the usual Lorentz force) is zero.
+      if (fieldSetup.chordFinder == nullptr) {
+        fieldSetup = buildMonopoleFieldSetup(mMagneticCharge / CLHEP::eplus, tpcDriftFieldMagnitude());
+      }
       installMonopoleTransport(pmanager, particle, fieldSetup);
 
       ++nAttached;

Generators/src/GeneratorFactory.cxx

@@ -17,6 +17,7 @@
 #include <Generators/BoxGenerator.h>
 #include <fairlogger/Logger.h>
 #include <SimConfig/SimConfig.h>
+#include <SimConfig/G4Params.h>
 #include <Generators/GeneratorFromFile.h>
 #include <Generators/GeneratorTParticle.h>
 #include <Generators/GeneratorTParticleParam.h>
@@ -86,7 +87,9 @@ void GeneratorFactory::setPrimaryGenerator(o2::conf::SimConfig const& conf, Fair
 
   /** generators **/
 
-  o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance());
+  // Monopole configurable mass added to the PDG database
+  o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(),
+                                       o2::conf::G4Params::Instance().monopoleMass);
   auto genconfig = conf.getGenerator();
 #if defined(GENERATORS_WITH_PYTHIA8) && defined(GENERATORS_WITH_HEPMC3)
   if (GeneratorHybridParam::Instance().switchExtToHybrid && (genconfig.compare("external") == 0 || genconfig.compare("extgen") == 0)) {

Steer/src/O2MCApplication.cxx

@@ -1753,12 +1753,16 @@ void addSpecialParticles()
   TVirtualMC::GetMC()->DefineParticle(-900000020, "AntiSexaquark", kPTUndefined, 2.0, 0.0, 4.35e+17, "Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, -2, kTRUE);
 
   // BSM Monopoles
+  // The transported mass has to match the one the generator used; see
+  // G4Params.monopoleMass and o2::sim::MonopoleMassDefaultGeV.
+  // To-do: find a way to define multiple masses for monopoles species
+  const double monopoleMass = o2::conf::G4Params::Instance().monopoleMass;
   // Symmetric monopoles: same electric and magnetic charge
-  TVirtualMC::GetMC()->DefineParticle(4110000, "Monopole_symm", kPTHadron, 100., 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
-  TVirtualMC::GetMC()->DefineParticle(-4110000, "AntiMonopole_symm", kPTHadron, 100., 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
+  TVirtualMC::GetMC()->DefineParticle(4110000, "Monopole_symm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
+  TVirtualMC::GetMC()->DefineParticle(-4110000, "AntiMonopole_symm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
   // Asymmetric monopoles: opposite electric and magnetic charge
-  TVirtualMC::GetMC()->DefineParticle(4120000, "Monopole_asymm", kPTHadron, 100., 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
-  TVirtualMC::GetMC()->DefineParticle(-4120000, "AntiMonopole_asymm", kPTHadron, 100., 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
+  TVirtualMC::GetMC()->DefineParticle(4120000, "Monopole_asymm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
+  TVirtualMC::GetMC()->DefineParticle(-4120000, "AntiMonopole_asymm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
 }
 
 void O2MCApplicationBase::AddParticles()

macro/o2sim.C

@@ -15,6 +15,7 @@
 #include <Generators/GeneratorFactory.h>
 #include <Generators/Generator.h>
 #include "SimulationDataFormat/O2DatabasePDG.h"
+#include "SimConfig/G4Params.h"
 #include "SimulationDataFormat/MCEventHeader.h"
 #include <SimConfig/SimConfig.h>
 #include <SimConfig/SimParams.h>
@@ -178,7 +179,11 @@ FairRunSim* o2sim_init(bool asservice, bool evalmat = false)
   run->Init();
 
   // add ALICE particles to TDatabasePDG singleton
-  o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance());
+  // The monopole mass has to be the one the transport uses; in the worker
+  // (asservice) the GeneratorFactory call is skipped, so this is the only
+  // place that registers it.
+  o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(),
+                                       o2::conf::G4Params::Instance().monopoleMass);
 
   long runStart = timestamp;
   {