Add DRTubesSDAction

Add DRTubesSDAction as in lopezzot/DREndcapTube v0.3. Comment out
RunAction and EventAction classes which are not needed for execution
inside k4geo. File name changed from DREndcapTubesSDAction to
DRTubesSDAction as this SDAction code will be used for DR barrel
detector using capillary tubes.

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)

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<DREndcapTubesSDData>::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<DREndcapTubesSDData>::defineCollections()
+{
+  std::string ROname = m_sensitive.readout().name();
+  m_collectionID = defineCollection<Geant4Calorimeter::Hit>(ROname + "ScinRight");
+  m_userData.collection_cher_right = defineCollection<Geant4Calorimeter::Hit>(ROname + "CherRight");
+  m_userData.collection_scin_left = defineCollection<Geant4Calorimeter::Hit>(ROname + "ScinLeft");
+  m_userData.collection_cher_left = defineCollection<Geant4Calorimeter::Hit>(ROname + "CherLeft");
+}
+
+// Function template specialization of Geant4SensitiveAction class.
+// Method that accesses the G4Step object at each track step.
+template<>
+bool Geant4SensitiveAction<DREndcapTubesSDData>::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<<it->first<<" "<<it->second<<std::endl;
+  // ---
+  // but this brute force method makes the simulaiton slower by more than two order of magnitudes
+  // (see https://github.com/AIDASoft/DD4hep/issues/1319).
+  // Therefore we use copynumbers instead of volIDs.
+
+#ifdef DREndcapTubesSDDebug
+  // Print out some info step-by-step in sensitive volumes
+  //
+  DRTubesSglHpr::PrintStepInfo(aStep);
+#endif
+
+  auto Edep = aStep->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<unsigned int>(aStep->GetPreStepPoint()->GetTouchable()->GetCopyNumber(3));
+  auto StaveID = static_cast<unsigned int>(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<G4OpBoundaryProcess*>(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<Geant4Calorimeter::Hit>(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<DREndcapTubesSDData> DREndcapTubesSDAction;
+}
+}  // namespace dd4hep
+DECLARE_GEANT4SENSITIVE(DREndcapTubesSDAction)
+
+//**************************************************************************