From 0fb2da00ee7e7b09db981ed308a60f19900e3086 Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Wed, 15 Jul 2026 16:11:37 +0200 Subject: [PATCH] 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()