Skip to content

Commit adf79c9

Browse files
authored
Integration of monopoles in O2 framework (#15602)
This commit adds the transport and detector response of magnetic monopoles to o2-sim when running with Geant4. Related to https://its.cern.ch/jira/browse/O2-6780 Until now, a monopole could neither be defined as a particle nor lose energy in the detectors, since it has no electric charge and the standard Geant4 physics does not ionise for it. - Four new particles, `Monopole_symm`/`Monopole_asymm` and their antiparticles (PDG codes ±4110000 and ±4120000), are registered in `O2DatabasePDG` and the Geant4 particle table. The helpers are in `MonopoleParticles.h`. - `O2MonopolePhysics` provides the monopole ionisation (`mplIoni`) and a dedicated transport process with the full dual Lorentz force. It does not replace the global chord finder or the stepper requested by the VMC macro. - ITS and TPC `ProcessHits` let neutral monopoles deposit energy. ITS gets a monopole response and a data format, TPC uses a Bethe-Bloch based ionisation. - `GeneratorFactory` takes the monopole mass so that generator and transport agree, and the monopole hook is combined with the fast-sim hook. Everything is steered by three new `G4Params` options: - `G4.monopole` (default `false`) switches the monopole physics on. - `G4.monopoleMagneticCharge` is the charge in units of the Dirac charge (default 1). - `G4.monopoleMass` is the mass in GeV (default 100) and has to agree with the mass in `TDatabasePDG`, otherwise the run is fatal. A non-positive magnetic charge or mass, or monopole physics requested but not installed, are also fatal. By default nothing changes for the user: the particles are only defined and the lookups in ITS and TPC are only done when `G4.monopole` is set. The changes were validated with Schwinger pair-produced monopoles and several values of the magnetic charge.
1 parent aa95339 commit adf79c9

12 files changed

Lines changed: 816 additions & 16 deletions

File tree

‎Common/SimConfig/include/SimConfig/G4Params.h‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414

1515
#include "CommonUtils/ConfigurableParam.h"
1616
#include "CommonUtils/ConfigurableParamHelper.h"
17+
#include "SimulationDataFormat/MonopoleParticles.h"
1718

1819
namespace o2
1920
{
@@ -80,6 +81,22 @@ struct G4Params : public o2::conf::ConfigurableParamHelper<G4Params> {
8081
bool g4scoring = false;
8182
bool g4fluenceweight = false;
8283

84+
// Enable magnetic-monopole ionisation as defined in Detectors/gconfig/O2MonopolePhysics.
85+
// The G4mplIonisation process (Ahlen stopping power) is attached to
86+
// the O2 monopole particles (PDG +-4110000 / +-4120000) on top of the chosen
87+
// reference physics list. Off by default
88+
bool monopole = false;
89+
// Monopole magnetic charge in units of the Dirac charge g_D = eplus/(2*alpha)
90+
// (~68.5 eplus). 1.0 corresponds to a single classic Dirac monopole.
91+
// The monopole is transported with the full dual Lorentz force
92+
// F = g*(B - v x E/c^2). The electric part acts inside the TPC field cage,
93+
// the only volume with an electric field, and its magnitude is the TPC's own
94+
// TPCGEMParam.ElectricField[0]
95+
float monopoleMagneticCharge = 1.f;
96+
// Monopole mass in GeV, preliminarly applied to all the monopoles species
97+
// To-do: allow different monopoles masses in the same run
98+
float monopoleMass = static_cast<float>(o2::sim::MonopoleMassDefaultGeV);
99+
83100
// Fast simulation. Empty fastSimModels (the default) disables the feature
84101
// entirely; see Detectors/gconfig/include/SimSetup/G4FastSimulation.h.
85102
std::string fastSimModels = ""; // comma-separated model names to activate
Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
/// \file MonopoleParticles.h
13+
/// \brief Identification of the BSM magnetic monopoles across the simulation
14+
/// @author M. Giacalone - September 2026
15+
///
16+
/// The monopole PDG codes are defined in O2MCApplication::AddParticles() and in
17+
/// O2DatabasePDG. They are collected here for better syncing inside the simulation code
18+
///
19+
/// To do: add dyons and other monopoles configurations when needed
20+
21+
#ifndef ALICEO2_SIMULATION_MONOPOLEPARTICLES_H_
22+
#define ALICEO2_SIMULATION_MONOPOLEPARTICLES_H_
23+
24+
namespace o2::sim
25+
{
26+
27+
/// Both species are electrically neutral for now, and the transport gives them the
28+
/// same magnetic charge, with its sign taken from the PDG sign.
29+
/// monopole intended to carry equal electric and magnetic charge
30+
constexpr int MonopolePdgSymm = 4110000;
31+
/// monopole intended to carry opposite electric and magnetic charge (not yet distinguished)
32+
constexpr int MonopolePdgAsymm = 4120000;
33+
34+
/// Default monopole mass in GeV.
35+
constexpr double MonopoleMassDefaultGeV = 100.;
36+
37+
/// true for the monopole and antimonopole of both species
38+
constexpr bool isMonopole(int pdg) noexcept
39+
{
40+
const int abspdg = pdg < 0 ? -pdg : pdg;
41+
return abspdg == MonopolePdgSymm || abspdg == MonopolePdgAsymm;
42+
}
43+
44+
} // namespace o2::sim
45+
46+
#endif

‎DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h‎

Lines changed: 38 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,9 +16,11 @@
1616
#ifndef O2_O2DATABASEPDG_H
1717
#define O2_O2DATABASEPDG_H
1818

19+
#include <cmath>
1920
#include <string>
2021
#include "TDatabasePDG.h"
2122
#include "TParticlePDG.h"
23+
#include "SimulationDataFormat/MonopoleParticles.h"
2224

2325
namespace o2
2426
{
@@ -45,9 +47,34 @@ class O2DatabasePDG
4547
}
4648

4749
// adds ALICE particles to a given TDatabasePDG instance
48-
static void addALICEParticles(TDatabasePDG* db = TDatabasePDG::Instance());
50+
// monopoleMass is the mass of all monopole species in GeV; it has to match
51+
// G4Params.monopoleMass or the generator and the transport will disagree
52+
// To fix: all quantities in the MCTrack dependent from GetMass() will get the default 100 GeV.
53+
// A possible fix would be to stop deriving the masses from the compiled-in table,
54+
// or a better solution should be found in the future
55+
static void addALICEParticles(TDatabasePDG* db = TDatabasePDG::Instance(),
56+
double monopoleMass = getMonopoleMass());
4957
static void addParticlesFromExternalFile(TDatabasePDG* db);
5058

59+
// mass of all monopole species in GeV used when they get registered; since TDatabasePDG
60+
// keeps the first registration of a PDG code, it has to be set before the first
61+
// Instance() or addALICEParticles() call, which hasMonopoleMass() can verify
62+
static double getMonopoleMass() { return sMonopoleMass; }
63+
static void setMonopoleMass(double monopoleMass) { sMonopoleMass = monopoleMass; }
64+
65+
// true if all monopole species are registered with the given mass in GeV; TDatabasePDG
66+
// keeps the first registration, so an earlier call with another mass makes this false
67+
static bool hasMonopoleMass(TDatabasePDG* db, double monopoleMass)
68+
{
69+
for (int pdg : {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm, o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm}) {
70+
const auto* particle = db->GetParticle(pdg);
71+
if (particle == nullptr || std::abs(particle->Mass() - monopoleMass) > 1.e-6 * monopoleMass) {
72+
return false;
73+
}
74+
}
75+
return true;
76+
}
77+
5178
// get particle's (if any) mass
5279
static Double_t MassImpl(TParticlePDG* particle, bool& success)
5380
{
@@ -79,13 +106,14 @@ class O2DatabasePDG
79106
O2DatabasePDG() = delete;
80107

81108
private:
109+
inline static double sMonopoleMass = o2::sim::MonopoleMassDefaultGeV;
82110
static constexpr int IONBASELOW{1000000000};
83111
static constexpr int IONBASEHIGH{1099999999};
84112
};
85113

86114
// by keeping this inline, we can use it in other parts of the code, for instance Framework or Analysis,
87115
// without needing to link against this library
88-
inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db)
116+
inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMass)
89117
{
90118
//
91119
// Add ALICE particles to the ROOT PDG data base
@@ -190,6 +218,14 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db)
190218
db->AddParticle("CHI2P_B2", " ", 10.269, kFALSE, 0.0, 0, "meson", 100555);
191219
db->AddParticle("UPSLON4S", " ", 10.580, kFALSE, 0.0, 0, "meson", 300553);
192220

221+
// BSM targeted inclusions
222+
// Monopoles with same electric and magnetic charge
223+
db->AddParticle("Monopole_symm", "Monopole_symm", monopoleMass, kTRUE, 0.0, 0, "BSM", 4110000);
224+
db->AddParticle("AntiMonopole_symm", "AntiMonopole_symm", monopoleMass, kTRUE, 0.0, 0, "BSM", -4110000);
225+
// Monopoles with opposite electric and magnetic charge
226+
db->AddParticle("Monopole_asymm", "Monopole_asymm", monopoleMass, kTRUE, 0.0, 0, "BSM", 4120000);
227+
db->AddParticle("AntiMonopole_asymm", "AntiMonopole_asymm", monopoleMass, kTRUE, 0.0, 0, "BSM", -4120000);
228+
193229
// IONS
194230
//
195231
// Done by default now from Pythia6 table

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

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@
2424

2525
#include "DetectorsBase/Stack.h"
2626
#include "SimulationDataFormat/TrackReference.h"
27+
#include "SimulationDataFormat/MonopoleParticles.h"
28+
#include "SimConfig/G4Params.h"
2729
#include "fairlogger/Logger.h" // for LOG, LOG_IF
2830

2931
// FairRoot includes
@@ -318,7 +320,15 @@ void Detector::InitializeO2Detector()
318320
Bool_t Detector::ProcessHits(FairVolume* vol)
319321
{
320322
// This method is called from the MC stepping
321-
if (!(fMC->TrackCharge())) {
323+
// Electrically neutral magnetic monopoles deposit energy in the
324+
// silicon through G4mplIonisation (Ahlen stopping power), so they must not be
325+
// rejected by the electric-charge gate. The PDG lookup only runs for neutral
326+
// particles in runs with monopole physics enabled.
327+
// To-do: handle dyons.
328+
static const bool sMonopoleIonisation = o2::conf::G4Params::Instance().monopole;
329+
const bool isNeutral = (fMC->TrackCharge() == 0);
330+
const bool isMonopole = isNeutral && sMonopoleIonisation && o2::sim::isMonopole(fMC->TrackPid());
331+
if (isNeutral && !isMonopole) {
322332
return kFALSE;
323333
}
324334

@@ -386,6 +396,10 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
386396
mTrackData.mHitStarted = true;
387397
}
388398
if (stopHit) {
399+
// A monopole that deposited nothing in the sensor leaves no hit
400+
if (isMonopole && mTrackData.mEnergyLoss <= 0.) {
401+
return kFALSE;
402+
}
389403
TLorentzVector positionStop;
390404
fMC->TrackPosition(positionStop);
391405
// Retrieve the indices with the volume path

‎Detectors/TPC/simulation/src/Detector.cxx‎

Lines changed: 30 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,8 @@
1717

1818
#include "DetectorsBase/Stack.h"
1919
#include "SimulationDataFormat/TrackReference.h"
20+
#include "SimulationDataFormat/MonopoleParticles.h"
21+
#include "SimConfig/G4Params.h"
2022

2123
#include "FairVolume.h" // for FairVolume
2224

@@ -120,11 +122,24 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
120122
/* This method is called from the MC stepping for the sensitive volume only */
121123
// LOG(info) << "tpc::ProcessHits";
122124
const double trackCharge = fMC->TrackCharge();
125+
// Magnetic monopoles have zero electric charge but ionise the gas through
126+
// their magnetic charge energy loss (G4mplIonisation), so they must not be
127+
// rejected by the electric-charge gate. The PDG lookup is only worth doing for neutral particle
128+
// To-do: add dyons case
129+
bool isMonopole = false;
123130
if (static_cast<int>(trackCharge) == 0) {
124-
125-
// set a very large step size for neutral particles
126-
fMC->SetMaxStep(1.e10);
127-
return kFALSE; // take only charged particles
131+
static const bool sMonopoleIonisation = o2::conf::G4Params::Instance().monopole;
132+
isMonopole = sMonopoleIonisation && o2::sim::isMonopole(fMC->TrackPid());
133+
if (!isMonopole) {
134+
// set a very large step size for neutral particles
135+
fMC->SetMaxStep(1.e10);
136+
return kFALSE; // take only charged particles
137+
}
138+
if (fMC->Edep() <= 0.) {
139+
// An ionising monopole depositing nothing on this step: there
140+
// is no hit to create, and the step limit is deliberately left alone
141+
return kFALSE;
142+
}
128143
}
129144

130145
// ===| SET THE LENGTH OF THE NEXT ENERGY LOSS STEP |=========================
@@ -198,7 +213,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
198213
// TODO: Add discussion about drawback
199214

200215
Int_t numberOfElectrons = 0;
201-
// I.H. - the type expected in addHit is short
216+
// I.H. - the type expected in addHit is float
202217

203218
// use Geant4 energy deposit directly for ionisation (Kr-83m calibration simulations)
204219
if (detParam.UseGeant4Edep) {
@@ -209,6 +224,16 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
209224
if (meanIon > 0.) {
210225
numberOfElectrons = static_cast<int>(gasParam.FanoFactorG4 * Gamma(meanIon / gasParam.FanoFactorG4));
211226
}
227+
} else if (isMonopole) {
228+
// ---| MONOPOLE IONISATION |--------
229+
// A magnetic monopole ionises the gas via G4mplIonisation (Ahlen stopping
230+
// power), which is not described by the electric-charge model (no e-charge)
231+
// Ionisation electrons calculated directly from the energy deposited in this step: Nel = Edep / Wion.
232+
// To check by TPC experts if this is actually the best way...
233+
numberOfElectrons = static_cast<int>(fMC->Edep() / static_cast<double>(gasParam.Wion));
234+
// The number of electrons is stored as a float in the HitGroup: maximum cap at
235+
// 2^24 (16777216) ==> largest integer a IEEE-754 float can represent exactly
236+
numberOfElectrons = TMath::Min(numberOfElectrons, 16777216);
212237
} else {
213238
// ---| Stepsize in cm |---
214239
const double stepSize = fMC->TrackStep();

‎Detectors/gconfig/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ if(TGeo2VecGeom_FOUND)
2828
endif()
2929

3030
set(G4SETUP_SOURCES src/G4Config.cxx src/G4RunConfiguration.cxx src/G4LocalFieldConstruction.cxx
31-
src/VecGeomNavigation.cxx)
31+
src/VecGeomNavigation.cxx src/O2MonopolePhysics.cxx)
3232
if(G4SETUP_WITH_VECGEOM)
3333
list(APPEND G4SETUP_SOURCES src/VecGeomG4Map.cxx src/VecGeomChecks.cxx src/VecGeomG4NavigatorBase.cxx
3434
src/VecGeomG4Navigator.cxx

‎Detectors/gconfig/g4Config.C‎

Lines changed: 16 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,7 @@ R__LOAD_LIBRARY(libgeant4vmc)
6161
#include "TG4RunConfiguration.h"
6262
#include "SimConfig/G4Params.h"
6363
#include "SimConfig/FluenceWeightCalculator.h"
64+
#include "SimSetup/O2MonopolePhysics.h"
6465
#include "SimConfig/G4ScoringMerger.h"
6566
#include "G4ScoringManager.hh"
6667
#include "G4VScoringMesh.hh"
@@ -124,11 +125,21 @@ void Config()
124125
LOG(fatal) << "Unsupported geometry navigation mode";
125126
}
126127

127-
// o2::g4config::G4RunConfiguration adds the fast-simulation hook and the local
128-
// magnetic fields; with neither configured it behaves like TG4RunConfiguration.
129-
auto runConfiguration = new o2::g4config::G4RunConfiguration(geomNavStr, physicsSetup,
130-
"stepLimiter+specialCuts",
131-
specialStacking, mtMode);
128+
TG4RunConfiguration* runConfiguration = nullptr;
129+
if (g4Params.monopole) {
130+
// Reference physics list + magnetic-monopole ionisation attached
131+
// The reference physics is unchanged.
132+
std::cout << "Monopole ionisation physics requested (G4.monopole=1)\n";
133+
runConfiguration = o2::g4config::createMonopoleRunConfiguration(
134+
geomNavStr, physicsSetup, "stepLimiter+specialCuts", specialStacking, mtMode,
135+
g4Params.monopoleMagneticCharge);
136+
} else {
137+
// o2::g4config::G4RunConfiguration adds the fast-simulation hook and the local
138+
// magnetic fields; with neither configured it behaves like TG4RunConfiguration.
139+
runConfiguration = new o2::g4config::G4RunConfiguration(geomNavStr, physicsSetup,
140+
"stepLimiter+specialCuts",
141+
specialStacking, mtMode);
142+
}
132143
if (g4Params.g4scoring) {
133144
runConfiguration->SetUseOfG4Scoring();
134145
if (g4Params.g4fluenceweight) {
Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
/// \author M+Giacalone - July 2026
13+
14+
#ifndef O2_G4CONFIG_O2MONOPOLEPHYSICS_H_
15+
#define O2_G4CONFIG_O2MONOPOLEPHYSICS_H_
16+
17+
#include <TString.h>
18+
19+
class TG4RunConfiguration;
20+
21+
namespace o2
22+
{
23+
namespace g4config
24+
{
25+
26+
/// Build a Geant4-VMC run configuration which, on top of the requested reference physics list,
27+
/// attaches the magnetic-monopole ionisation process (G4mplIonisation) and transport to
28+
/// the monopole particles already defined by O2 (PDG +-4110000 / +-4120000).
29+
///
30+
/// Only used when G4Params.monopole == true;
31+
/// Otherwise o2::g4config::G4RunConfiguration is used (see g4Config.C)
32+
///
33+
/// \param userGeometry VMC geometry-navigation string (as for TG4RunConfiguration)
34+
/// \param physicsList the reference physics-list selection string
35+
/// \param specialProcess the VMC special-process selection string
36+
/// \param specialStacking VMC special-stacking flag
37+
/// \param mtApplication multithreading flag
38+
/// \param magneticChargeDirac magnetic charge of the monopole in units of the
39+
/// Dirac charge g_D = eplus/(2*alpha) (~68.5 eplus);
40+
/// 1.0 (default) reproduces the classic single Dirac monopole
41+
TG4RunConfiguration* createMonopoleRunConfiguration(const TString& userGeometry,
42+
const TString& physicsList,
43+
const TString& specialProcess,
44+
bool specialStacking,
45+
bool mtApplication,
46+
double magneticChargeDirac);
47+
48+
} // namespace g4config
49+
} // namespace o2
50+
51+
#endif // O2_G4CONFIG_O2MONOPOLEPHYSICS_H_

0 commit comments

Comments
 (0)