From 0fb2da00ee7e7b09db981ed308a60f19900e3086 Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Wed, 15 Jul 2026 16:11:37 +0200 Subject: [PATCH 01/11] 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() From ad1ae639fc2b05846fc07e5f98bab52cf52a9a2c Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Tue, 21 Jul 2026 14:17:55 +0200 Subject: [PATCH 02/11] Include full Lorentz dual force for monopoles --- Common/SimConfig/include/SimConfig/G4Params.h | 4 + Detectors/TPC/simulation/src/Detector.cxx | 10 +- Detectors/gconfig/src/O2MonopolePhysics.cxx | 243 ++++++++++++++++++ 3 files changed, 252 insertions(+), 5 deletions(-) diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index be71960427589..b84afc9c33b50 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -62,6 +62,10 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { 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/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 09ee899e324c8..47fe376074909 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -198,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(); @@ -219,9 +219,9 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // 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); + // 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; @@ -244,7 +244,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol) 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 + nel_step = TMath::Min(nel_step, gasParam.MaxElePerStep); // 300 electrons corresponds to 10 keV numberOfElectrons += nel_step; } } diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 4fba37a067678..d0815b854d22d 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -29,6 +29,11 @@ /// Geant4 libraries (libG4processes) #include "SimSetup/O2MonopolePhysics.h" +#include "CommonUtils/ConfigurableParam.h" + +#include // needed to instantiate getValueAs<> + +#include #include @@ -43,10 +48,24 @@ #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + #include #include #include +#include #include namespace o2 @@ -61,8 +80,213 @@ namespace // +-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. @@ -126,6 +350,14 @@ class O2MonopolePhysics : public G4VUserPhysicsList // 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 @@ -134,6 +366,17 @@ class O2MonopolePhysics : public G4VUserPhysicsList 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: From 399d3ac86555ca5faa79734279300fb29fa1991c Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Wed, 2 Sep 2026 10:36:00 +0200 Subject: [PATCH 03/11] Add ITS monopoles response + data format for monopoles --- .../SimulationDataFormat/MonopoleParticles.h | 41 +++++++++++++++++++ .../ITSMFT/ITS/simulation/src/Detector.cxx | 13 +++++- Detectors/TPC/simulation/src/Detector.cxx | 3 +- Detectors/gconfig/src/O2MonopolePhysics.cxx | 10 ++--- 4 files changed, 59 insertions(+), 8 deletions(-) create mode 100644 DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h new file mode 100644 index 0000000000000..4f618e933e11b --- /dev/null +++ b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h @@ -0,0 +1,41 @@ +// 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; + +/// 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 diff --git a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx index 63d7a8ad8dfa2..b965ff50219e4 100644 --- a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx +++ b/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,11 @@ 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 + const bool isMonopole = o2::sim::isMonopole(fMC->TrackPid()); + if (!(fMC->TrackCharge()) && !isMonopole) { return kFALSE; } @@ -386,6 +391,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 diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 47fe376074909..800847f450a63 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -17,6 +17,7 @@ #include "DetectorsBase/Stack.h" #include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include "FairVolume.h" // for FairVolume @@ -115,7 +116,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // 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); + const bool isMonopole = o2::sim::isMonopole(trackPdg); if (static_cast(trackCharge) == 0) { // Fall through only for monopoles when ionisation is enabled. // The behaviour for the other neutral particles remains as before. diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index d0815b854d22d..a61aee58b8cdb 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -30,6 +30,7 @@ #include "SimSetup/O2MonopolePhysics.h" #include "CommonUtils/ConfigurableParam.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include // needed to instantiate getValueAs<> @@ -79,13 +80,10 @@ namespace // O2MCApplication::AddParticles() and O2DatabasePDG: // +-4110000 : "symmetric" monopoles // +-4120000 : "asymmetric" monopoles -constexpr std::array gMonopolePDGs = {4110000, -4110000, 4120000, -4120000}; +constexpr std::array gMonopolePDGs = {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm, + o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm}; -inline bool isMonopolePDG(int pdg) -{ - const int abspdg = std::abs(pdg); - return abspdg == 4110000 || abspdg == 4120000; -} +inline bool isMonopolePDG(int pdg) { return o2::sim::isMonopole(pdg); } /// 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 From 3ea9362a80a78cd4c01eb8a30efd7375b1a29f20 Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Fri, 11 Sep 2026 10:32:26 +0200 Subject: [PATCH 04/11] Combine monopole + fastsim hook --- Detectors/gconfig/src/O2MonopolePhysics.cxx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index a61aee58b8cdb..70ad9812ff301 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -40,6 +40,7 @@ #include "TG4RunConfiguration.h" #include "TG4ComposedPhysicsList.h" +#include "FastSim/G4FastSimulation.h" #include #include @@ -384,20 +385,22 @@ class O2MonopolePhysics : public G4VUserPhysicsList //____________________________________________________________________________ /// TG4RunConfiguration that appends O2MonopolePhysics to the composed physics /// list which VMC builds for the requested reference list. -class O2G4RunConfiguration : public TG4RunConfiguration +/// Derives from o2::fastsim::G4RunConfiguration (rather than TG4RunConfiguration +/// directly) so that monopole ionisation and the G4 fast-simulation can be combined +class O2G4RunConfiguration : public o2::fastsim::G4RunConfiguration { 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), + : o2::fastsim::G4RunConfiguration(userGeometry, physicsList, specialProcess, specialStacking, mtApplication), mMagneticCharge(magneticChargeEplus) { } G4VUserPhysicsList* CreatePhysicsList() override { - G4VUserPhysicsList* physicsList = TG4RunConfiguration::CreatePhysicsList(); + G4VUserPhysicsList* physicsList = o2::fastsim::G4RunConfiguration::CreatePhysicsList(); if (auto* composed = dynamic_cast(physicsList)) { composed->AddPhysicsList(new O2MonopolePhysics(mMagneticCharge)); LOG(info) << "O2G4RunConfiguration: monopole ionisation physics registered " From 1d62ef1e42684746f2d8bff3a480988931840faa Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Tue, 22 Sep 2026 15:08:36 +0200 Subject: [PATCH 05/11] Implemented Sandro comments + give monopoles own transport process Previous version was configuring O2MonopolePhysics globally which replaced the chord finder on the run G4FieldManager and also called G4Transportation::EnableMa gneticMoment(true), which are botbh global switched. Since NystromRK4 is incompatible with the monopole physics, the G4ClassicalRK4 was used, which I think is wrong, as I wanted to keep the original stepper the VMC macro asks for. Now the monopoles chord finder is built and not installed. In addition some minor things that were raised by Sandro's review and found while rechecking the code: - Now SetMaxstep(1.e10) is not called for monopoles, as that limit belongs to the dri ft volume rather than the track - TPC and ITS ProcessHits look up the monopole PDG code for neutral tracks only (reduction in calls) - Throw fatal in case Monopole Physics is requested but not installed properly The commit, together with the rest has been validated by running with Schwinger pair monopoles + various gD options. --- .../ITSMFT/ITS/simulation/src/Detector.cxx | 9 +- Detectors/TPC/simulation/src/Detector.cxx | 17 +- Detectors/gconfig/src/O2MonopolePhysics.cxx | 162 +++++++++++++++--- 3 files changed, 154 insertions(+), 34 deletions(-) diff --git a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx index 86dee749fd807..e5fffe20f6fa2 100644 --- a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx +++ b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx @@ -321,9 +321,12 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // This method is called from the MC stepping // 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 - const bool isMonopole = o2::sim::isMonopole(fMC->TrackPid()); - if (!(fMC->TrackCharge()) && !isMonopole) { + // 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; } diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 1ae560f26f3b7..8dc577a3ef509 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -122,17 +122,22 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // 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 = o2::sim::isMonopole(trackPdg); + // 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(trackCharge) == 0) { - // Fall through only for monopoles when ionisation is enabled. - // The behaviour for the other neutral particles remains as before. - if (!isMonopole || fMC->Edep() <= 0.) { + isMonopole = 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.) { + // The monopole deposits nothing when no ionisation process is attached to + // it (G4.monopole=0), so no hit to make. + return kFALSE; + } } // ===| SET THE LENGTH OF THE NEXT ENERGY LOSS STEP |========================= diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 70ad9812ff301..6d220e7825f83 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -60,6 +60,7 @@ #include #include #include +#include #include #include @@ -234,8 +235,20 @@ class O2MonopoleEquation : public G4EquationOfMotion }; //____________________________________________________________________________ -/// Make the global field integrate O2MonopoleEquation instead of the default -/// electric-charge-only equation. +/// The field manager of the run, together with a second chord finder that +/// integrates O2MonopoleEquation +struct MonopoleFieldSetup { + G4FieldManager* fieldManager = nullptr; + G4ChordFinder* chordFinder = nullptr; +}; + +//____________________________________________________________________________ +/// Build a chord finder that integrates O2MonopoleEquation instead of the +/// default electric-charge-only equation. +/// +/// It is deliberately NOT installed on the field manager here. The default +/// chord finder is what Geant4-VMC configured for this run (NystromRK4 unless +/// the macro says otherwise) and every ordinary particle keeps it /// /// SetUserEquationOfMotion() in Geant4-VMC is not usable here: it /// registers the object with TG4GeometryManager, and the field integrator is @@ -244,20 +257,22 @@ class O2MonopoleEquation : public G4EquationOfMotion /// /// \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) +inline MonopoleFieldSetup buildMonopoleFieldSetup(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; + LOG(fatal) << "O2MonopolePhysics: no G4FieldManager, the monopole equation of motion " + "cannot be installed; rerun with G4.monopole=0 if that is what you want"; + 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; + LOG(fatal) << "O2MonopolePhysics: no magnetic field attached to the field manager, the " + "monopole equation of motion cannot be installed; rerun with G4.monopole=0 " + "if that is what you want"; + return {}; } auto* equation = new O2MonopoleEquation(magneticField, magneticChargeEplusUnits, tpcDriftFieldGeant4); @@ -273,17 +288,111 @@ inline void installMonopoleFieldIntegrator(double magneticChargeEplusUnits, doub 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)), " + LOG(info) << "O2MonopolePhysics: monopole equation of motion built (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 = " + LOG(info) << "O2MonopolePhysics: monopole equation of motion built (F = g*B), magnetic charge = " << magneticChargeEplusUnits << " eplus (TPC drift field coupling disabled)"; } + return {fieldManager, chordFinder}; +} + +//____________________________________________________________________________ +/// Transportation for the monopole species only. +/// +/// These must be true for the monopoles and false for every other particle: +/// +/// - the field manager has to use the monopole chord finder, otherwise O2MonopoleEquation +/// might not be evaluated at all (such as it happens with the default NystromRK4) +/// - G4Transportation has to consider the magnetic moment +/// +/// Ordinary tracks keep the stepper the run was configured with, and neutral particles that +/// happen to carry a magnetic moment (neutrons above all) keep their +/// straight-line transport. +/// +/// This mirrors G4MonopoleTransportation from the Geant4 monopole example, but +/// by derives directly from G4Transportation, so that it follows the Geant4 versions +class O2MonopoleTransportation : public G4Transportation +{ + public: + O2MonopoleTransportation(G4FieldManager* fieldManager, G4ChordFinder* monopoleChordFinder) + : G4Transportation(0), mFieldManager(fieldManager), mMonopoleChordFinder(monopoleChordFinder) + { + } + + G4double AlongStepGetPhysicalInteractionLength(const G4Track& track, G4double previousStepSize, + G4double currentMinimumStep, G4double& currentSafety, + G4GPILSelection* selection) override + { + const G4bool previousMoment = G4Transportation::EnableMagneticMoment(true); + auto* previousChordFinder = mFieldManager->GetChordFinder(); + mFieldManager->SetChordFinder(mMonopoleChordFinder); + + const G4double length = G4Transportation::AlongStepGetPhysicalInteractionLength( + track, previousStepSize, currentMinimumStep, currentSafety, selection); + + mFieldManager->SetChordFinder(previousChordFinder); + G4Transportation::EnableMagneticMoment(previousMoment); + return length; + } + + private: + G4FieldManager* mFieldManager; ///< field manager of the run, not owned + G4ChordFinder* mMonopoleChordFinder; ///< installed only for the duration of a monopole step, not owned +}; + +//____________________________________________________________________________ +/// Replace the transport of one monopole species with +/// O2MonopoleTransportation, keeping it first in the DoIt vectors exactly as +/// G4VUserPhysicsList::AddTransportation() left it. +inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4ParticleDefinition* particle, + const MonopoleFieldSetup& fieldSetup) +{ + G4VProcess* existing = nullptr; + G4ProcessVector* plist = pmanager->GetProcessList(); + for (G4int ip = 0; ip < static_cast(plist->size()); ++ip) { + if (dynamic_cast((*plist)[ip]) != nullptr) { + existing = (*plist)[ip]; + break; + } + } + if (existing == nullptr) { + LOG(fatal) << "O2MonopolePhysics: " << particle->GetParticleName() + << " has no transportation process to replace; the monopole could not be " + "coupled to the field"; + return false; + } + if (dynamic_cast(existing) != nullptr) { + // Parallel worlds are in use. O2MonopoleTransportation derives from plain + // G4Transportation, so swapping it in would drop the parallel-world + // navigation; refuse rather than silently mis-navigate. + LOG(fatal) << "O2MonopolePhysics: " << particle->GetParticleName() + << " uses G4CoupledTransportation (parallel worlds); the monopole transportation " + "does not support that. Rerun with G4.monopole=0 or without parallel worlds"; + return false; + } + + // One G4Transportation instance is shared by every particle + // (G4VUserPhysicsList::AddTransportation creates a single one), so it is + // detached from this particle only and must not be deleted. + pmanager->RemoveProcess(existing); + + auto* transportation = new O2MonopoleTransportation(fieldSetup.fieldManager, fieldSetup.chordFinder); + pmanager->AddProcess(transportation); + // Transportation has to be first in the DoIt vectors + // + // Geant4 prints "Set Ordering First is invoked twice for Transportation to + // " (ProcMan113, JustWarning) once per call here, because + // G4ProcessManager latches isSetOrderingFirstInvoked and RemoveProcess() does + // not clear it. The insertion is performed before that check and is correct; + // the warning in the stdout is expected and harmless. + pmanager->SetProcessOrderingToFirst(transportation, idxAlongStep); + pmanager->SetProcessOrderingToFirst(transportation, idxPostStep); + return true; } //____________________________________________________________________________ @@ -310,6 +419,13 @@ class O2MonopolePhysics : public G4VUserPhysicsList void ConstructProcess() override { auto* table = G4ParticleTable::GetParticleTable(); + + // Built once and shared by all monopoles species: O2MonopoleEquation reads the + // sign of the magnetic charge off the track, so one chord finder serves + // monopoles and anti-monopoles + const MonopoleFieldSetup fieldSetup = + buildMonopoleFieldSetup(mMagneticCharge / CLHEP::eplus, tpcDriftFieldMagnitude()); + int nAttached = 0; for (int pdg : gMonopolePDGs) { auto* particle = table->FindParticle(pdg); @@ -357,25 +473,20 @@ class O2MonopolePhysics : public G4VUserPhysicsList if (particle->GetPDGMagneticMoment() == 0.) { particle->SetPDGMagneticMoment(gMonopoleFieldGateMoment); } + + // 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. + installMonopoleTransport(pmanager, particle, fieldSetup); + ++nAttached; - LOG(info) << "O2MonopolePhysics: attached G4mplIonisation to " + LOG(info) << "O2MonopolePhysics: attached G4mplIonisation and monopole transport 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: @@ -406,8 +517,9 @@ class O2G4RunConfiguration : public o2::fastsim::G4RunConfiguration 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"; + LOG(fatal) << "O2G4RunConfiguration: physics list is not a TG4ComposedPhysicsList, " + "monopole ionisation cannot be enabled; rerun with G4.monopole=0 if that " + "is what you want"; } return physicsList; } From 26f3cee935d09843131d091bd9bca402d85797dd Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Tue, 22 Sep 2026 18:38:15 +0200 Subject: [PATCH 06/11] Make monopoles mass configurable via parameter --- Common/SimConfig/include/SimConfig/G4Params.h | 4 ++++ .../SimulationDataFormat/MonopoleParticles.h | 3 +++ .../SimulationDataFormat/O2DatabasePDG.h | 19 +++++++++++++------ Detectors/TPC/simulation/src/Detector.cxx | 8 +++++--- Detectors/gconfig/src/O2MonopolePhysics.cxx | 13 ++++++++----- Generators/src/GeneratorFactory.cxx | 5 ++++- Steer/src/O2MCApplication.cxx | 12 ++++++++---- macro/o2sim.C | 7 ++++++- 8 files changed, 51 insertions(+), 20 deletions(-) diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 54dcf861f3251..97b31b3e3c5e3 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -14,6 +14,7 @@ #include "CommonUtils/ConfigurableParam.h" #include "CommonUtils/ConfigurableParamHelper.h" +#include "SimulationDataFormat/MonopoleParticles.h" namespace o2 { @@ -67,6 +68,9 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { // 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(o2::sim::MonopoleMassDefaultGeV); // Fast simulation. Empty fastSimModels (the default) disables the feature // entirely; see Detectors/gconfig/include/SimSetup/G4FastSimulation.h. diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h index 4f618e933e11b..93b9f1f67f915 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h +++ b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h @@ -29,6 +29,9 @@ 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 { diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index a24a36134becb..2c9b966bd098a 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -19,6 +19,7 @@ #include #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 @@ -191,11 +198,11 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db) // 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); + 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", 100., kTRUE, 0.0, 0, "BSM", 4120000); - db->AddParticle("AntiMonopole_asymm", "AntiMonopole_asymm", 100., kTRUE, 0.0, 0, "BSM", -4120000); + 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 // diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 8dc577a3ef509..8359e820a03a3 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -18,6 +18,7 @@ #include "DetectorsBase/Stack.h" #include "SimulationDataFormat/TrackReference.h" #include "SimulationDataFormat/MonopoleParticles.h" +#include "SimConfig/G4Params.h" #include "FairVolume.h" // for FairVolume @@ -127,15 +128,16 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // To-do: add dyons case bool isMonopole = false; if (static_cast(trackCharge) == 0) { - isMonopole = o2::sim::isMonopole(fMC->TrackPid()); + 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.) { - // The monopole deposits nothing when no ionisation process is attached to - // it (G4.monopole=0), so no hit to make. + // An ionising monopole depositing nothing on this step: there + // is no hit to create, and the step limit is deliberately left alone return kFALSE; } } diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 6d220e7825f83..86c5a4cd36f83 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -420,11 +420,11 @@ class O2MonopolePhysics : public G4VUserPhysicsList { auto* table = G4ParticleTable::GetParticleTable(); - // Built once and shared by all monopoles species: O2MonopoleEquation reads the - // sign of the magnetic charge off the track, so one chord finder serves - // monopoles and anti-monopoles - const MonopoleFieldSetup fieldSetup = - buildMonopoleFieldSetup(mMagneticCharge / CLHEP::eplus, tpcDriftFieldMagnitude()); + // Built on the first monopole species actually found, then shared by the rest: + // O2MonopoleEquation reads the sign of the magnetic charge off the track, so + // one chord finder serves monopoles and anti-monopoles + // No monopoles in the particles table == no monopole ionisation attached + MonopoleFieldSetup fieldSetup; int nAttached = 0; for (int pdg : gMonopolePDGs) { @@ -477,6 +477,9 @@ class O2MonopolePhysics : public G4VUserPhysicsList // 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. + if (fieldSetup.chordFinder == nullptr) { + fieldSetup = buildMonopoleFieldSetup(mMagneticCharge / CLHEP::eplus, tpcDriftFieldMagnitude()); + } installMonopoleTransport(pmanager, particle, fieldSetup); ++nAttached; diff --git a/Generators/src/GeneratorFactory.cxx b/Generators/src/GeneratorFactory.cxx index f62e8dfddb849..3c3644bb4c3bb 100644 --- a/Generators/src/GeneratorFactory.cxx +++ b/Generators/src/GeneratorFactory.cxx @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -86,7 +87,9 @@ void GeneratorFactory::setPrimaryGenerator(o2::conf::SimConfig const& conf, Fair /** generators **/ - o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance()); + // Monopole configurable mass added to the PDG database + o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(), + o2::conf::G4Params::Instance().monopoleMass); auto genconfig = conf.getGenerator(); #if defined(GENERATORS_WITH_PYTHIA8) && defined(GENERATORS_WITH_HEPMC3) if (GeneratorHybridParam::Instance().switchExtToHybrid && (genconfig.compare("external") == 0 || genconfig.compare("extgen") == 0)) { diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index c8c339af506ce..c805e0f1011d6 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/Steer/src/O2MCApplication.cxx @@ -1753,12 +1753,16 @@ void addSpecialParticles() 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 + // The transported mass has to match the one the generator used; see + // G4Params.monopoleMass and o2::sim::MonopoleMassDefaultGeV. + // To-do: find a way to define multiple masses for monopoles species + const double monopoleMass = o2::conf::G4Params::Instance().monopoleMass; // 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); + TVirtualMC::GetMC()->DefineParticle(4110000, "Monopole_symm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); + TVirtualMC::GetMC()->DefineParticle(-4110000, "AntiMonopole_symm", kPTHadron, monopoleMass, 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); + TVirtualMC::GetMC()->DefineParticle(4120000, "Monopole_asymm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); + TVirtualMC::GetMC()->DefineParticle(-4120000, "AntiMonopole_asymm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); } void O2MCApplicationBase::AddParticles() diff --git a/macro/o2sim.C b/macro/o2sim.C index fd44de6e93b54..b88251705873b 100644 --- a/macro/o2sim.C +++ b/macro/o2sim.C @@ -15,6 +15,7 @@ #include #include #include "SimulationDataFormat/O2DatabasePDG.h" +#include "SimConfig/G4Params.h" #include "SimulationDataFormat/MCEventHeader.h" #include #include @@ -178,7 +179,11 @@ FairRunSim* o2sim_init(bool asservice, bool evalmat = false) run->Init(); // add ALICE particles to TDatabasePDG singleton - o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance()); + // The monopole mass has to be the one the transport uses; in the worker + // (asservice) the GeneratorFactory call is skipped, so this is the only + // place that registers it. + o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(), + o2::conf::G4Params::Instance().monopoleMass); long runStart = timestamp; { From 0a43242877b7238b892ab8b767bd441139023fc4 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 24 Sep 2026 17:34:12 +0200 Subject: [PATCH 07/11] Monopoles: changes from further review - Register mplIoni through G4PhysicsListHelper. - Move the TPC Bethe-Bloch computation back into the branch that uses it. - Run the ITS monopole PDG lookup only when G4.monopole is set. - Define the monopole particles only when G4.monopole is set. - Make non-positive monopole magnetic charge or mass fatal. - Make a mismatch between G4.monopoleMass and TDatabasePDG fatal. - Warn when the TPC drift field is zero. - Align the species comment with the code. Co-Authored-By: Claude Opus 5.5 --- .../SimulationDataFormat/MonopoleParticles.h | 6 +++-- .../SimulationDataFormat/O2DatabasePDG.h | 14 +++++++++++ .../ITSMFT/ITS/simulation/src/Detector.cxx | 12 +++++----- Detectors/TPC/simulation/src/Detector.cxx | 24 +++++++++---------- Detectors/gconfig/src/O2MonopolePhysics.cxx | 22 +++++++++++++---- Generators/src/GeneratorFactory.cxx | 7 ++++-- Steer/src/O2MCApplication.cxx | 9 ++++++- macro/o2sim.C | 7 ++++-- 8 files changed, 72 insertions(+), 29 deletions(-) diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h index 93b9f1f67f915..c1f546fd5e569 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h +++ b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h @@ -24,9 +24,11 @@ namespace o2::sim { -/// monopole carrying equal electric and magnetic charge +/// 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 carrying opposite electric and magnetic charge +/// monopole intended to carry opposite electric and magnetic charge (not yet distinguished) constexpr int MonopolePdgAsymm = 4120000; /// Default monopole mass in GeV. diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index 2c9b966bd098a..ca1c2f4c2aa20 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -16,6 +16,7 @@ #ifndef O2_O2DATABASEPDG_H #define O2_O2DATABASEPDG_H +#include #include #include "TDatabasePDG.h" #include "TParticlePDG.h" @@ -54,6 +55,19 @@ class O2DatabasePDG 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) { diff --git a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx index e5fffe20f6fa2..6f4f1fa5fe028 100644 --- a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx +++ b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx @@ -25,6 +25,7 @@ #include "DetectorsBase/Stack.h" #include "SimulationDataFormat/TrackReference.h" #include "SimulationDataFormat/MonopoleParticles.h" +#include "SimConfig/G4Params.h" #include "fairlogger/Logger.h" // for LOG, LOG_IF // FairRoot includes @@ -321,11 +322,12 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // This method is called from the MC stepping // 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. + // rejected by the electric-charge gate. The PDG lookup only runs for neutral + // particles in runs with monopole physics enabled. // To-do: handle dyons. + static const bool sMonopoleIonisation = o2::conf::G4Params::Instance().monopole; const bool isNeutral = (fMC->TrackCharge() == 0); - const bool isMonopole = isNeutral && o2::sim::isMonopole(fMC->TrackPid()); + const bool isMonopole = isNeutral && sMonopoleIonisation && o2::sim::isMonopole(fMC->TrackPid()); if (isNeutral && !isMonopole) { return kFALSE; } @@ -394,9 +396,7 @@ 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 + // A monopole that deposited nothing in the sensor leaves no hit if (isMonopole && mTrackData.mEnergyLoss <= 0.) { return kFALSE; } diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index 8359e820a03a3..7e81863d3bff4 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -216,18 +216,6 @@ Bool_t Detector::ProcessHits(FairVolume* vol) Int_t numberOfElectrons = 0; // 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(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) { // We have multiple collisions and add fluctuations: smear nel using @@ -248,6 +236,18 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // 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(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(fMC->GetRandom()->Poisson(meanNcoll)); diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 86c5a4cd36f83..dc1605743d82e 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -45,6 +45,7 @@ #include #include #include +#include #include #include #include @@ -92,7 +93,8 @@ inline bool isMonopolePDG(int pdg) { return o2::sim::isMonopole(pdg); } /// so that G4Transportation queries the field for it constexpr G4double gMonopoleFieldGateMoment = 1.0e-20; -// Extent of the TPC drift region +// Extent of the TPC field cage, i.e. the region with the drift field (smaller than +// the TPC_Drift volume, which extends beyond the endcaps and the field-cage rods) constexpr G4double gTPCFieldCageRMin = 83.5 * CLHEP::cm; constexpr G4double gTPCFieldCageRMax = 254.5 * CLHEP::cm; constexpr G4double gTPCFieldCageZMax = 249.525 * CLHEP::cm; @@ -107,6 +109,11 @@ inline double tpcDriftFieldMagnitude() } try { const double valueKVPerCm = o2::conf::ConfigurableParam::getValueAs("TPCGEMParam.ElectricField[0]"); + if (valueKVPerCm <= 0.) { + LOG(warn) << "O2MonopolePhysics: TPCGEMParam.ElectricField[0] = " << valueKVPerCm + << " kV/cm; the monopole is not coupled to a drift field"; + return 0.; + } return valueKVPerCm * CLHEP::kilovolt / CLHEP::cm; } catch (...) { LOG(warn) << "O2MonopolePhysics: the TPC is in the geometry but TPCGEMParam is not " @@ -159,6 +166,9 @@ class O2MonopoleEquation : public G4EquationOfMotion // 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. + // TODO: both species are electrically neutral today, so their magnetic sign is + // the same function of the PDG sign; once they carry electric charge the + // asymmetric species (4120000) needs the opposite relative sign. double signedMagneticCharge = 0.; int pdg = 0; if (const G4Track* track = currentTrack()) { @@ -462,9 +472,8 @@ class O2MonopolePhysics : public G4VUserPhysicsList 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); + // Ordered by G4PhysicsListHelper from the process subtype, as in the Geant4 monopole example + G4PhysicsListHelper::GetPhysicsListHelper()->RegisterProcess(new G4mplIonisation(mMagneticCharge), particle); // 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 @@ -542,6 +551,11 @@ TG4RunConfiguration* createMonopoleRunConfiguration(const TString& userGeometry, // 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. + // G4mplIonisation silently replaces a zero charge with one Dirac charge, while the + // equation of motion would use the zero, so non-positive values are rejected + if (!(magneticChargeDirac > 0.)) { + LOG(fatal) << "O2MonopolePhysics: G4.monopoleMagneticCharge must be positive, got " << magneticChargeDirac; + } const double gDirac = CLHEP::eplus / (2.0 * CLHEP::fine_structure_const); const double magneticChargeEplus = magneticChargeDirac * gDirac; LOG(info) << "Monopole ionisation enabled: magnetic charge = " << magneticChargeDirac diff --git a/Generators/src/GeneratorFactory.cxx b/Generators/src/GeneratorFactory.cxx index 3c3644bb4c3bb..1b6adad30052e 100644 --- a/Generators/src/GeneratorFactory.cxx +++ b/Generators/src/GeneratorFactory.cxx @@ -88,8 +88,11 @@ void GeneratorFactory::setPrimaryGenerator(o2::conf::SimConfig const& conf, Fair /** generators **/ // Monopole configurable mass added to the PDG database - o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(), - o2::conf::G4Params::Instance().monopoleMass); + const double monopoleMass = o2::conf::G4Params::Instance().monopoleMass; + o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(), monopoleMass); + if (!o2::O2DatabasePDG::hasMonopoleMass(TDatabasePDG::Instance(), monopoleMass)) { + LOG(fatal) << "Monopoles were registered in TDatabasePDG with a mass other than G4.monopoleMass = " << monopoleMass << " GeV"; + } auto genconfig = conf.getGenerator(); #if defined(GENERATORS_WITH_PYTHIA8) && defined(GENERATORS_WITH_HEPMC3) if (GeneratorHybridParam::Instance().switchExtToHybrid && (genconfig.compare("external") == 0 || genconfig.compare("extgen") == 0)) { diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index c805e0f1011d6..2c3647f9cd163 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/Steer/src/O2MCApplication.cxx @@ -1752,11 +1752,18 @@ void addSpecialParticles() 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 + // BSM Monopoles, defined only when their physics is enabled (G4.monopole=1), so that a + // monopole primary without it is reported by the engine instead of crossing the detector unseen. // The transported mass has to match the one the generator used; see // G4Params.monopoleMass and o2::sim::MonopoleMassDefaultGeV. // To-do: find a way to define multiple masses for monopoles species + if (!o2::conf::G4Params::Instance().monopole) { + return; + } const double monopoleMass = o2::conf::G4Params::Instance().monopoleMass; + if (!(monopoleMass > 0.)) { + LOG(fatal) << "G4.monopoleMass must be positive, got " << monopoleMass; + } // Symmetric monopoles: same electric and magnetic charge TVirtualMC::GetMC()->DefineParticle(4110000, "Monopole_symm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); TVirtualMC::GetMC()->DefineParticle(-4110000, "AntiMonopole_symm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); diff --git a/macro/o2sim.C b/macro/o2sim.C index b88251705873b..c99cd012d4d73 100644 --- a/macro/o2sim.C +++ b/macro/o2sim.C @@ -182,8 +182,11 @@ FairRunSim* o2sim_init(bool asservice, bool evalmat = false) // The monopole mass has to be the one the transport uses; in the worker // (asservice) the GeneratorFactory call is skipped, so this is the only // place that registers it. - o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(), - o2::conf::G4Params::Instance().monopoleMass); + const double monopoleMass = o2::conf::G4Params::Instance().monopoleMass; + o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance(), monopoleMass); + if (!o2::O2DatabasePDG::hasMonopoleMass(TDatabasePDG::Instance(), monopoleMass)) { + LOG(fatal) << "Monopoles were registered in TDatabasePDG with a mass other than G4.monopoleMass = " << monopoleMass << " GeV"; + } long runStart = timestamp; { From e74a2e0ea1eea43d2da9995877ade8717e5e1b1f Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Wed, 7 Oct 2026 12:37:46 +0200 Subject: [PATCH 08/11] Implement Ivana code review + clean up code + minor fixes --- .../include/SimSetup/O2MonopolePhysics.h | 6 +- Detectors/gconfig/src/O2MonopolePhysics.cxx | 162 ++++++++---------- Steer/src/O2MCApplication.cxx | 10 +- 3 files changed, 84 insertions(+), 94 deletions(-) diff --git a/Detectors/gconfig/include/SimSetup/O2MonopolePhysics.h b/Detectors/gconfig/include/SimSetup/O2MonopolePhysics.h index 46f67623d2191..7e8508552c3e7 100644 --- a/Detectors/gconfig/include/SimSetup/O2MonopolePhysics.h +++ b/Detectors/gconfig/include/SimSetup/O2MonopolePhysics.h @@ -23,12 +23,12 @@ namespace o2 namespace g4config { -/// Build a Geant4-VMC run configuration which attaches the magnetic-monopole ionisation process -/// (G4mplIonisation) to the requested reference physics list. +/// Build a Geant4-VMC run configuration which, on top of the requested reference physics list, +/// attaches the magnetic-monopole ionisation process (G4mplIonisation) and transport to /// the monopole particles already defined by O2 (PDG +-4110000 / +-4120000). /// /// Only used when G4Params.monopole == true; -/// By default, the standard TG4RunConfiguration is created +/// Otherwise o2::g4config::G4RunConfiguration is used (see g4Config.C) /// /// \param userGeometry VMC geometry-navigation string (as for TG4RunConfiguration) /// \param physicsList the reference physics-list selection string diff --git a/Detectors/gconfig/src/O2MonopolePhysics.cxx b/Detectors/gconfig/src/O2MonopolePhysics.cxx index 8f1eb72dda631..1bb6eeaf97852 100644 --- a/Detectors/gconfig/src/O2MonopolePhysics.cxx +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -38,7 +38,6 @@ #include -#include "TG4RunConfiguration.h" #include "TG4ComposedPhysicsList.h" #include "SimSetup/G4RunConfiguration.h" @@ -48,19 +47,16 @@ #include #include #include -#include #include #include #include #include #include -#include #include -#include #include +#include #include -#include #include #include #include @@ -70,7 +66,6 @@ #include #include -#include namespace o2 { @@ -86,8 +81,6 @@ namespace constexpr std::array gMonopolePDGs = {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm, o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm}; -inline bool isMonopolePDG(int pdg) { return o2::sim::isMonopole(pdg); } - /// 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 @@ -157,27 +150,16 @@ class O2MonopoleEquation : public G4EquationOfMotion { } + /// Sign of the magnetic charge of the monopole about to be transported (+1 or -1), + /// set by O2MonopoleTransportation before each step + void setMagneticChargeSign(int sign) { mMagneticChargeSign = sign; } + 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. - // TODO: both species are electrically neutral today, so their magnetic sign is - // the same function of the PDG sign; once they carry electric charge the - // asymmetric species (4120000) needs the opposite relative sign. - 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; + mMagCharge = CLHEP::eplus * mMagneticChargeSign * mMagneticChargeEplus * CLHEP::c_light; } void EvaluateRhsGivenB(const G4double y[], const G4double B[3], G4double dydx[]) const override @@ -205,8 +187,7 @@ class O2MonopoleEquation : public G4EquationOfMotion // 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.) { + if (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 @@ -227,17 +208,8 @@ class O2MonopoleEquation : public G4EquationOfMotion } 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) + int mMagneticChargeSign = 1; ///< sign of the magnetic charge of the current monopole G4double mTPCDriftField = 0.; ///< TPC drift field, Geant4 units; 0 disables the coupling G4double mElCharge = 0.; G4double mMagCharge = 0.; @@ -245,20 +217,23 @@ class O2MonopoleEquation : public G4EquationOfMotion }; //____________________________________________________________________________ -/// The field manager of the run, together with a second chord finder that -/// integrates O2MonopoleEquation +/// The chord finder integrating O2MonopoleEquation, and that equation struct MonopoleFieldSetup { - G4FieldManager* fieldManager = nullptr; G4ChordFinder* chordFinder = nullptr; + O2MonopoleEquation* equation = nullptr; ///< to pass the sign of the magnetic charge at each step }; //____________________________________________________________________________ /// Build a chord finder that integrates O2MonopoleEquation instead of the -/// default electric-charge-only equation. +/// default electric-charge-only equation, in the field of the run. /// -/// It is deliberately NOT installed on the field manager here. The default -/// chord finder is what Geant4-VMC configured for this run (NystromRK4 unless -/// the macro says otherwise) and every ordinary particle keeps it +/// It is deliberately NOT installed on a field manager here. The default +/// chord finder is what Geant4-VMC configured for the field (NystromRK4 unless +/// the macro says otherwise) and every ordinary particle keeps it. +/// O2MonopoleTransportation installs it only for the duration of a monopole step, on the +/// field manager of the current volume: the global one, or one of the local ones created for +/// the volumes with their own field parameters (G4LocalFieldConstruction). The local field +/// managers carry the same field as the global one, so one chord finder serves them all. /// /// SetUserEquationOfMotion() in Geant4-VMC is not usable here: it /// registers the object with TG4GeometryManager, and the field integrator is @@ -308,7 +283,7 @@ inline MonopoleFieldSetup buildMonopoleFieldSetup(double magneticChargeEplusUnit LOG(info) << "O2MonopolePhysics: monopole equation of motion built (F = g*B), magnetic charge = " << magneticChargeEplusUnits << " eplus (TPC drift field coupling disabled)"; } - return {fieldManager, chordFinder}; + return {chordFinder, equation}; } //____________________________________________________________________________ @@ -316,50 +291,82 @@ inline MonopoleFieldSetup buildMonopoleFieldSetup(double magneticChargeEplusUnit /// /// These must be true for the monopoles and false for every other particle: /// -/// - the field manager has to use the monopole chord finder, otherwise O2MonopoleEquation -/// might not be evaluated at all (such as it happens with the default NystromRK4) +/// - the field manager of the current volume (the global one, or a local one) has to use the +/// monopole chord finder, otherwise O2MonopoleEquation might not be evaluated at all +/// (such as it happens with the default NystromRK4) /// - G4Transportation has to consider the magnetic moment +/// - the field manager has to declare that the field changes the energy: the g*B force does +/// work on the monopole. Otherwise G4Transportation would reset the kinetic energy at the +/// end of every step to the one at its start, keeping only the change of direction /// /// Ordinary tracks keep the stepper the run was configured with, and neutral particles that /// happen to carry a magnetic moment (neutrons above all) keep their /// straight-line transport. /// -/// This mirrors G4MonopoleTransportation from the Geant4 monopole example, but -/// by derives directly from G4Transportation, so that it follows the Geant4 versions +/// This mirrors G4MonopoleTransportation from the Geant4 monopole example, which takes the +/// kinetic energy and the time of flight from the integrated track, but derives directly +/// from G4Transportation, so that it follows the Geant4 versions class O2MonopoleTransportation : public G4Transportation { public: - O2MonopoleTransportation(G4FieldManager* fieldManager, G4ChordFinder* monopoleChordFinder) - : G4Transportation(0), mFieldManager(fieldManager), mMonopoleChordFinder(monopoleChordFinder) + O2MonopoleTransportation(const MonopoleFieldSetup& fieldSetup) + : G4Transportation(0), mChordFinder(fieldSetup.chordFinder), mEquation(fieldSetup.equation) { } + void StartTracking(G4Track* track) override + { + G4Transportation::StartTracking(track); + // G4Transportation::StartTracking() clears the integration state only of the chord finders + // installed on the field managers. Without this the step estimate of the previous monopole + // would carry over, and the result would depend on which events the worker simulated before + mChordFinder->ResetStepEstimate(); + } + G4double AlongStepGetPhysicalInteractionLength(const G4Track& track, G4double previousStepSize, G4double currentMinimumStep, G4double& currentSafety, G4GPILSelection* selection) override { + // the field manager the propagator uses in this volume, as G4Transportation finds it + auto* fieldManager = GetPropagatorInField()->FindAndSetFieldManager(track.GetVolume()); + if (fieldManager->GetDetectorField() == nullptr) { + // volume without field: straight-line transport + return G4Transportation::AlongStepGetPhysicalInteractionLength(track, previousStepSize, currentMinimumStep, + currentSafety, selection); + } + // The sign of the magnetic charge follows the PDG sign, so monopole and anti-monopole + // are pushed in opposite directions along B. + // TODO: both species are electrically neutral currently, so their magnetic sign is + // the same function of the PDG sign; once they carry electric charge the + // asymmetric species (4120000) needs the opposite relative sign. + mEquation->setMagneticChargeSign(track.GetDefinition()->GetPDGEncoding() > 0 ? 1 : -1); + const G4bool previousMoment = G4Transportation::EnableMagneticMoment(true); - auto* previousChordFinder = mFieldManager->GetChordFinder(); - mFieldManager->SetChordFinder(mMonopoleChordFinder); + auto* previousChordFinder = fieldManager->GetChordFinder(); + const G4bool previousChangesEnergy = fieldManager->DoesFieldChangeEnergy(); + fieldManager->SetChordFinder(mChordFinder); + // the end kinetic energy and time of flight are then taken from the integrated track + fieldManager->SetFieldChangesEnergy(true); const G4double length = G4Transportation::AlongStepGetPhysicalInteractionLength( track, previousStepSize, currentMinimumStep, currentSafety, selection); - mFieldManager->SetChordFinder(previousChordFinder); + fieldManager->SetFieldChangesEnergy(previousChangesEnergy); + fieldManager->SetChordFinder(previousChordFinder); G4Transportation::EnableMagneticMoment(previousMoment); return length; } private: - G4FieldManager* mFieldManager; ///< field manager of the run, not owned - G4ChordFinder* mMonopoleChordFinder; ///< installed only for the duration of a monopole step, not owned + G4ChordFinder* mChordFinder; ///< installed only for the duration of a monopole step, not owned + O2MonopoleEquation* mEquation; ///< the equation integrated by mChordFinder, not owned }; //____________________________________________________________________________ /// Replace the transport of one monopole species with /// O2MonopoleTransportation, keeping it first in the DoIt vectors exactly as /// G4VUserPhysicsList::AddTransportation() left it. -inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4ParticleDefinition* particle, +inline void installMonopoleTransport(G4ProcessManager* pmanager, const G4ParticleDefinition* particle, const MonopoleFieldSetup& fieldSetup) { G4VProcess* existing = nullptr; @@ -374,7 +381,7 @@ inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4Particl LOG(fatal) << "O2MonopolePhysics: " << particle->GetParticleName() << " has no transportation process to replace; the monopole could not be " "coupled to the field"; - return false; + return; } if (dynamic_cast(existing) != nullptr) { // Parallel worlds are in use. O2MonopoleTransportation derives from plain @@ -383,7 +390,7 @@ inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4Particl LOG(fatal) << "O2MonopolePhysics: " << particle->GetParticleName() << " uses G4CoupledTransportation (parallel worlds); the monopole transportation " "does not support that. Rerun with G4.monopole=0 or without parallel worlds"; - return false; + return; } // One G4Transportation instance is shared by every particle @@ -391,7 +398,7 @@ inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4Particl // detached from this particle only and must not be deleted. pmanager->RemoveProcess(existing); - auto* transportation = new O2MonopoleTransportation(fieldSetup.fieldManager, fieldSetup.chordFinder); + auto* transportation = new O2MonopoleTransportation(fieldSetup); pmanager->AddProcess(transportation); // Transportation has to be first in the DoIt vectors // @@ -402,17 +409,17 @@ inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4Particl // the warning in the stdout is expected and harmless. pmanager->SetProcessOrderingToFirst(transportation, idxAlongStep); pmanager->SetProcessOrderingToFirst(transportation, idxPostStep); - return true; } //____________________________________________________________________________ -/// Minimal physics list whose only job is to attach the magnetic-monopole -/// ionisation process to the already-defined O2 monopole particles. +/// Minimal physics list attaching the magnetic-monopole ionisation and transport +/// 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. +/// any particle (they are defined by O2MCApplication::AddParticles()). For each +/// monopole it adds mplIoni and replaces the transportation of the reference +/// physics list with O2MonopoleTransportation. class O2MonopolePhysics : public G4VUserPhysicsList { public: @@ -431,7 +438,7 @@ class O2MonopolePhysics : public G4VUserPhysicsList auto* table = G4ParticleTable::GetParticleTable(); // Built on the first monopole species actually found, then shared by the rest: - // O2MonopoleEquation reads the sign of the magnetic charge off the track, so + // O2MonopoleTransportation sets the sign of the magnetic charge at each step, so // one chord finder serves monopoles and anti-monopoles // No monopoles in the particles table == no monopole ionisation attached MonopoleFieldSetup fieldSetup; @@ -451,28 +458,9 @@ class O2MonopolePhysics : public G4VUserPhysicsList 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); - } - } - // Ordered by G4PhysicsListHelper from the process subtype, as in the Geant4 monopole example + // The only ionisation of the monopoles: they are defined as kPTUndefined, for which + // Geant4-VMC attaches no process. Ordered by G4PhysicsListHelper from the process + // subtype, as in the Geant4 monopole example G4PhysicsListHelper::GetPhysicsListHelper()->RegisterProcess(new G4mplIonisation(mMagneticCharge), particle); // G4Transportation only looks up the field when the particle has a diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index 9230451c473fe..82fdda4eace58 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/Steer/src/O2MCApplication.cxx @@ -1764,12 +1764,14 @@ void addSpecialParticles() if (!(monopoleMass > 0.)) { LOG(fatal) << "G4.monopoleMass must be positive, got " << monopoleMass; } + // kPTUndefined: Geant4-VMC attaches no process to them (a kPTHadron would get hIoni and hmsc); + // their ionisation and transport are attached by O2MonopolePhysics // Symmetric monopoles: same electric and magnetic charge - TVirtualMC::GetMC()->DefineParticle(4110000, "Monopole_symm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); - TVirtualMC::GetMC()->DefineParticle(-4110000, "AntiMonopole_symm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); + TVirtualMC::GetMC()->DefineParticle(4110000, "Monopole_symm", kPTUndefined, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); + TVirtualMC::GetMC()->DefineParticle(-4110000, "AntiMonopole_symm", kPTUndefined, monopoleMass, 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, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); - TVirtualMC::GetMC()->DefineParticle(-4120000, "AntiMonopole_asymm", kPTHadron, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); + TVirtualMC::GetMC()->DefineParticle(4120000, "Monopole_asymm", kPTUndefined, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); + TVirtualMC::GetMC()->DefineParticle(-4120000, "AntiMonopole_asymm", kPTUndefined, monopoleMass, 0.0, 1e10, "BSM", 0.0, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE); } void O2MCApplicationBase::AddParticles() From dafdd1909e3f54a4f67a4adab0c452e792d54909 Mon Sep 17 00:00:00 2001 From: tubagundem Date: Tue, 13 Oct 2026 00:07:53 +0200 Subject: [PATCH 09/11] Store the monopole mass in the MCEventHeader and use it when reading kinematics The monopole mass of a simulation (G4.monopoleMass) was only known to the simulation itself. Consumers of the output (digitizer, reconstruction MC labels, QC, analysis) took the compiled-in default of 100 GeV through O2DatabasePDG, so MCTrack::GetMass(), GetEnergy(), GetRapidity() and Get4Momentum() were wrong for monopoles of any other mass. - O2HitMerger: apply the G4.* config key-values (they were not applied in the merger process) and write the mass to the MCEventHeader under o2::sim::MonopoleMassHeaderKey when G4.monopole is set. - MCKinematicsReader: take the mass from the header when headers or tracks are loaded and hand it to O2DatabasePDG; warn if kinematics of different masses are read in the same process. - O2DatabasePDG: setMonopoleMass()/getMonopoleMass(); Mass() returns the set mass for the monopole PDG codes. Unset by default, so the simulation itself keeps using TDatabasePDG. Limitations: entries already registered in a TDatabasePDG keep their mass (ROOT has no setter); G4 parameters given through --configFile are not applied in the merger; files produced before this change have no key. Co-Authored-By: Claude Sonnet 5.5 --- .../SimulationDataFormat/MonopoleParticles.h | 4 ++ .../SimulationDataFormat/O2DatabasePDG.h | 60 +++++++++++++++---- Steer/src/MCKinematicsReader.cxx | 37 ++++++++++++ run/O2HitMerger.cxx | 26 ++++++++ 4 files changed, 114 insertions(+), 13 deletions(-) diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h index c1f546fd5e569..ca3b0d8ff0979 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h +++ b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h @@ -34,6 +34,10 @@ constexpr int MonopolePdgAsymm = 4120000; /// Default monopole mass in GeV. constexpr double MonopoleMassDefaultGeV = 100.; +/// Key under which the monopole mass (GeV) used in a simulation is stored in the MCEventHeader, +/// so that readers of the kinematics can use the mass the events were produced with. +constexpr const char* MonopoleMassHeaderKey = "monopoleMass"; + /// true for the monopole and antimonopole of both species constexpr bool isMonopole(int pdg) noexcept { diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index 2a4743c8c87d3..c6302f6dec523 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -16,6 +16,7 @@ #ifndef O2_O2DATABASEPDG_H #define O2_O2DATABASEPDG_H +#include #include #include #include "TDatabasePDG.h" @@ -39,7 +40,8 @@ class O2DatabasePDG static bool initialized = false; // initialize this --> adds particles to TDatabasePDG; auto db = TDatabasePDG::Instance(); if (!initialized) { - addALICEParticles(db); + // uses the monopole mass of the simulation if a reader has already set it + addALICEParticles(db, getMonopoleMass()); initialized = true; } return db; @@ -68,6 +70,26 @@ class O2DatabasePDG return true; } + // Sets the mass (GeV) of the monopole species for consumers of simulated data, which do not + // have access to G4Params.monopoleMass. It is taken from the MCEventHeader of the + // simulation (key o2::sim::MonopoleMassHeaderKey) by MCKinematicsReader. + // It only affects Mass() (hence MCTrack::GetMass() / GetEnergy() / GetRapidity()) and + // Instance() if called before the first Instance(): particles already registered in a + // TDatabasePDG keep their mass, since ROOT offers no way to change it afterwards. + // Not used (and not needed) in the simulation itself, where TDatabasePDG is registered + // with G4Params.monopoleMass. + static void setMonopoleMass(double massGeV) { monopoleMassOverride().store(massGeV); } + + // mass set by setMonopoleMass(), or the default if it was never set + static double getMonopoleMass() + { + const double mass = monopoleMassOverride().load(); + return mass > 0. ? mass : o2::sim::MonopoleMassDefaultGeV; + } + + // true if a monopole mass was set by setMonopoleMass() + static bool hasMonopoleMassOverride() { return monopoleMassOverride().load() > 0.; } + // get particle's (if any) mass static Double_t MassImpl(TParticlePDG* particle, bool& success) { @@ -82,6 +104,11 @@ class O2DatabasePDG // determine particle to get mass for based on PDG static Double_t Mass(int pdg, bool& success, TDatabasePDG* db = O2DatabasePDG::Instance()) { + if (hasMonopoleMassOverride() && o2::sim::isMonopole(pdg)) { + // the table may have been filled with the default mass before the simulation mass was known + success = true; + return getMonopoleMass(); + } if (pdg < IONBASELOW || pdg > IONBASEHIGH) { // not an ion, return immediately return MassImpl(db->GetParticle(pdg), success); @@ -99,6 +126,13 @@ class O2DatabasePDG O2DatabasePDG() = delete; private: + // negative: not set (the TDatabasePDG content is used) + static std::atomic& monopoleMassOverride() + { + static std::atomic mass{-1.}; + return mass; + } + static constexpr int IONBASELOW{1000000000}; static constexpr int IONBASEHIGH{1099999999}; }; @@ -161,7 +195,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa db->AddParticle("A_0+ ", " ", 0.9960, kFALSE, 0.0, +3, "meson", 9000211); db->AddParticle("A_0- ", " ", 0.9960, kFALSE, 0.0, -3, "meson", -9000211); - //db->AddParticle("F0P0 ", " ", 0.9960, kFALSE, 0.0, 0, "meson", 9010221); + // db->AddParticle("F0P0 ", " ", 0.9960, kFALSE, 0.0, 0, "meson", 9010221); db->AddParticle("KDL_2+ ", " ", 1.773, kFALSE, 0.0, +3, "meson", 10325); db->AddParticle("KDL_2- ", " ", 1.773, kFALSE, 0.0, -3, "meson", -10325); @@ -261,7 +295,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa db->AddParticle("FeedbackPhoton", "FeedbackPhoton", 0, kFALSE, 0, 0, "Special", kspe + 51); - //Hyper nuclei and exotica + // Hyper nuclei and exotica ionCode = 1010010030; if (!db->GetParticle(ionCode)) { db->AddParticle("HyperTriton", "HyperTriton", 2.991134, kFALSE, @@ -274,25 +308,25 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa 2.5e-15, 3, "Ion", ionCode); } - //hyper hydrogen 4 ground state + // hyper hydrogen 4 ground state ionCode = 1010010040; if (!db->GetParticle(ionCode)) { db->AddParticle("Hyperhydrog4", "Hyperhydrog4", 3.922434, kFALSE, 2.5e-15, 3, "Ion", ionCode); } - //anti hyper hydrogen 4 ground state + // anti hyper hydrogen 4 ground state ionCode = -1010010040; if (!db->GetParticle(ionCode)) { db->AddParticle("AntiHyperhydrog4", "AntiHyperhydrog4", 3.922434, kFALSE, 2.5e-15, 3, "Ion", ionCode); } - //hyper helium 4 ground state + // hyper helium 4 ground state ionCode = 1010020040; if (!db->GetParticle(ionCode)) { db->AddParticle("Hyperhelium4", "Hyperhelium4", 3.921728, kFALSE, 2.5e-15, 6, "Ion", ionCode); } - //anti hyper helium 4 ground state + // anti hyper helium 4 ground state ionCode = -1010020040; if (!db->GetParticle(ionCode)) { db->AddParticle("AntiHyperhelium4", "AntiHyperhelium4", 3.921728, kFALSE, @@ -642,9 +676,9 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa 0.040, -3, "Resonance", -ionCode); } - //Additional Hidden Strangeness Pentaquarks + // Additional Hidden Strangeness Pentaquarks - //Ps +/- + // Ps +/- ionCode = 9322132; if (!db->GetParticle(ionCode)) { db->AddParticle("Ps_Plus_1870", "Ps_Plus_1870", 1.870, kFALSE, @@ -682,7 +716,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa 0.10, -3, "Resonance", -ionCode); } - //Ps0 + // Ps0 ionCode = 9322131; if (!db->GetParticle(ionCode)) { db->AddParticle("Ps_0_1870", "Ps_0_1870", 1.870, kFALSE, @@ -786,14 +820,14 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa // d*(2380) - dibaryon resonance - ionCode = 900010020; //Arbitrary choice - as deuteron but with leading 9 instead of 10 + ionCode = 900010020; // Arbitrary choice - as deuteron but with leading 9 instead of 10 if (!db->GetParticle(ionCode)) { db->AddParticle("d*_2380", "d*_2380", 2.38, kFALSE, 0.070, 3, "Resonance", ionCode); } db->AddAntiParticle("d*_2380_bar", -ionCode); - //Sexaquark (uuddss): compact, neutral and stable hypothetical bound state (arxiv.org/abs/1708.08951) + // Sexaquark (uuddss): compact, neutral and stable hypothetical bound state (arxiv.org/abs/1708.08951) ionCode = 900000020; if (!db->GetParticle(ionCode)) { db->AddParticle("Sexaquark", "Sexaquark", 2.0, kTRUE, 0.0, 0, "Special", ionCode); @@ -826,4 +860,4 @@ inline void O2DatabasePDG::addParticlesFromExternalFile(TDatabasePDG* db) } // namespace o2 -#endif //O2_O2DATABASEPDG_H +#endif // O2_O2DATABASEPDG_H diff --git a/Steer/src/MCKinematicsReader.cxx b/Steer/src/MCKinematicsReader.cxx index 21024dba78368..08a400354494e 100644 --- a/Steer/src/MCKinematicsReader.cxx +++ b/Steer/src/MCKinematicsReader.cxx @@ -13,6 +13,8 @@ #include "Steer/MCKinematicsReader.h" #include "SimulationDataFormat/MCEventHeader.h" #include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/MonopoleParticles.h" +#include "SimulationDataFormat/O2DatabasePDG.h" #include #include #include @@ -82,6 +84,29 @@ void MCKinematicsReader::initTracksForSource(int source) const } } +namespace +{ +// Makes the monopole mass used by the simulation, which is stored in the event header, the +// one of O2DatabasePDG; otherwise MCTrack::GetMass() (and GetEnergy(), GetRapidity(), ...) would +// give the default mass for monopoles of a run with another G4.monopoleMass. +void applyMonopoleMass(o2::dataformats::MCEventHeader const& header) +{ + if (!header.hasInfo(o2::sim::MonopoleMassHeaderKey)) { + return; + } + bool valid = false; + const double mass = header.getInfo(o2::sim::MonopoleMassHeaderKey, valid); + if (!valid || !(mass > 0.)) { + return; + } + if (o2::O2DatabasePDG::hasMonopoleMassOverride() && std::abs(o2::O2DatabasePDG::getMonopoleMass() - mass) > 1.e-6 * mass) { + LOG(warn) << "Kinematics with different monopole masses (" << o2::O2DatabasePDG::getMonopoleMass() << " and " << mass + << " GeV) are read in the same process; monopole masses will be wrong for one of them"; + } + o2::O2DatabasePDG::setMonopoleMass(mass); +} +} // namespace + void MCKinematicsReader::loadTracksForSourceAndEvent(int source, int event) const { auto chain = mInputChains[source]; @@ -95,6 +120,17 @@ void MCKinematicsReader::loadTracksForSourceAndEvent(int source, int event) cons mTracks[source][event] = new std::vector; *mTracks[source][event] = *loadtracks; delete loadtracks; + // the monopole mass has to be known before the tracks are used + if (auto hbr = chain->GetBranch("MCEventHeader.")) { + o2::dataformats::MCEventHeader* header = nullptr; + hbr->SetAddress(&header); + hbr->GetEntry(event); + if (header) { + applyMonopoleMass(*header); + } + delete header; + hbr->ResetAddress(); + } } } } @@ -120,6 +156,7 @@ void MCKinematicsReader::loadHeadersForSource(int source) const for (int event = 0; event < br->GetEntries(); ++event) { br->GetEntry(event); mHeaders[source][event] = *header; + applyMonopoleMass(mHeaders[source][event]); } delete header; header = nullptr; diff --git a/run/O2HitMerger.cxx b/run/O2HitMerger.cxx index b79c18bd5ca7a..fc9739c75b976 100644 --- a/run/O2HitMerger.cxx +++ b/run/O2HitMerger.cxx @@ -41,6 +41,9 @@ #include "PrimaryServerState.h" #include +#include +#include +#include #include #include #include @@ -147,6 +150,23 @@ void O2HitMerger::InitTask() // we didn't manage to get a configuration --> better to fail LOG(fatal) << "No configuration received. Aborting"; } + // The config key-values are not applied to the ConfigurableParams in this process. The merger needs the G4 ones + // to record the monopole mass in the event headers (G4.monopole, G4.monopoleMass), so apply those. + // Note: G4 parameters given via a config file (--configFile) are not applied here. + { + std::string g4Tokens; + std::stringstream tokens(o2::conf::SimConfig::Instance().getKeyValueString()); + std::string token; + while (std::getline(tokens, token, ';')) { + const auto first = token.find_first_not_of(" \t"); + if (first != std::string::npos && token.compare(first, 3, "G4.") == 0) { + g4Tokens += token.substr(first) + ";"; + } + } + if (!g4Tokens.empty()) { + o2::conf::ConfigurableParam::updateFromString(g4Tokens); + } + } mAsService = o2::conf::SimConfig::Instance().asService(); mForwardKine = o2::conf::SimConfig::Instance().forwardKine(); mWriteToDisc = o2::conf::SimConfig::Instance().writeToDisc(); @@ -694,6 +714,12 @@ bool O2HitMerger::mergeAndFlushData() } } + // record the mass used for the monopoles, such that readers of the kinematics + // (MCTrack::GetMass() and derived quantities) do not fall back to the default + if (eventheader && o2::conf::G4Params::Instance().monopole) { + eventheader->putInfo(o2::sim::MonopoleMassHeaderKey, o2::conf::G4Params::Instance().monopoleMass); + } + // now see which events can be discarded in any case due to no hits if (confref.isFilterOutNoHitEvents()) { if (eventheader && eventheader->getMCEventStats().getNHits() == 0) { From d4421b0cfcfc771567ae04af3cac32f17845c8fe Mon Sep 17 00:00:00 2001 From: tgundem Date: Tue, 13 Oct 2026 21:11:59 +0200 Subject: [PATCH 10/11] Revert "Store the monopole mass in the MCEventHeader and use it when reading kinematics" This reverts commit dafdd1909e3f54a4f67a4adab0c452e792d54909. --- .../SimulationDataFormat/MonopoleParticles.h | 4 -- .../SimulationDataFormat/O2DatabasePDG.h | 60 ++++--------------- Steer/src/MCKinematicsReader.cxx | 37 ------------ run/O2HitMerger.cxx | 26 -------- 4 files changed, 13 insertions(+), 114 deletions(-) diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h index ca3b0d8ff0979..c1f546fd5e569 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h +++ b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h @@ -34,10 +34,6 @@ constexpr int MonopolePdgAsymm = 4120000; /// Default monopole mass in GeV. constexpr double MonopoleMassDefaultGeV = 100.; -/// Key under which the monopole mass (GeV) used in a simulation is stored in the MCEventHeader, -/// so that readers of the kinematics can use the mass the events were produced with. -constexpr const char* MonopoleMassHeaderKey = "monopoleMass"; - /// true for the monopole and antimonopole of both species constexpr bool isMonopole(int pdg) noexcept { diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index c6302f6dec523..2a4743c8c87d3 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -16,7 +16,6 @@ #ifndef O2_O2DATABASEPDG_H #define O2_O2DATABASEPDG_H -#include #include #include #include "TDatabasePDG.h" @@ -40,8 +39,7 @@ class O2DatabasePDG static bool initialized = false; // initialize this --> adds particles to TDatabasePDG; auto db = TDatabasePDG::Instance(); if (!initialized) { - // uses the monopole mass of the simulation if a reader has already set it - addALICEParticles(db, getMonopoleMass()); + addALICEParticles(db); initialized = true; } return db; @@ -70,26 +68,6 @@ class O2DatabasePDG return true; } - // Sets the mass (GeV) of the monopole species for consumers of simulated data, which do not - // have access to G4Params.monopoleMass. It is taken from the MCEventHeader of the - // simulation (key o2::sim::MonopoleMassHeaderKey) by MCKinematicsReader. - // It only affects Mass() (hence MCTrack::GetMass() / GetEnergy() / GetRapidity()) and - // Instance() if called before the first Instance(): particles already registered in a - // TDatabasePDG keep their mass, since ROOT offers no way to change it afterwards. - // Not used (and not needed) in the simulation itself, where TDatabasePDG is registered - // with G4Params.monopoleMass. - static void setMonopoleMass(double massGeV) { monopoleMassOverride().store(massGeV); } - - // mass set by setMonopoleMass(), or the default if it was never set - static double getMonopoleMass() - { - const double mass = monopoleMassOverride().load(); - return mass > 0. ? mass : o2::sim::MonopoleMassDefaultGeV; - } - - // true if a monopole mass was set by setMonopoleMass() - static bool hasMonopoleMassOverride() { return monopoleMassOverride().load() > 0.; } - // get particle's (if any) mass static Double_t MassImpl(TParticlePDG* particle, bool& success) { @@ -104,11 +82,6 @@ class O2DatabasePDG // determine particle to get mass for based on PDG static Double_t Mass(int pdg, bool& success, TDatabasePDG* db = O2DatabasePDG::Instance()) { - if (hasMonopoleMassOverride() && o2::sim::isMonopole(pdg)) { - // the table may have been filled with the default mass before the simulation mass was known - success = true; - return getMonopoleMass(); - } if (pdg < IONBASELOW || pdg > IONBASEHIGH) { // not an ion, return immediately return MassImpl(db->GetParticle(pdg), success); @@ -126,13 +99,6 @@ class O2DatabasePDG O2DatabasePDG() = delete; private: - // negative: not set (the TDatabasePDG content is used) - static std::atomic& monopoleMassOverride() - { - static std::atomic mass{-1.}; - return mass; - } - static constexpr int IONBASELOW{1000000000}; static constexpr int IONBASEHIGH{1099999999}; }; @@ -195,7 +161,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa db->AddParticle("A_0+ ", " ", 0.9960, kFALSE, 0.0, +3, "meson", 9000211); db->AddParticle("A_0- ", " ", 0.9960, kFALSE, 0.0, -3, "meson", -9000211); - // db->AddParticle("F0P0 ", " ", 0.9960, kFALSE, 0.0, 0, "meson", 9010221); + //db->AddParticle("F0P0 ", " ", 0.9960, kFALSE, 0.0, 0, "meson", 9010221); db->AddParticle("KDL_2+ ", " ", 1.773, kFALSE, 0.0, +3, "meson", 10325); db->AddParticle("KDL_2- ", " ", 1.773, kFALSE, 0.0, -3, "meson", -10325); @@ -295,7 +261,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa db->AddParticle("FeedbackPhoton", "FeedbackPhoton", 0, kFALSE, 0, 0, "Special", kspe + 51); - // Hyper nuclei and exotica + //Hyper nuclei and exotica ionCode = 1010010030; if (!db->GetParticle(ionCode)) { db->AddParticle("HyperTriton", "HyperTriton", 2.991134, kFALSE, @@ -308,25 +274,25 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa 2.5e-15, 3, "Ion", ionCode); } - // hyper hydrogen 4 ground state + //hyper hydrogen 4 ground state ionCode = 1010010040; if (!db->GetParticle(ionCode)) { db->AddParticle("Hyperhydrog4", "Hyperhydrog4", 3.922434, kFALSE, 2.5e-15, 3, "Ion", ionCode); } - // anti hyper hydrogen 4 ground state + //anti hyper hydrogen 4 ground state ionCode = -1010010040; if (!db->GetParticle(ionCode)) { db->AddParticle("AntiHyperhydrog4", "AntiHyperhydrog4", 3.922434, kFALSE, 2.5e-15, 3, "Ion", ionCode); } - // hyper helium 4 ground state + //hyper helium 4 ground state ionCode = 1010020040; if (!db->GetParticle(ionCode)) { db->AddParticle("Hyperhelium4", "Hyperhelium4", 3.921728, kFALSE, 2.5e-15, 6, "Ion", ionCode); } - // anti hyper helium 4 ground state + //anti hyper helium 4 ground state ionCode = -1010020040; if (!db->GetParticle(ionCode)) { db->AddParticle("AntiHyperhelium4", "AntiHyperhelium4", 3.921728, kFALSE, @@ -676,9 +642,9 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa 0.040, -3, "Resonance", -ionCode); } - // Additional Hidden Strangeness Pentaquarks + //Additional Hidden Strangeness Pentaquarks - // Ps +/- + //Ps +/- ionCode = 9322132; if (!db->GetParticle(ionCode)) { db->AddParticle("Ps_Plus_1870", "Ps_Plus_1870", 1.870, kFALSE, @@ -716,7 +682,7 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa 0.10, -3, "Resonance", -ionCode); } - // Ps0 + //Ps0 ionCode = 9322131; if (!db->GetParticle(ionCode)) { db->AddParticle("Ps_0_1870", "Ps_0_1870", 1.870, kFALSE, @@ -820,14 +786,14 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMa // d*(2380) - dibaryon resonance - ionCode = 900010020; // Arbitrary choice - as deuteron but with leading 9 instead of 10 + ionCode = 900010020; //Arbitrary choice - as deuteron but with leading 9 instead of 10 if (!db->GetParticle(ionCode)) { db->AddParticle("d*_2380", "d*_2380", 2.38, kFALSE, 0.070, 3, "Resonance", ionCode); } db->AddAntiParticle("d*_2380_bar", -ionCode); - // Sexaquark (uuddss): compact, neutral and stable hypothetical bound state (arxiv.org/abs/1708.08951) + //Sexaquark (uuddss): compact, neutral and stable hypothetical bound state (arxiv.org/abs/1708.08951) ionCode = 900000020; if (!db->GetParticle(ionCode)) { db->AddParticle("Sexaquark", "Sexaquark", 2.0, kTRUE, 0.0, 0, "Special", ionCode); @@ -860,4 +826,4 @@ inline void O2DatabasePDG::addParticlesFromExternalFile(TDatabasePDG* db) } // namespace o2 -#endif // O2_O2DATABASEPDG_H +#endif //O2_O2DATABASEPDG_H diff --git a/Steer/src/MCKinematicsReader.cxx b/Steer/src/MCKinematicsReader.cxx index 08a400354494e..21024dba78368 100644 --- a/Steer/src/MCKinematicsReader.cxx +++ b/Steer/src/MCKinematicsReader.cxx @@ -13,8 +13,6 @@ #include "Steer/MCKinematicsReader.h" #include "SimulationDataFormat/MCEventHeader.h" #include "SimulationDataFormat/TrackReference.h" -#include "SimulationDataFormat/MonopoleParticles.h" -#include "SimulationDataFormat/O2DatabasePDG.h" #include #include #include @@ -84,29 +82,6 @@ void MCKinematicsReader::initTracksForSource(int source) const } } -namespace -{ -// Makes the monopole mass used by the simulation, which is stored in the event header, the -// one of O2DatabasePDG; otherwise MCTrack::GetMass() (and GetEnergy(), GetRapidity(), ...) would -// give the default mass for monopoles of a run with another G4.monopoleMass. -void applyMonopoleMass(o2::dataformats::MCEventHeader const& header) -{ - if (!header.hasInfo(o2::sim::MonopoleMassHeaderKey)) { - return; - } - bool valid = false; - const double mass = header.getInfo(o2::sim::MonopoleMassHeaderKey, valid); - if (!valid || !(mass > 0.)) { - return; - } - if (o2::O2DatabasePDG::hasMonopoleMassOverride() && std::abs(o2::O2DatabasePDG::getMonopoleMass() - mass) > 1.e-6 * mass) { - LOG(warn) << "Kinematics with different monopole masses (" << o2::O2DatabasePDG::getMonopoleMass() << " and " << mass - << " GeV) are read in the same process; monopole masses will be wrong for one of them"; - } - o2::O2DatabasePDG::setMonopoleMass(mass); -} -} // namespace - void MCKinematicsReader::loadTracksForSourceAndEvent(int source, int event) const { auto chain = mInputChains[source]; @@ -120,17 +95,6 @@ void MCKinematicsReader::loadTracksForSourceAndEvent(int source, int event) cons mTracks[source][event] = new std::vector; *mTracks[source][event] = *loadtracks; delete loadtracks; - // the monopole mass has to be known before the tracks are used - if (auto hbr = chain->GetBranch("MCEventHeader.")) { - o2::dataformats::MCEventHeader* header = nullptr; - hbr->SetAddress(&header); - hbr->GetEntry(event); - if (header) { - applyMonopoleMass(*header); - } - delete header; - hbr->ResetAddress(); - } } } } @@ -156,7 +120,6 @@ void MCKinematicsReader::loadHeadersForSource(int source) const for (int event = 0; event < br->GetEntries(); ++event) { br->GetEntry(event); mHeaders[source][event] = *header; - applyMonopoleMass(mHeaders[source][event]); } delete header; header = nullptr; diff --git a/run/O2HitMerger.cxx b/run/O2HitMerger.cxx index fc9739c75b976..b79c18bd5ca7a 100644 --- a/run/O2HitMerger.cxx +++ b/run/O2HitMerger.cxx @@ -41,9 +41,6 @@ #include "PrimaryServerState.h" #include -#include -#include -#include #include #include #include @@ -150,23 +147,6 @@ void O2HitMerger::InitTask() // we didn't manage to get a configuration --> better to fail LOG(fatal) << "No configuration received. Aborting"; } - // The config key-values are not applied to the ConfigurableParams in this process. The merger needs the G4 ones - // to record the monopole mass in the event headers (G4.monopole, G4.monopoleMass), so apply those. - // Note: G4 parameters given via a config file (--configFile) are not applied here. - { - std::string g4Tokens; - std::stringstream tokens(o2::conf::SimConfig::Instance().getKeyValueString()); - std::string token; - while (std::getline(tokens, token, ';')) { - const auto first = token.find_first_not_of(" \t"); - if (first != std::string::npos && token.compare(first, 3, "G4.") == 0) { - g4Tokens += token.substr(first) + ";"; - } - } - if (!g4Tokens.empty()) { - o2::conf::ConfigurableParam::updateFromString(g4Tokens); - } - } mAsService = o2::conf::SimConfig::Instance().asService(); mForwardKine = o2::conf::SimConfig::Instance().forwardKine(); mWriteToDisc = o2::conf::SimConfig::Instance().writeToDisc(); @@ -714,12 +694,6 @@ bool O2HitMerger::mergeAndFlushData() } } - // record the mass used for the monopoles, such that readers of the kinematics - // (MCTrack::GetMass() and derived quantities) do not fall back to the default - if (eventheader && o2::conf::G4Params::Instance().monopole) { - eventheader->putInfo(o2::sim::MonopoleMassHeaderKey, o2::conf::G4Params::Instance().monopoleMass); - } - // now see which events can be discarded in any case due to no hits if (confref.isFilterOutNoHitEvents()) { if (eventheader && eventheader->getMCEventStats().getNHits() == 0) { From 2f1280379d9f91f169a3db7fe2c0e48ba6911669 Mon Sep 17 00:00:00 2001 From: tubagundem Date: Thu, 15 Oct 2026 15:05:21 +0200 Subject: [PATCH 11/11] Extending track references --- .../include/CPVSimulation/Detector.h | 1 + Detectors/CPV/simulation/src/Detector.cxx | 9 + .../include/FDDSimulation/Detector.h | 1 + Detectors/FIT/FDD/simulation/src/Detector.cxx | 15 +- .../include/FT0Simulation/Detector.h | 5 +- Detectors/FIT/FT0/simulation/src/Detector.cxx | 225 +++++++++--------- .../include/FV0Simulation/Detector.h | 1 + Detectors/FIT/FV0/simulation/src/Detector.cxx | 9 + .../include/MFTSimulation/Detector.h | 1 + .../ITSMFT/MFT/simulation/src/Detector.cxx | 9 + Detectors/PHOS/simulation/src/Detector.cxx | 1 + .../include/TOFSimulation/Detector.h | 5 +- Detectors/TOF/simulation/src/Detector.cxx | 11 +- 13 files changed, 177 insertions(+), 116 deletions(-) diff --git a/Detectors/CPV/simulation/include/CPVSimulation/Detector.h b/Detectors/CPV/simulation/include/CPVSimulation/Detector.h index c9ece3ba482ec..3bd41f89c7ea2 100644 --- a/Detectors/CPV/simulation/include/CPVSimulation/Detector.h +++ b/Detectors/CPV/simulation/include/CPVSimulation/Detector.h @@ -157,6 +157,7 @@ class Detector : public o2::base::DetImpl // Simulation std::vector* mHits = nullptr; //! Collection of CPV hits + int mTrackRefLastTrackID = -1; //! last track that received a TrackReference in the current event template friend class o2::base::DetImpl; diff --git a/Detectors/CPV/simulation/src/Detector.cxx b/Detectors/CPV/simulation/src/Detector.cxx index ee6ceef3de1e3..c2f07a157ceec 100644 --- a/Detectors/CPV/simulation/src/Detector.cxx +++ b/Detectors/CPV/simulation/src/Detector.cxx @@ -30,6 +30,7 @@ #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Stack.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include #include @@ -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(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); diff --git a/Detectors/FIT/FDD/simulation/include/FDDSimulation/Detector.h b/Detectors/FIT/FDD/simulation/include/FDDSimulation/Detector.h index 485197b30d0bb..363c05dd4fec4 100644 --- a/Detectors/FIT/FDD/simulation/include/FDDSimulation/Detector.h +++ b/Detectors/FIT/FDD/simulation/include/FDDSimulation/Detector.h @@ -71,6 +71,7 @@ class Detector : public o2::base::DetImpl /// Container for data points std::vector* mHits = nullptr; + int mTrackRefLastTrackID = -1; //! last track that received a TrackReference in the current event Detector& operator=(const Detector&); diff --git a/Detectors/FIT/FDD/simulation/src/Detector.cxx b/Detectors/FIT/FDD/simulation/src/Detector.cxx index 3bb945ddd1bba..14e92379d433b 100644 --- a/Detectors/FIT/FDD/simulation/src/Detector.cxx +++ b/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(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 "<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; } diff --git a/Detectors/FIT/FT0/simulation/include/FT0Simulation/Detector.h b/Detectors/FIT/FT0/simulation/include/FT0Simulation/Detector.h index ebfe6e4a1ba3a..a35af48ba7334 100644 --- a/Detectors/FIT/FT0/simulation/include/FT0Simulation/Detector.h +++ b/Detectors/FIT/FT0/simulation/include/FT0Simulation/Detector.h @@ -149,8 +149,8 @@ class Detector : public o2::base::DetImpl std::vector 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 /// Container for data points std::vector* mHits = nullptr; + int mTrackRefLastTrackID = -1; //! last track that received a TrackReference in the current event // Define volume IDs int mREGVolID = -1; diff --git a/Detectors/FIT/FT0/simulation/src/Detector.cxx b/Detectors/FIT/FT0/simulation/src/Detector.cxx index 3b7bff8d51633..134dd252971a2 100644 --- a/Detectors/FIT/FT0/simulation/src/Detector.cxx +++ b/Detectors/FIT/FT0/simulation/src/Detector.cxx @@ -37,6 +37,7 @@ #include "FT0Base/Geometry.h" #include "FT0Simulation/Detector.h" #include "DetectorsBase/Stack.h" +#include "SimulationDataFormat/MonopoleParticles.h" using namespace o2::ft0; using o2::ft0::Geometry; @@ -134,7 +135,7 @@ void Detector::ConstructGeometry() TGeoCombiTrans* comCable[nCellsC]; TString nameCom; - //Additional elements for the C-side frame + // Additional elements for the C-side frame TGeoCombiTrans* plateCom[nCellsC]; TGeoMedium* Al = gGeoManager->GetMedium("FT0_Aluminium$"); TGeoCompositeShape* plateCompositeShape = new TGeoCompositeShape("plateCompositeShape", cPlateShapeString().c_str()); @@ -161,26 +162,26 @@ void Detector::ConstructGeometry() TGeoVolume* cables = SetCablesSize(itr); LOG(debug) << " C " << mPosModuleCx[ic] << " " << mPosModuleCy[ic]; // cables->Print(); - //Additional shift (+0.1) introduced to cable planes so they don't overlap the C-side frame + // Additional shift (+0.1) introduced to cable planes so they don't overlap the C-side frame comCable[ic] = new TGeoCombiTrans(mPosModuleCx[ic], mPosModuleCy[ic], mPosModuleCz[ic] + mInStart[2] + 0.2 + 0.1, rot[ic]); TGeoHMatrix hmCable = *comCable[ic]; TGeoHMatrix* phCable = new TGeoHMatrix(hmCable); stlinC->AddNode(cables, itr, comCable[ic]); - //C-side frame elements - module plates + // C-side frame elements - module plates plateCom[ic] = new TGeoCombiTrans(mPosModuleCx[ic], mPosModuleCy[ic], (mPosModuleCz[ic] - 3), rot[ic]); TGeoHMatrix hmPlate = *plateCom[ic]; TGeoHMatrix* phPlate = new TGeoHMatrix(hmPlate); stlinC->AddNode(plateVol, itr, phPlate); } - //Add C-side frame + // Add C-side frame stlinC->AddNode(constructFrameCGeometry(), nCellsA + nCellsC + 1); TGeoVolume* alice = gGeoManager->GetVolume("barrel"); - //Add A-side detector - alice->AddNode(stlinA, 1, new TGeoTranslation(0, 30., zdetA + 0.63)); //offset to avoid overlap with FV0 + // Add A-side detector + alice->AddNode(stlinA, 1, new TGeoTranslation(0, 30., zdetA + 0.63)); // offset to avoid overlap with FV0 - //Add C-side detector + // Add C-side detector TGeoRotation* rotC = new TGeoRotation("rotC", 90., 0., 90., 90., 180., 0.); alice->AddNode(stlinC, 1, new TGeoCombiTrans(0., 30., -zdetC, rotC)); @@ -214,7 +215,7 @@ void Detector::SetOneMCP(TGeoVolume* ins) const double kWrapClearance = 1.e-4; // cm Float_t ptopref[3] = {static_cast(ptop[0] + 2 * prfv[0] + kWrapClearance), static_cast(ptop[1] + 2 * prfh[1] + kWrapClearance), 1.}; - Float_t pmcp[3] = {2.949, 2.949, 0.66}; // MCP + Float_t pmcp[3] = {2.949, 2.949, 0.66}; // MCP Float_t pmcpinner[3] = {2.749, 2.749, 0.1}; Float_t pmcpbase[3] = {2.949, 2.949, 0.675}; Float_t pmcpside[3] = {0.15, 2.949, 0.65}; @@ -431,18 +432,18 @@ void Detector::addAlignableVolumes() const } } -//Construction of FT0-A support structure -//The frame is constructed by first building a block of Aluminum as a mother volume from which details can -//be subtracted. First, 6 boxe shapes are subtracted from around the edges and 2 from the center of the frame -//to create the fin shapes at the edges of the detector and the cross shape at the center of the detector. -//These boxe shapes are then subtracted again but from positions reflected in both x and y which is -//reflects the symmetry of the frame. Then a loop is used to subtract out the PMT sockets which are also -//box shapes from the positions given in the Geometry.cxx file. In the loop, after a socket is subtracted, -//either an inner or an outer plate group is placed inside of it. Inner and outer plate groups are -//both composed of a cover plate and a cable plate with fiber heads subtracted from the cable plate. -//The orientation of these holes differs between the inner and outer plates making them distinct. -//Contributors: Joe Crowley (2019-20), Jason Pruitt (2020-21), Sam Christensen (2021-22), -//and Jennifer Klay (2019-22) from Cal Poly SLO. +// Construction of FT0-A support structure +// The frame is constructed by first building a block of Aluminum as a mother volume from which details can +// be subtracted. First, 6 boxe shapes are subtracted from around the edges and 2 from the center of the frame +// to create the fin shapes at the edges of the detector and the cross shape at the center of the detector. +// These boxe shapes are then subtracted again but from positions reflected in both x and y which is +// reflects the symmetry of the frame. Then a loop is used to subtract out the PMT sockets which are also +// box shapes from the positions given in the Geometry.cxx file. In the loop, after a socket is subtracted, +// either an inner or an outer plate group is placed inside of it. Inner and outer plate groups are +// both composed of a cover plate and a cable plate with fiber heads subtracted from the cable plate. +// The orientation of these holes differs between the inner and outer plates making them distinct. +// Contributors: Joe Crowley (2019-20), Jason Pruitt (2020-21), Sam Christensen (2021-22), +// and Jennifer Klay (2019-22) from Cal Poly SLO. TGeoVolume* Detector::constructFrameAGeometry() { // define the media @@ -452,20 +453,20 @@ TGeoVolume* Detector::constructFrameAGeometry() // make a volume assembly for the frame TGeoVolumeAssembly* FT0_Frame = new TGeoVolumeAssembly("FT0_Frame"); - //Define the block of aluminum that forms the frame - Double_t blockdX = 37.1; //slightly larger in x - Double_t blockdY = 36.85; //than y - Double_t blockdZ = 6.95; //thickness of frame and back plates + // Define the block of aluminum that forms the frame + Double_t blockdX = 37.1; // slightly larger in x + Double_t blockdY = 36.85; // than y + Double_t blockdZ = 6.95; // thickness of frame and back plates TGeoBBox* block = new TGeoBBox("block", blockdX / 2, blockdY / 2, blockdZ / 2); - //To form the outer frame shape with fins that bolt it to the FV0, remove - //aluminum in six chunks (boxes) from two sides, then reflect these to remove - //from the other sides. As viewed from the back side of the detector, count - //clockwise from bottom left for numbering. - Double_t box1dX = 1.57; //narrower - Double_t box1dY = 6.55; //than it is tall - Double_t box1PosX = -(blockdX / 2 - box1dX / 2); //placement on the frame block - Double_t box1PosY = 0; //starts at the middle + // To form the outer frame shape with fins that bolt it to the FV0, remove + // aluminum in six chunks (boxes) from two sides, then reflect these to remove + // from the other sides. As viewed from the back side of the detector, count + // clockwise from bottom left for numbering. + Double_t box1dX = 1.57; // narrower + Double_t box1dY = 6.55; // than it is tall + Double_t box1PosX = -(blockdX / 2 - box1dX / 2); // placement on the frame block + Double_t box1PosY = 0; // starts at the middle TGeoBBox* box1 = new TGeoBBox("box1", box1dX / 2, box1dY / 2, blockdZ / 2); TGeoTranslation* box1Tr1 = new TGeoTranslation("box1Tr1", box1PosX, box1PosY, 0); box1Tr1->RegisterYourself(); @@ -482,14 +483,14 @@ TGeoVolume* Detector::constructFrameAGeometry() TGeoTranslation* box2Tr2 = new TGeoTranslation("box2Tr2", -box2PosX, -box2PosY, 0); box2Tr2->RegisterYourself(); - //Box 3 is shallower than the others to preserve the aluminum fin where the - //FT0 is bolted to the FV0 + // Box 3 is shallower than the others to preserve the aluminum fin where the + // FT0 is bolted to the FV0 Double_t box3dX = 12.7; Double_t box3dY = 3; Double_t box3dZ = 5.45; Double_t box3PosX = -(blockdX / 2 - box2dX - box3dZ / 2); Double_t box3PosY = blockdY / 2 - box3dY / 2; - Double_t box3PosZ = blockdZ / 2 - box3dZ / 2; //subtract from the back, leaving fin on the front + Double_t box3PosZ = blockdZ / 2 - box3dZ / 2; // subtract from the back, leaving fin on the front TGeoBBox* box3 = new TGeoBBox("box3", box3dX / 2, box3dY / 2, box3dZ / 2); TGeoTranslation* box3Tr1 = new TGeoTranslation("box3Tr1", box3PosX, box3PosY, box3PosZ); box3Tr1->RegisterYourself(); @@ -516,25 +517,25 @@ TGeoVolume* Detector::constructFrameAGeometry() TGeoTranslation* box5Tr2 = new TGeoTranslation("box5Tr2", -box5PosX, -box5PosY, 0); box5Tr2->RegisterYourself(); - //Similar to box 3, box 6 is shallower in z to leave aluminum for the fin that - //bolts FT0 to FV0 + // Similar to box 3, box 6 is shallower in z to leave aluminum for the fin that + // bolts FT0 to FV0 Double_t box6dX = 2.9; Double_t box6dY = 12.2; Double_t box6dZ = 5.45; Double_t box6PosX = blockdX / 2 - box6dX / 2; Double_t box6PosY = blockdY / 2 - box5dY - box6dY / 2; - Double_t box6PosZ = blockdZ / 2 - box6dZ / 2; //subtract from the back, leaving fin at the front + Double_t box6PosZ = blockdZ / 2 - box6dZ / 2; // subtract from the back, leaving fin at the front TGeoBBox* box6 = new TGeoBBox("box6", box6dX / 2, box6dY / 2, box6dZ / 2); TGeoTranslation* box6Tr1 = new TGeoTranslation("box6Tr1", box6PosX, box6PosY, box6PosZ); box6Tr1->RegisterYourself(); TGeoTranslation* box6Tr2 = new TGeoTranslation("box6Tr2", -box6PosX, -box6PosY, box6PosZ); box6Tr2->RegisterYourself(); - //The central hole that accommodates the beam pipe is not the same on all four sides - //so we define two rectangular boxes - one vertical and one horizontal - and copy/rotate them - //to remove the aluminum in a "+" shape at the center - // cbox is a central rectangle - Double_t cbox1dX = 7.175; //horizontal center box + // The central hole that accommodates the beam pipe is not the same on all four sides + // so we define two rectangular boxes - one vertical and one horizontal - and copy/rotate them + // to remove the aluminum in a "+" shape at the center + // cbox is a central rectangle + Double_t cbox1dX = 7.175; // horizontal center box Double_t cbox1dY = 5.5; Double_t cbox1Xoffset = 14.425; Double_t cbox1PosX = -(blockdX / 2 - cbox1Xoffset - cbox1dX / 2); @@ -545,7 +546,7 @@ TGeoVolume* Detector::constructFrameAGeometry() TGeoTranslation* cbox1Tr2 = new TGeoTranslation("cbox1Tr2", -cbox1PosX, -cbox1PosY, 0); cbox1Tr2->RegisterYourself(); - Double_t cbox2dX = 5.75; //vertical center box + Double_t cbox2dX = 5.75; // vertical center box Double_t cbox2dY = 6.575; Double_t cbox2Yoffset = 14.425; @@ -557,11 +558,11 @@ TGeoVolume* Detector::constructFrameAGeometry() TGeoTranslation* cbox2Tr2 = new TGeoTranslation("cbox2Tr2", -cbox2PosX, -cbox2PosY, 0); cbox2Tr2->RegisterYourself(); - //The two L-shaped pieces that form the frame have a small 1mm gap between them, - //where they come together. As viewed from the back, the gaps are on the upper - //left and lower right, so that for the center column of modules, the upper two - //are shifted slightly to the right (as viewed from the back) and the lower two - //are shifted slightly to the left (as viewed from the back) + // The two L-shaped pieces that form the frame have a small 1mm gap between them, + // where they come together. As viewed from the back, the gaps are on the upper + // left and lower right, so that for the center column of modules, the upper two + // are shifted slightly to the right (as viewed from the back) and the lower two + // are shifted slightly to the left (as viewed from the back) Double_t gapBoxdX = 0.1; Double_t gapBoxdY = blockdY / 2; Double_t gapPosX = -(sPmtSide / 2 + sEps + gapBoxdX / 2); @@ -572,8 +573,8 @@ TGeoVolume* Detector::constructFrameAGeometry() TGeoTranslation* gapBoxTr2 = new TGeoTranslation("gapBoxTr2", -gapPosX, -gapPosY, 0); gapBoxTr2->RegisterYourself(); - //Create a string to define the complete frame object shape - //Start from the aluminum block then subtract the boxes + // Create a string to define the complete frame object shape + // Start from the aluminum block then subtract the boxes std::string frameACompositeString = "block "; frameACompositeString += "- box1:box1Tr1 - box1:box1Tr2 "; frameACompositeString += "- box2:box2Tr1 - box2:box2Tr2 "; @@ -585,33 +586,33 @@ TGeoVolume* Detector::constructFrameAGeometry() frameACompositeString += "- cbox2:cbox2Tr1 - cbox2:cbox2Tr2 "; frameACompositeString += "- gapBox:gapBoxTr1 - gapBox:gapBoxTr2"; - //The next section defines the objects that form the sockets in the - //frame for the sensitive elements and the individual cover plates - //at the front of the detector which include the optical fiber - //heads that permit the LED pulser light to reach the quartz radiator - //surfaces of each module + // The next section defines the objects that form the sockets in the + // frame for the sensitive elements and the individual cover plates + // at the front of the detector which include the optical fiber + // heads that permit the LED pulser light to reach the quartz radiator + // surfaces of each module - //There are two fiber head configurations, called "inner" and "outer" - //with different locations and angles of the fiber heads. - Double_t coverPlatedZ = 0.2; //top cover thickness - Double_t fiberPlatedZ = 0.5; //fiberhead plate is underneath + // There are two fiber head configurations, called "inner" and "outer" + // with different locations and angles of the fiber heads. + Double_t coverPlatedZ = 0.2; // top cover thickness + Double_t fiberPlatedZ = 0.5; // fiberhead plate is underneath - //Each fiber is guided to a small rectangular opening in the plate - Double_t opticalFiberHeaddY = 0.52; //narrow side - Double_t opticalFiberHeaddX = 1.142; //long side + // Each fiber is guided to a small rectangular opening in the plate + Double_t opticalFiberHeaddY = 0.52; // narrow side + Double_t opticalFiberHeaddX = 1.142; // long side - //The "top" two fiber heads are positioned at slightly different - //locations than the bottom two, which are also rotated + // The "top" two fiber heads are positioned at slightly different + // locations than the bottom two, which are also rotated //"Outer" fiberhead placements Double_t fh1TopPosX = -1.6; Double_t fh1TopPosY = 1.325; Double_t fh1BotPosX = 1.555; Double_t fh1BotPosY = 1.249; - Double_t fh1BotAngle = 16; //degrees + Double_t fh1BotAngle = 16; // degrees //"Inner" fiberhead placements - //All of these are placed at an angle + // All of these are placed at an angle Double_t fh2TopPosX = -1.563; Double_t fh2TopPosY = 1.4625; Double_t fh2TopAngle = 60; @@ -645,7 +646,7 @@ TGeoVolume* Detector::constructFrameAGeometry() TGeoCombiTrans* fh2BotTr2 = new TGeoCombiTrans("fh2BotTr2", fh2BotPosX, -fh2BotPosY, 0, new TGeoRotation("fh2BotRot2", fh2BotAngle, 0, 0)); fh2BotTr2->RegisterYourself(); - //Create a string that defines the plate group for the outer plates + // Create a string that defines the plate group for the outer plates std::string outerPlateGroupString = "fiberPlate "; outerPlateGroupString += "- opticalFiberHead:fh1TopTr1 "; outerPlateGroupString += "- opticalFiberHead:fh1TopTr2 "; @@ -653,10 +654,10 @@ TGeoVolume* Detector::constructFrameAGeometry() outerPlateGroupString += "- opticalFiberHead:fh1BotTr2 "; outerPlateGroupString += "+ coverPlate:coverPlateTr"; - //Create the composite shape for the outer plates + // Create the composite shape for the outer plates TGeoCompositeShape* outerPlateGroup = new TGeoCompositeShape("outerPlateGroup", outerPlateGroupString.c_str()); - //Create a string that defines the plate group for the inner plates + // Create a string that defines the plate group for the inner plates std::string innerPlateGroupString = "fiberPlate "; innerPlateGroupString += "- opticalFiberHead:fh2TopTr1 "; innerPlateGroupString += "- opticalFiberHead:fh2TopTr2 "; @@ -664,38 +665,38 @@ TGeoVolume* Detector::constructFrameAGeometry() innerPlateGroupString += "- opticalFiberHead:fh2BotTr2 "; innerPlateGroupString += "+ coverPlate:coverPlateTr"; - //Create the composite shape for the inner plates + // Create the composite shape for the inner plates TGeoCompositeShape* innerPlateGroup = new TGeoCompositeShape("innerPlateGroup", innerPlateGroupString.c_str()); - //The sockets that are cut out of the aluminum block for the senitive elements - //to fit into are offset slightly in z to leave a thin plate of aluminum at the - //back - the back plate covers + // The sockets that are cut out of the aluminum block for the senitive elements + // to fit into are offset slightly in z to leave a thin plate of aluminum at the + // back - the back plate covers Double_t backPlanedZ = 0.25; Double_t socketdZ = blockdZ - backPlanedZ; - //Define the socket volume as a box of vacuum + // Define the socket volume as a box of vacuum TGeoVolume* socket = gGeoManager->MakeBox("Socket", Vacuum, sPmtSide / 2 + sEps, sPmtSide / 2 + sEps, socketdZ / 2); - //Define the orientation angles of the plate groups that will cover - //the sockets holding the sensitive elements + // Define the orientation angles of the plate groups that will cover + // the sockets holding the sensitive elements Double_t rotAngle[Geometry::NCellsA] = {0, 0, -90, -90, -90, 0, 0, -90, -90, -90, 0, 0, 180, 180, 90, 90, 90, 180, 180, 90, 90, 90, 180, 180}; - //Define the socket and plate group translations + // Define the socket and plate group translations TGeoTranslation* trSocket[Geometry::NCellsA]; TString nameTrSocket; TGeoCombiTrans* trPlateGroup[Geometry::NCellsA]; TString nameTrPlateGroup; TString namePGRot; - //Loop over the number of modules, subtracting the sockets and adding back in the - //plate groups at the position of each module + // Loop over the number of modules, subtracting the sockets and adding back in the + // plate groups at the position of each module for (Int_t itr = 0; itr < Geometry::NCellsA; itr++) { nameTrSocket = Form("trSocket%i", itr + 1); float z = -backPlanedZ / 4.0; trSocket[itr] = new TGeoTranslation(nameTrSocket.Data(), mPosModuleAx[itr], mPosModuleAy[itr], z); trSocket[itr]->RegisterYourself(); - frameACompositeString += "- Socket:"; //subtract it from the aluminum block - frameACompositeString += nameTrSocket.Data(); //at its corresponding location + frameACompositeString += "- Socket:"; // subtract it from the aluminum block + frameACompositeString += nameTrSocket.Data(); // at its corresponding location nameTrPlateGroup = Form("trPlateGroup%i", itr + 1); namePGRot = Form("pgRot%i", itr + 1); @@ -704,29 +705,29 @@ TGeoVolume* Detector::constructFrameAGeometry() trPlateGroup[itr]->RegisterYourself(); if (itr == 0 || itr == 2 || itr == 3 || itr == 4 || itr == 5 || itr == 10 || itr == 13 || itr == 18 || itr == 19 || itr == 20 || itr == 21 || itr == 23) { - frameACompositeString += " + outerPlateGroup:"; //add the outer plate group back on to these modules + frameACompositeString += " + outerPlateGroup:"; // add the outer plate group back on to these modules frameACompositeString += nameTrPlateGroup.Data(); frameACompositeString += " "; } else { - frameACompositeString += " + innerPlateGroup:"; //or add the inner plate group back on to all other modules + frameACompositeString += " + innerPlateGroup:"; // or add the inner plate group back on to all other modules frameACompositeString += nameTrPlateGroup.Data(); frameACompositeString += " "; } } - //Finally, define the A side frame object from the complete composite shape defined above + // Finally, define the A side frame object from the complete composite shape defined above TGeoVolume* frameA = new TGeoVolume("frameA", new TGeoCompositeShape("frameA", frameACompositeString.c_str()), Al); - //Add the frame object to the mother volume + // Add the frame object to the mother volume FT0_Frame->AddNode(frameA, 1); return FT0_Frame; } -//C-side Support Structure -//This code was written by Jason Pruitt and Sam Christensen of Cal Poly in 2021 -//They followed a similar method as for the A-side frame but had to account for the -//spherical geometry of the C-side. +// C-side Support Structure +// This code was written by Jason Pruitt and Sam Christensen of Cal Poly in 2021 +// They followed a similar method as for the A-side frame but had to account for the +// spherical geometry of the C-side. TGeoVolume* Detector::constructFrameCGeometry() { @@ -842,7 +843,7 @@ TGeoVolume* Detector::constructFrameCGeometry() TGeoCombiTrans* comCorners2 = new TGeoCombiTrans("comCorners2", xCorn2, yCorn2, zCorn2, rotCorners); comCorners2->RegisterYourself(); - //Create a string that defines the composite shape + // Create a string that defines the composite shape std::string shellString = ""; shellString += "sphere1"; // start with spherical shell - this will be reflected shellString += "- sphere2:rotTr1"; // copy and combitrans a subtraction @@ -852,7 +853,7 @@ TGeoVolume* Detector::constructFrameCGeometry() shellString += "- framecornerBox:comCorners"; // subtract corners shellString += "- framecornerBox:comCorners2"; // - //Create string that defines the back plate composite shape + // Create string that defines the back plate composite shape std::string backPlateString = ""; backPlateString += "sphere3"; backPlateString += "- sphere4:rotTr1"; @@ -862,8 +863,8 @@ TGeoVolume* Detector::constructFrameCGeometry() backPlateString += "- framecornerBox:comCorners"; backPlateString += "- framecornerBox:comCorners2"; - //These could be set up to use the values in the geometry file after some - //investigation of subtle differences... + // These could be set up to use the values in the geometry file after some + // investigation of subtle differences... static constexpr Double_t xi[NCellsC] = {-15.038271418735729, 15.038271418735729, -15.003757581112167, 15.003757581112167, -9.02690018974363, 9.02690018974363, -9.026897413747076, 9.026897413747076, -9.026896531935773, 9.026896531935773, -3.0004568618531313, 3.0004568618531313, -3.0270795197907225, 3.0270795197907225, 3.0003978432927543, -3.0003978432927543, 3.0270569670429572, -3.0270569670429572, 9.026750365564254, -9.026750365564254, 9.026837450695885, -9.026837450695885, 9.026849243816981, -9.026849243816981, 15.038129472387304, -15.038129472387304, 15.003621961057961, -15.003621961057961}; static constexpr Double_t yi[NCellsC] = {3.1599494336464455, -3.1599494336464455, 9.165191680982874, -9.165191680982874, 3.1383331772537426, -3.1383331772537426, 9.165226363918643, -9.165226363918643, 15.141616002932361, -15.141616002932361, 9.16517861649866, -9.16517861649866, 15.188854859073416, -15.188854859073416, 9.165053319552113, -9.165053319552113, 15.188703787345304, -15.188703787345304, 3.138263189805292, -3.138263189805292, 9.165104089644917, -9.165104089644917, 15.141494417823818, -15.141494417823818, 3.1599158563428644, -3.1599158563428644, 9.165116302773846, -9.165116302773846}; @@ -1027,31 +1028,31 @@ TGeoVolume* Detector::constructFrameCGeometry() std::string Detector::cPlateShapeString() { - Double_t prismHeight = 0.3895; //height of vertical edge of square prism part of base - Double_t prismSide = 5.9; //width and length of square prism part of base - Double_t radCurve = 81.9469; //radius of curvature of top part of base + Double_t prismHeight = 0.3895; // height of vertical edge of square prism part of base + Double_t prismSide = 5.9; // width and length of square prism part of base + Double_t radCurve = 81.9469; // radius of curvature of top part of base Double_t delHeight = radCurve * (1.0 - TMath::Sqrt(1.0 - 0.5 * TMath::Power(prismSide / radCurve, 2.0))); - //height from top of square prism to center of curved top surface of base + // height from top of square prism to center of curved top surface of base - Double_t heightBase = prismHeight + delHeight; //from center of bottom to center of top - Double_t sliceSide = 5.3; //side lengths of slice's flat top + Double_t heightBase = prismHeight + delHeight; // from center of bottom to center of top + Double_t sliceSide = 5.3; // side lengths of slice's flat top Double_t heightBaseBox = 2 * heightBase; Double_t totalHeight = 0.5; Double_t sliceHeight = 0.5 - heightBase; - //cable dimensions and distances + // cable dimensions and distances Double_t cableHoleWidth = 0.3503; Double_t cableHoleLength = 0.9003; - Double_t cableHoleDepth = 1; //really big just to punch a hole + Double_t cableHoleDepth = 1; // really big just to punch a hole - //sholes denotes "straight holes" and rholes denote "rotated holes" - //all distances measured from edges of slice - //up and down sholes + // sholes denotes "straight holes" and rholes denote "rotated holes" + // all distances measured from edges of slice + // up and down sholes Double_t sHolesBottomEdge = 1.585; Double_t sHolesTopEdge = 0.515; Double_t sHolesAvgTopBottom = (sHolesBottomEdge + sHolesTopEdge) / 2.0; - Double_t sHolesUpFromCenter = ((sliceSide / 2.0) - sHolesAvgTopBottom); //amount up in x the sholes need to move - //left and right sholes + Double_t sHolesUpFromCenter = ((sliceSide / 2.0) - sHolesAvgTopBottom); // amount up in x the sholes need to move + // left and right sholes Double_t sHolesFarEdge = 1.585; Double_t sHolesNearEdge = 1.065; Double_t sHolesAvgNearFar = (sHolesFarEdge + sHolesNearEdge) / 2.0; @@ -1115,7 +1116,7 @@ std::string Detector::cPlateShapeString() return plateString; } -//End Support structure code +// End Support structure code //////////////////////////////////////////// Bool_t Detector::ProcessHits(FairVolume* v) @@ -1123,6 +1124,13 @@ Bool_t Detector::ProcessHits(FairVolume* v) Int_t copy; Int_t volID = fMC->CurrentVolID(copy); + if (fMC->IsTrackEntering() && !fMC->IsNewTrack() && (fMC->TrackCharge() != 0 || o2::sim::isMonopole(fMC->TrackPid())) && volID == mTOPVolID) { // one reference per charged particle (or monopole), first entry from outside + auto o2stack = static_cast(fMC->GetStack()); + if (o2stack->GetCurrentTrackNumber() != mTrackRefLastTrackID) { + mTrackRefLastTrackID = o2stack->GetCurrentTrackNumber(); + o2stack->addTrackReference(o2::TrackReference(*fMC, GetDetId())); + } + } TVirtualMCStack* stack = fMC->GetStack(); Int_t quadrant, mcp; @@ -1203,6 +1211,7 @@ void Detector::Reset() if (!o2::utils::ShmManager::Instance().isOperational()) { mHits->clear(); } + mTrackRefLastTrackID = -1; } void Detector::CreateMaterials() diff --git a/Detectors/FIT/FV0/simulation/include/FV0Simulation/Detector.h b/Detectors/FIT/FV0/simulation/include/FV0Simulation/Detector.h index f9126061bf214..fdff96b6eafa8 100644 --- a/Detectors/FIT/FV0/simulation/include/FV0Simulation/Detector.h +++ b/Detectors/FIT/FV0/simulation/include/FV0Simulation/Detector.h @@ -107,6 +107,7 @@ class Detector : public o2::base::DetImpl private: /// Container for hits std::vector* mHits = nullptr; + int mTrackRefLastTrackID = -1; //! last track that received a TrackReference in the current event /// Geometry pointer Geometry* mGeometry = nullptr; //! diff --git a/Detectors/FIT/FV0/simulation/src/Detector.cxx b/Detectors/FIT/FV0/simulation/src/Detector.cxx index 07eb9053bf3b8..87b0b63607499 100644 --- a/Detectors/FIT/FV0/simulation/src/Detector.cxx +++ b/Detectors/FIT/FV0/simulation/src/Detector.cxx @@ -24,6 +24,7 @@ #include "FV0Simulation/Detector.h" #include "Framework/Logger.h" #include "DetectorsBase/Stack.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include "FV0Base/Geometry.h" using namespace o2::fv0; @@ -80,6 +81,13 @@ Bool_t Detector::ProcessHits(FairVolume* v) { // This method is called from the MC stepping + 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(fMC->GetStack()); + if (o2stack->GetCurrentTrackNumber() != mTrackRefLastTrackID) { + mTrackRefLastTrackID = o2stack->GetCurrentTrackNumber(); + o2stack->addTrackReference(o2::TrackReference(*fMC, GetDetId())); + } + } // Track only charged particles and photons bool isPhotonTrack = false; Int_t particlePdg = fMC->TrackPid(); @@ -150,6 +158,7 @@ void Detector::Reset() if (!o2::utils::ShmManager::Instance().isOperational()) { mHits->clear(); } + mTrackRefLastTrackID = -1; } void Detector::EndOfEvent() diff --git a/Detectors/ITSMFT/MFT/simulation/include/MFTSimulation/Detector.h b/Detectors/ITSMFT/MFT/simulation/include/MFTSimulation/Detector.h index 1e0a58fc80217..3c44705f747c0 100644 --- a/Detectors/ITSMFT/MFT/simulation/include/MFTSimulation/Detector.h +++ b/Detectors/ITSMFT/MFT/simulation/include/MFTSimulation/Detector.h @@ -148,6 +148,7 @@ class Detector : public o2::base::DetImpl private: /// Container for hit data std::vector* mHits; + int mTrackRefLastTrackID = -1; //! last track that received a TrackReference in the current event Detector(const Detector&); Detector& operator=(const Detector&); diff --git a/Detectors/ITSMFT/MFT/simulation/src/Detector.cxx b/Detectors/ITSMFT/MFT/simulation/src/Detector.cxx index 7a5243fde8063..c25b4b168f862 100644 --- a/Detectors/ITSMFT/MFT/simulation/src/Detector.cxx +++ b/Detectors/ITSMFT/MFT/simulation/src/Detector.cxx @@ -25,6 +25,7 @@ #include "Field/MagneticField.h" #include "DetectorsBase/Stack.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include #include "FairRootManager.h" @@ -96,6 +97,13 @@ Bool_t Detector::ProcessHits(FairVolume* vol) { // This method is called from the MC stepping + 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(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; @@ -744,4 +752,5 @@ void Detector::Reset() if (!o2::utils::ShmManager::Instance().isOperational()) { mHits->clear(); } + mTrackRefLastTrackID = -1; } diff --git a/Detectors/PHOS/simulation/src/Detector.cxx b/Detectors/PHOS/simulation/src/Detector.cxx index a8ff99b3805f1..086c4dd600a3d 100644 --- a/Detectors/PHOS/simulation/src/Detector.cxx +++ b/Detectors/PHOS/simulation/src/Detector.cxx @@ -163,6 +163,7 @@ Bool_t Detector::ProcessHits(FairVolume* v) if (isNewPartile) { // mark track to be kept by stack stack->addHit(GetDetId()); + stack->addTrackReference(o2::TrackReference(*fMC, GetDetId())); // PHOS entry (new superparent) } Double_t lostenergy = fMC->Edep(); diff --git a/Detectors/TOF/simulation/include/TOFSimulation/Detector.h b/Detectors/TOF/simulation/include/TOFSimulation/Detector.h index 1aa0efb6ef61c..1dfd730efe716 100644 --- a/Detectors/TOF/simulation/include/TOFSimulation/Detector.h +++ b/Detectors/TOF/simulation/include/TOFSimulation/Detector.h @@ -158,8 +158,9 @@ class Detector : public o2::base::DetImpl Int_t mEventNr; // event count Int_t mTOFSectors[o2::tof::Geo::NSECTORS]; - Bool_t mTOFHoles; // flag to allow for holes in front of the PHOS - Int_t mLastChannelID = -1; ///< Last channel seen by the hit + Bool_t mTOFHoles; // flag to allow for holes in front of the PHOS + Int_t mLastChannelID = -1; ///< Last channel seen by the hit + Int_t mTrackRefLastTrackID = -1; //! last track that received a TrackReference in the current event /// container for data points std::vector* mHits; //! diff --git a/Detectors/TOF/simulation/src/Detector.cxx b/Detectors/TOF/simulation/src/Detector.cxx index a9af40e208d81..b87645d373913 100644 --- a/Detectors/TOF/simulation/src/Detector.cxx +++ b/Detectors/TOF/simulation/src/Detector.cxx @@ -27,6 +27,7 @@ #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/MaterialManager.h" #include "DetectorsBase/Stack.h" +#include "SimulationDataFormat/MonopoleParticles.h" #include #include @@ -75,6 +76,13 @@ void Detector::InitializeO2Detector() Bool_t Detector::ProcessHits(FairVolume* v) { // This method is called from the MC stepping for the sensitive volume only + 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(fMC->GetStack()); + if (o2stack->GetCurrentTrackNumber() != mTrackRefLastTrackID) { + mTrackRefLastTrackID = o2stack->GetCurrentTrackNumber(); + o2stack->addTrackReference(o2::TrackReference(*fMC, GetDetId())); + } + } if (static_cast(fMC->TrackCharge()) == 0) { // set a very large step size for neutral particles return kFALSE; // take only charged particles @@ -130,6 +138,7 @@ void Detector::Reset() mHits->clear(); } mLastChannelID = -1; + mTrackRefLastTrackID = -1; } void Detector::CreateMaterials() @@ -1179,7 +1188,7 @@ void Detector::createBackZone(Double_t xtof, Double_t ytof, Double_t zlenA) cons if (mTOFHoles) { TVirtualMC::GetMC()->Gsvolu("FAIB", "BOX ", getMediumID(kAir), par, 3); // Air } - TVirtualMC::GetMC()->Gsvolu("FAIC", "BOX ", getMediumID(kAir), par, 3); // Air + TVirtualMC::GetMC()->Gsvolu("FAIC", "BOX ", getMediumID(kAir), par, 3); // Air Double_t feaParam[3] = {Geo::FEAPARAMETERS[0], Geo::FEAPARAMETERS[1], Geo::FEAPARAMETERS[2]}; Double_t feaRoof1[3] = {Geo::ROOF1PARAMETERS[0], Geo::ROOF1PARAMETERS[1], Geo::ROOF1PARAMETERS[2]};