From 0a43242877b7238b892ab8b767bd441139023fc4 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 24 Sep 2026 17:34:12 +0200 Subject: [PATCH] Monopoles: changes from further review - Register mplIoni through G4PhysicsListHelper. - Move the TPC Bethe-Bloch computation back into the branch that uses it. - Run the ITS monopole PDG lookup only when G4.monopole is set. - Define the monopole particles only when G4.monopole is set. - Make non-positive monopole magnetic charge or mass fatal. - Make a mismatch between G4.monopoleMass and TDatabasePDG fatal. - Warn when the TPC drift field is zero. - Align the species comment with the code. Co-Authored-By: Claude Opus 5.5 --- .../SimulationDataFormat/MonopoleParticles.h | 6 +++-- .../SimulationDataFormat/O2DatabasePDG.h | 14 +++++++++++ .../ITSMFT/ITS/simulation/src/Detector.cxx | 12 +++++----- Detectors/TPC/simulation/src/Detector.cxx | 24 +++++++++---------- Detectors/gconfig/src/O2MonopolePhysics.cxx | 22 +++++++++++++---- Generators/src/GeneratorFactory.cxx | 7 ++++-- Steer/src/O2MCApplication.cxx | 9 ++++++- macro/o2sim.C | 7 ++++-- 8 files changed, 72 insertions(+), 29 deletions(-) diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h index 93b9f1f67f915..c1f546fd5e569 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h +++ b/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. diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index 2c9b966bd098a..ca1c2f4c2aa20 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -16,6 +16,7 @@ #ifndef O2_O2DATABASEPDG_H #define O2_O2DATABASEPDG_H +#include #include #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) { diff --git a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx index e5fffe20f6fa2..6f4f1fa5fe028 100644 --- a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx +++ b/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; } diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 8359e820a03a3..7e81863d3bff4 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/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(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(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(fMC->GetRandom()->Poisson(meanNcoll)); diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 86c5a4cd36f83..dc1605743d82e 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -45,6 +45,7 @@ #include #include #include +#include #include #include #include @@ -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("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 diff --git a/Generators/src/GeneratorFactory.cxx b/Generators/src/GeneratorFactory.cxx index 3c3644bb4c3bb..1b6adad30052e 100644 --- a/Generators/src/GeneratorFactory.cxx +++ b/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)) { diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index c805e0f1011d6..2c3647f9cd163 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/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); diff --git a/macro/o2sim.C b/macro/o2sim.C index b88251705873b..c99cd012d4d73 100644 --- a/macro/o2sim.C +++ b/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; {