diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
index d7c2949b4416a5d5e62809a620c3cb33f1b4d274..36783972a794fab389b648b8b36f772a3cf7ff06 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
@@ -8,6 +8,7 @@
#include "DD4hep/Detector.h"
#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
#include "DDRec/Vector3D.h"
#include "GaudiKernel/INTupleSvc.h"
@@ -109,9 +110,10 @@ StatusCode DCHDigiAlg::execute()
id_hits_map[id] = vhit ;
}
}
-
- m_n_sim = 0;
- m_n_digi = 0 ;
+ if(m_WriteAna){
+ m_n_sim = 0;
+ m_n_digi = 0 ;
+ }
for(std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> >::iterator iter = id_hits_map.begin(); iter != id_hits_map.end(); iter++)
{
unsigned long long wcellid = iter->first;
@@ -133,8 +135,10 @@ StatusCode DCHDigiAlg::execute()
TVector3 Wstart(0,0,0);
TVector3 Wend (0,0,0);
m_segmentation->cellposition(wcellid, Wstart, Wend);
- Wstart = 10*Wstart;// from DD4HEP cm to mm
- Wend = 10*Wend ;
+ float dd4hep_mm = dd4hep::mm;
+ //std::cout<<"dd4hep_mm="<<dd4hep_mm<<std::endl;
+ Wstart =(1/dd4hep_mm)* Wstart;// from DD4HEP cm to mm
+ Wend =(1/dd4hep_mm)* Wend ;
//std::cout<<"wcellid="<<wcellid<<",chamber="<<chamber<<",layer="<<layer<<",cellID="<<cellID<<",s_x="<<Wstart.x()<<",s_y="<<Wstart.y()<<",s_z="<<Wstart.z()<<",E_x="<<Wend.x()<<",E_y="<<Wend.y()<<",E_z="<<Wend.z()<<std::endl;
TVector3 denominator = (Wend-Wstart) ;
diff --git a/Examples/options/tut_detsim_digi_truthTracker_SDT_dedx.py b/Examples/options/tut_detsim_digi_truthTracker_SDT_dedx.py
new file mode 100644
index 0000000000000000000000000000000000000000..926428b1e5ce27bfad74d01255c5878d86fa4f12
--- /dev/null
+++ b/Examples/options/tut_detsim_digi_truthTracker_SDT_dedx.py
@@ -0,0 +1,218 @@
+#!/usr/bin/env python
+
+import os
+print(os.environ["DD4HEP_LIBRARY_PATH"])
+import sys
+# sys.exit(0)
+
+from Gaudi.Configuration import *
+
+##############################################################################
+# Random Number Svc
+##############################################################################
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Geant4
+rndmengine.SetSingleton = True
+rndmengine.Seeds = [42]
+
+# rndmgensvc = RndmGenSvc("RndmGenSvc")
+# rndmgensvc.Engine = rndmengine.name()
+
+
+##############################################################################
+# Event Data Svc
+##############################################################################
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc")
+
+
+##############################################################################
+# Geometry Svc
+##############################################################################
+
+# geometry_option = "CepC_v4-onlyTracker.xml"
+geometry_option = "det.xml"
+
+if not os.getenv("DETDRIFTCHAMBERROOT"):
+ print("Can't find the geometry. Please setup envvar DETCEPCV4ROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETDRIFTCHAMBERROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+#geosvc.compact = "/workfs/higgs/fangwx/fork_DedxAlg_20201217/CEPCSW/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01.xml"
+#geosvc.compact = "/workfs/higgs/fangwx/fork_DedxAlg_20201217/CEPCSW/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01_DCH.xml"
+##############################################################################
+# Physics Generator
+##############################################################################
+from Configurables import GenAlgo
+from Configurables import GtGunTool
+from Configurables import StdHepRdr
+from Configurables import SLCIORdr
+from Configurables import HepMCRdr
+from Configurables import GenPrinter
+
+gun = GtGunTool("GtGunTool")
+# gun.Particles = ["pi+"]
+# gun.EnergyMins = [100.] # GeV
+# gun.EnergyMaxs = [100.] # GeV
+
+
+# gun.PositionXs = [100.] # mm
+# gun.PositionYs = [100.] # mm
+# gun.PositionZs = [0.] # mm
+
+gun.Particles = ["K-"]
+
+gun.EnergyMins = [0.1] # GeV
+gun.EnergyMaxs = [20] # GeV
+
+gun.ThetaMins = [90] # rad; 45deg
+gun.ThetaMaxs = [90] # rad; 45deg
+
+gun.PhiMins = [90] # rad; 0deg
+gun.PhiMaxs = [90] # rad; 360deg
+#gun.PhiMins = [0] # rad; 0deg
+#gun.PhiMaxs = [360] # rad; 360deg
+
+# stdheprdr = StdHepRdr("StdHepRdr")
+# stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
+
+# lciordr = SLCIORdr("SLCIORdr")
+# lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio"
+
+# hepmcrdr = HepMCRdr("HepMCRdr")
+# hepmcrdr.Input = "example_UsingIterators.txt"
+
+genprinter = GenPrinter("GenPrinter")
+
+genalg = GenAlgo("GenAlgo")
+genalg.GenTools = ["GtGunTool"]
+# genalg.GenTools = ["StdHepRdr"]
+# genalg.GenTools = ["StdHepRdr", "GenPrinter"]
+# genalg.GenTools = ["SLCIORdr", "GenPrinter"]
+# genalg.GenTools = ["HepMCRdr", "GenPrinter"]
+
+##############################################################################
+# Detector Simulation
+##############################################################################
+from Configurables import DetSimSvc
+
+detsimsvc = DetSimSvc("DetSimSvc")
+
+# from Configurables import ExampleAnaElemTool
+# example_anatool = ExampleAnaElemTool("ExampleAnaElemTool")
+
+from Configurables import DetSimAlg
+
+detsimalg = DetSimAlg("DetSimAlg")
+#detsimalg.RunMacs = ["Examples/options/noDecay.mac"]
+#detsimalg.RunCmds = ["Examples/options/noDecay.mac"]
+
+if int(os.environ.get("VIS", 0)):
+ detsimalg.VisMacs = ["vis.mac"]
+
+detsimalg.RunCmds = [
+# "/tracking/verbose 1",
+]
+detsimalg.AnaElems = [
+ # example_anatool.name()
+ # "ExampleAnaElemTool"
+ "Edm4hepWriterAnaElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+
+from Configurables import AnExampleDetElemTool
+example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
+
+from Configurables import CalorimeterSensDetTool
+from Configurables import DriftChamberSensDetTool
+
+calo_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
+driftchamber_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+
+# dedxoption = "DummyDedxSimTool"
+dedxoption = "BetheBlochEquationDedxSimTool"
+
+driftchamber_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+ dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+ dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+ dedx_simtool.material_Z = 7 # approximate to Air
+ dedx_simtool.material_A = 14
+ dedx_simtool.scale = 1
+ dedx_simtool.resolution = 0
+
+##############################################################################
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = ["MyTuples DATAFILE='Dedx_ana.root' OPT='NEW' TYP='ROOT'"]
+
+##############################################################################
+# DCHDigiAlg
+##############################################################################
+from Configurables import DCHDigiAlg
+dCHDigiAlg = DCHDigiAlg("DCHDigiAlg")
+dCHDigiAlg.readout = "DriftChamberHitsCollection"
+dCHDigiAlg.drift_velocity = 40#um/ns
+dCHDigiAlg.mom_threshold = 0 #GeV
+dCHDigiAlg.SimDCHitCollection = "DriftChamberHitsCollection"
+dCHDigiAlg.DigiDCHitCollection = "DigiDCHitsCollection"
+dCHDigiAlg.AssociationCollection = "DCHAssociationCollection"
+dCHDigiAlg.WriteAna = False
+
+##############################################################################
+# TruthTrackerAlg
+##############################################################################
+from Configurables import TruthTrackerAlg
+truthTrackerAlg = TruthTrackerAlg("TruthTrackerAlg")
+truthTrackerAlg.MCParticle= "MCParticle"
+truthTrackerAlg.DigiDCHitCollection="DigiDCHitsCollection"
+truthTrackerAlg.DCHitAssociationCollection="DCHAssociationCollection"
+truthTrackerAlg.DCTrackCollection="DCTrackCollection"
+truthTrackerAlg.debug = 1
+
+##############################################################################
+# DedxAlg
+##############################################################################
+from Configurables import RecDCHDedxAlg
+dedxAlg = RecDCHDedxAlg("RecDCHDedxAlg")
+dedxAlg.debug = True
+dedxAlg.method = 1
+dedxAlg.sampling_option = "BetheBlochEquationDedxSimTool"
+dedxAlg.dedx_resolution = 0.05
+dedxAlg.DCTrackCollection = "DCTrackCollection"
+dedxAlg.DCHitAssociationCollection = "DCHAssociationCollection"
+dedxAlg.WriteAna = True
+##############################################################################
+# POD I/O
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "DedxAlg.root"
+out.outputCommands = ["keep *"]
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, truthTrackerAlg, dedxAlg],
+ EvtSel = 'NONE',
+ EvtMax = 10,
+ ExtSvc = [rndmengine, dsvc, geosvc],
+ HistogramPersistency = "ROOT",
+ OutputLevel=INFO
+)
diff --git a/Reconstruction/DCHDedx/CMakeLists.txt b/Reconstruction/DCHDedx/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3167363915e71288cd35e089f94e250e7c49da8f
--- /dev/null
+++ b/Reconstruction/DCHDedx/CMakeLists.txt
@@ -0,0 +1,26 @@
+gaudi_subdir(DCHDedx v0r0)
+
+find_package(EDM4HEP REQUIRED )
+find_package(DD4hep COMPONENTS DDCore DDRec REQUIRED)
+
+
+gaudi_depends_on_subdirs(
+ k4FWCore
+ Simulation/DetSimInterface
+ Detector/DetInterface
+ Detector/DetSegmentation
+)
+
+set(srcs
+ src/*.cpp
+)
+
+# Modules
+gaudi_add_module(DCHDedx ${srcs}
+ INCLUDE_DIRS GaudiKernel
+ LINK_LIBRARIES GaudiAlgLib GaudiKernel
+ DetSegmentation
+ -Wl,--no-as-needed
+ EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+ -Wl,--as-needed
+)
diff --git a/Reconstruction/DCHDedx/src/RecDCHDedxAlg.cpp b/Reconstruction/DCHDedx/src/RecDCHDedxAlg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9025bb88bf618cd4ffd5ec8cab6a06fc321af0a9
--- /dev/null
+++ b/Reconstruction/DCHDedx/src/RecDCHDedxAlg.cpp
@@ -0,0 +1,187 @@
+#include "RecDCHDedxAlg.h"
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoParticleAssociationCollection.h"
+#include "edm4hep/ReconstructedParticleCollection.h"
+#include "DDSegmentation/Segmentation.h"
+#include "DetSegmentation/GridDriftChamber.h"
+#include "DetInterface/IGeomSvc.h"
+//#include "edm4hep/SimCalorimeterHitCollection.h"
+//#include "edm4hep/CaloHitContributionCollection.h"
+
+#include "DD4hep/Detector.h"
+#include "DD4hep/IDDescriptor.h"
+#include "DD4hep/Plugins.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "CLHEP/Random/RandGauss.h"
+
+DECLARE_COMPONENT(RecDCHDedxAlg)
+
+RecDCHDedxAlg::RecDCHDedxAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc)
+{
+ declareProperty("DCHitAssociationCollection", m_dcHitAssociationCol, "Handle of association collection");
+ declareProperty("DCTrackCollection", m_dcTrackCol,"Handle of input Track collection");
+}
+
+StatusCode RecDCHDedxAlg::initialize()
+{
+ if(m_WriteAna){
+
+ NTuplePtr nt( ntupleSvc(), "MyTuples/Dedx_evt" );
+ if ( nt ) m_tuple = nt;
+ else {
+ m_tuple = ntupleSvc()->book( "MyTuples/Dedx_evt", CLID_ColumnWiseTuple, "Dedx" );
+ if ( m_tuple ) {
+ m_tuple->addItem( "N_track" , m_n_track, 0, 1000 ).ignore();
+ m_tuple->addItem( "track_dedx" , m_n_track, m_track_dedx ).ignore();
+ m_tuple->addItem( "track_dedx_BB" , m_n_track, m_track_dedx_BB ).ignore();
+ m_tuple->addItem( "track_px" , m_n_track, m_track_px ).ignore();
+ m_tuple->addItem( "track_py" , m_n_track, m_track_py ).ignore();
+ m_tuple->addItem( "track_pz" , m_n_track, m_track_pz ).ignore();
+ m_tuple->addItem( "track_mass" , m_n_track, m_track_mass ).ignore();
+ m_tuple->addItem( "track_pid" , m_n_track, m_track_pid ).ignore();
+ }
+ else { // did not manage to book the N tuple....
+ error() << " Cannot book N-tuple:" << long( m_tuple ) << endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ }
+
+ return GaudiAlgorithm::initialize();
+}
+
+
+StatusCode RecDCHDedxAlg::execute()
+{
+
+ if(m_WriteAna){
+ m_n_track = 0;
+ }
+ if(m_method==1){// get the primary mc particle of the reconstructed track and do dedx sampling using beta-gamma curve
+
+ 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 dedx_ori = tmp_track.getDEdx();
+ double dedx = 0.0;
+ dedx = m_dedx_simtool->dedx(tmp_mc);
+ float dedx_smear = CLHEP::RandGauss::shoot(m_scale, m_resolution);
+ if(m_debug) std::cout<<"ori dedx="<<dedx_ori<<", dedx from sampling ="<<dedx<<",smear="<<dedx_smear<<",final="<<dedx*dedx_smear<<",scale="<<m_scale<<",res="<<m_resolution<<std::endl;
+ ptrkCol->at(j).setDEdx(dedx*dedx_smear);
+ if(m_WriteAna){
+ m_track_dedx[m_n_track]= dedx*dedx_smear;
+ m_track_dedx_BB[m_n_track]= dedx;
+ 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 1
+ if(m_WriteAna){
+ StatusCode status = m_tuple->write();
+ if ( status.isFailure() ) {
+ error() << " Cannot fill N-tuple:" << long( m_tuple ) << endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode RecDCHDedxAlg::finalize()
+{
+ return GaudiAlgorithm::finalize();
+}
+
+double RecDCHDedxAlg::cal_dedx_bitrunc(float truncate, std::vector<double> phlist, int & usedhit )
+{
+ sort(phlist.begin(),phlist.end());
+ int nsampl = (int)( phlist.size()*truncate );
+ int smpl = (int)(phlist.size()*(truncate+0.05));
+ int min_cut = (int)( phlist.size()*0.05 + 0.5 );
+ double qSum = 0;
+ unsigned i = 0;
+ for(std::vector<double>::iterator ql= phlist.begin();ql!=phlist.end();ql++)
+ {
+ i++;
+ if(i<= smpl && i>=min_cut ) qSum += (*ql);
+ }
+ double trunc=-99;
+ usedhit = smpl-min_cut+1;
+ if(usedhit>0) trunc=qSum/usedhit;
+
+ return trunc;
+}
+/*
+double RecDCHDedxAlg::BetheBlochEquationDedx(const edm4hep::MCParticle& mcp)
+{
+ int z = mcp.getCharge();
+ if (z == 0) return 0;
+ float m_I = m_material_Z*10; // Approximate
+
+ double M = mcp.getMass();
+ double gammabeta=sqrt(mcp.getMomentum()[0]*mcp.getMomentum()[0]+mcp.getMomentum()[1]*mcp.getMomentum()[1]+mcp.getMomentum()[2]*mcp.getMomentum()[2])/M;
+ if(gammabeta<0.01)return 0;//too low momentum
+ float beta = gammabeta/sqrt(1.0+pow(gammabeta,2));
+ float gamma = gammabeta/beta;
+ M = M*pow(10,9);//to eV
+ float Tmax = 2*m_me*pow(gammabeta,2)/(1+(2*gamma*m_me/M)+pow(m_me/M,2));
+ float dedx = m_K*pow(z,2)*m_material_Z*(0.5*log(2*m_me*pow(gammabeta,2)*Tmax/pow(m_I,2))-pow(beta,2))/(m_material_A*pow(beta,2));
+ dedx = dedx*CLHEP::RandGauss::shoot(m_scale, m_resolution);
+ return dedx; // (CLHEP::MeV/CLHEP::cm) / (CLHEP::g/CLHEP::cm3)
+}
+*/
diff --git a/Reconstruction/DCHDedx/src/RecDCHDedxAlg.h b/Reconstruction/DCHDedx/src/RecDCHDedxAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2a5b9f58775ffe47675741184a908eaa460b76e
--- /dev/null
+++ b/Reconstruction/DCHDedx/src/RecDCHDedxAlg.h
@@ -0,0 +1,75 @@
+#ifndef RecDCHDedxAlg_h
+#define RecDCHDedxAlg_h
+
+#include "edm4hep/MCParticleCollection.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/NTuple.h"
+#include "k4FWCore/DataHandle.h"
+#include "DD4hep/Fields.h"
+#include "DetSimInterface/IDedxSimTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+class IGeomSvc;
+namespace dd4hep {
+ class Detector;
+ namespace DDSegmentation{
+ class GridDriftChamber;
+ class BitFieldCoder;
+ }
+}
+namespace edm4hep {
+ class MCParticleCollection;
+ class TrackerHitCollection;
+ class TrackCollection;
+ class MCRecoTrackerAssociationCollection;
+ class ReconstructedParticleCollection;
+ class MCRecoParticleAssociationCollection;
+}
+
+class RecDCHDedxAlg: public GaudiAlgorithm
+{
+ public:
+ RecDCHDedxAlg(const std::string& name, ISvcLocator* svcLoc);
+
+ virtual StatusCode initialize();
+ virtual StatusCode execute();
+ virtual StatusCode finalize();
+ double cal_dedx_bitrunc(float truncate, std::vector<double> phlist, int & usedhit );
+
+ private:
+ ToolHandle<IDedxSimTool> m_dedx_simtool;
+ Gaudi::Property<std::string> m_dedx_sim_option{ this, "sampling_option", ""};
+
+ //reader
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_dcHitAssociationCol{ "DCHitAssociationCollection",Gaudi::DataHandle::Reader, this};
+ //writer
+ DataHandle<edm4hep::TrackCollection> m_dcTrackCol{"DCTrackCollection", Gaudi::DataHandle::Writer, this};
+
+ Gaudi::Property<int> m_debug{ this, "debug", false};
+ Gaudi::Property<int> m_method{ this, "method", 1};
+ Gaudi::Property<float> m_truncate{ this, "truncate", 0.7};
+ Gaudi::Property<bool> m_WriteAna{ this, "WriteAna", false};
+ Gaudi::Property<float> m_scale{this, "dedx_scale", 1};
+ Gaudi::Property<float> m_resolution{this, "dedx_resolution", 0.};
+
+ NTuple::Tuple* m_tuple = nullptr ;
+ NTuple::Item<long> m_n_track;
+ NTuple::Item<long> m_hit;
+ NTuple::Array<double> m_track_dedx;
+ NTuple::Array<double> m_track_dedx_BB;
+ NTuple::Array<double> m_track_px;
+ NTuple::Array<double> m_track_py;
+ NTuple::Array<double> m_track_pz;
+ NTuple::Array<double> m_track_mass;
+ NTuple::Array<int> m_track_pid;
+ NTuple::Array<double> m_hit_x ;
+ NTuple::Array<double> m_hit_y ;
+ NTuple::Array<double> m_hit_z ;
+ NTuple::Array<double> m_hit_dedx;
+ NTuple::Item<double> m_mc_px;
+ NTuple::Item<double> m_mc_py;
+ NTuple::Item<double> m_mc_pz;
+
+};
+
+#endif
diff --git a/Simulation/DetSimDedx/CMakeLists.txt b/Simulation/DetSimDedx/CMakeLists.txt
index 69e83034f8d6224856e93a44fed8e9badc955c7c..858ae456b31863aac77a6e3feb779bb9456338e6 100644
--- a/Simulation/DetSimDedx/CMakeLists.txt
+++ b/Simulation/DetSimDedx/CMakeLists.txt
@@ -8,6 +8,8 @@ gaudi_depends_on_subdirs(
find_package(Geant4 REQUIRED ui_all vis_all)
include(${Geant4_USE_FILE})
find_package(DD4hep COMPONENTS DDG4 REQUIRED)
+find_package(EDM4HEP REQUIRED )
+include_directories(${EDM4HEP_INCLUDE_DIR})
set(DetSimDedx_srcs
src/DummyDedxSimTool.cpp
@@ -20,5 +22,7 @@ gaudi_add_module(DetSimDedx ${DetSimDedx_srcs}
DD4hep
${DD4hep_COMPONENT_LIBRARIES}
GaudiKernel
+ -Wl,--no-as-needed
+ EDM4HEP::edm4hep EDM4HEP::edm4hepDict
)
diff --git a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
index f8dab312c8eb867116cb537732976803a47a3575..724b9b7ce8c4033ed3b1721578c71e378a242aa6 100644
--- a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
+++ b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
@@ -28,8 +28,26 @@ double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep)
float Tmax = 2*m_me*pow(gammabeta,2)/(1+(2*gamma*m_me/M)+pow(m_me/M,2));
float dedx = m_K*pow(z,2)*material_Z*(0.5*log(2*m_me*pow(gammabeta,2)*Tmax/pow(m_I,2))-pow(beta,2))/(material_A*pow(beta,2));
dedx = dedx*CLHEP::RandGauss::shoot(m_scale, m_resolution);
- double Dedx = dedx * (CLHEP::MeV/CLHEP::cm) / (CLHEP::g/CLHEP::cm3) ;
- return Dedx*material_density;
+ double Dedx = dedx * (CLHEP::MeV/CLHEP::cm) ;
+ return Dedx;
+}
+
+double BetheBlochEquationDedxSimTool::dedx(const edm4hep::MCParticle& mcp) {
+
+ int z = mcp.getCharge();
+ if (z == 0) return 0;
+ float m_I = m_material_Z*10; // Approximate
+
+ double M = mcp.getMass();
+ double gammabeta=sqrt(mcp.getMomentum()[0]*mcp.getMomentum()[0]+mcp.getMomentum()[1]*mcp.getMomentum()[1]+mcp.getMomentum()[2]*mcp.getMomentum()[2])/M;
+ if(gammabeta<0.01)return 0;//too low momentum
+ float beta = gammabeta/sqrt(1.0+pow(gammabeta,2));
+ float gamma = gammabeta/beta;
+ M = M*pow(10,9);//to eV
+ float Tmax = 2*m_me*pow(gammabeta,2)/(1+(2*gamma*m_me/M)+pow(m_me/M,2));
+ float dedx = m_K*pow(z,2)*m_material_Z*(0.5*log(2*m_me*pow(gammabeta,2)*Tmax/pow(m_I,2))-pow(beta,2))/(m_material_A*pow(beta,2));
+ dedx = dedx*CLHEP::RandGauss::shoot(m_scale, m_resolution);
+ return dedx; // (CLHEP::MeV/CLHEP::cm) / (CLHEP::g/CLHEP::cm3)
}
StatusCode BetheBlochEquationDedxSimTool::initialize()
diff --git a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h
index dfb671e3d9ddb692bcc4ea5eb63f2d5a9f06266d..1979a526d5682eed9fef9d267605160c2d293f46 100644
--- a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h
+++ b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h
@@ -3,6 +3,7 @@
#include "DetSimInterface/IDedxSimTool.h"
#include <GaudiKernel/AlgTool.h>
+#include "edm4hep/MCParticle.h"
class BetheBlochEquationDedxSimTool: public extends<AlgTool, IDedxSimTool> {
public:
@@ -11,6 +12,7 @@ class BetheBlochEquationDedxSimTool: public extends<AlgTool, IDedxSimTool> {
StatusCode initialize() override;
StatusCode finalize() override;
double dedx(const G4Step* aStep) override;
+ double dedx(const edm4hep::MCParticle& mc) override;
private:
@@ -18,7 +20,7 @@ class BetheBlochEquationDedxSimTool: public extends<AlgTool, IDedxSimTool> {
Gaudi::Property<float> m_material_A{this, "material_A", 4};
Gaudi::Property<float> m_material_density{this, "material_density", 0.000178};//g/cm^3
Gaudi::Property<float> m_scale{this, "scale", 1};
- Gaudi::Property<float> m_resolution{this, "resolution", 0.01};
+ Gaudi::Property<float> m_resolution{this, "resolution", 0};
float m_me;// Here me is the electron rest mass
float m_K; // K was set as a constant.
float m_I; // Mean excitation energy
diff --git a/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp b/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp
index 291b8c205eafd630a9e45f9be64a99cbd0827c87..60212b7134f9021b05c1aa2f097e87ee7f4ec614 100644
--- a/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp
+++ b/Simulation/DetSimDedx/src/DummyDedxSimTool.cpp
@@ -22,3 +22,6 @@ double DummyDedxSimTool::dedx(const G4Step* aStep) {
return result;
}
+double DummyDedxSimTool::dedx(const edm4hep::MCParticle& mc) {
+ return -1;
+}
diff --git a/Simulation/DetSimDedx/src/DummyDedxSimTool.h b/Simulation/DetSimDedx/src/DummyDedxSimTool.h
index ce7f3e3e4a49430a2f9736e03c4538c2cb87878f..879ddaaf4f98a67a8b8c39deaf34555bd33b46f9 100644
--- a/Simulation/DetSimDedx/src/DummyDedxSimTool.h
+++ b/Simulation/DetSimDedx/src/DummyDedxSimTool.h
@@ -3,6 +3,7 @@
#include "GaudiKernel/AlgTool.h"
#include "DetSimInterface/IDedxSimTool.h"
+#include "edm4hep/MCParticle.h"
class DummyDedxSimTool: public extends<AlgTool, IDedxSimTool> {
@@ -15,6 +16,7 @@ public:
/// Overriding dedx tool
double dedx(const G4Step* aStep) override;
+ double dedx(const edm4hep::MCParticle& mc) override;
};
diff --git a/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h b/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h
index f47ceb0211c3a61553adca623e443f25e46c5038..331d2c458772656c3b7b92c1a4a9d55d2c8cc3fe 100644
--- a/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h
+++ b/Simulation/DetSimInterface/DetSimInterface/IDedxSimTool.h
@@ -15,6 +15,9 @@
#include "GaudiKernel/IAlgTool.h"
class G4Step;
+namespace edm4hep{
+ class MCParticle;
+}
class IDedxSimTool: virtual public IAlgTool {
public:
@@ -23,6 +26,7 @@ public:
virtual ~IDedxSimTool() {}
virtual double dedx(const G4Step* aStep) = 0;
+ virtual double dedx(const edm4hep::MCParticle& mc) = 0;
};