Monopoles: changes from further review

DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h

@@ -24,9 +24,11 @@
 namespace o2::sim
 {
 
-/// monopole carrying equal electric and magnetic charge
+/// Both species are electrically neutral for now, and the transport gives them the
+/// same magnetic charge, with its sign taken from the PDG sign.
+/// monopole intended to carry equal electric and magnetic charge
 constexpr int MonopolePdgSymm = 4110000;
-/// monopole carrying opposite electric and magnetic charge
+/// monopole intended to carry opposite electric and magnetic charge (not yet distinguished)
 constexpr int MonopolePdgAsymm = 4120000;
 
 /// Default monopole mass in GeV.

DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h

@@ -16,6 +16,7 @@
 #ifndef O2_O2DATABASEPDG_H
 #define O2_O2DATABASEPDG_H
 
+#include <cmath>
 #include <string>
 #include "TDatabasePDG.h"
 #include "TParticlePDG.h"
@@ -54,6 +55,19 @@ class O2DatabasePDG
                                 double monopoleMass = o2::sim::MonopoleMassDefaultGeV);
   static void addParticlesFromExternalFile(TDatabasePDG* db);
 
+  // true if all monopole species are registered with the given mass in GeV; TDatabasePDG
+  // keeps the first registration, so an earlier call with another mass makes this false
+  static bool hasMonopoleMass(TDatabasePDG* db, double monopoleMass)
+  {
+    for (int pdg : {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm, o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm}) {
+      const auto* particle = db->GetParticle(pdg);
+      if (particle == nullptr || std::abs(particle->Mass() - monopoleMass) > 1.e-6 * monopoleMass) {
+        return false;
+      }
+    }
+    return true;
+  }
+
   // get particle's (if any) mass
   static Double_t MassImpl(TParticlePDG* particle, bool& success)
   {

Detectors/ITSMFT/ITS/simulation/src/Detector.cxx

@@ -25,6 +25,7 @@
 #include "DetectorsBase/Stack.h"
 #include "SimulationDataFormat/TrackReference.h"
 #include "SimulationDataFormat/MonopoleParticles.h"
+#include "SimConfig/G4Params.h"
 #include "fairlogger/Logger.h" // for LOG, LOG_IF
 
 // FairRoot includes
@@ -321,11 +322,12 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
   // This method is called from the MC stepping
   // Electrically neutral magnetic monopoles deposit energy in the
   // silicon through G4mplIonisation (Ahlen stopping power), so they must not be
-  // rejected by the electric-charge gate. PDG lookup
-  // never runs for ordinary charged production.
+  // rejected by the electric-charge gate. The PDG lookup only runs for neutral
+  // particles in runs with monopole physics enabled.
   // To-do: handle dyons.
+  static const bool sMonopoleIonisation = o2::conf::G4Params::Instance().monopole;
   const bool isNeutral = (fMC->TrackCharge() == 0);
-  const bool isMonopole = isNeutral && o2::sim::isMonopole(fMC->TrackPid());
+  const bool isMonopole = isNeutral && sMonopoleIonisation && o2::sim::isMonopole(fMC->TrackPid());
   if (isNeutral && !isMonopole) {
     return kFALSE;
   }
@@ -394,9 +396,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
     mTrackData.mHitStarted = true;
   }
   if (stopHit) {
-    // A monopole reaches this point even when no ionisation process is attached to
-    // it (G4.monopole=0), in which case it crosses the sensor depositing nothing.
-    // Storing such empty hits would only inflate the hit file, so they are skipped
+    // A monopole that deposited nothing in the sensor leaves no hit
     if (isMonopole && mTrackData.mEnergyLoss <= 0.) {
       return kFALSE;
     }

Detectors/TPC/simulation/src/Detector.cxx

@@ -216,18 +216,6 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
   Int_t numberOfElectrons = 0;
   // I.H. - the type expected in addHit is float
 
-  // ---| Stepsize in cm |---
-  const double stepSize = fMC->TrackStep();
-
-  double betaGamma = momentum.P() / fMC->TrackMass();
-  betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg
-
-  // ---| number of primary ionisations per cm |---
-  const double primaryElectronsPerCM =
-    gasParam.Nprim * BetheBlochAleph(static_cast<float>(betaGamma), gasParam.BetheBlochParam[0],
-                                     gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2],
-                                     gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]);
-
   // use Geant4 energy deposit directly for ionisation (Kr-83m calibration simulations)
   if (detParam.UseGeant4Edep) {
     // We have multiple collisions and add fluctuations: smear nel using
@@ -248,6 +236,18 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
     // 2^24 (16777216) ==> largest integer a IEEE-754 float can represent exactly
     numberOfElectrons = TMath::Min(numberOfElectrons, 16777216);
   } else {
+    // ---| Stepsize in cm |---
+    const double stepSize = fMC->TrackStep();
+
+    double betaGamma = momentum.P() / fMC->TrackMass();
+    betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg
+
+    // ---| number of primary ionisations per cm |---
+    const double primaryElectronsPerCM =
+      gasParam.Nprim * BetheBlochAleph(static_cast<float>(betaGamma), gasParam.BetheBlochParam[0],
+                                       gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2],
+                                       gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]);
+
     // ---| mean number of collisions and random for this event |---
     const double meanNcoll = stepSize * trackCharge * trackCharge * primaryElectronsPerCM;
     const int nColl = static_cast<int>(fMC->GetRandom()->Poisson(meanNcoll));

