diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 2a333a39e4242..b84afc9c33b50 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -54,6 +54,20 @@ 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. + // 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/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 8c9c45b6c13a8..47fe376074909 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 |========================= @@ -191,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(); @@ -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, gasParam.MaxElePerStep); // 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 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; + 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, gasParam.MaxElePerStep); // 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..d0815b854d22d --- /dev/null +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -0,0 +1,438 @@ +// 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 "CommonUtils/ConfigurableParam.h" + +#include // needed to instantiate getValueAs<> + +#include + +#include + +#include "TG4RunConfiguration.h" +#include "TG4ComposedPhysicsList.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#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}; + +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. +/// +/// 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); + + // 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 + << "), 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: + 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()