diff --git a/Examples/options/tut_detsim_digi_truthTracker_SDT_dedx.py b/Examples/options/tut_detsim_digi_truthTracker_SDT_dedx.py
index 926428b1e5ce27bfad74d01255c5878d86fa4f12..892572df60aac3d093d145135c7872e4506e9409 100644
--- a/Examples/options/tut_detsim_digi_truthTracker_SDT_dedx.py
+++ b/Examples/options/tut_detsim_digi_truthTracker_SDT_dedx.py
@@ -187,12 +187,16 @@ truthTrackerAlg.debug = 1
##############################################################################
# DedxAlg
##############################################################################
+from Configurables import GFDndxSimTool
+dndx_simtool = GFDndxSimTool("GFDndxSimTool")
from Configurables import RecDCHDedxAlg
dedxAlg = RecDCHDedxAlg("RecDCHDedxAlg")
dedxAlg.debug = True
-dedxAlg.method = 1
-dedxAlg.sampling_option = "BetheBlochEquationDedxSimTool"
-dedxAlg.dedx_resolution = 0.05
+#dedxAlg.method = 1
+#dedxAlg.sampling_option = "BetheBlochEquationDedxSimTool"
+dedxAlg.method = 2
+dedxAlg.sampling_option = "GFDndxSimTool"
+dedxAlg.dedx_resolution = 0.0
dedxAlg.DCTrackCollection = "DCTrackCollection"
dedxAlg.DCHitAssociationCollection = "DCHAssociationCollection"
dedxAlg.WriteAna = True
diff --git a/Reconstruction/DCHDedx/src/RecDCHDedxAlg.cpp b/Reconstruction/DCHDedx/src/RecDCHDedxAlg.cpp
index 9025bb88bf618cd4ffd5ec8cab6a06fc321af0a9..eaac6ffe1be6efd69df095a9500ba33562ca3c0d 100644
--- a/Reconstruction/DCHDedx/src/RecDCHDedxAlg.cpp
+++ b/Reconstruction/DCHDedx/src/RecDCHDedxAlg.cpp
@@ -131,6 +131,74 @@ StatusCode RecDCHDedxAlg::execute()
}
}
}//method 1
+ if(m_method==2){// get the primary mc particle of the reconstructed track and do dndx sampling according the mc beta-gamma
+
+ m_dedx_simtool = ToolHandle<IDedxSimTool>(m_dedx_sim_option.value());
+ if (!m_dedx_simtool) {
+ error() << "Failed to find dedx sampling tool :"<<m_dedx_sim_option.value()<< endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ const edm4hep::TrackCollection* trkCol=nullptr;
+ trkCol=m_dcTrackCol.get();
+ if(nullptr==trkCol){
+ throw ("Error: TrackCollection not found.");
+ return StatusCode::SUCCESS;
+ }
+ edm4hep::TrackCollection* ptrkCol = const_cast<edm4hep::TrackCollection*>(trkCol);
+
+ const edm4hep::MCRecoTrackerAssociationCollection* assoCol=nullptr;
+ assoCol = m_dcHitAssociationCol.get();
+ if(nullptr==assoCol){
+ throw ("Error: MCRecoTrackerAssociationCollection not found.");
+ return StatusCode::SUCCESS;
+ }
+ for(unsigned j=0; j<ptrkCol->size(); j++){
+ std::map< edm4hep::ConstMCParticle, int > map_mc_count;
+ auto tmp_track = ptrkCol->at(j);
+ for(unsigned k=0; k< tmp_track.trackerHits_size(); k++){
+ for(unsigned z=0; z< assoCol->size(); z++){
+ if(assoCol->at(z).getRec().id() != tmp_track.getTrackerHits(k).id()) continue;
+
+ auto tmp_mc = assoCol->at(z).getSim().getMCParticle();
+ if(tmp_mc.parents_size()!=0) continue;
+ if(map_mc_count.find(tmp_mc) != map_mc_count.end()) map_mc_count[tmp_mc] += 1;
+ else map_mc_count[tmp_mc] = 1;
+
+ }
+ }
+ edm4hep::MCParticle tmp_mc;
+ int max_cout = 0;
+ for(auto iter = map_mc_count.begin(); iter != map_mc_count.end(); iter++ ){
+ if(iter->second > max_cout){
+ max_cout = iter->second;
+ //tmp_mc = iter->first;
+ tmp_mc.setMass(iter->first.getMass());
+ tmp_mc.setCharge(iter->first.getCharge());
+ tmp_mc.setMomentum(iter->first.getMomentum());
+ if(m_debug){
+ std::cout<<"mass="<<iter->first.getMass()<<",charge="<<iter->first.getCharge()<<",mom_x"<<iter->first.getMomentum()[0]<<",mom_y"<<iter->first.getMomentum()[1]<<",mom_z"<<iter->first.getMomentum()[2]<<std::endl;
+ }
+ }
+ }
+ double dndx = 0.0;
+ double betagamma = tmp_mc.getMass() !=0 ? sqrt(tmp_mc.getMomentum()[0]*tmp_mc.getMomentum()[0] + tmp_mc.getMomentum()[1]*tmp_mc.getMomentum()[1] + tmp_mc.getMomentum()[2]*tmp_mc.getMomentum()[2] )/tmp_mc.getMass() : -1 ;
+ dndx = m_dedx_simtool->dndx(betagamma);
+ float dndx_smear = CLHEP::RandGauss::shoot(m_scale, m_resolution);
+ if(m_debug) std::cout<<"dndx from sampling ="<<dndx<<",smear="<<dndx_smear<<",final="<<dndx*dndx_smear<<",scale="<<m_scale<<",res="<<m_resolution<<std::endl;
+ ptrkCol->at(j).setDEdx(dndx*dndx_smear);//update the dndx
+ if(m_WriteAna){
+ m_track_dedx[m_n_track]= dndx*dndx_smear;
+ m_track_dedx_BB[m_n_track]= dndx;
+ m_track_px [m_n_track]= tmp_mc.getMomentum()[0];
+ m_track_py [m_n_track]= tmp_mc.getMomentum()[1];
+ m_track_pz [m_n_track]= tmp_mc.getMomentum()[2];
+ m_track_mass [m_n_track]= tmp_mc.getMass();
+ m_track_pid [m_n_track]= tmp_mc.getPDG();
+ m_n_track ++;
+ }
+ }
+ }//method 2
if(m_WriteAna){
StatusCode status = m_tuple->write();
if ( status.isFailure() ) {
diff --git a/Simulation/DetSimDedx/CMakeLists.txt b/Simulation/DetSimDedx/CMakeLists.txt
index 858ae456b31863aac77a6e3feb779bb9456338e6..e90963eb2b7a45a1a451fd8d76c047bd3fb6c813 100644
--- a/Simulation/DetSimDedx/CMakeLists.txt
+++ b/Simulation/DetSimDedx/CMakeLists.txt
@@ -14,6 +14,7 @@ include_directories(${EDM4HEP_INCLUDE_DIR})
set(DetSimDedx_srcs
src/DummyDedxSimTool.cpp
src/BetheBlochEquationDedxSimTool.cpp
+ src/GFDndxSimTool.cpp
)
gaudi_add_module(DetSimDedx ${DetSimDedx_srcs}
diff --git a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
index 724b9b7ce8c4033ed3b1721578c71e378a242aa6..4c5ff269eaa4f0bc30bbe935f4ebd43b00f32ad8 100644
--- a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
+++ b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
@@ -49,6 +49,9 @@ double BetheBlochEquationDedxSimTool::dedx(const edm4hep::MCParticle& mcp) {
dedx = dedx*CLHEP::RandGauss::shoot(m_scale, m_resolution);
return dedx; // (CLHEP::MeV/CLHEP::cm) / (CLHEP::g/CLHEP::cm3)
}
+double BetheBlochEquationDedxSimTool::dndx(double betagamma) {
+ return -1;
+}
StatusCode BetheBlochEquationDedxSimTool::initialize()
{
diff --git a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h
index 1979a526d5682eed9fef9d267605160c2d293f46..a4d08870ccfe3634772001d6d4f6136c3c7217f8 100644
--- a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h
+++ b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h
@@ -13,6 +13,7 @@ class BetheBlochEquationDedxSimTool: public extends<AlgTool, IDedxSimTool> {
StatusCode finalize() override;
double dedx(const G4Step* aStep) override;
double dedx(const edm4hep::MCParticle& mc) override;
+ double dndx(double betagamma) override;
private:
diff --git a/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp b/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp
index 60212b7134f9021b05c1aa2f097e87ee7f4ec614..a5c3244ce4f292dac27a2249d0d770f94c15040f 100644
--- a/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp
+++ b/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp
@@ -25,3 +25,6 @@ double DummyDedxSimTool::dedx(const G4Step* aStep) {
double DummyDedxSimTool::dedx(const edm4hep::MCParticle& mc) {
return -1;
}
+double DummyDedxSimTool::dndx(double betagamma) {
+ return -1;
+}
diff --git a/Simulation/DetSimDedx/src/DummyDedxSimTool.h b/Simulation/DetSimDedx/src/DummyDedxSimTool.h
index 879ddaaf4f98a67a8b8c39deaf34555bd33b46f9..58ac1cf63b05518b229b5b3e0fb329767e64d0e7 100644
--- a/Simulation/DetSimDedx/src/DummyDedxSimTool.h
+++ b/Simulation/DetSimDedx/src/DummyDedxSimTool.h
@@ -17,6 +17,7 @@ public:
/// Overriding dedx tool
double dedx(const G4Step* aStep) override;
double dedx(const edm4hep::MCParticle& mc) override;
+ double dndx(double betagamma) override;
};
diff --git a/Simulation/DetSimDedx/src/GFDndxSimTool.cpp b/Simulation/DetSimDedx/src/GFDndxSimTool.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..96592f7fcbe689df63895ecf1b2043baa0a50be7
--- /dev/null
+++ b/Simulation/DetSimDedx/src/GFDndxSimTool.cpp
@@ -0,0 +1,28 @@
+#include "GFDndxSimTool.h"
+
+#include "G4Step.hh"
+
+DECLARE_COMPONENT(GFDndxSimTool);
+
+StatusCode GFDndxSimTool::initialize() {
+ StatusCode sc;
+
+ return sc;
+}
+
+StatusCode GFDndxSimTool::finalize() {
+ StatusCode sc;
+
+ return sc;
+}
+
+double GFDndxSimTool::dedx(const G4Step* aStep) {
+ double result = aStep->GetTotalEnergyDeposit();
+ return result;
+}
+double GFDndxSimTool::dedx(const edm4hep::MCParticle& mc) {
+ return -1;
+}
+double GFDndxSimTool::dndx(double betagamma) {
+ return -999;
+}
diff --git a/Simulation/DetSimDedx/src/GFDndxSimTool.h b/Simulation/DetSimDedx/src/GFDndxSimTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..3fbc2a8375544a3cdf90563c138f17f455795635
--- /dev/null
+++ b/Simulation/DetSimDedx/src/GFDndxSimTool.h
@@ -0,0 +1,24 @@
+#ifndef GFDndxSimTool_h
+#define GFDndxSimTool_h
+
+#include "GaudiKernel/AlgTool.h"
+#include "DetSimInterface/IDedxSimTool.h"
+#include "edm4hep/MCParticle.h"
+
+class GFDndxSimTool: public extends<AlgTool, IDedxSimTool> {
+
+public:
+ using extends::extends;
+
+ /// Overriding initialize and finalize
+ StatusCode initialize() override;
+ StatusCode finalize() override;
+
+ /// Overriding dedx tool
+ double dedx(const G4Step* aStep) override;
+ double dedx(const edm4hep::MCParticle& mc) override;
+ double dndx(double betagamma) override;
+
+};
+
+#endif
diff --git a/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h b/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h
index 331d2c458772656c3b7b92c1a4a9d55d2c8cc3fe..8c804500a7b4322e1998a04a7a8a43f46ba3f499 100644
--- a/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h
+++ b/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h
@@ -27,6 +27,7 @@ public:
virtual double dedx(const G4Step* aStep) = 0;
virtual double dedx(const edm4hep::MCParticle& mc) = 0;
+ virtual double dndx(double betagamma) = 0;
};