Detectors/gconfig/src/O2MonopolePhysics.cxx

@@ -45,6 +45,7 @@
 #include <G4VUserPhysicsList.hh>
 #include <G4ParticleTable.hh>
 #include <G4ParticleDefinition.hh>
+#include <G4PhysicsListHelper.hh>
 #include <G4ProcessManager.hh>
 #include <G4ProcessVector.hh>
 #include <G4VEnergyLossProcess.hh>
@@ -92,7 +93,8 @@ inline bool isMonopolePDG(int pdg) { return o2::sim::isMonopole(pdg); }
 /// so that G4Transportation queries the field for it
 constexpr G4double gMonopoleFieldGateMoment = 1.0e-20;
 
-// Extent of the TPC drift region
+// Extent of the TPC field cage, i.e. the region with the drift field (smaller than
+// the TPC_Drift volume, which extends beyond the endcaps and the field-cage rods)
 constexpr G4double gTPCFieldCageRMin = 83.5 * CLHEP::cm;
 constexpr G4double gTPCFieldCageRMax = 254.5 * CLHEP::cm;
 constexpr G4double gTPCFieldCageZMax = 249.525 * CLHEP::cm;
@@ -107,6 +109,11 @@ inline double tpcDriftFieldMagnitude()
   }
   try {
     const double valueKVPerCm = o2::conf::ConfigurableParam::getValueAs<float>("TPCGEMParam.ElectricField[0]");
+    if (valueKVPerCm <= 0.) {
+      LOG(warn) << "O2MonopolePhysics: TPCGEMParam.ElectricField[0] = " << valueKVPerCm
+                << " kV/cm; the monopole is not coupled to a drift field";
+      return 0.;
+    }
     return valueKVPerCm * CLHEP::kilovolt / CLHEP::cm;
   } catch (...) {
     LOG(warn) << "O2MonopolePhysics: the TPC is in the geometry but TPCGEMParam is not "
@@ -159,6 +166,9 @@ class O2MonopoleEquation : public G4EquationOfMotion
     // Non-zero only while an O2 monopole is being transported. The sign follows
     // the PDG sign, so monopole and anti-monopole are pushed in opposite
     // directions along B.
+    // TODO: both species are electrically neutral today, so their magnetic sign is
+    // the same function of the PDG sign; once they carry electric charge the
+    // asymmetric species (4120000) needs the opposite relative sign.
     double signedMagneticCharge = 0.;
     int pdg = 0;
     if (const G4Track* track = currentTrack()) {
@@ -462,9 +472,8 @@ class O2MonopolePhysics : public G4VUserPhysicsList
           pmanager->RemoveProcess(proc);
         }
       }
