TPC MC monopoles and extended track references

Common/SimConfig/include/SimConfig/G4Params.h

@@ -14,6 +14,7 @@
 
 #include "CommonUtils/ConfigurableParam.h"
 #include "CommonUtils/ConfigurableParamHelper.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
 
 namespace o2
 {
@@ -80,6 +81,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,46 @@
+// 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
+{
+
+/// Both species are electrically neutral for now, and the transport gives them the
+/// same magnetic charge, with its sign taken from the PDG sign.
+/// monopole intended to carry equal electric and magnetic charge
+constexpr int MonopolePdgSymm = 4110000;
+/// monopole intended to carry opposite electric and magnetic charge (not yet distinguished)
+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

@@ -16,9 +16,11 @@
 #ifndef O2_O2DATABASEPDG_H
 #define O2_O2DATABASEPDG_H
 
+#include <cmath>
 #include <string>
 #include "TDatabasePDG.h"
 #include "TParticlePDG.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
 
 namespace o2
 {
@@ -44,9 +46,28 @@ 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);
 
+  // true if all monopole species are registered with the given mass in GeV; TDatabasePDG
+  // keeps the first registration, so an earlier call with another mass makes this false
+  static bool hasMonopoleMass(TDatabasePDG* db, double monopoleMass)
+  {
+    for (int pdg : {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm, o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm}) {
+      const auto* particle = db->GetParticle(pdg);
+      if (particle == nullptr || std::abs(particle->Mass() - monopoleMass) > 1.e-6 * monopoleMass) {
+        return false;
+      }
+    }
+    return true;
+  }
+
   // get particle's (if any) mass
   static Double_t MassImpl(TParticlePDG* particle, bool& success)
   {
@@ -84,7 +105,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 +210,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/CPV/simulation/include/CPVSimulation/Detector.h

@@ -157,6 +157,7 @@ class Detector : public o2::base::DetImpl<Detector>
 
   // Simulation
   std::vector<Hit>* mHits = nullptr; //! Collection of CPV hits
+  int mTrackRefLastTrackID = -1;     //! last track that received a TrackReference in the current event
 
   template <typename Det>
   friend class o2::base::DetImpl;

Detectors/CPV/simulation/src/Detector.cxx

@@ -30,6 +30,7 @@
 
 #include "DetectorsBase/GeometryManager.h"
 #include "DetectorsBase/Stack.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
 
 #include <boost/algorithm/string/predicate.hpp>
 #include <boost/range/irange.hpp>
@@ -106,6 +107,7 @@ void Detector::Reset()
   if (!o2::utils::ShmManager::Instance().isOperational()) {
     mHits->clear();
   }
+  mTrackRefLastTrackID = -1;
 }
 
 void Detector::Register() { FairRootManager::Instance()->RegisterAny(addNameTo("Hit").data(), mHits, kTRUE); }
@@ -123,6 +125,13 @@ Bool_t Detector::ProcessHits(FairVolume* v)
   // ported from AliRoot 2019
   // ------------------------------------------------------------------------
 
+  if (fMC->IsTrackEntering() && !fMC->IsNewTrack() && (fMC->TrackCharge() != 0 || o2::sim::isMonopole(fMC->TrackPid()))) { // one reference per charged particle (or monopole), first entry from outside
+    auto o2stack = static_cast<o2::data::Stack*>(fMC->GetStack());
+    if (o2stack->GetCurrentTrackNumber() != mTrackRefLastTrackID) {
+      mTrackRefLastTrackID = o2stack->GetCurrentTrackNumber();
+      o2stack->addTrackReference(o2::TrackReference(*fMC, GetDetId()));
+    }
+  }
   // analyze only charged
   if (fMC->TrackCharge() == 0) {
     fMC->SetMaxStep(1.e10);

Detectors/FIT/FDD/simulation/include/FDDSimulation/Detector.h

@@ -71,6 +71,7 @@ class Detector : public o2::base::DetImpl<Detector>
 
   /// Container for data points
   std::vector<o2::fdd::Hit>* mHits = nullptr;
+  int mTrackRefLastTrackID = -1; //! last track that received a TrackReference in the current event
 
   Detector& operator=(const Detector&);
 

Detectors/FIT/FDD/simulation/src/Detector.cxx

@@ -21,6 +21,7 @@
 #include "FDDSimulation/Detector.h"
 
 #include "DetectorsBase/Stack.h"
+#include "SimulationDataFormat/MonopoleParticles.h"
 #include "Field/MagneticField.h"
 
 #include "TVirtualMC.h"
@@ -103,6 +104,13 @@ void Detector::InitializeO2Detector()
 Bool_t Detector::ProcessHits(FairVolume* vol)
 {
   // This is step manager
+  if (fMC->IsTrackEntering() && !fMC->IsNewTrack() && (fMC->TrackCharge() != 0 || o2::sim::isMonopole(fMC->TrackPid()))) { // one reference per charged particle (or monopole), first entry from outside
+    auto o2stack = static_cast<o2::data::Stack*>(fMC->GetStack());
+    if (o2stack->GetCurrentTrackNumber() != mTrackRefLastTrackID) {
+      mTrackRefLastTrackID = o2stack->GetCurrentTrackNumber();
+      o2stack->addTrackReference(o2::TrackReference(*fMC, GetDetId()));
+    }
+  }
   // do not track neutral particles
   if (!(fMC->TrackCharge())) {
     return kFALSE;
@@ -138,7 +146,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
   // FDA 3  =  12-15
 
   Float_t fFDDLightYield(12.8e6);
-  //BC420 yield is 0.64 of antracene which is 20k photons/MeV = 12800/MeV = 12.8e6/GeV
+  // BC420 yield is 0.64 of antracene which is 20k photons/MeV = 12800/MeV = 12.8e6/GeV
 
   Float_t destep_ad = fMC->Edep();
   Int_t nPhotonsInStep_ad = Int_t(fFDDLightYield * destep_ad);
@@ -164,7 +172,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
   if (fMC->IsTrackExiting() || fMC->IsTrackStop() || fMC->IsTrackDisappeared()) {
     Int_t trackID = fMC->GetStack()->GetCurrentTrackNumber();
 
-    Float_t time = fMC->TrackTime() * 1.0e9; //time from seconds to ns
+    Float_t time = fMC->TrackTime() * 1.0e9; // time from seconds to ns
 
     addHit(trackID, ADsector, vPos, time, eloss_ad, nPhotons_ad);
     return kTRUE;
@@ -176,7 +184,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
 //_____________________________________________________________________________
 Hit* Detector::addHit(int trackID, unsigned short detID, const TVector3& Pos, double Time, double eLoss, int nPhot)
 {
-  //LOG(info) << "FDD hit "<<trackID<<" "<<detID<<" "<<Time<<" "<<eLoss<<" "<<nPhot;
+  // LOG(info) << "FDD hit "<<trackID<<" "<<detID<<" "<<Time<<" "<<eLoss<<" "<<nPhot;
   mHits->emplace_back(trackID, detID, Pos, Time, eLoss, nPhot);
   auto stack = (o2::data::Stack*)fMC->GetStack();
   stack->addHit(GetDetId());
@@ -373,4 +381,5 @@ void Detector::Reset()
   if (!o2::utils::ShmManager::Instance().isOperational()) {
     mHits->clear();
   }
+  mTrackRefLastTrackID = -1;
 }

Detectors/FIT/FT0/simulation/include/FT0Simulation/Detector.h

@@ -149,8 +149,8 @@ class Detector : public o2::base::DetImpl<Detector>
   std::vector<Double_t> mRindexFrontWindow;
 
   // Define the aluminium frame for the detector
-  TGeoVolume* constructFrameAGeometry(); //A-side
-  TGeoVolume* constructFrameCGeometry(); //C-side
+  TGeoVolume* constructFrameAGeometry(); // A-side
+  TGeoVolume* constructFrameCGeometry(); // C-side
   std::string cPlateShapeString();
 
   // BEGIN: Support structure constants
@@ -179,6 +179,7 @@ class Detector : public o2::base::DetImpl<Detector>
 
   /// Container for data points
   std::vector<o2::ft0::HitType>* mHits = nullptr;
+  int mTrackRefLastTrackID = -1; //! last track that received a TrackReference in the current event
 
   // Define volume IDs
   int mREGVolID = -1;