Integration of monopoles in O2 framework

Common/SimConfig/include/SimConfig/G4Params.h

@@ -14,6 +14,7 @@
 
 #include "CommonUtils/ConfigurableParam.h"
 #include "CommonUtils/ConfigurableParamHelper.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
 
 namespace o2
 {
@@ -55,6 +56,22 @@ struct G4Params : public o2::conf::ConfigurableParamHelper<G4Params> {
   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;
+  // 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<float>(o2::sim::MonopoleMassDefaultGeV);
+
   // Fast simulation. Empty fastSimModels (the default) disables the feature
   // entirely; see Detectors/gconfig/include/SimSetup/G4FastSimulation.h.
   std::string fastSimModels = "";   // comma-separated model names to activate

DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h

@@ -0,0 +1,44 @@
+// 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;
+
+/// Default monopole mass in GeV.
+constexpr double MonopoleMassDefaultGeV = 100.;
+
+/// 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

DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h

@@ -19,6 +19,7 @@
 #include <string>
 #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
@@ -189,6 +196,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", 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", monopoleMass, kTRUE, 0.0, 0, "BSM", 4120000);
+  db->AddParticle("AntiMonopole_asymm", "AntiMonopole_asymm", monopoleMass, kTRUE, 0.0, 0, "BSM", -4120000);
+
   // IONS
   //
   // Done by default now from Pythia6 table

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,14 @@ 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. 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;
   }
 
@@ -386,6 +394,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

Detectors/TPC/simulation/src/Detector.cxx

@@ -17,6 +17,8 @@
 
 #include "DetectorsBase/Stack.h"
 #include "SimulationDataFormat/TrackReference.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
+#include "SimConfig/G4Params.h"
 
 #include "FairVolume.h" // for FairVolume
 
@@ -120,11 +122,24 @@ 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), 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<int>(trackCharge) == 0) {
-
-    // set a very large step size for neutral particles
-    fMC->SetMaxStep(1.e10);
-    return kFALSE; // take only charged particles
+    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.) {
+      // An ionising monopole depositing nothing on this step: there
+      // is no hit to create, and the step limit is deliberately left alone
+      return kFALSE;
+    }
   }
 
   // ===| SET THE LENGTH OF THE NEXT ENERGY LOSS STEP |=========================
@@ -199,7 +214,19 @@ 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();
+
+  double betaGamma = momentum.P() / fMC->TrackMass();
+  betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg
+
+  // ---| number of primary ionisations per cm |---
+  const double primaryElectronsPerCM =
+    gasParam.Nprim * BetheBlochAleph(static_cast<float>(betaGamma), gasParam.BetheBlochParam[0],
+                                     gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2],
+                                     gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]);
 
   // use Geant4 energy deposit directly for ionisation (Kr-83m calibration simulations)
   if (detParam.UseGeant4Edep) {
@@ -210,19 +237,17 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
     if (meanIon > 0.) {
       numberOfElectrons = static_cast<int>(gasParam.FanoFactorG4 * Gamma(meanIon / gasParam.FanoFactorG4));
     }
+  } else 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<int>(fMC->Edep() / static_cast<double>(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 {
-    // ---| Stepsize in cm |---
-    const double stepSize = fMC->TrackStep();
-
-    double betaGamma = momentum.P() / fMC->TrackMass();
-    betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg
-
-    // ---| number of primary ionisations per cm |---
-    const double primaryElectronsPerCM =
-      gasParam.Nprim * BetheBlochAleph(static_cast<float>(betaGamma), gasParam.BetheBlochParam[0],
-                                       gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2],
-                                       gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]);
-
     // ---| mean number of collisions and random for this event |---
     const double meanNcoll = stepSize * trackCharge * trackCharge * primaryElectronsPerCM;
     const int nColl = static_cast<int>(fMC->GetRandom()->Poisson(meanNcoll));

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 O2::FastSim
 )
 

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"
 #include "SimConfig/G4ScoringMerger.h"
 #include "G4ScoringManager.hh"
 #include "G4VScoringMesh.hh"
@@ -119,12 +120,22 @@ void Config()
     LOG(fatal) << "Unsupported geometry navigation mode";
   }
 
-  // o2::fastsim::G4RunConfiguration differs from TG4RunConfiguration only in
-  // providing the fast-simulation hook; with G4.fastSimModels empty it behaves
-  // identically.
-  auto runConfiguration = new o2::fastsim::G4RunConfiguration(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 {
+    // o2::fastsim::G4RunConfiguration differs from TG4RunConfiguration only in
+    // providing the fast-simulation hook; with G4.fastSimModels empty it behaves
+    // identically.
+    runConfiguration = new o2::fastsim::G4RunConfiguration(geomNavStr, physicsSetup,
+                                                           "stepLimiter+specialCuts",
+                                                           specialStacking, mtMode);
+  }
   if (g4Params.g4scoring) {
     runConfiguration->SetUseOfG4Scoring();
     if (g4Params.g4fluenceweight) {

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 <TString.h>
+
+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_