From 0fb2da00ee7e7b09db981ed308a60f19900e3086 Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Wed, 15 Jul 2026 16:11:37 +0200 Subject: [PATCH 1/6] Integration of monopoles in O2 framework --- Common/SimConfig/include/SimConfig/G4Params.h | 10 + .../SimulationDataFormat/O2DatabasePDG.h | 8 + Detectors/TPC/simulation/src/Detector.cxx | 73 ++++--- Detectors/gconfig/CMakeLists.txt | 2 +- Detectors/gconfig/g4Config.C | 15 +- .../include/SimSetup/O2MonopolePhysics.h | 51 +++++ Detectors/gconfig/src/O2MonopolePhysics.cxx | 195 ++++++++++++++++++ Steer/src/O2MCApplication.cxx | 8 + 8 files changed, 332 insertions(+), 30 deletions(-) create mode 100644 Detectors/gconfig/include/SimSetup/O2MonopolePhysics.h create mode 100644 Detectors/gconfig/src/O2MonopolePhysics.cxx diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 2a333a39e4242..be71960427589 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -54,6 +54,16 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { bool g4scoring = false; bool g4fluenceweight = false; + + // Enable magnetic-monopole ionisation as defined in Detectors/gconfig/O2MonopolePhysics. + // The G4mplIonisation process (Ahlen stopping power) is attached to + // the O2 monopole particles (PDG +-4110000 / +-4120000) on top of the chosen + // reference physics list. Off by default + bool monopole = false; + // Monopole magnetic charge in units of the Dirac charge g_D = eplus/(2*alpha) + // (~68.5 eplus). 1.0 corresponds to a single classic Dirac monopole. + float monopoleMagneticCharge = 1.f; + O2ParamDef(G4Params, "G4"); }; diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index 9111f548bb9d4..52cb5173fc3f0 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -189,6 +189,14 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db) db->AddParticle("CHI2P_B2", " ", 10.269, kFALSE, 0.0, 0, "meson", 100555); db->AddParticle("UPSLON4S", " ", 10.580, kFALSE, 0.0, 0, "meson", 300553); + // 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); + // 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); + // IONS // // Done by default now from Pythia6 table diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 1a7c0fc25802b..09ee899e324c8 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -112,11 +112,18 @@ Bool_t Detector::ProcessHits(FairVolume* vol) /* This method is called from the MC stepping for the sensitive volume only */ // LOG(info) << "tpc::ProcessHits"; const double trackCharge = fMC->TrackCharge(); + // Magnetic monopoles have zero electric charge but ionise the gas through + // their magnetic charge energy loss (G4mplIonisation). + const int trackPdg = fMC->TrackPid(); + const bool isMonopole = (TMath::Abs(trackPdg) == 4110000 || TMath::Abs(trackPdg) == 4120000); if (static_cast(trackCharge) == 0) { - - // set a very large step size for neutral particles - fMC->SetMaxStep(1.e10); - return kFALSE; // take only charged particles + // Fall through only for monopoles when ionisation is enabled. + // The behaviour for the other neutral particles remains as before. + if (!isMonopole || fMC->Edep() <= 0.) { + // set a very large step size for neutral particles + fMC->SetMaxStep(1.e10); + return kFALSE; // take only charged particles + } } // ===| SET THE LENGTH OF THE NEXT ENERGY LOSS STEP |========================= @@ -205,29 +212,41 @@ Bool_t Detector::ProcessHits(FairVolume* vol) 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)); - - // Variables needed to generate random powerlaw distributed energy loss - const double alpha_p1 = 1. - gasParam.Exp; // NA49/G3 value - const double oneOverAlpha_p1 = 1. / alpha_p1; - const double eMin = gasParam.Ipot; - const double eMax = gasParam.Eend; - const double kMin = TMath::Power(eMin, alpha_p1); - const double kMax = TMath::Power(eMax, alpha_p1); - const double wIon = gasParam.Wion; - - for (Int_t n = 0; n < nColl; n++) { - // Use GEANT3 / NA49 expression: - // P(eDep) ~ k * edep^-gasParam.getExp() - // eMin(~I) < eDep < eMax(300 electrons) - // k fixed so that Int_Emin^EMax P(Edep) = 1. - const double rndm = fMC->GetRandom()->Rndm(); - const double eDep = TMath::Power((kMax - kMin) * rndm + kMin, oneOverAlpha_p1); - int nel_step = static_cast(((eDep - eMin) / wIon) + 1); - nel_step = TMath::Min(nel_step, 300); // 300 electrons corresponds to 10 keV - numberOfElectrons += nel_step; + if (isMonopole) { + // ---| MONOPOLE IONISATION |-------- + // A magnetic monopole ionises the gas via G4mplIonisation (Ahlen stopping + // power), which is not described by the electric-charge model (no e-charge) + // Ionisation electrons calculated directly from the energy deposited in this step: Nel = Edep / Wion. + // To check by TPC experts if this is actually the best way... + numberOfElectrons = static_cast(fMC->Edep() / static_cast(gasParam.Wion)); + // The number of electrons is stored as a short in the HitGroup; cap it to + // avoid overflow for the very high monopole dE/dx. + numberOfElectrons = TMath::Min(numberOfElectrons, 32000); + } else { + // ---| 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)); + + // Variables needed to generate random powerlaw distributed energy loss + const double alpha_p1 = 1. - gasParam.Exp; // NA49/G3 value + const double oneOverAlpha_p1 = 1. / alpha_p1; + const double eMin = gasParam.Ipot; + const double eMax = gasParam.Eend; + const double kMin = TMath::Power(eMin, alpha_p1); + const double kMax = TMath::Power(eMax, alpha_p1); + const double wIon = gasParam.Wion; + + for (Int_t n = 0; n < nColl; n++) { + // Use GEANT3 / NA49 expression: + // P(eDep) ~ k * edep^-gasParam.getExp() + // eMin(~I) < eDep < eMax(300 electrons) + // k fixed so that Int_Emin^EMax P(Edep) = 1. + const double rndm = fMC->GetRandom()->Rndm(); + const double eDep = TMath::Power((kMax - kMin) * rndm + kMin, oneOverAlpha_p1); + int nel_step = static_cast(((eDep - eMin) / wIon) + 1); + nel_step = TMath::Min(nel_step, 300); // 300 electrons corresponds to 10 keV + numberOfElectrons += nel_step; + } } // LOG(info) << "tpc::AddHit" << FairLogger::endl << "Eloss: " diff --git a/Detectors/gconfig/CMakeLists.txt b/Detectors/gconfig/CMakeLists.txt index 444f125b7cbb7..f06936c4c6bc9 100644 --- a/Detectors/gconfig/CMakeLists.txt +++ b/Detectors/gconfig/CMakeLists.txt @@ -15,7 +15,7 @@ o2_add_library(G3Setup ) o2_add_library(G4Setup - SOURCES src/G4Config.cxx + SOURCES src/G4Config.cxx src/O2MonopolePhysics.cxx PUBLIC_LINK_LIBRARIES MC::Geant4VMC MC::Geant4 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup ) diff --git a/Detectors/gconfig/g4Config.C b/Detectors/gconfig/g4Config.C index 16374cf9fd4a3..625ad2b93762b 100644 --- a/Detectors/gconfig/g4Config.C +++ b/Detectors/gconfig/g4Config.C @@ -61,6 +61,7 @@ R__LOAD_LIBRARY(libgeant4vmc) #include "TG4RunConfiguration.h" #include "SimConfig/G4Params.h" #include "SimConfig/FluenceWeightCalculator.h" +#include "SimSetup/O2MonopolePhysics.h" #endif #include "commonConfig.C" @@ -114,8 +115,18 @@ void Config() LOG(fatal) << "Unsupported geometry navigation mode"; } - auto runConfiguration = new TG4RunConfiguration(geomNavStr, physicsSetup, "stepLimiter+specialCuts", - specialStacking, mtMode); + TG4RunConfiguration* runConfiguration = nullptr; + if (g4Params.monopole) { + // Reference physics list + magnetic-monopole ionisation attached + // The reference physics is unchanged. + std::cout << "Monopole ionisation physics requested (G4.monopole=1)\n"; + runConfiguration = o2::g4config::createMonopoleRunConfiguration( + geomNavStr, physicsSetup, "stepLimiter+specialCuts", specialStacking, mtMode, + g4Params.monopoleMagneticCharge); + } else { + runConfiguration = new TG4RunConfiguration(geomNavStr, physicsSetup, "stepLimiter+specialCuts", + specialStacking, mtMode); + } if (g4Params.g4scoring) { runConfiguration->SetUseOfG4Scoring(); if (g4Params.g4fluenceweight) { diff --git a/Detectors/gconfig/include/SimSetup/O2MonopolePhysics.h b/Detectors/gconfig/include/SimSetup/O2MonopolePhysics.h new file mode 100644 index 0000000000000..46f67623d2191 --- /dev/null +++ b/Detectors/gconfig/include/SimSetup/O2MonopolePhysics.h @@ -0,0 +1,51 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \author M+Giacalone - July 2026 + +#ifndef O2_G4CONFIG_O2MONOPOLEPHYSICS_H_ +#define O2_G4CONFIG_O2MONOPOLEPHYSICS_H_ + +#include + +class TG4RunConfiguration; + +namespace o2 +{ +namespace g4config +{ + +/// Build a Geant4-VMC run configuration which attaches the magnetic-monopole ionisation process +/// (G4mplIonisation) to the requested reference physics list. +/// the monopole particles already defined by O2 (PDG +-4110000 / +-4120000). +/// +/// Only used when G4Params.monopole == true; +/// By default, the standard TG4RunConfiguration is created +/// +/// \param userGeometry VMC geometry-navigation string (as for TG4RunConfiguration) +/// \param physicsList the reference physics-list selection string +/// \param specialProcess the VMC special-process selection string +/// \param specialStacking VMC special-stacking flag +/// \param mtApplication multithreading flag +/// \param magneticChargeDirac magnetic charge of the monopole in units of the +/// Dirac charge g_D = eplus/(2*alpha) (~68.5 eplus); +/// 1.0 (default) reproduces the classic single Dirac monopole +TG4RunConfiguration* createMonopoleRunConfiguration(const TString& userGeometry, + const TString& physicsList, + const TString& specialProcess, + bool specialStacking, + bool mtApplication, + double magneticChargeDirac); + +} // namespace g4config +} // namespace o2 + +#endif // O2_G4CONFIG_O2MONOPOLEPHYSICS_H_ diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx new file mode 100644 index 0000000000000..4fba37a067678 --- /dev/null +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -0,0 +1,195 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file O2MonopolePhysics.cxx +/// \brief Opt-in magnetic-monopole ionisation physics for the O2 Geant4 engine. +/// \author M+Giacalone - July 2026 +/// +/// O2 defines the BSM monopole particles (PDG +-4110000 / +-4120000) via +/// O2MCApplication::AddParticles(), so they are transported by Geant4, but no +/// energy-loss process is attached to them by any of the stock reference +/// physics lists. +/// +/// This file adds - the magnetic-monopole ionisation process G4mplIonisation +/// (Ahlen stopping power, Rev. Mod. Phys. 52 (1980) 121) to those particles, when requested. +/// The implementation takes Geant4 `examples/extended/exoticphysics/monopole` example (and the +/// equivalent geant4_vmc "monopole" physics builder). The main difference is that +/// here the process is attached to the pre-existing O2 monopole particles (with correct PDG ID() +/// rather than to a freshly created G4Monopole +/// +/// G4mplIonisation and G4mplIonisationWithDeltaModel are taken directly from the standard +/// Geant4 libraries (libG4processes) + +#include "SimSetup/O2MonopolePhysics.h" + +#include + +#include "TG4RunConfiguration.h" +#include "TG4ComposedPhysicsList.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace o2 +{ +namespace g4config +{ + +namespace +{ +// PDG codes of the O2 monopole species, as defined in +// O2MCApplication::AddParticles() and O2DatabasePDG: +// +-4110000 : "symmetric" monopoles +// +-4120000 : "asymmetric" monopoles +constexpr std::array gMonopolePDGs = {4110000, -4110000, 4120000, -4120000}; +} // namespace + +//____________________________________________________________________________ +/// Minimal physics list whose only job is to attach the magnetic-monopole +/// ionisation process to the already-defined O2 monopole particles. +/// +/// It is implemented as a G4VUserPhysicsList so that it can be handed to VMC's +/// TG4ComposedPhysicsList::AddPhysicsList(). It intentionally does NOT create +/// any particle and does NOT add transportation (both are handled by the +/// reference physics list), it only adds one extra process. +class O2MonopolePhysics : public G4VUserPhysicsList +{ + public: + /// \param magneticChargeEplus monopole magnetic charge in Geant4 internal + /// (eplus) units + explicit O2MonopolePhysics(double magneticChargeEplus) : mMagneticCharge(magneticChargeEplus) {} + + // Particles are already defined by O2MCApplication::AddParticles(). + void ConstructParticle() override {} + + // Cuts are handled by the reference physics list. + void SetCuts() override {} + + void ConstructProcess() override + { + auto* table = G4ParticleTable::GetParticleTable(); + int nAttached = 0; + for (int pdg : gMonopolePDGs) { + auto* particle = table->FindParticle(pdg); + if (particle == nullptr) { + // Not necessarily an error: only the species actually produced by the + // generator strictly need this. Report and continue. + LOG(info) << "O2MonopolePhysics: monopole PDG " << pdg + << " not present in the particle table, skipping"; + continue; + } + auto* pmanager = particle->GetProcessManager(); + if (pmanager == nullptr) { + LOG(warning) << "O2MonopolePhysics: no process manager for PDG " << pdg << ", skipping"; + continue; + } + // The reference EM physics list attaches the standard hadron ionisation + // (hIoni) to the monopole, which is wrong and makes G4LossTableManager crash + // while merging the two dE/dx tables. + // Any pre-existing energy-loss process is removed so that mplIoni is the + // monopole single, correct ionisation process. This is taken from Geant4 + // exoticphysics/monopole example + { + std::vector toRemove; + G4ProcessVector* plist = pmanager->GetProcessList(); + for (G4int ip = 0; ip < static_cast(plist->size()); ++ip) { + G4VProcess* proc = (*plist)[ip]; + if (dynamic_cast(proc) != nullptr) { + toRemove.push_back(proc); + } + } + for (auto* proc : toRemove) { + LOG(info) << "O2MonopolePhysics: removing pre-existing energy-loss process " + << proc->GetProcessName() << " from " << particle->GetParticleName(); + 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); + ++nAttached; + LOG(info) << "O2MonopolePhysics: attached G4mplIonisation to " + << particle->GetParticleName() << " (PDG " << pdg + << "), magnetic charge = " << mMagneticCharge / CLHEP::eplus << " eplus"; + } + if (nAttached == 0) { + LOG(warning) << "O2MonopolePhysics: no monopole particle found; no ionisation attached"; + } + } + + private: + double mMagneticCharge; ///< magnetic charge in Geant4 internal (eplus) units +}; + +//____________________________________________________________________________ +/// TG4RunConfiguration that appends O2MonopolePhysics to the composed physics +/// list which VMC builds for the requested reference list. +class O2G4RunConfiguration : public TG4RunConfiguration +{ + public: + O2G4RunConfiguration(const TString& userGeometry, const TString& physicsList, + const TString& specialProcess, Bool_t specialStacking, + Bool_t mtApplication, double magneticChargeEplus) + : TG4RunConfiguration(userGeometry, physicsList, specialProcess, specialStacking, mtApplication), + mMagneticCharge(magneticChargeEplus) + { + } + + G4VUserPhysicsList* CreatePhysicsList() override + { + G4VUserPhysicsList* physicsList = TG4RunConfiguration::CreatePhysicsList(); + if (auto* composed = dynamic_cast(physicsList)) { + composed->AddPhysicsList(new O2MonopolePhysics(mMagneticCharge)); + LOG(info) << "O2G4RunConfiguration: monopole ionisation physics registered " + "on the composed physics list"; + } else { + LOG(error) << "O2G4RunConfiguration: physics list is not a TG4ComposedPhysicsList, " + "monopole ionisation could NOT be enabled"; + } + return physicsList; + } + + private: + double mMagneticCharge; ///< magnetic charge in Geant4 internal (eplus) units +}; + +//____________________________________________________________________________ +TG4RunConfiguration* createMonopoleRunConfiguration(const TString& userGeometry, + const TString& physicsList, + const TString& specialProcess, + bool specialStacking, + bool mtApplication, + double magneticChargeDirac) +{ + // 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. + const double gDirac = CLHEP::eplus / (2.0 * CLHEP::fine_structure_const); + const double magneticChargeEplus = magneticChargeDirac * gDirac; + LOG(info) << "Monopole ionisation enabled: magnetic charge = " << magneticChargeDirac + << " Dirac charge(s) = " << magneticChargeEplus / CLHEP::eplus << " eplus"; + return new O2G4RunConfiguration(userGeometry, physicsList, specialProcess, + specialStacking, mtApplication, magneticChargeEplus); +} + +} // namespace g4config +} // namespace o2 diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index 1e3f925042d01..5577db119daa1 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/Steer/src/O2MCApplication.cxx @@ -1602,6 +1602,14 @@ void addSpecialParticles() //Sexaquark (uuddss): compact, neutral and stable hypothetical bound state (arxiv.org/abs/1708.08951) 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 + // 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); + // 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); } void O2MCApplicationBase::AddParticles() From ad1ae639fc2b05846fc07e5f98bab52cf52a9a2c Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Tue, 21 Jul 2026 14:17:55 +0200 Subject: [PATCH 2/6] Include full Lorentz dual force for monopoles --- Common/SimConfig/include/SimConfig/G4Params.h | 4 + Detectors/TPC/simulation/src/Detector.cxx | 10 +- Detectors/gconfig/src/O2MonopolePhysics.cxx | 243 ++++++++++++++++++ 3 files changed, 252 insertions(+), 5 deletions(-) diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index be71960427589..b84afc9c33b50 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -62,6 +62,10 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { bool monopole = false; // Monopole magnetic charge in units of the Dirac charge g_D = eplus/(2*alpha) // (~68.5 eplus). 1.0 corresponds to a single classic Dirac monopole. + // The monopole is transported with the full dual Lorentz force + // F = g*(B - v x E/c^2). The electric part acts inside the TPC field cage, + // the only volume with an electric field, and its magnitude is the TPC's own + // TPCGEMParam.ElectricField[0] float monopoleMagneticCharge = 1.f; O2ParamDef(G4Params, "G4"); diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 09ee899e324c8..47fe376074909 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -198,7 +198,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // TODO: Add discussion about drawback Int_t numberOfElectrons = 0; - // I.H. - the type expected in addHit is short + // I.H. - the type expected in addHit is float // ---| Stepsize in cm |--- const double stepSize = fMC->TrackStep(); @@ -219,9 +219,9 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // Ionisation electrons calculated directly from the energy deposited in this step: Nel = Edep / Wion. // To check by TPC experts if this is actually the best way... numberOfElectrons = static_cast(fMC->Edep() / static_cast(gasParam.Wion)); - // The number of electrons is stored as a short in the HitGroup; cap it to - // avoid overflow for the very high monopole dE/dx. - numberOfElectrons = TMath::Min(numberOfElectrons, 32000); + // The number of electrons is stored as a float in the HitGroup: maximum cap at + // 2^24 (16777216) ==> largest integer a IEEE-754 float can represent exactly + numberOfElectrons = TMath::Min(numberOfElectrons, 16777216); } else { // ---| mean number of collisions and random for this event |--- const double meanNcoll = stepSize * trackCharge * trackCharge * primaryElectronsPerCM; @@ -244,7 +244,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol) const double rndm = fMC->GetRandom()->Rndm(); const double eDep = TMath::Power((kMax - kMin) * rndm + kMin, oneOverAlpha_p1); int nel_step = static_cast(((eDep - eMin) / wIon) + 1); - nel_step = TMath::Min(nel_step, 300); // 300 electrons corresponds to 10 keV + nel_step = TMath::Min(nel_step, gasParam.MaxElePerStep); // 300 electrons corresponds to 10 keV numberOfElectrons += nel_step; } } diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 4fba37a067678..d0815b854d22d 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -29,6 +29,11 @@ /// Geant4 libraries (libG4processes) #include "SimSetup/O2MonopolePhysics.h" +#include "CommonUtils/ConfigurableParam.h" + +#include // needed to instantiate getValueAs<> + +#include #include @@ -43,10 +48,24 @@ #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + #include #include #include +#include #include namespace o2 @@ -61,8 +80,213 @@ namespace // +-4110000 : "symmetric" monopoles // +-4120000 : "asymmetric" monopoles constexpr std::array gMonopolePDGs = {4110000, -4110000, 4120000, -4120000}; + +inline bool isMonopolePDG(int pdg) +{ + const int abspdg = std::abs(pdg); + return abspdg == 4110000 || abspdg == 4120000; +} + +/// GEANT4 only queries the magnetic field when a particle has non-zero electric +/// charge or non-zero magnetic moment (μ) and the monopole has neither. So this is a workaround +/// so that G4Transportation queries the field for it +constexpr G4double gMonopoleFieldGateMoment = 1.0e-20; + +// Extent of the TPC drift region +constexpr G4double gTPCFieldCageRMin = 83.5 * CLHEP::cm; +constexpr G4double gTPCFieldCageRMax = 254.5 * CLHEP::cm; +constexpr G4double gTPCFieldCageZMax = 249.525 * CLHEP::cm; + +/// Magnitude of the TPC drift field, in Geant4 units; 0 when there is no TPC. +inline double tpcDriftFieldMagnitude() +{ + if (gGeoManager == nullptr || gGeoManager->GetVolume("TPC_Drift") == nullptr) { + LOG(info) << "O2MonopolePhysics: no TPC in the geometry of this run, " + "the monopole is not coupled to a drift field"; + return 0.; + } + try { + const double valueKVPerCm = o2::conf::ConfigurableParam::getValueAs("TPCGEMParam.ElectricField[0]"); + return valueKVPerCm * CLHEP::kilovolt / CLHEP::cm; + } catch (...) { + LOG(warn) << "O2MonopolePhysics: the TPC is in the geometry but TPCGEMParam is not " + "registered; no drift-field coupling for the monopole"; + return 0.; + } +} + +/// Electric field of the TPC drift region at a specific position. +/// Returns false outside the field cage. +inline bool tpcDriftField(const G4double position[3], G4double driftField, G4double E[3]) +{ + const G4double z = position[2]; + if (std::fabs(z) >= gTPCFieldCageZMax) { + return false; + } + const G4double r2 = position[0] * position[0] + position[1] * position[1]; + if (r2 < gTPCFieldCageRMin * gTPCFieldCageRMin || r2 > gTPCFieldCageRMax * gTPCFieldCageRMax) { + return false; + } + E[0] = 0.; + E[1] = 0.; + E[2] = (z >= 0.) ? -driftField : driftField; + return true; +} } // namespace +//____________________________________________________________________________ +/// Equation of motion implementing the dual Lorentz force of a magnetic charge. +class O2MonopoleEquation : public G4EquationOfMotion +{ + public: + /// \param field the magnetic field to integrate in + /// \param magneticChargeEplusUnits monopole magnetic charge expressed in eplus + /// units (one Dirac charge = 1/(2*alpha) ~ 68.5) + /// \param tpcDriftFieldGeant4 TPC drift field in Geant4 units + O2MonopoleEquation(G4Field* field, double magneticChargeEplusUnits, double tpcDriftFieldGeant4) + : G4EquationOfMotion(field), + mMagneticChargeEplus(magneticChargeEplusUnits), + mTPCDriftField(tpcDriftFieldGeant4) + { + } + + void SetChargeMomentumMass(G4ChargeState particleChargeState, + G4double /*momentum*/, G4double particleMass) override + { + mElCharge = CLHEP::eplus * particleChargeState.GetCharge() * CLHEP::c_light; + mMassCof = particleMass * particleMass; + + // 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. + double signedMagneticCharge = 0.; + int pdg = 0; + if (const G4Track* track = currentTrack()) { + pdg = track->GetDefinition()->GetPDGEncoding(); + if (isMonopolePDG(pdg)) { + signedMagneticCharge = (pdg > 0) ? mMagneticChargeEplus : -mMagneticChargeEplus; + } + } + mMagCharge = CLHEP::eplus * signedMagneticCharge * CLHEP::c_light; + } + + void EvaluateRhsGivenB(const G4double y[], const G4double B[3], G4double dydx[]) const override + { + const G4double pSquared = y[3] * y[3] + y[4] * y[4] + y[5] * y[5]; + if (pSquared <= 0.) { + for (int i = 0; i < 8; ++i) { + dydx[i] = 0.; + } + return; + } + const G4double energy = std::sqrt(pSquared + mMassCof); + const G4double pModuleInverse = 1.0 / std::sqrt(pSquared); + const G4double cofEl = mElCharge * pModuleInverse; + const G4double cofMag = mMagCharge * energy * pModuleInverse; + + dydx[0] = y[3] * pModuleInverse; + dydx[1] = y[4] * pModuleInverse; + dydx[2] = y[5] * pModuleInverse; + + // magnetic charge -> force along B; electric charge -> the usual v x B + dydx[3] = cofMag * B[0] + cofEl * (y[4] * B[2] - y[5] * B[1]); + dydx[4] = cofMag * B[1] + cofEl * (y[5] * B[0] - y[3] * B[2]); + dydx[5] = cofMag * B[2] + cofEl * (y[3] * B[1] - y[4] * B[0]); + + // Coupling of the magnetic charge to an electric field: the full dual + // Lorentz force is F = g*(B - v x E/c^2). + // mMagCharge is zero for every non-monopole, so nothing else is affected. + if (mMagCharge != 0. && mTPCDriftField > 0.) { + G4double E[3] = {0., 0., 0.}; + if (tpcDriftField(y, mTPCDriftField, E)) { + // Relative to cofMag this carries 1/(c*energy), so the term is of order + // beta*E/(c*B) compared with the g*B term (similar to G4MonopoleEq, which uses + // d(p)/ds = g*(c*energy*B - p x E)/(p*c)). + const G4double cofMagE = mMagCharge * pModuleInverse / CLHEP::c_light; + const G4double dEx = cofMagE * (y[4] * E[2] - y[5] * E[1]); + const G4double dEy = cofMagE * (y[5] * E[0] - y[3] * E[2]); + const G4double dEz = cofMagE * (y[3] * E[1] - y[4] * E[0]); + dydx[3] -= dEx; + dydx[4] -= dEy; + dydx[5] -= dEz; + } + } + + dydx[6] = 0.; // not used + dydx[7] = energy * pModuleInverse / CLHEP::c_light; // inverse velocity + } + + private: + static const G4Track* currentTrack() + { + auto* eventManager = G4EventManager::GetEventManager(); + if (eventManager == nullptr) { + return nullptr; + } + auto* trackingManager = eventManager->GetTrackingManager(); + return trackingManager != nullptr ? trackingManager->GetTrack() : nullptr; + } + + double mMagneticChargeEplus; ///< |g| in eplus units (1 g_D = 1/(2*alpha) ~ 68.5) + G4double mTPCDriftField = 0.; ///< TPC drift field, Geant4 units; 0 disables the coupling + G4double mElCharge = 0.; + G4double mMagCharge = 0.; + G4double mMassCof = 0.; +}; + +//____________________________________________________________________________ +/// Make the global field integrate O2MonopoleEquation instead of the default +/// electric-charge-only equation. +/// +/// SetUserEquationOfMotion() in Geant4-VMC is not usable here: it +/// registers the object with TG4GeometryManager, and the field integrator is +/// built before that registration is consulted, so the equation is never +/// actually called. +/// +/// \param magneticChargeEplusUnits monopole magnetic charge in eplus units +/// \param tpcDriftFieldGeant4 TPC drift field in Geant4 units (0 = disabled) +inline void installMonopoleFieldIntegrator(double magneticChargeEplusUnits, double tpcDriftFieldGeant4) +{ + auto* transportationManager = G4TransportationManager::GetTransportationManager(); + auto* fieldManager = transportationManager != nullptr ? transportationManager->GetFieldManager() : nullptr; + if (fieldManager == nullptr) { + LOG(error) << "O2MonopolePhysics: no G4FieldManager, monopole equation of motion NOT installed"; + return; + } + auto* magneticField = + const_cast(dynamic_cast(fieldManager->GetDetectorField())); + if (magneticField == nullptr) { + LOG(error) << "O2MonopolePhysics: no magnetic field attached to the field manager, " + "monopole equation of motion NOT installed"; + return; + } + + auto* equation = new O2MonopoleEquation(magneticField, magneticChargeEplusUnits, tpcDriftFieldGeant4); + // A generic (non-helix) integrator is required: the helix steppers hard-code + // the constant-curvature motion of an electric charge. 8 variables so that the + // time-of-flight component is integrated too. + auto* stepper = new G4ClassicalRK4(equation, 8); + + // keep the accuracy Geant4-VMC configured for this field + auto* previous = fieldManager->GetChordFinder(); + const G4double stepMinimum = 1.0e-2 * CLHEP::mm; + auto* chordFinder = new G4ChordFinder(magneticField, stepMinimum, stepper); + if (previous != nullptr) { + chordFinder->SetDeltaChord(previous->GetDeltaChord()); + } + fieldManager->SetChordFinder(chordFinder); + + if (tpcDriftFieldGeant4 > 0.) { + LOG(info) << "O2MonopolePhysics: monopole equation of motion installed (F = g*(B - v x E/c^2)), " + "magnetic charge = " + << magneticChargeEplusUnits << " eplus, TPC drift field = " + << tpcDriftFieldGeant4 / (CLHEP::volt / CLHEP::cm) << " V/cm"; + } else { + LOG(info) << "O2MonopolePhysics: monopole equation of motion installed (F = g*B), magnetic charge = " + << magneticChargeEplusUnits << " eplus (TPC drift field coupling disabled)"; + } +} + //____________________________________________________________________________ /// Minimal physics list whose only job is to attach the magnetic-monopole /// ionisation process to the already-defined O2 monopole particles. @@ -126,6 +350,14 @@ class O2MonopolePhysics : public G4VUserPhysicsList // 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); + + // 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 + // neither, so the field would never be queried. Hence a gate magnetic moment is set + // to allow the query to happen. + if (particle->GetPDGMagneticMoment() == 0.) { + particle->SetPDGMagneticMoment(gMonopoleFieldGateMoment); + } ++nAttached; LOG(info) << "O2MonopolePhysics: attached G4mplIonisation to " << particle->GetParticleName() << " (PDG " << pdg @@ -134,6 +366,17 @@ class O2MonopolePhysics : public G4VUserPhysicsList if (nAttached == 0) { LOG(warning) << "O2MonopolePhysics: no monopole particle found; no ionisation attached"; } + + // This static switch is what makes G4Transportation consider the μ of the monopole + // It is global, so electrically neutral particles that already carry a momentum (neutrons) are + // now propagated through the field as well; O2MonopoleEquation gives them + // exactly zero force, so their trajectories are unchanged. + G4Transportation::EnableMagneticMoment(true); + + // 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. + installMonopoleFieldIntegrator(mMagneticCharge / CLHEP::eplus, tpcDriftFieldMagnitude()); } private: From 399d3ac86555ca5faa79734279300fb29fa1991c Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Wed, 2 Sep 2026 10:36:00 +0200 Subject: [PATCH 3/6] Add ITS monopoles response + data format for monopoles --- .../SimulationDataFormat/MonopoleParticles.h | 41 +++++++++++++++++++ .../ITSMFT/ITS/simulation/src/Detector.cxx | 13 +++++- Detectors/TPC/simulation/src/Detector.cxx | 3 +- Detectors/gconfig/src/O2MonopolePhysics.cxx | 10 ++--- 4 files changed, 59 insertions(+), 8 deletions(-) create mode 100644 DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h new file mode 100644 index 0000000000000..4f618e933e11b --- /dev/null +++ b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h @@ -0,0 +1,41 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file MonopoleParticles.h +/// \brief Identification of the BSM magnetic monopoles across the simulation +/// @author M. Giacalone - September 2026 +/// +/// The monopole PDG codes are defined in O2MCApplication::AddParticles() and in +/// O2DatabasePDG. They are collected here for better syncing inside the simulation code +/// +/// To do: add dyons and other monopoles configurations when needed + +#ifndef ALICEO2_SIMULATION_MONOPOLEPARTICLES_H_ +#define ALICEO2_SIMULATION_MONOPOLEPARTICLES_H_ + +namespace o2::sim +{ + +/// monopole carrying equal electric and magnetic charge +constexpr int MonopolePdgSymm = 4110000; +/// monopole carrying opposite electric and magnetic charge +constexpr int MonopolePdgAsymm = 4120000; + +/// true for the monopole and antimonopole of both species +constexpr bool isMonopole(int pdg) noexcept +{ + const int abspdg = pdg < 0 ? -pdg : pdg; + return abspdg == MonopolePdgSymm || abspdg == MonopolePdgAsymm; +} + +} // namespace o2::sim + +#endif diff --git a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx index 63d7a8ad8dfa2..b965ff50219e4 100644 --- a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx +++ b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx @@ -24,6 +24,7 @@ #include "DetectorsBase/Stack.h" #include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include "fairlogger/Logger.h" // for LOG, LOG_IF // FairRoot includes @@ -318,7 +319,11 @@ void Detector::InitializeO2Detector() Bool_t Detector::ProcessHits(FairVolume* vol) { // This method is called from the MC stepping - if (!(fMC->TrackCharge())) { + // 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 + const bool isMonopole = o2::sim::isMonopole(fMC->TrackPid()); + if (!(fMC->TrackCharge()) && !isMonopole) { return kFALSE; } @@ -386,6 +391,12 @@ 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 + if (isMonopole && mTrackData.mEnergyLoss <= 0.) { + return kFALSE; + } TLorentzVector positionStop; fMC->TrackPosition(positionStop); // Retrieve the indices with the volume path diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 47fe376074909..800847f450a63 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -17,6 +17,7 @@ #include "DetectorsBase/Stack.h" #include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include "FairVolume.h" // for FairVolume @@ -115,7 +116,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // Magnetic monopoles have zero electric charge but ionise the gas through // their magnetic charge energy loss (G4mplIonisation). const int trackPdg = fMC->TrackPid(); - const bool isMonopole = (TMath::Abs(trackPdg) == 4110000 || TMath::Abs(trackPdg) == 4120000); + const bool isMonopole = o2::sim::isMonopole(trackPdg); if (static_cast(trackCharge) == 0) { // Fall through only for monopoles when ionisation is enabled. // The behaviour for the other neutral particles remains as before. diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index d0815b854d22d..a61aee58b8cdb 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -30,6 +30,7 @@ #include "SimSetup/O2MonopolePhysics.h" #include "CommonUtils/ConfigurableParam.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include // needed to instantiate getValueAs<> @@ -79,13 +80,10 @@ namespace // O2MCApplication::AddParticles() and O2DatabasePDG: // +-4110000 : "symmetric" monopoles // +-4120000 : "asymmetric" monopoles -constexpr std::array gMonopolePDGs = {4110000, -4110000, 4120000, -4120000}; +constexpr std::array gMonopolePDGs = {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm, + o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm}; -inline bool isMonopolePDG(int pdg) -{ - const int abspdg = std::abs(pdg); - return abspdg == 4110000 || abspdg == 4120000; -} +inline bool isMonopolePDG(int pdg) { return o2::sim::isMonopole(pdg); } /// GEANT4 only queries the magnetic field when a particle has non-zero electric /// charge or non-zero magnetic moment (μ) and the monopole has neither. So this is a workaround From 3ea9362a80a78cd4c01eb8a30efd7375b1a29f20 Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Fri, 11 Sep 2026 10:32:26 +0200 Subject: [PATCH 4/6] Combine monopole + fastsim hook --- Detectors/gconfig/src/O2MonopolePhysics.cxx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index a61aee58b8cdb..70ad9812ff301 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -40,6 +40,7 @@ #include "TG4RunConfiguration.h" #include "TG4ComposedPhysicsList.h" +#include "FastSim/G4FastSimulation.h" #include #include @@ -384,20 +385,22 @@ class O2MonopolePhysics : public G4VUserPhysicsList //____________________________________________________________________________ /// TG4RunConfiguration that appends O2MonopolePhysics to the composed physics /// list which VMC builds for the requested reference list. -class O2G4RunConfiguration : public TG4RunConfiguration +/// Derives from o2::fastsim::G4RunConfiguration (rather than TG4RunConfiguration +/// directly) so that monopole ionisation and the G4 fast-simulation can be combined +class O2G4RunConfiguration : public o2::fastsim::G4RunConfiguration { public: O2G4RunConfiguration(const TString& userGeometry, const TString& physicsList, const TString& specialProcess, Bool_t specialStacking, Bool_t mtApplication, double magneticChargeEplus) - : TG4RunConfiguration(userGeometry, physicsList, specialProcess, specialStacking, mtApplication), + : o2::fastsim::G4RunConfiguration(userGeometry, physicsList, specialProcess, specialStacking, mtApplication), mMagneticCharge(magneticChargeEplus) { } G4VUserPhysicsList* CreatePhysicsList() override { - G4VUserPhysicsList* physicsList = TG4RunConfiguration::CreatePhysicsList(); + G4VUserPhysicsList* physicsList = o2::fastsim::G4RunConfiguration::CreatePhysicsList(); if (auto* composed = dynamic_cast(physicsList)) { composed->AddPhysicsList(new O2MonopolePhysics(mMagneticCharge)); LOG(info) << "O2G4RunConfiguration: monopole ionisation physics registered " From 1d62ef1e42684746f2d8bff3a480988931840faa Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Tue, 22 Sep 2026 15:08:36 +0200 Subject: [PATCH 5/6] Implemented Sandro comments + give monopoles own transport process Previous version was configuring O2MonopolePhysics globally which replaced the chord finder on the run G4FieldManager and also called G4Transportation::EnableMa gneticMoment(true), which are botbh global switched. Since NystromRK4 is incompatible with the monopole physics, the G4ClassicalRK4 was used, which I think is wrong, as I wanted to keep the original stepper the VMC macro asks for. Now the monopoles chord finder is built and not installed. In addition some minor things that were raised by Sandro's review and found while rechecking the code: - Now SetMaxstep(1.e10) is not called for monopoles, as that limit belongs to the dri ft volume rather than the track - TPC and ITS ProcessHits look up the monopole PDG code for neutral tracks only (reduction in calls) - Throw fatal in case Monopole Physics is requested but not installed properly The commit, together with the rest has been validated by running with Schwinger pair monopoles + various gD options. --- .../ITSMFT/ITS/simulation/src/Detector.cxx | 9 +- Detectors/TPC/simulation/src/Detector.cxx | 17 +- Detectors/gconfig/src/O2MonopolePhysics.cxx | 162 +++++++++++++++--- 3 files changed, 154 insertions(+), 34 deletions(-) diff --git a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx index 86dee749fd807..e5fffe20f6fa2 100644 --- a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx +++ b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx @@ -321,9 +321,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 - const bool isMonopole = o2::sim::isMonopole(fMC->TrackPid()); - if (!(fMC->TrackCharge()) && !isMonopole) { + // rejected by the electric-charge gate. PDG lookup + // never runs for ordinary charged production. + // To-do: handle dyons. + const bool isNeutral = (fMC->TrackCharge() == 0); + const bool isMonopole = isNeutral && o2::sim::isMonopole(fMC->TrackPid()); + if (isNeutral && !isMonopole) { return kFALSE; } diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 1ae560f26f3b7..8dc577a3ef509 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -122,17 +122,22 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // LOG(info) << "tpc::ProcessHits"; const double trackCharge = fMC->TrackCharge(); // Magnetic monopoles have zero electric charge but ionise the gas through - // their magnetic charge energy loss (G4mplIonisation). - const int trackPdg = fMC->TrackPid(); - const bool isMonopole = o2::sim::isMonopole(trackPdg); + // their magnetic charge energy loss (G4mplIonisation), so they must not be + // rejected by the electric-charge gate. The PDG lookup is only worth doing for neutral particle + // To-do: add dyons case + bool isMonopole = false; if (static_cast(trackCharge) == 0) { - // Fall through only for monopoles when ionisation is enabled. - // The behaviour for the other neutral particles remains as before. - if (!isMonopole || fMC->Edep() <= 0.) { + isMonopole = 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. + return kFALSE; + } } // ===| SET THE LENGTH OF THE NEXT ENERGY LOSS STEP |========================= diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 70ad9812ff301..6d220e7825f83 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -60,6 +60,7 @@ #include #include #include +#include #include #include @@ -234,8 +235,20 @@ class O2MonopoleEquation : public G4EquationOfMotion }; //____________________________________________________________________________ -/// Make the global field integrate O2MonopoleEquation instead of the default -/// electric-charge-only equation. +/// The field manager of the run, together with a second chord finder that +/// integrates O2MonopoleEquation +struct MonopoleFieldSetup { + G4FieldManager* fieldManager = nullptr; + G4ChordFinder* chordFinder = nullptr; +}; + +//____________________________________________________________________________ +/// Build a chord finder that integrates O2MonopoleEquation instead of the +/// default electric-charge-only equation. +/// +/// It is deliberately NOT installed on the field manager here. The default +/// chord finder is what Geant4-VMC configured for this run (NystromRK4 unless +/// the macro says otherwise) and every ordinary particle keeps it /// /// SetUserEquationOfMotion() in Geant4-VMC is not usable here: it /// registers the object with TG4GeometryManager, and the field integrator is @@ -244,20 +257,22 @@ class O2MonopoleEquation : public G4EquationOfMotion /// /// \param magneticChargeEplusUnits monopole magnetic charge in eplus units /// \param tpcDriftFieldGeant4 TPC drift field in Geant4 units (0 = disabled) -inline void installMonopoleFieldIntegrator(double magneticChargeEplusUnits, double tpcDriftFieldGeant4) +inline MonopoleFieldSetup buildMonopoleFieldSetup(double magneticChargeEplusUnits, double tpcDriftFieldGeant4) { auto* transportationManager = G4TransportationManager::GetTransportationManager(); auto* fieldManager = transportationManager != nullptr ? transportationManager->GetFieldManager() : nullptr; if (fieldManager == nullptr) { - LOG(error) << "O2MonopolePhysics: no G4FieldManager, monopole equation of motion NOT installed"; - return; + LOG(fatal) << "O2MonopolePhysics: no G4FieldManager, the monopole equation of motion " + "cannot be installed; rerun with G4.monopole=0 if that is what you want"; + return {}; } auto* magneticField = const_cast(dynamic_cast(fieldManager->GetDetectorField())); if (magneticField == nullptr) { - LOG(error) << "O2MonopolePhysics: no magnetic field attached to the field manager, " - "monopole equation of motion NOT installed"; - return; + LOG(fatal) << "O2MonopolePhysics: no magnetic field attached to the field manager, the " + "monopole equation of motion cannot be installed; rerun with G4.monopole=0 " + "if that is what you want"; + return {}; } auto* equation = new O2MonopoleEquation(magneticField, magneticChargeEplusUnits, tpcDriftFieldGeant4); @@ -273,17 +288,111 @@ inline void installMonopoleFieldIntegrator(double magneticChargeEplusUnits, doub if (previous != nullptr) { chordFinder->SetDeltaChord(previous->GetDeltaChord()); } - fieldManager->SetChordFinder(chordFinder); if (tpcDriftFieldGeant4 > 0.) { - LOG(info) << "O2MonopolePhysics: monopole equation of motion installed (F = g*(B - v x E/c^2)), " + LOG(info) << "O2MonopolePhysics: monopole equation of motion built (F = g*(B - v x E/c^2)), " "magnetic charge = " << magneticChargeEplusUnits << " eplus, TPC drift field = " << tpcDriftFieldGeant4 / (CLHEP::volt / CLHEP::cm) << " V/cm"; } else { - LOG(info) << "O2MonopolePhysics: monopole equation of motion installed (F = g*B), magnetic charge = " + LOG(info) << "O2MonopolePhysics: monopole equation of motion built (F = g*B), magnetic charge = " << magneticChargeEplusUnits << " eplus (TPC drift field coupling disabled)"; } + return {fieldManager, chordFinder}; +} + +//____________________________________________________________________________ +/// Transportation for the monopole species only. +/// +/// These must be true for the monopoles and false for every other particle: +/// +/// - the field manager has to use the monopole chord finder, otherwise O2MonopoleEquation +/// might not be evaluated at all (such as it happens with the default NystromRK4) +/// - G4Transportation has to consider the magnetic moment +/// +/// Ordinary tracks keep the stepper the run was configured with, and neutral particles that +/// happen to carry a magnetic moment (neutrons above all) keep their +/// straight-line transport. +/// +/// This mirrors G4MonopoleTransportation from the Geant4 monopole example, but +/// by derives directly from G4Transportation, so that it follows the Geant4 versions +class O2MonopoleTransportation : public G4Transportation +{ + public: + O2MonopoleTransportation(G4FieldManager* fieldManager, G4ChordFinder* monopoleChordFinder) + : G4Transportation(0), mFieldManager(fieldManager), mMonopoleChordFinder(monopoleChordFinder) + { + } + + G4double AlongStepGetPhysicalInteractionLength(const G4Track& track, G4double previousStepSize, + G4double currentMinimumStep, G4double& currentSafety, + G4GPILSelection* selection) override + { + const G4bool previousMoment = G4Transportation::EnableMagneticMoment(true); + auto* previousChordFinder = mFieldManager->GetChordFinder(); + mFieldManager->SetChordFinder(mMonopoleChordFinder); + + const G4double length = G4Transportation::AlongStepGetPhysicalInteractionLength( + track, previousStepSize, currentMinimumStep, currentSafety, selection); + + mFieldManager->SetChordFinder(previousChordFinder); + G4Transportation::EnableMagneticMoment(previousMoment); + return length; + } + + private: + G4FieldManager* mFieldManager; ///< field manager of the run, not owned + G4ChordFinder* mMonopoleChordFinder; ///< installed only for the duration of a monopole step, not owned +}; + +//____________________________________________________________________________ +/// Replace the transport of one monopole species with +/// O2MonopoleTransportation, keeping it first in the DoIt vectors exactly as +/// G4VUserPhysicsList::AddTransportation() left it. +inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4ParticleDefinition* particle, + const MonopoleFieldSetup& fieldSetup) +{ + G4VProcess* existing = nullptr; + G4ProcessVector* plist = pmanager->GetProcessList(); + for (G4int ip = 0; ip < static_cast(plist->size()); ++ip) { + if (dynamic_cast((*plist)[ip]) != nullptr) { + existing = (*plist)[ip]; + break; + } + } + if (existing == nullptr) { + LOG(fatal) << "O2MonopolePhysics: " << particle->GetParticleName() + << " has no transportation process to replace; the monopole could not be " + "coupled to the field"; + return false; + } + if (dynamic_cast(existing) != nullptr) { + // Parallel worlds are in use. O2MonopoleTransportation derives from plain + // G4Transportation, so swapping it in would drop the parallel-world + // navigation; refuse rather than silently mis-navigate. + LOG(fatal) << "O2MonopolePhysics: " << particle->GetParticleName() + << " uses G4CoupledTransportation (parallel worlds); the monopole transportation " + "does not support that. Rerun with G4.monopole=0 or without parallel worlds"; + return false; + } + + // One G4Transportation instance is shared by every particle + // (G4VUserPhysicsList::AddTransportation creates a single one), so it is + // detached from this particle only and must not be deleted. + pmanager->RemoveProcess(existing); + + auto* transportation = new O2MonopoleTransportation(fieldSetup.fieldManager, fieldSetup.chordFinder); + pmanager->AddProcess(transportation); + // Transportation has to be first in the DoIt vectors + // + // Geant4 prints "Set Ordering First is invoked twice for Transportation to + // " (ProcMan113, JustWarning) once per call here, because + // G4ProcessManager latches isSetOrderingFirstInvoked and RemoveProcess() does + // not clear it. The insertion is performed before that check and is correct; + // the warning in the stdout is expected and harmless. + pmanager->SetProcessOrderingToFirst(transportation, idxAlongStep); + pmanager->SetProcessOrderingToFirst(transportation, idxPostStep); + return true; } //____________________________________________________________________________ @@ -310,6 +419,13 @@ class O2MonopolePhysics : public G4VUserPhysicsList void ConstructProcess() override { 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()); + int nAttached = 0; for (int pdg : gMonopolePDGs) { auto* particle = table->FindParticle(pdg); @@ -357,25 +473,20 @@ class O2MonopolePhysics : public G4VUserPhysicsList if (particle->GetPDGMagneticMoment() == 0.) { particle->SetPDGMagneticMoment(gMonopoleFieldGateMoment); } + + // 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. + installMonopoleTransport(pmanager, particle, fieldSetup); + ++nAttached; - LOG(info) << "O2MonopolePhysics: attached G4mplIonisation to " + LOG(info) << "O2MonopolePhysics: attached G4mplIonisation and monopole transport to " << particle->GetParticleName() << " (PDG " << pdg << "), magnetic charge = " << mMagneticCharge / CLHEP::eplus << " eplus"; } if (nAttached == 0) { LOG(warning) << "O2MonopolePhysics: no monopole particle found; no ionisation attached"; } - - // This static switch is what makes G4Transportation consider the μ of the monopole - // It is global, so electrically neutral particles that already carry a momentum (neutrons) are - // now propagated through the field as well; O2MonopoleEquation gives them - // exactly zero force, so their trajectories are unchanged. - G4Transportation::EnableMagneticMoment(true); - - // 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. - installMonopoleFieldIntegrator(mMagneticCharge / CLHEP::eplus, tpcDriftFieldMagnitude()); } private: @@ -406,8 +517,9 @@ class O2G4RunConfiguration : public o2::fastsim::G4RunConfiguration LOG(info) << "O2G4RunConfiguration: monopole ionisation physics registered " "on the composed physics list"; } else { - LOG(error) << "O2G4RunConfiguration: physics list is not a TG4ComposedPhysicsList, " - "monopole ionisation could NOT be enabled"; + LOG(fatal) << "O2G4RunConfiguration: physics list is not a TG4ComposedPhysicsList, " + "monopole ionisation cannot be enabled; rerun with G4.monopole=0 if that " + "is what you want"; } return physicsList; } From 26f3cee935d09843131d091bd9bca402d85797dd Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Tue, 22 Sep 2026 18:38:15 +0200 Subject: [PATCH 6/6] Make monopoles mass configurable via parameter --- Common/SimConfig/include/SimConfig/G4Params.h | 4 ++++ .../SimulationDataFormat/MonopoleParticles.h | 3 +++ .../SimulationDataFormat/O2DatabasePDG.h | 19 +++++++++++++------ Detectors/TPC/simulation/src/Detector.cxx | 8 +++++--- Detectors/gconfig/src/O2MonopolePhysics.cxx | 13 ++++++++----- Generators/src/GeneratorFactory.cxx | 5 ++++- Steer/src/O2MCApplication.cxx | 12 ++++++++---- macro/o2sim.C | 7 ++++++- 8 files changed, 51 insertions(+), 20 deletions(-) diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 54dcf861f3251..97b31b3e3c5e3 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/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 { // 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(o2::sim::MonopoleMassDefaultGeV); // Fast simulation. Empty fastSimModels (the default) disables the feature // entirely; see Detectors/gconfig/include/SimSetup/G4FastSimulation.h. diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h index 4f618e933e11b..93b9f1f67f915 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h +++ b/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 { diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index a24a36134becb..2c9b966bd098a 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -19,6 +19,7 @@ #include #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 // diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 8dc577a3ef509..8359e820a03a3 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/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(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; } } diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 6d220e7825f83..86c5a4cd36f83 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/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; diff --git a/Generators/src/GeneratorFactory.cxx b/Generators/src/GeneratorFactory.cxx index f62e8dfddb849..3c3644bb4c3bb 100644 --- a/Generators/src/GeneratorFactory.cxx +++ b/Generators/src/GeneratorFactory.cxx @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -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)) { diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index c8c339af506ce..c805e0f1011d6 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/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() diff --git a/macro/o2sim.C b/macro/o2sim.C index fd44de6e93b54..b88251705873b 100644 --- a/macro/o2sim.C +++ b/macro/o2sim.C @@ -15,6 +15,7 @@ #include #include #include "SimulationDataFormat/O2DatabasePDG.h" +#include "SimConfig/G4Params.h" #include "SimulationDataFormat/MCEventHeader.h" #include #include @@ -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; {