From 15068d41a1f1e43894d9b8921c2056b4051d2245 Mon Sep 17 00:00:00 2001
From: Chengdong Fu <fucd@ihep.ac.cn>
Date: Tue, 17 Nov 2020 15:05:30 +0800
Subject: [PATCH] add special TPC SD migrated from Marlin processor
---
.../DetSimGeom/src/AnExampleDetElemTool.cpp | 11 +
.../DetSimGeom/src/AnExampleDetElemTool.h | 1 +
Simulation/DetSimSD/CMakeLists.txt | 3 +
.../DetSimSD/DetSimSD/DDG4SensitiveDetector.h | 4 +-
.../DetSimSD/src/DDG4SensitiveDetector.cpp | 4 +-
.../src/TimeProjectionChamberSensDetTool.cpp | 57 +++
.../src/TimeProjectionChamberSensDetTool.h | 40 +++
...TimeProjectionChamberSensitiveDetector.cpp | 331 ++++++++++++++++++
.../TimeProjectionChamberSensitiveDetector.h | 68 ++++
9 files changed, 515 insertions(+), 4 deletions(-)
create mode 100644 Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.cpp
create mode 100644 Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.h
create mode 100644 Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.cpp
create mode 100644 Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.h
diff --git a/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp b/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp
index fffebf79..001675db 100644
--- a/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp
+++ b/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp
@@ -120,6 +120,16 @@ AnExampleDetElemTool::ConstructSDandField() {
warning() << "DriftChamberSensDetTool is not found. " << endmsg;
}
}
+ else if (nam == "TPC") {
+ m_tpc_sdtool = ToolHandle<ISensDetTool>("TimeProjectionChamberSensDetTool");
+ if (m_tpc_sdtool) {
+ info() << "Find the TimeProjectionChamberSensDetTool" << endmsg;
+ g4sd = m_tpc_sdtool->createSD(nam);
+ }
+ else {
+ warning() << "TimeProjectionChamberSensDetTool is not found, and default tracker SD will be used" << endmsg;
+ }
+ }
}
}
@@ -199,6 +209,7 @@ AnExampleDetElemTool::initialize() {
m_calo_sdtool = ToolHandle<ISensDetTool>("CalorimeterSensDetTool");
m_driftchamber_sdtool = ToolHandle<ISensDetTool>("DriftChamberSensDetTool");
+ m_tpc_sdtool = ToolHandle<ISensDetTool>("TimeProjectionChamberSensDetTool");
return sc;
}
diff --git a/Simulation/DetSimGeom/src/AnExampleDetElemTool.h b/Simulation/DetSimGeom/src/AnExampleDetElemTool.h
index 31fcfc0e..908de073 100644
--- a/Simulation/DetSimGeom/src/AnExampleDetElemTool.h
+++ b/Simulation/DetSimGeom/src/AnExampleDetElemTool.h
@@ -34,6 +34,7 @@ private:
SmartIF<IGeomSvc> m_geosvc;
ToolHandle<ISensDetTool> m_calo_sdtool;
ToolHandle<ISensDetTool> m_driftchamber_sdtool;
+ ToolHandle<ISensDetTool> m_tpc_sdtool;
};
#endif
diff --git a/Simulation/DetSimSD/CMakeLists.txt b/Simulation/DetSimSD/CMakeLists.txt
index 1828c1e1..33d00e18 100644
--- a/Simulation/DetSimSD/CMakeLists.txt
+++ b/Simulation/DetSimSD/CMakeLists.txt
@@ -19,6 +19,9 @@ set(DetSimSD_srcs
src/DriftChamberSensDetTool.cpp
src/DriftChamberSensitiveDetector.cpp
+
+ src/TimeProjectionChamberSensDetTool.cpp
+ src/TimeProjectionChamberSensitiveDetector.cpp
)
gaudi_add_module(DetSimSD ${DetSimSD_srcs}
diff --git a/Simulation/DetSimSD/DetSimSD/DDG4SensitiveDetector.h b/Simulation/DetSimSD/DetSimSD/DDG4SensitiveDetector.h
index f7806f7d..cee8826e 100644
--- a/Simulation/DetSimSD/DetSimSD/DDG4SensitiveDetector.h
+++ b/Simulation/DetSimSD/DetSimSD/DDG4SensitiveDetector.h
@@ -45,12 +45,12 @@ public:
/// Returns the volumeID of the sensitive volume corresponding to the step -
/// combining the VolIDS of the complete geometry path (Geant4TouchableHistory)
// from the current sensitive volume to the world volume
- virtual long long getVolumeID(G4Step* step);
+ virtual long long getVolumeID(const G4Step* step);
/// Returns the volumeID of the sensitive volume corresponding to the step -
/// combining the VolIDS of the complete geometry path (Geant4TouchableHistory)
// from the current sensitive volume to the world volume
- virtual long long getCellID(G4Step* step);
+ virtual long long getCellID(const G4Step* step);
protected:
diff --git a/Simulation/DetSimSD/src/DDG4SensitiveDetector.cpp b/Simulation/DetSimSD/src/DDG4SensitiveDetector.cpp
index 25ce06c0..9d12dde3 100644
--- a/Simulation/DetSimSD/src/DDG4SensitiveDetector.cpp
+++ b/Simulation/DetSimSD/src/DDG4SensitiveDetector.cpp
@@ -48,7 +48,7 @@ DDG4SensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE) {
}
long long
-DDG4SensitiveDetector::getVolumeID(G4Step* aStep) {
+DDG4SensitiveDetector::getVolumeID(const G4Step* aStep) {
dd4hep::sim::Geant4StepHandler step(aStep);
dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager();
@@ -58,7 +58,7 @@ DDG4SensitiveDetector::getVolumeID(G4Step* aStep) {
}
long long
-DDG4SensitiveDetector::getCellID(G4Step* step) {
+DDG4SensitiveDetector::getCellID(const G4Step* step) {
dd4hep::sim::Geant4StepHandler h(step);
dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager();
dd4hep::VolumeID volID = volMgr.volumeID(h.preTouchable());
diff --git a/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.cpp b/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.cpp
new file mode 100644
index 00000000..511b2b18
--- /dev/null
+++ b/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.cpp
@@ -0,0 +1,57 @@
+#include "TimeProjectionChamberSensDetTool.h"
+
+#include "G4VSensitiveDetector.hh"
+
+#include "DD4hep/Detector.h"
+
+#include "TimeProjectionChamberSensitiveDetector.h"
+
+#include "CLHEP/Units/SystemOfUnits.h"
+
+DECLARE_COMPONENT(TimeProjectionChamberSensDetTool);
+
+StatusCode TimeProjectionChamberSensDetTool::initialize() {
+ StatusCode sc;
+ debug() << "initialize() " << endmsg;
+ debug() << "TypeOption = " << m_sdTypeOption.value() << endmsg;
+
+ m_geosvc = service<IGeomSvc>("GeomSvc");
+ if (!m_geosvc) {
+ error() << "Failed to find GeomSvc." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ return AlgTool::initialize();
+}
+
+StatusCode TimeProjectionChamberSensDetTool::finalize() {
+ StatusCode sc;
+
+ return sc;
+}
+
+G4VSensitiveDetector* TimeProjectionChamberSensDetTool::createSD(const std::string& name) {
+ debug() << "createSD" << endmsg;
+
+ dd4hep::Detector* dd4hep_geo = m_geosvc->lcdd();
+
+ G4VSensitiveDetector* sd = nullptr;
+ if (name == "TPC") {
+ if(m_sdTypeOption==1/*FIXME: const defined by CEPC whole description*/){
+ TimeProjectionChamberSensitiveDetector* tpcsd = new TimeProjectionChamberSensitiveDetector(name, *dd4hep_geo);
+ // switch units to Geant4's, Geant4 uses CLHEP units in fact
+ tpcsd->setThreshold(m_threshold/dd4hep::eV*CLHEP::eV); //FIXME: let SD get threshold from dd4hep_geo object
+ tpcsd->setSameStepLimit(m_sameStepLimit);
+ tpcsd->setWriteMCTruthForLowPtHits(m_writeMCTruthForLowPtHits);
+ tpcsd->setLowPtCut(m_lowPtCut/dd4hep::MeV*CLHEP::MeV);
+ tpcsd->setLowPtMaxHitSeparation(m_lowPtMaxHitSeparation/dd4hep::mm*CLHEP::mm);
+ sd = tpcsd;
+ info() << "TPC will use TimeProjectionChamberSensitiveDetector" << endmsg;
+ }
+ else{
+ info() << "TPC will use default Geant4TrackerHit SD " << endmsg;
+ }
+ }
+
+ return sd;
+}
diff --git a/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.h b/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.h
new file mode 100644
index 00000000..89a8593e
--- /dev/null
+++ b/Simulation/DetSimSD/src/TimeProjectionChamberSensDetTool.h
@@ -0,0 +1,40 @@
+#ifndef TimeProjectionChamberSensDetTool_h
+#define TimeProjectionChamberSensDetTool_h
+
+/*
+ * TimeProjectionChamberSensDetTool is used to create Time Projection Chamber SD.
+ */
+
+#include "GaudiKernel/AlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "DetSimInterface/ISensDetTool.h"
+#include "DetInterface/IGeomSvc.h"
+
+#include "DD4hep/DD4hepUnits.h"
+
+class TimeProjectionChamberSensDetTool: public extends<AlgTool, ISensDetTool> {
+
+ public:
+
+ using extends::extends;
+
+ /// Overriding initialize and finalize
+ StatusCode initialize() override;
+ StatusCode finalize() override;
+
+ /// Override ISensDetTool
+ virtual G4VSensitiveDetector* createSD(const std::string& name) override;
+
+private:
+
+ // in order to initialize SD, we need to get the lcdd()
+ SmartIF<IGeomSvc> m_geosvc;
+ Gaudi::Property<int> m_sdTypeOption{this, "TypeOption", 0};
+ Gaudi::Property<double> m_threshold{this, "ThresholdEnergyDeposit", 32*dd4hep::eV};
+ Gaudi::Property<double> m_lowPtCut{this, "LowPtCut", 10*dd4hep::MeV};
+ Gaudi::Property<double> m_lowPtMaxHitSeparation{this, "LowPtMaxHitSeparation", 5*dd4hep::mm};
+ Gaudi::Property<bool> m_sameStepLimit{this, "SameStepLimit", true};
+ Gaudi::Property<bool> m_writeMCTruthForLowPtHits{this, "WriteMCTruthForLowPtHits", false};
+};
+
+#endif
diff --git a/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.cpp b/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.cpp
new file mode 100644
index 00000000..fdb1ff09
--- /dev/null
+++ b/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.cpp
@@ -0,0 +1,331 @@
+#include "TimeProjectionChamberSensitiveDetector.h"
+
+#include "G4Step.hh"
+#include "G4VProcess.hh"
+//#include "G4SteppingManager.hh"
+#include "G4SDManager.hh"
+//#include "UserTrackInformation.hh"
+#include "DD4hep/DD4hepUnits.h"
+
+TimeProjectionChamberSensitiveDetector::TimeProjectionChamberSensitiveDetector(const std::string& name,
+ dd4hep::Detector& description)
+ : DDG4SensitiveDetector(name, description){
+
+ G4String CollName1=name+"Collection";
+ collectionName.insert(CollName1);
+
+ G4String CollName2=name+"SpacePointCollection";
+ collectionName.insert(CollName2);
+
+ G4String CollName3=name+"LowPtCollection";
+ collectionName.insert(CollName3);
+}
+
+void TimeProjectionChamberSensitiveDetector::Initialize(G4HCofThisEvent* HCE){
+ m_hc = new HitCollection(GetName(), collectionName[0]);
+ int HCID = G4SDManager::GetSDMpointer()->GetCollectionID(m_hc);
+ HCE->AddHitsCollection( HCID, m_hc );
+
+ m_spaceHC = new HitCollection(GetName(), collectionName[1]);
+ int HCIDSpace = G4SDManager::GetSDMpointer()->GetCollectionID(m_spaceHC);
+ HCE->AddHitsCollection( HCIDSpace, m_spaceHC );
+
+ m_lowPtHC = new HitCollection(GetName(), collectionName[2]);
+ int HCIDLowPt = G4SDManager::GetSDMpointer()->GetCollectionID(m_lowPtHC);
+ HCE->AddHitsCollection( HCIDLowPt, m_lowPtHC );
+
+ CumulativeNumSteps=0;
+
+ dEInPadRow = 0.0;
+ CumulativeEnergyDeposit = 0.0;
+ globalTimeAtPadRingCentre=0.0;
+ pathLengthInPadRow=0.0;
+ CumulativePathLength=0.0;
+ CrossingOfPadRingCentre.SetXYZ(0,0,0);// = dd4hep::Position(0,0,0);
+ MomentumAtPadRingCentre.SetXYZ(0,0,0);// = dd4hep::sim::Momentum(0,0,0);
+
+ CumulativeMeanPosition.SetXYZ(0,0,0);// = dd4hep::Position(0,0,0);
+ CumulativeMeanMomentum.SetXYZ(0,0,0);// = dd4hep::sim::Momentum(0,0,0);
+
+ PreviousPostStepPosition.SetXYZ(0,0,0);// = dd4hep::Position(0,0,0);
+ CurrentPDGEncoding=0;
+ CurrentTrackID=-1;
+ CurrentGlobalTime=0.0;
+ CurrentCopyNumber=0;
+
+ /*
+ if( m_lowPtStepLimit ) {
+ std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtStepLimit ON" << std::endl;
+ }
+ else {
+ std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtStepLimit OFF" << std::endl;
+ }
+
+ if( m_writeMCParticleForLowPtHits ) {
+ std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtStoreMCPForHits ON" << std::endl;
+ }
+ else {
+ std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtStoreMCPForHits OFF" << std::endl;
+ }
+
+ std::cout << "TimeProjectionChamberSensitiveDetector: Threshold for Energy Deposit = " << m_thresholdEnergyDeposit / CLHEP::eV << " eV" << G4endl;
+ std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPtCut = " << m_lowPtCut / CLHEP::MeV << " MeV "<< G4endl;
+ std::cout << "TimeProjectionChamberSensitiveDetector: TPCLowPt Max hit separation "<< m_lowPtMaxHitSeparation / CLHEP::mm << " mm" << G4endl;
+ */
+}
+
+G4bool TimeProjectionChamberSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*){
+ // FIXME:
+ // (i) in the following algorithm if a particle "appears" within a pad-ring half and
+ // leaves without passing through the middle of the pad-ring it will not create a hit in
+ // this ring
+ // (ii) a particle that crosses the boundry between two pad-ring halves will have the hit
+ // placed on this surface at the last crossing point, and will be assinged the total energy
+ // deposited in the whole pad-ring. This is a possible source of bias for the hit
+
+ G4TouchableHandle touchPost = step->GetPostStepPoint()->GetTouchableHandle();
+ G4TouchableHandle touchPre = step->GetPreStepPoint()->GetTouchableHandle();
+ dd4hep::sim::Geant4StepHandler h(step);
+
+ if (fabs(h.trackDef()->GetPDGCharge()) < 0.01) return true;
+
+ const dd4hep::Position PrePosition = h.prePos();
+ const dd4hep::Position PostPosition = h.postPos();
+ const dd4hep::sim::Momentum Momentum = h.postMom(); //use post step momentum, Mokka history
+
+ float ptSQRD = Momentum.X()*Momentum.X()+Momentum.Y()*Momentum.Y();
+ if( ptSQRD >= (m_lowPtCut*m_lowPtCut) ){
+ // Step finishes at a geometric boundry
+ if(step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary){
+ //if(h.postStepStatus() == "GeomBoundary"){
+ // step within the same pair of upper and lower pad ring halves
+ if(int(touchPre->GetCopyNumber()/2)==int(touchPost->GetCopyNumber()/2)) {
+ //this step must have ended on the boundry between these two pad ring halfs
+ //record the tracks coordinates at this position
+ //and return
+
+ CrossingOfPadRingCentre.SetXYZ(PostPosition.X(),PostPosition.Y(),PostPosition.Z());
+ MomentumAtPadRingCentre = Momentum;
+ dEInPadRow += h.deposit();
+ globalTimeAtPadRingCentre = h.trkTime();
+ pathLengthInPadRow += h.stepLength();
+
+ // G4cout << "step must have ended on the boundry between these two pad ring halfs" << G4endl;
+ // G4cout << "CrossingOfPadRingCentre = "
+ // << sqrt( CrossingOfPadRingCentre[0]*CrossingOfPadRingCentre[0]
+ // +CrossingOfPadRingCentre[1]*CrossingOfPadRingCentre[1] )
+ // << G4endl;
+
+ return true;
+ }
+ else if(!(CrossingOfPadRingCentre.X()==0.0 && CrossingOfPadRingCentre.Y()==0.0 && CrossingOfPadRingCentre.Z()==0.0)) {
+ // the above IF statment is to catch the case where the particle "appears" in this pad-row half volume and
+ // leaves with out crossing the pad-ring centre, as mentioned above
+
+ //it is leaving the pad ring couplet
+ //write out a hit
+ //make sure particle is added to MC list
+ //and return
+ // G4cout << "step must be leaving the pad ring couplet" << G4endl;
+ // G4cout << "write out hit at "
+ // << sqrt( CrossingOfPadRingCentre[0]*CrossingOfPadRingCentre[0]
+ // +CrossingOfPadRingCentre[1]*CrossingOfPadRingCentre[1] )
+ // << " " << "dEdx = " << step->GetTotalEnergyDeposit()+dEInPadRow
+ // << " " << "step length = " << step->GetStepLength()+pathLengthInPadRow
+ // << G4endl;
+
+ double dE = h.deposit()+dEInPadRow;
+ if ( dE > m_thresholdEnergyDeposit || m_trackingPhysicsListELossOn == false ) {
+ // needed for delta electrons: all particles causing hits have to be saved
+ // TODO: set track flag to save
+
+ dd4hep::sim::Geant4TrackerHit* hit = new dd4hep::sim::Geant4TrackerHit(h.trkID(),h.trkPdgID(),dE,globalTimeAtPadRingCentre);
+ hit->cellID = touchPre->GetCopyNumber()/2+1;//getCellID(step);
+ hit->energyDeposit = dE;
+ hit->position = CrossingOfPadRingCentre;
+ hit->momentum = MomentumAtPadRingCentre;
+ hit->length = h.stepLength()+pathLengthInPadRow;
+ m_hc->insert(hit);
+ }
+
+ // zero cumulative variables
+ dEInPadRow = 0.0;
+ globalTimeAtPadRingCentre = 0.0;
+ pathLengthInPadRow = 0.0;
+ CrossingOfPadRingCentre.SetXYZ(0.,0.,0.);
+ MomentumAtPadRingCentre.SetXYZ(0.,0.,0.);
+ return true;
+ }
+ }
+ //case for which the step remains within geometric volume
+ //FIXME: need and another IF case to catch particles which Stop within the padring
+ else if(step->GetPostStepPoint()->GetStepStatus() != fGeomBoundary){
+ //else if(h.postStepStatus() != "GeomBoundary") {
+ //the step is not limited by the step length
+ if(step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName()!="StepLimiter"){
+ // if(particle not stoped){
+ //add the dEdx and return
+ // G4cout << "Step ended by Physics Process: Add dEdx and carry on" << G4endl;
+ dEInPadRow += h.deposit();
+ pathLengthInPadRow += h.stepLength();
+ return true;
+ //}
+ //else{
+ // write out the hit and clear counters
+ //}
+ }
+ else {
+ //double position_x = PostPosition[0];
+ //double position_y = PostPosition[1];
+ //double position_z = PostPosition[2];
+
+ const dd4hep::sim::Momentum PreMomentum = h.preMom();
+ const dd4hep::sim::Momentum PostMomentum = h.postMom();
+
+ //double momentum_x = PostMomentum[0];
+ //double momentum_y = PostMomentum[1];
+ //double momentum_z = PostMomentum[2];
+
+
+ // G4cout << "step must have been stopped by the step limiter" << G4endl;
+ // G4cout << "write out hit at "
+ // << sqrt( position_x*position_x
+ // +position_y*position_y )
+ // << " " << "dEdx = " << step->GetTotalEnergyDeposit()
+ // << " " << "step length = " << step->GetStepLength()
+ // << G4endl;
+
+ // write out step limited hit
+ // these are just space point hits so do not save the dE, which is set to ZERO
+ // if ( step->GetTotalEnergyDeposit() > fThresholdEnergyDeposit )
+
+ // needed for delta electrons: all particles causing hits have to be saved in the LCIO file
+ // TODO: set track flag to be saved
+ //UserTrackInformation *info = (UserTrackInformation*)(step->GetTrack()->GetUserInformation());
+ //if (info) info->GetTheTrackSummary()->SetToBeSaved();
+
+ dd4hep::sim::Geant4TrackerHit* hit = new dd4hep::sim::Geant4TrackerHit(h.trkID(),h.trkPdgID(),0.0,h.trkTime());
+ hit->cellID = touchPre->GetCopyNumber()/2+1;//getCellID(step);
+ hit->energyDeposit = 0.0;
+ hit->position = PostPosition;
+ hit->momentum = PostMomentum;
+ hit->length = h.stepLength();
+ m_spaceHC->insert(hit);
+
+ // add dE and pathlegth and return
+ dEInPadRow += h.deposit();
+ pathLengthInPadRow += h.stepLength();
+ return true;
+ }
+ }
+ }
+ else if(m_lowPtStepLimit) { // low pt tracks will be treated differently as their step length is limited by the special low pt steplimiter
+ if ( ( PreviousPostStepPosition - PrePosition ).R() > 1.0e-6 * dd4hep::mm ) {
+ // This step does not continue the previous path. Deposit the energy and begin a new Pt hit.
+ if (CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
+ DepositLowPtHit();
+ }
+ else {
+ // reset the cumulative variables if the hit has not been deposited.
+ // The previous track has ended and the cumulated energy left at the end
+ // was not enough to ionize
+ //G4cout << "reset due to new track , discarding " << CumulativeEnergyDeposit / eV << " eV" << std::endl;
+ ResetCumulativeVariables();
+ }
+ }
+ else {
+ //G4cout << "continuing track" << endl;
+ }
+
+ CumulateLowPtStep(h);
+
+ // check whether to deposit the hit
+ if( ( CumulativePathLength > m_lowPtMaxHitSeparation ) ) {
+ // hit is deposited because the step limit is reached and there is enough energy
+ // to ionize
+ if ( CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
+ DepositLowPtHit();
+ }
+ //else {
+ //G4cout << "not deposited, energy is " << CumulativeEnergyDeposit/eV << " eV" << std::endl;
+ //}
+ }
+ else { // even if the step lenth has not been reached the hit might
+ // be deposited because the particle track ends
+ if ( h.post->GetKineticEnergy() == 0 ) {
+ // only deposit the hit if the energy is high enough
+ if (CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
+ DepositLowPtHit();
+ }
+ else { // energy is not enoug to ionize.
+ // However, the track has ended and the energy is discarded and not added to the next step
+ //G4cout << "reset due to end of track, discarding " << CumulativeEnergyDeposit/eV << " eV" << std::endl;
+ ResetCumulativeVariables();
+ }
+ }
+ }
+
+ PreviousPostStepPosition = h.postPos();
+
+ return true;
+ }
+
+ return true;
+}
+
+void TimeProjectionChamberSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE){
+ // There might be one low Pt hit which has not been added to the collection.
+ if (CumulativeEnergyDeposit > m_thresholdEnergyDeposit) {
+ DepositLowPtHit();
+ }
+}
+
+void TimeProjectionChamberSensitiveDetector::DepositLowPtHit(){
+ dd4hep::sim::Geant4TrackerHit* hit = new dd4hep::sim::Geant4TrackerHit(CurrentTrackID,CurrentPDGEncoding,CumulativeEnergyDeposit,CurrentGlobalTime);
+ hit->cellID = CurrentCopyNumber;
+ hit->energyDeposit = CumulativeEnergyDeposit; // FIXME: also use the dedx
+ hit->position = CumulativeMeanPosition;
+ hit->momentum = CumulativeMeanMomentum;
+ hit->length = CumulativePathLength;
+ m_lowPtHC->insert(hit);
+ // reset the cumulative variables after positioning the hit
+ ResetCumulativeVariables();
+}
+
+void TimeProjectionChamberSensitiveDetector::ResetCumulativeVariables(){
+ CumulativeMeanPosition.SetXYZ(0.,0.,0.);// = dd4hep::Position(0.0,0.0,0.0);
+ CumulativeMeanMomentum.SetXYZ(0.,0.,0.);// = dd4hep::sim::Momentum(0.0,0.0,0.0);
+ CumulativeNumSteps = 0;
+ CumulativeEnergyDeposit = 0;
+ CumulativePathLength = 0;
+}
+
+void TimeProjectionChamberSensitiveDetector::CumulateLowPtStep(dd4hep::sim::Geant4StepHandler& h){
+ //dd4hep::Position prePos = h.prePos();
+ //dd4hep::Position postPos = h.postPos();
+ //dd4hep::Position direction = postPos - prePos;
+ dd4hep::Position meanPosition = h.avgPosition();//0.5*(prePos+postPos);
+ //double hit_len = direction.R();
+ dd4hep::Position meanMomentum = 0.5*(h.preMom() + h.postMom());
+
+ ++CumulativeNumSteps;
+ CumulativeMeanPosition = ( (CumulativeMeanPosition*(CumulativeNumSteps-1)) + meanPosition ) / CumulativeNumSteps;
+ CumulativeMeanMomentum = ( (CumulativeMeanMomentum*(CumulativeNumSteps-1)) + meanMomentum ) / CumulativeNumSteps;
+ CumulativeEnergyDeposit += h.deposit();
+ CumulativePathLength += h.stepLength();
+ CurrentPDGEncoding = h.trkPdgID();
+ CurrentTrackID = h.trkID();
+ CurrentGlobalTime = h.trkTime();
+ CurrentCopyNumber = h.preTouchable()->GetCopyNumber()/2+1;//getCellID( h.step );
+
+ // needed for delta electrons: all particles causing hits have to be saved in the LCIO file
+ if ( CumulativeEnergyDeposit > m_thresholdEnergyDeposit ) {
+ // This low Pt hit will eventually be saved, so set the flag to store the particle
+
+ // writing the MC Particles can be turned on or off for the lowPt particles
+ if(m_writeMCParticleForLowPtHits){
+ // TODO: here set write flag for this track
+
+ }
+ }
+}
diff --git a/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.h b/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.h
new file mode 100644
index 00000000..bcb91b43
--- /dev/null
+++ b/Simulation/DetSimSD/src/TimeProjectionChamberSensitiveDetector.h
@@ -0,0 +1,68 @@
+// *********************************************************
+// Implemented from Mokka
+// *********************************************************
+//
+// $Id: TimeProjectionChamberSensitiveDetector.hh,v 1.1 2009/05/13 08:22:57 steve Exp $
+
+#ifndef TimeProjectionChamberSensitiveDetector_h
+#define TimeProjectionChamberSensitiveDetector_h
+
+#include "DetSimSD/DDG4SensitiveDetector.h"
+#include "DDG4/Defs.h"
+
+class TimeProjectionChamberSensitiveDetector: public DDG4SensitiveDetector {
+ public:
+ TimeProjectionChamberSensitiveDetector(const std::string& name, dd4hep::Detector& description);
+
+ void Initialize(G4HCofThisEvent* HCE);
+ G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
+ void EndOfEvent(G4HCofThisEvent* HCE);
+
+ void setThreshold(double e){m_thresholdEnergyDeposit=e;};
+ void setSameStepLimit(bool flag){m_lowPtStepLimit=(!flag);};
+ void setWriteMCTruthForLowPtHits(bool flag){m_writeMCParticleForLowPtHits=flag;};
+ void setLowPtCut(double cut){m_lowPtCut=cut;};
+ void setLowPtMaxHitSeparation(double length){m_lowPtMaxHitSeparation=length;};
+
+ /// helper function to avoid code duplication, writes a low Pt hit to the collection
+ void DepositLowPtHit();
+
+ /// helper function to avoid code duplication, resets all cumulative variables
+ void ResetCumulativeVariables();
+
+ /// helper function to avoid code duplication,
+ /// adds energy, track length and momentum of a low pt step to the cumulative variables
+ void CumulateLowPtStep(dd4hep::sim::Geant4StepHandler& h);
+
+ protected:
+
+ double m_thresholdEnergyDeposit = 0;
+ bool m_lowPtStepLimit = false;
+ bool m_writeMCParticleForLowPtHits = false;
+ double m_lowPtCut = 0;
+ double m_lowPtMaxHitSeparation = 0;
+ bool m_trackingPhysicsListELossOn = true;
+
+ HitCollection* m_hc = nullptr;
+ HitCollection* m_spaceHC = nullptr;
+ HitCollection* m_lowPtHC = nullptr;
+
+ dd4hep::Position CrossingOfPadRingCentre;
+ dd4hep::sim::Momentum MomentumAtPadRingCentre;
+ double dEInPadRow;
+ double globalTimeAtPadRingCentre;
+ double pathLengthInPadRow;
+ double CumulativePathLength;
+ double CumulativeEnergyDeposit;
+ dd4hep::Position CumulativeMeanPosition;
+ dd4hep::sim::Momentum CumulativeMeanMomentum;
+ int CumulativeNumSteps;
+
+ dd4hep::Position PreviousPostStepPosition; //< the end point of the previous step
+ int CurrentPDGEncoding; //< the PDG encoding of the particle causing the cumulative energy deposit
+ int CurrentTrackID; //< the TrackID of the particle causing the cumulative energy deposit
+ double CurrentGlobalTime; ///< the global time of the track causing the cumulative energy deposit
+ int CurrentCopyNumber; ///< copy number of the preStepPoint's TouchableHandle for the cumulative energy deposit
+
+};
+#endif
--
GitLab