-      // Ordering (AtRest, AlongStep, PostStep) = (-1, 1, 1) as in the Geant4
-      // monopole example: continuous-and-discrete energy loss, not active at rest.
-      pmanager->AddProcess(new G4mplIonisation(mMagneticCharge), -1, 1, 1);
+      // Ordered by G4PhysicsListHelper from the process subtype, as in the Geant4 monopole example
+      G4PhysicsListHelper::GetPhysicsListHelper()->RegisterProcess(new G4mplIonisation(mMagneticCharge), particle);
 
       // G4Transportation only looks up the field when the particle has a
       // non-zero electric charge or a non-zero magnetic moment (μ) and a monopole has
@@ -542,6 +551,11 @@ TG4RunConfiguration* createMonopoleRunConfiguration(const TString& userGeometry,
   // One Dirac charge g_D = eplus / (2*alpha) ~= 68.5 eplus (Dirac quantisation).
   // magneticChargeDirac counts Dirac charges (1.0 == classic single monopole),
   // matching the "g_1" convention used by the O2 monopole generator input.
+  // G4mplIonisation silently replaces a zero charge with one Dirac charge, while the
+  // equation of motion would use the zero, so non-positive values are rejected
+  if (!(magneticChargeDirac > 0.)) {
+    LOG(fatal) << "O2MonopolePhysics: G4.monopoleMagneticCharge must be positive, got " << magneticChargeDirac;
+  }
   const double gDirac = CLHEP::eplus / (2.0 * CLHEP::fine_structure_const);
   const double magneticChargeEplus = magneticChargeDirac * gDirac;
   LOG(info) << "Monopole ionisation enabled: magnetic charge = " << magneticChargeDirac

Generators/src/GeneratorFactory.cxx

@@ -88,8 +88,11 @@ void GeneratorFactory::setPrimaryGenerator(o2::conf::SimConfig const& conf, Fair
   /** generators **/
 
   // Monopole configurable mass added to the PDG database
-  o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(),
-                                       o2::conf::G4Params::Instance().monopoleMass);
+  const double monopoleMass = o2::conf::G4Params::Instance().monopoleMass;
+  o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(), monopoleMass);
+  if (!o2::O2DatabasePDG::hasMonopoleMass(TDatabasePDG::Instance(), monopoleMass)) {
+    LOG(fatal) << "Monopoles were registered in TDatabasePDG with a mass other than G4.monopoleMass = " << monopoleMass << " GeV";
+  }
   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

@@ -1752,11 +1752,18 @@ void addSpecialParticles()
   TVirtualMC::GetMC()->DefineParticle(900000020, "Sexaquark", kPTUndefined, 2.0, 0.0, 4.35e+17, "Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kTRUE);
   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
+  // BSM Monopoles, defined only when their physics is enabled (G4.monopole=1), so that a
+  // monopole primary without it is reported by the engine instead of crossing the detector unseen.
   // 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
+  if (!o2::conf::G4Params::Instance().monopole) {
+    return;
+  }
   const double monopoleMass = o2::conf::G4Params::Instance().monopoleMass;
+  if (!(monopoleMass > 0.)) {
+    LOG(fatal) << "G4.monopoleMass must be positive, got " << monopoleMass;
+  }
   // Symmetric monopoles: same electric and magnetic charge
   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);

macro/o2sim.C

@@ -182,8 +182,11 @@ FairRunSim* o2sim_init(bool asservice, bool evalmat = false)
   // 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);
+  const double monopoleMass = o2::conf::G4Params::Instance().monopoleMass;
+  o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(), monopoleMass);
+  if (!o2::O2DatabasePDG::hasMonopoleMass(TDatabasePDG::Instance(), monopoleMass)) {
+    LOG(fatal) << "Monopoles were registered in TDatabasePDG with a mass other than G4.monopoleMass = " << monopoleMass << " GeV";
+  }
 
   long runStart = timestamp;
   {