diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 25e06801560c9..216f6f8d71cf2 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 { @@ -80,6 +81,22 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { bool g4scoring = false; bool g4fluenceweight = false; + // Enable magnetic-monopole ionisation as defined in Detectors/gconfig/O2MonopolePhysics. + // The G4mplIonisation process (Ahlen stopping power) is attached to + // the O2 monopole particles (PDG +-4110000 / +-4120000) on top of the chosen + // reference physics list. Off by default + bool monopole = false; + // Monopole magnetic charge in units of the Dirac charge g_D = eplus/(2*alpha) + // (~68.5 eplus). 1.0 corresponds to a single classic Dirac monopole. + // The monopole is transported with the full dual Lorentz force + // F = g*(B - v x E/c^2). The electric part acts inside the TPC field cage, + // the only volume with an electric field, and its magnitude is the TPC's own + // TPCGEMParam.ElectricField[0] + float monopoleMagneticCharge = 1.f; + // 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. std::string fastSimModels = ""; // comma-separated model names to activate diff --git a/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h new file mode 100644 index 0000000000000..c1f546fd5e569 --- /dev/null +++ b/DataFormats/simulation/include/SimulationDataFormat/MonopoleParticles.h @@ -0,0 +1,46 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file MonopoleParticles.h +/// \brief Identification of the BSM magnetic monopoles across the simulation +/// @author M. Giacalone - September 2026 +/// +/// The monopole PDG codes are defined in O2MCApplication::AddParticles() and in +/// O2DatabasePDG. They are collected here for better syncing inside the simulation code +/// +/// To do: add dyons and other monopoles configurations when needed + +#ifndef ALICEO2_SIMULATION_MONOPOLEPARTICLES_H_ +#define ALICEO2_SIMULATION_MONOPOLEPARTICLES_H_ + +namespace o2::sim +{ + +/// Both species are electrically neutral for now, and the transport gives them the +/// same magnetic charge, with its sign taken from the PDG sign. +/// monopole intended to carry equal electric and magnetic charge +constexpr int MonopolePdgSymm = 4110000; +/// monopole intended to carry opposite electric and magnetic charge (not yet distinguished) +constexpr int MonopolePdgAsymm = 4120000; + +/// Default monopole mass in GeV. +constexpr double MonopoleMassDefaultGeV = 100.; + +/// true for the monopole and antimonopole of both species +constexpr bool isMonopole(int pdg) noexcept +{ + const int abspdg = pdg < 0 ? -pdg : pdg; + return abspdg == MonopolePdgSymm || abspdg == MonopolePdgAsymm; +} + +} // namespace o2::sim + +#endif diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index bd2f81f6a0338..2a4743c8c87d3 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -16,9 +16,11 @@ #ifndef O2_O2DATABASEPDG_H #define O2_O2DATABASEPDG_H +#include #include #include "TDatabasePDG.h" #include "TParticlePDG.h" +#include "SimulationDataFormat/MonopoleParticles.h" namespace o2 { @@ -44,9 +46,28 @@ class O2DatabasePDG } // adds ALICE particles to a given TDatabasePDG instance - static void addALICEParticles(TDatabasePDG* db = TDatabasePDG::Instance()); + // monopoleMass is the mass of all monopole species in GeV; it has to match + // G4Params.monopoleMass or the generator and the transport will disagree + // To fix: all quantities in the MCTrack dependent from GetMass() will get the default 100 GeV. + // A possible fix would be to stop deriving the masses from the compiled-in table, + // or a better solution should be found in the future + static void addALICEParticles(TDatabasePDG* db = TDatabasePDG::Instance(), + double monopoleMass = o2::sim::MonopoleMassDefaultGeV); static void addParticlesFromExternalFile(TDatabasePDG* db); + // true if all monopole species are registered with the given mass in GeV; TDatabasePDG + // keeps the first registration, so an earlier call with another mass makes this false + static bool hasMonopoleMass(TDatabasePDG* db, double monopoleMass) + { + for (int pdg : {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm, o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm}) { + const auto* particle = db->GetParticle(pdg); + if (particle == nullptr || std::abs(particle->Mass() - monopoleMass) > 1.e-6 * monopoleMass) { + return false; + } + } + return true; + } + // get particle's (if any) mass static Double_t MassImpl(TParticlePDG* particle, bool& success) { @@ -84,7 +105,7 @@ class O2DatabasePDG // by keeping this inline, we can use it in other parts of the code, for instance Framework or Analysis, // without needing to link against this library -inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db) +inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db, double monopoleMass) { // // Add ALICE particles to the ROOT PDG data base @@ -189,6 +210,14 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db) db->AddParticle("CHI2P_B2", " ", 10.269, kFALSE, 0.0, 0, "meson", 100555); db->AddParticle("UPSLON4S", " ", 10.580, kFALSE, 0.0, 0, "meson", 300553); + // BSM targeted inclusions + // Monopoles with same electric and magnetic charge + db->AddParticle("Monopole_symm", "Monopole_symm", monopoleMass, kTRUE, 0.0, 0, "BSM", 4110000); + db->AddParticle("AntiMonopole_symm", "AntiMonopole_symm", monopoleMass, kTRUE, 0.0, 0, "BSM", -4110000); + // Monopoles with opposite electric and magnetic charge + db->AddParticle("Monopole_asymm", "Monopole_asymm", monopoleMass, kTRUE, 0.0, 0, "BSM", 4120000); + db->AddParticle("AntiMonopole_asymm", "AntiMonopole_asymm", monopoleMass, kTRUE, 0.0, 0, "BSM", -4120000); + // IONS // // Done by default now from Pythia6 table 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/ITS/simulation/src/Detector.cxx b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx index c01657d4a7c39..6f4f1fa5fe028 100644 --- a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx +++ b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx @@ -24,6 +24,8 @@ #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 @@ -318,7 +320,15 @@ 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. 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 && sMonopoleIonisation && o2::sim::isMonopole(fMC->TrackPid()); + if (isNeutral && !isMonopole) { return kFALSE; } @@ -386,6 +396,10 @@ Bool_t Detector::ProcessHits(FairVolume* vol) mTrackData.mHitStarted = true; } if (stopHit) { + // A monopole that deposited nothing in the sensor leaves no hit + if (isMonopole && mTrackData.mEnergyLoss <= 0.) { + return kFALSE; + } TLorentzVector positionStop; fMC->TrackPosition(positionStop); // Retrieve the indices with the volume path 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]}; diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index da75d8f81f466..28e8eff25cfc0 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -17,6 +17,8 @@ #include "DetectorsBase/Stack.h" #include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/MonopoleParticles.h" +#include "SimConfig/G4Params.h" #include "FairVolume.h" // for FairVolume @@ -120,11 +122,24 @@ Bool_t Detector::ProcessHits(FairVolume* vol) /* This method is called from the MC stepping for the sensitive volume only */ // LOG(info) << "tpc::ProcessHits"; const double trackCharge = fMC->TrackCharge(); + // Magnetic monopoles have zero electric charge but ionise the gas through + // their magnetic charge energy loss (G4mplIonisation), so they must not be + // rejected by the electric-charge gate. The PDG lookup is only worth doing for neutral particle + // To-do: add dyons case + bool isMonopole = false; if (static_cast(trackCharge) == 0) { - - // set a very large step size for neutral particles - fMC->SetMaxStep(1.e10); - return kFALSE; // take only charged particles + static const bool sMonopoleIonisation = o2::conf::G4Params::Instance().monopole; + isMonopole = sMonopoleIonisation && o2::sim::isMonopole(fMC->TrackPid()); + if (!isMonopole) { + // set a very large step size for neutral particles + fMC->SetMaxStep(1.e10); + return kFALSE; // take only charged particles + } + if (fMC->Edep() <= 0.) { + // An ionising monopole depositing nothing on this step: there + // is no hit to create, and the step limit is deliberately left alone + return kFALSE; + } } // ===| SET THE LENGTH OF THE NEXT ENERGY LOSS STEP |========================= @@ -198,7 +213,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 // use Geant4 energy deposit directly for ionisation (Kr-83m calibration simulations) if (detParam.UseGeant4Edep) { @@ -209,6 +224,16 @@ Bool_t Detector::ProcessHits(FairVolume* vol) if (meanIon > 0.) { numberOfElectrons = static_cast(gasParam.FanoFactorG4 * Gamma(meanIon / gasParam.FanoFactorG4)); } + } else if (isMonopole) { + // ---| MONOPOLE IONISATION |-------- + // A magnetic monopole ionises the gas via G4mplIonisation (Ahlen stopping + // power), which is not described by the electric-charge model (no e-charge) + // Ionisation electrons calculated directly from the energy deposited in this step: Nel = Edep / Wion. + // To check by TPC experts if this is actually the best way... + numberOfElectrons = static_cast(fMC->Edep() / static_cast(gasParam.Wion)); + // The number of electrons is stored as a float in the HitGroup: maximum cap at + // 2^24 (16777216) ==> largest integer a IEEE-754 float can represent exactly + numberOfElectrons = TMath::Min(numberOfElectrons, 16777216); } else { // ---| Stepsize in cm |--- const double stepSize = fMC->TrackStep(); diff --git a/Detectors/gconfig/CMakeLists.txt b/Detectors/gconfig/CMakeLists.txt index d6fffad875401..a4d207d2e3543 100644 --- a/Detectors/gconfig/CMakeLists.txt +++ b/Detectors/gconfig/CMakeLists.txt @@ -28,7 +28,7 @@ if(TGeo2VecGeom_FOUND) endif() set(G4SETUP_SOURCES src/G4Config.cxx src/G4RunConfiguration.cxx src/G4LocalFieldConstruction.cxx - src/VecGeomNavigation.cxx) + src/VecGeomNavigation.cxx src/O2MonopolePhysics.cxx) if(G4SETUP_WITH_VECGEOM) list(APPEND G4SETUP_SOURCES src/VecGeomG4Map.cxx src/VecGeomChecks.cxx src/VecGeomG4NavigatorBase.cxx src/VecGeomG4Navigator.cxx diff --git a/Detectors/gconfig/g4Config.C b/Detectors/gconfig/g4Config.C index 3cc2133961ab4..43d49a7a4796d 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" #include "SimConfig/G4ScoringMerger.h" #include "G4ScoringManager.hh" #include "G4VScoringMesh.hh" @@ -124,11 +125,21 @@ void Config() LOG(fatal) << "Unsupported geometry navigation mode"; } - // o2::g4config::G4RunConfiguration adds the fast-simulation hook and the local - // magnetic fields; with neither configured it behaves like TG4RunConfiguration. - auto runConfiguration = new o2::g4config::G4RunConfiguration(geomNavStr, physicsSetup, - "stepLimiter+specialCuts", - specialStacking, mtMode); + TG4RunConfiguration* runConfiguration = nullptr; + if (g4Params.monopole) { + // Reference physics list + magnetic-monopole ionisation attached + // The reference physics is unchanged. + std::cout << "Monopole ionisation physics requested (G4.monopole=1)\n"; + runConfiguration = o2::g4config::createMonopoleRunConfiguration( + geomNavStr, physicsSetup, "stepLimiter+specialCuts", specialStacking, mtMode, + g4Params.monopoleMagneticCharge); + } else { + // o2::g4config::G4RunConfiguration adds the fast-simulation hook and the local + // magnetic fields; with neither configured it behaves like TG4RunConfiguration. + runConfiguration = new o2::g4config::G4RunConfiguration(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..7e8508552c3e7 --- /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, 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; +/// 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 +/// \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..1bb6eeaf97852 --- /dev/null +++ b/Detectors/gconfig/src/O2MonopolePhysics.cxx @@ -0,0 +1,556 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file O2MonopolePhysics.cxx +/// \brief Opt-in magnetic-monopole ionisation physics for the O2 Geant4 engine. +/// \author M+Giacalone - July 2026 +/// +/// O2 defines the BSM monopole particles (PDG +-4110000 / +-4120000) via +/// O2MCApplication::AddParticles(), so they are transported by Geant4, but no +/// energy-loss process is attached to them by any of the stock reference +/// physics lists. +/// +/// This file adds - the magnetic-monopole ionisation process G4mplIonisation +/// (Ahlen stopping power, Rev. Mod. Phys. 52 (1980) 121) to those particles, when requested. +/// The implementation takes Geant4 `examples/extended/exoticphysics/monopole` example (and the +/// equivalent geant4_vmc "monopole" physics builder). The main difference is that +/// here the process is attached to the pre-existing O2 monopole particles (with correct PDG ID() +/// rather than to a freshly created G4Monopole +/// +/// G4mplIonisation and G4mplIonisationWithDeltaModel are taken directly from the standard +/// Geant4 libraries (libG4processes) + +#include "SimSetup/O2MonopolePhysics.h" +#include "CommonUtils/ConfigurableParam.h" +#include "SimulationDataFormat/MonopoleParticles.h" + +#include // needed to instantiate getValueAs<> + +#include + +#include + +#include "TG4ComposedPhysicsList.h" +#include "SimSetup/G4RunConfiguration.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace o2 +{ +namespace g4config +{ + +namespace +{ +// PDG codes of the O2 monopole species, as defined in +// O2MCApplication::AddParticles() and O2DatabasePDG: +// +-4110000 : "symmetric" monopoles +// +-4120000 : "asymmetric" monopoles +constexpr std::array gMonopolePDGs = {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm, + o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm}; + +/// 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 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; + +/// 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]"); + 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 " + "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) + { + } + + /// 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; + mMagCharge = CLHEP::eplus * mMagneticChargeSign * mMagneticChargeEplus * 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). + 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 + // 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: + 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.; + G4double mMassCof = 0.; +}; + +//____________________________________________________________________________ +/// The chord finder integrating O2MonopoleEquation, and that equation +struct MonopoleFieldSetup { + 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, in the field of the run. +/// +/// 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 +/// 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 MonopoleFieldSetup buildMonopoleFieldSetup(double magneticChargeEplusUnits, double tpcDriftFieldGeant4) +{ + auto* transportationManager = G4TransportationManager::GetTransportationManager(); + auto* fieldManager = transportationManager != nullptr ? transportationManager->GetFieldManager() : nullptr; + if (fieldManager == nullptr) { + 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(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); + // 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()); + } + + if (tpcDriftFieldGeant4 > 0.) { + 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 built (F = g*B), magnetic charge = " + << magneticChargeEplusUnits << " eplus (TPC drift field coupling disabled)"; + } + return {chordFinder, equation}; +} + +//____________________________________________________________________________ +/// Transportation for the monopole species only. +/// +/// These must be true for the monopoles and false for every other particle: +/// +/// - 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, 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(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 = 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); + + fieldManager->SetFieldChangesEnergy(previousChangesEnergy); + fieldManager->SetChordFinder(previousChordFinder); + G4Transportation::EnableMagneticMoment(previousMoment); + return length; + } + + private: + 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 void 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; + } + 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; + } + + // 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); + 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); +} + +//____________________________________________________________________________ +/// 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 (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: + /// \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(); + + // Built on the first monopole species actually found, then shared by the rest: + // 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; + + 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 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 + // 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); + } + + // 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; + 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"; + } + } + + 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. +/// Derives from o2::g4config::G4RunConfiguration (rather than TG4RunConfiguration +/// directly) so that monopole runs keep the G4 fast-simulation hook and the local magnetic fields +class O2G4RunConfiguration : public o2::g4config::G4RunConfiguration +{ + public: + O2G4RunConfiguration(const TString& userGeometry, const TString& physicsList, + const TString& specialProcess, Bool_t specialStacking, + Bool_t mtApplication, double magneticChargeEplus) + : o2::g4config::G4RunConfiguration(userGeometry, physicsList, specialProcess, specialStacking, mtApplication), + mMagneticCharge(magneticChargeEplus) + { + } + + G4VUserPhysicsList* CreatePhysicsList() override + { + G4VUserPhysicsList* physicsList = o2::g4config::G4RunConfiguration::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(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; + } + + 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. + // 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 + << " Dirac charge(s) = " << magneticChargeEplus / CLHEP::eplus << " eplus"; + return new O2G4RunConfiguration(userGeometry, physicsList, specialProcess, + specialStacking, mtApplication, magneticChargeEplus); +} + +} // namespace g4config +} // namespace o2 diff --git a/Generators/src/GeneratorFactory.cxx b/Generators/src/GeneratorFactory.cxx index 85077a7a5db21..7c59d70a7c864 100644 --- a/Generators/src/GeneratorFactory.cxx +++ b/Generators/src/GeneratorFactory.cxx @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -140,7 +141,12 @@ void GeneratorFactory::setPrimaryGenerator(o2::conf::SimConfig const& conf, Fair /** generators **/ - o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance()); + // Monopole configurable mass added to the PDG database + 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) std::string hyperloopExtHybridConfigFile; // set when IS_HYPERLOOP is defined diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index a834b4a21c98d..82fdda4eace58 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/Steer/src/O2MCApplication.cxx @@ -1751,6 +1751,27 @@ 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, 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; + } + // 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", 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", 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() diff --git a/macro/o2sim.C b/macro/o2sim.C index fd44de6e93b54..c99cd012d4d73 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,14 @@ 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. + 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; {