diff --git a/CMakeLists.txt b/CMakeLists.txt index 811779559..e6f765ad1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -117,6 +117,8 @@ file(GLOB G4sources ./plugins/Geant4Output2EDM4hep_DRC.cpp ./plugins/DRCaloFastSimModel.cpp ./plugins/DRCaloFastSimModel.h + ./plugins/DRTubesSDAction.hh + ./plugins/DRTubesSDAction.cpp ) if(DD4HEP_USE_PYROOT) diff --git a/plugins/DRTubesSDAction.cpp b/plugins/DRTubesSDAction.cpp new file mode 100644 index 000000000..ba238f15d --- /dev/null +++ b/plugins/DRTubesSDAction.cpp @@ -0,0 +1,291 @@ +//************************************************************************** +// \file DRTubesSDAction.cpp +// \brief: Implementation of Geant4SensitiveAction template class for +// dual-readout-tubes calorimeters +// \author: Lorenzo Pezzotti (CERN) @lopezzot +// \start date: 11 August 2024 +//************************************************************************** + +// Includers from DD4HEP +#include "DD4hep/Segmentations.h" +#include "DD4hep/Version.h" +#include "DDG4/Factories.h" +#include "DDG4/Geant4EventAction.h" +#include "DDG4/Geant4GeneratorAction.h" +#include "DDG4/Geant4Mapping.h" +#include "DDG4/Geant4RunAction.h" +#include "DDG4/Geant4SensDetAction.inl" + +// Includers from Geant4 +#include "G4OpBoundaryProcess.hh" +#include "G4OpticalPhoton.hh" +#include "G4Poisson.hh" +#include "G4ProcessManager.hh" +#include "G4Types.hh" +#include "G4UserSteppingAction.hh" +#include "G4VProcess.hh" +#include "globals.hh" + +// Includers from project files +#include "DRTubesSglHpr.hh" + +// #define DREndcapTubesSDDebug + +namespace dd4hep +{ +namespace sim +{ +class DREndcapTubesSDData +{ + // Constructor and destructor + // + public: + DREndcapTubesSDData() = default; + ~DREndcapTubesSDData() = default; + + // Methods + // + public: + //void beginRun(const G4Run* run) + //{ + //fRunAction = new DREndcapTubesRunAction; + //fEvtAction = new DREndcapTubesEvtAction; + //fRunAction->BeginOfRunAction(run); + //} + + //void endRun(const G4Run* run) { fRunAction->EndOfRunAction(run); } + + //void beginEvent(const G4Event* event) { fEvtAction->BeginOfEventAction(event); } + + //void endEvent(const G4Event* event) { fEvtAction->EndOfEventAction(event); } + + // Fields + // + public: + //DREndcapTubesRunAction* fRunAction; + //DREndcapTubesEvtAction* fEvtAction; + Geant4Sensitive* sensitive{}; + int collection_cher_right; + int collection_cher_left; + int collection_scin_left; +}; +} // namespace sim +} // namespace dd4hep + +namespace dd4hep +{ +namespace sim +{ + +// Function template specialization of Geant4SensitiveAction class. +// Define actions +template<> +void Geant4SensitiveAction::initialize() +{ + //eventAction().callAtBegin(&m_userData, &DREndcapTubesSDData::beginEvent); + //eventAction().callAtEnd(&m_userData, &DREndcapTubesSDData::endEvent); + + //runAction().callAtBegin(&m_userData, &DREndcapTubesSDData::beginRun); + //runAction().callAtEnd(&m_userData, &DREndcapTubesSDData::endRun); + + m_userData.sensitive = this; +} + +// Function template specialization of Geant4SensitiveAction class. +// Define collections created by this sensitivie action object +template<> +void Geant4SensitiveAction::defineCollections() +{ + std::string ROname = m_sensitive.readout().name(); + m_collectionID = defineCollection(ROname + "ScinRight"); + m_userData.collection_cher_right = defineCollection(ROname + "CherRight"); + m_userData.collection_scin_left = defineCollection(ROname + "ScinLeft"); + m_userData.collection_cher_left = defineCollection(ROname + "CherLeft"); +} + +// Function template specialization of Geant4SensitiveAction class. +// Method that accesses the G4Step object at each track step. +template<> +bool Geant4SensitiveAction::process(const G4Step* aStep, + G4TouchableHistory* /*history*/) +{ + // NOTE: Here we do manipulation of the signal in each fiber (Scintillating and Cherenkov) + // to compute the calorimeter signal and populate the corresponding hit. + // Sensitive volumes are associated in the steering file using the DD4hep regexSensitiveDetector + // and matching the substring DRETS. Usually DD4hep volIDs are retrieved with something like this + // --- + // dd4hep::BitFieldCoder + // decoder("stave:10,tower:6,air:1,col:16,row:16,clad:1,core:1,cherenkov:1"); + // auto VolID = volumeID(aStep); auto CherenkovID = decoder.get(VolID,"cherenkov"); + // --- + // However using regexSensitiveDetector does not populate the cache that allows using the + // volumeID() method. One can still access volIDs in this sensitive detector action with something + // like this + // --- + // G4VPhysicalVolume* PhysVol = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(); + // Geant4Mapping& Mapping = Geant4Mapping::instance(); + // PlacedVolume PlacedVol = Mapping.placement(PhysVol); + // const PlacedVolumeExtension::VolIDs& TestIDs = PlacedVol.volIDs(); + // auto it = TestIDs.find("name"); + // std::cout<first<<" "<second<GetTotalEnergyDeposit(); + auto cpNo = aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(); + // The second bit of the CopyNumber corresponds to the "core" entry: + // 1 if the step is in the fiber core (S or C) and 0 if it is + // in the fiber cladding (C only) + unsigned int CoreID = (cpNo & 0b10) >> 1; // take CpNo 2nd bit + bool IsCherClad = (CoreID == 0); + // The first bit of the CopyNumber corresponds to the "cherenkov" entry + // 1 for C fibers and 0 for S fibers + unsigned int CherenkovID = cpNo & 0b1; + bool IsScin = (CherenkovID == 0); + // Skip this step if edep is 0 and it is a scintillating fiber + if (IsScin && Edep == 0.) return true; + // If it is a track inside the cherenkov CLADDING skip this step, + // if it is an optical photon kill it first + if (IsCherClad) { + if (aStep->GetTrack()->GetParticleDefinition() == G4OpticalPhoton::Definition()) { + aStep->GetTrack()->SetTrackStatus(fStopAndKill); + } + return true; + } + + // Now we are inside fibers' core volume (either Scintillating or Cherenkov) + + // We recreate the TubeID from the tube copynumber: + // fist16 bits for the columnID and second 16 bits for the rowID + auto TubeID = aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(2); + unsigned int ColumnID = TubeID >> 16; + unsigned int RawID = TubeID & 0xFFFF; + auto TowerID = static_cast(aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(3)); + auto StaveID = static_cast(aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(4)); + + VolumeID VolID = 0; // recreate the 64-bit VolumeID + BitFieldCoder bc("stave:10,tower:6,air:1,col:16,row:16,clad:1,core:1,cherenkov:1"); + bc.set(VolID, "stave" , StaveID); + bc.set(VolID, "tower" , TowerID); + bc.set(VolID, "air", 0); + bc.set(VolID, "col", ColumnID); + bc.set(VolID, "row", RawID); + bc.set(VolID, "clad", 1); + bc.set(VolID, "core", CoreID); + bc.set(VolID, "cherenkov", CherenkovID); + + bool IsRight = (aStep->GetPreStepPoint()->GetPosition().z() > 0.); + + // We now calculate the signal in S and C fiber according to the step contribution + // + G4double steplength = aStep->GetStepLength(); + G4int signalhit = 0; + Geant4HitCollection* coll = (IsRight && IsScin) ? collection(m_collectionID) + : (IsRight && !IsScin) ? collection(m_userData.collection_cher_right) + : (!IsRight && IsScin) ? collection(m_userData.collection_scin_left) + : collection(m_userData.collection_cher_left); + + if (IsScin) { // it is a scintillating fiber + + //m_userData.fEvtAction->AddEdepScin(Edep); + + if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() == 0 || steplength == 0.) { + return true; // not ionizing particle + } + G4double distance_to_sipm = DRTubesSglHpr::GetDistanceToSiPM(aStep); + signalhit = + DRTubesSglHpr::SmearSSignal(DRTubesSglHpr::ApplyBirks(Edep, steplength)); + signalhit = DRTubesSglHpr::AttenuateSSignal(signalhit, distance_to_sipm); + if (signalhit == 0) return true; + //m_userData.fEvtAction->AddSglScin(signalhit); + } // end of scintillating fibre sigal calculation + + else { // it is a Cherenkov fiber + // save mc truth info in analysismanager auxiliary outputfile + //m_userData.fEvtAction->AddEdepCher(Edep); + // calculate the signal in terms of Cherenkov photo-electrons + if (aStep->GetTrack()->GetParticleDefinition() == G4OpticalPhoton::Definition()) { + G4OpBoundaryProcessStatus theStatus = Undefined; + G4ProcessManager* OpManager = G4OpticalPhoton::OpticalPhoton()->GetProcessManager(); + + if (OpManager) { + G4int MAXofPostStepLoops = OpManager->GetPostStepProcessVector()->entries(); + G4ProcessVector* fPostStepDoItVector = OpManager->GetPostStepProcessVector(typeDoIt); + + for (G4int i = 0; i < MAXofPostStepLoops; i++) { + G4VProcess* fCurrentProcess = (*fPostStepDoItVector)[i]; + G4OpBoundaryProcess* fOpProcess = dynamic_cast(fCurrentProcess); + if (fOpProcess) { + theStatus = fOpProcess->GetStatus(); + break; + } + } + } + + switch (theStatus) { + case TotalInternalReflection: { + // Kill Cherenkov photons inside fibers travelling towards the inner tip + if (!DRTubesSglHpr::IsReflectedForward(aStep)) return true; + G4double distance_to_sipm = DRTubesSglHpr::GetDistanceToSiPM(aStep); + G4int c_signal = DRTubesSglHpr::SmearCSignal(); + signalhit = DRTubesSglHpr::AttenuateCSignal(c_signal, distance_to_sipm); + if (signalhit == 0) return true; + // save mc truth info in analysismanager auxiliary outputfile + //m_userData.fEvtAction->AddSglCher(signalhit); + aStep->GetTrack()->SetTrackStatus(fStopAndKill); + break; + } + default: + aStep->GetTrack()->SetTrackStatus(fStopAndKill); + return true; + } // end of swich cases + } // end of optical photon + else { + return true; + } + } // end of Cherenkov fiber + + // We are going to create an hit per each fiber with a signal above 0 + // Each fiber is identified with a unique volID + // + Geant4Calorimeter::Hit* hit = coll->findByKey(VolID); // the hit + if (!hit) { // if the hit does not exist yet, create it + hit = new Geant4Calorimeter::Hit(); + hit->cellID = VolID; // this should be assigned only once + G4ThreeVector FiberVec = DRTubesSglHpr::CalculateFiberPosition(aStep); + Position FiberPos(FiberVec.x(), FiberVec.y(), FiberVec.z()); + hit->position = FiberPos; // this should be assigned only once + // Note, when the hit is saved in edm4hep format the energyDeposit is + // divided by 1000, i.e. it translates from MeV (Geant4 unit) to GeV (DD4hep unit). + // Here I am using this field to save photo-electrons, so I multiply it by 1000 + hit->energyDeposit = signalhit * 1000; + coll->add(VolID, hit); // add the hit to the hit collection + } + else { // if the hit exists already, increment its fields + hit->energyDeposit += signalhit * 1000; + } + return true; +} // end of Geant4SensitiveAction::process() method specialization + +} // namespace sim +} // namespace dd4hep + +//--- Factory declaration +namespace dd4hep +{ +namespace sim +{ +typedef Geant4SensitiveAction DREndcapTubesSDAction; +} +} // namespace dd4hep +DECLARE_GEANT4SENSITIVE(DREndcapTubesSDAction) + +//************************************************************************** diff --git a/plugins/DRTubesSglHpr.hh b/plugins/DRTubesSglHpr.hh new file mode 100644 index 000000000..6baa3fc10 --- /dev/null +++ b/plugins/DRTubesSglHpr.hh @@ -0,0 +1,178 @@ +//************************************************************************** +// \file DRTubesSglHpr.hh +// \brief: definition of DRTubesSglHpr class +// \author: Lorenzo Pezzotti (CERN) @lopezzot +// \start date: 13 August 2024 +//************************************************************************** + +// Header-only utility class to store methods needed when computing +// the DREndcapTubes signals in scintillating and Cherenkov fibers + +#ifndef DRTubesSglHpr_h +# define DRTubesSglHpr_h 1 + +// Includers from Geant4 +// +# include "G4Tubs.hh" +# include "globals.hh" + +# include "CLHEP/Units/SystemOfUnits.h" + +class DRTubesSglHpr +{ + public: + DRTubesSglHpr() = delete; + + private: + // Fields + // + static constexpr G4double fk_B = 0.126; // Birks constant + static constexpr G4double fSAttenuationLength = 1000.0 * CLHEP::m; + static constexpr G4double fCAttenuationLength = 1000.0 * CLHEP::m; + + public: + // Methods + // + // Apply Briks Law + static constexpr G4double ApplyBirks(const G4double& de, const G4double& steplength) + { + return (de / steplength) / (1 + fk_B * (de / steplength)) * steplength; + } + + // Smear S signal according to Poissonian fluctuations and light yield + static G4int SmearSSignal(const G4double& satde) { return G4Poisson(satde * 9.5); } + + // Smear C signal according to Poissonian fluctuations and light yield + static G4int SmearCSignal() { return G4Poisson(0.177); } + + // Calculate distance from step in fiber to SiPM + inline static G4double GetDistanceToSiPM(const G4Step* step, bool prestep); + static G4double GetDistanceToSiPM(const G4Step* step) { return GetDistanceToSiPM(step, true); } + + // Calculate how many photons survived after light attenuation + inline static G4int AttenuateHelper(const G4int& signal, const G4double& distance, + const G4double& attenuation_length); + + // Attenuate light in fibers + static G4int AttenuateSSignal(const G4int& signal, const G4double& distance) + { + return AttenuateHelper(signal, distance, fSAttenuationLength); + } + + static G4int AttenuateCSignal(const G4int& signal, const G4double& distance) + { + return AttenuateHelper(signal, distance, fCAttenuationLength); + } + + inline static G4ThreeVector CalculateFiberPosition(const G4Step* step); + + // Check if photon is travelling towards SiPM + inline static bool IsReflectedForward(const G4Step* step); + + // Print step info for debugging purposes + inline static void PrintStepInfo(const G4Step* aStep); +}; + +inline G4double DRTubesSglHpr::GetDistanceToSiPM(const G4Step* step, bool prestep) +{ + // Get the pre-step point + const G4StepPoint* StepPoint = prestep ? step->GetPreStepPoint() : step->GetPostStepPoint(); + // Get the global position of the step point + G4ThreeVector globalPos = StepPoint->GetPosition(); + // Get the local position of the step point in the current volume's coordinate system + G4ThreeVector localPos = + StepPoint->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(globalPos); + + // Get the logical volume of the current step + G4LogicalVolume* currentVolume = StepPoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume(); + // Get the solid associated with the logical volume + G4Tubs* solid = dynamic_cast(currentVolume->GetSolid()); + // Get the dimensions of the solid (size of the volume) + G4double size = solid->GetZHalfLength(); + + G4double distance_to_sipm = size - localPos.z(); + return distance_to_sipm; +} + +inline G4int DRTubesSglHpr::AttenuateHelper(const G4int& signal, const G4double& distance, + const G4double& attenuation_length) +{ + double probability_of_survival = exp(-distance / attenuation_length); + + G4int survived_photons = 0; + for (int i = 0; i < signal; i++) { + // Simulate drawing between 0 and 1 with probability x of getting 1 + if (G4UniformRand() <= probability_of_survival) survived_photons++; + } + return survived_photons; +} + +inline G4ThreeVector DRTubesSglHpr::CalculateFiberPosition(const G4Step* step) +{ + G4TouchableHandle theTouchable = step->GetPreStepPoint()->GetTouchableHandle(); + G4ThreeVector origin(0., 0., 0.); + G4ThreeVector zdir(0., 0., 1.); + G4ThreeVector vectPos = + theTouchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(origin); + G4ThreeVector direction = + theTouchable->GetHistory()->GetTopTransform().Inverse().TransformAxis(zdir); + + // Get the logical volume of the current step + G4LogicalVolume* currentVolume = + step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume(); + // Get the solid associated with the logical volume + G4Tubs* solid = dynamic_cast(currentVolume->GetSolid()); + // Get the dimensions of the solid (size of the volume) + G4double size = solid->GetZHalfLength(); + G4double lengthfiber = size * 2.; + G4ThreeVector Halffibervect = direction * lengthfiber / 2; + // Fibre tip position + G4ThreeVector vectPostip = vectPos - Halffibervect; + // SiPM position + // G4ThreeVector SiPMvecPos = vectPos + Halffibervect; + + return vectPostip; +} + +// Check if photon is travelling towards SiPM +// If a (Cherenkov) photon is internally reflected in fibers +// it might travel towards the SiPM or the inner finer tip. +// As we do not consider reflections at the inner fiber tip +// photons travelling backwards should be killed. +inline bool DRTubesSglHpr::IsReflectedForward(const G4Step* step) +{ + double PreStepDistance = GetDistanceToSiPM(step, true); + double PostStepDistance = GetDistanceToSiPM(step, false); + bool IsReflectedForward = (PostStepDistance < PreStepDistance) ? true : false; + return IsReflectedForward; +} + +inline void DRTubesSglHpr::PrintStepInfo(const G4Step* aStep) +{ + std::cout << "-------------------------------" << std::endl; + std::cout << "--> DREndcapTubes: track info: " << std::endl; + std::cout << "----> Track #: " << aStep->GetTrack()->GetTrackID() << " " + << "Step #: " << aStep->GetTrack()->GetCurrentStepNumber() << " " + << "Volume: " << aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() + << " " << std::endl; + std::cout << "--> DREndcapTubes:: position info(mm): " << std::endl; + std::cout << "----> x: " << aStep->GetPreStepPoint()->GetPosition().x() + << " y: " << aStep->GetPreStepPoint()->GetPosition().y() + << " z: " << aStep->GetPreStepPoint()->GetPosition().z() << std::endl; + std::cout << "--> DREndcapTubes: particle info: " << std::endl; + std::cout << "----> Particle " << aStep->GetTrack()->GetParticleDefinition()->GetParticleName() + << " " + << "Dep(MeV) " << aStep->GetTotalEnergyDeposit() << " " + << "Mat " << aStep->GetPreStepPoint()->GetMaterial()->GetName() << " " + << "Vol " << aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() + << " " + << "CpNo " << aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber() << " " + << "CpNo1 " << aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(1) << " " + << "CpNo2 " << aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(2) << " " + << "CpNo3 " << aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(3) << " " + << "CpNo4 " << aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(4) << std::endl; +} + +#endif // DRTubesSglHpr_h + +//**************************************************************************