Merge branch 'dndx' into 'master'

dN/dx for TPC and DCH

See merge request cepc/CEPCSW!33

Detector/DetCRD/scripts/TDR_o1_v01/tracking.py

@@ -39,6 +39,10 @@ gearsvc = GearSvc("GearSvc")
 from Configurables import TrackSystemSvc
 tracksystemsvc = TrackSystemSvc("TrackSystemSvc")
 
+from Configurables import SimplePIDSvc
+pidsvc = SimplePIDSvc("SimplePIDSvc")
+pidsvc.ParFile = "/cefs/higgs/zhaog/dndx_files/dNdx_TPC.root"
+
 from Configurables import PodioInput
 podioinput = PodioInput("PodioReader", collections=[
 #    "EventHeader",
@@ -207,10 +211,15 @@ full.SETHitToTrackDistance = 5.
 full.MinChi2ProbForSiliconTracks = 0
 #full.OutputLevel = DEBUG
 
+from Configurables import TPCDndxAlg
+tpc_dndx = TPCDndxAlg("TPCDndxAlg")
+tpc_dndx.Method = "Simple"
+
 from Configurables import TrackParticleRelationAlg
 tpr = TrackParticleRelationAlg("Track2Particle")
 tpr.MCParticleCollection = "MCParticle"
-tpr.TrackList = ["CompleteTracks"]
+tpr.TrackList = ["CompleteTracks", "ClupatraTracks"]
+tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
 #tpr.OutputLevel = DEBUG
 
 # output
@@ -222,10 +231,10 @@ out.outputCommands = ["keep *"]
 # ApplicationMgr
 from Configurables import ApplicationMgr
 mgr = ApplicationMgr(
-    TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, tracking, forward, subset, clupatra, full, tpr, out],
+    TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, out],
     EvtSel = 'NONE',
     EvtMax = 5,
-    ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc],
+    ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc],
     HistogramPersistency = 'ROOT',
     OutputLevel = ERROR
 )

Reconstruction/CMakeLists.txt

@@ -5,4 +5,5 @@ add_subdirectory(SiliconTracking)
 add_subdirectory(Tracking)
 add_subdirectory(RecGenfitAlg)
 add_subdirectory(RecAssociationMaker)
+add_subdirectory(ParticleID)
 add_subdirectory(TofRecAlg)

Reconstruction/ParticleID/CMakeLists.txt

+# gaudi_add_header_only_library(ParticleIDLib)
+
+# Modules
+gaudi_add_module(ParticleID
+                 SOURCES src/TPCDndxAlg.cpp
+                         src/DCHDndxAlg.cpp
+                 LINK SimplePIDSvc
+                      Gaudi::GaudiAlgLib
+                      Gaudi::GaudiKernel
+                      DataHelperLib
+                      DetSegmentation
+                      DetInterface
+                      ${GSL_LIBRARIES}
+                      ${GEAR_LIBRARIES}
+                      ${LCIO_LIBRARIES}
+                      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+                      k4FWCore::k4FWCore
+)
+
+target_include_directories(ParticleID PUBLIC
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include
+  $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
+install(TARGETS ParticleID
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)

Reconstruction/ParticleID/src/DCHDndxAlg.cpp

+#include "DCHDndxAlg.h"
+#include "DataHelper/HelixClass.h"
+#include "DataHelper/Navigation.h"
+#include "DetInterface/IGeomSvc.h"
+#include "UTIL/ILDConf.h"
+
+#include "DD4hep/Detector.h"
+#include "edm4hep/Hypothesis.h"
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/Quantity.h"
+#include "edm4hep/Vector3d.h"
+#include "lcio.h"
+
+#include <cmath>
+#include <fstream>
+
+using namespace edm4hep;
+
+DECLARE_COMPONENT(DCHDndxAlg)
+
+DCHDndxAlg::DCHDndxAlg(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
+    // Input
+    declareProperty("SDTRecTrackCollection", _trackCol, "handler of the input track collection");
+    declareProperty("SDTRecTrackCollectionParticleAssociation", _trkParAssCol, "handler of the input track particle association collection");
+
+    // Output
+    declareProperty("DndxTracks", _dndxCol, "handler of the collection of dN/dx tracks");
+}
+
+StatusCode DCHDndxAlg::initialize() {
+    info() << "Initilize DndxAlg ..." << endmsg;
+
+    if (m_method == "Simple") {
+        m_pid_svc = service("SimplePIDSvc");
+    }
+    else {
+        m_pid_svc = nullptr;
+    }
+
+    m_geom_svc = service("GeomSvc");
+
+    return GaudiAlgorithm::initialize();
+}
+
+StatusCode DCHDndxAlg::execute() {
+    info() << "Dndx reconstruction started" << endmsg;
+    // static std::ofstream check_file("log.txt");
+    // static int ievent = 0;
+    // check_file << "--------                 Event " << ievent++ << " --------" << std::endl;
+
+    const edm4hep::TrackCollection* trkcol = nullptr;
+    const edm4hep::MCRecoTrackParticleAssociationCollection* trkparasscol = nullptr;
+    try {
+        trkcol = _trackCol.get();
+        trkparasscol = _trkParAssCol.get();
+    }
+    catch(...) {
+        //
+    }
+
+    if (trkcol == nullptr || trkparasscol == nullptr) {
+        return StatusCode::SUCCESS;
+    }
+
+    RecDqdxCollection* outCol = _dndxCol.createAndPut();
+
+    // check_file << "before track loop" << std::endl;
+
+    for (std::size_t i = 0; i < trkcol->size(); i++) {
+        // check_file << "  Track " << i << std::endl;
+        Track trk(trkcol->at(i));
+        
+        /// MC truth information
+        int pdgid = -999;
+        double p_truth = -999.;
+        
+        float max_weight = -999.;
+        int max_weight_idx = -1;
+        int ass_idx = 0;
+        for (auto ass : *trkparasscol) {
+            if (ass.getRec() == trk) {
+                float weight = ass.getWeight();
+                if (weight > max_weight) {
+                    max_weight = weight;
+                    max_weight_idx = ass_idx;
+                }
+            }
+            ass_idx++;
+        }
+        if (max_weight_idx < 0) continue;
+        pdgid = trkparasscol->at(max_weight_idx).getSim().getPDG();
+        double px = trkparasscol->at(max_weight_idx).getSim().getMomentum()[0];
+        double py = trkparasscol->at(max_weight_idx).getSim().getMomentum()[1];
+        double pz = trkparasscol->at(max_weight_idx).getSim().getMomentum()[2];
+        p_truth = sqrt(px*px + py*py + pz*pz);
+
+        // check_file << "before track par" << std::endl;
+
+        /// Track parameters
+        podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
+        if (hitcol.size() == 0) return StatusCode::SUCCESS;
+
+        // check_file << "before get track state" << std::endl;
+        int ihit_first = -1;
+        int ihit_last = -1;
+        TrackState trk_par;
+        for (auto it = trk.trackStates_begin(); it != trk.trackStates_end(); it++) {
+            if (it->location == 0) {
+                trk_par = *it;
+                break;
+            }
+        }
+        if (trk_par.tanLambda > 0.1) {
+            getFirstAndLastHitsByZ(hitcol, ihit_first, ihit_last);
+        }
+        else {
+            // check_file << "before get first and last by r" << std::endl;
+            getFirstAndLastHitsByRadius(hitcol, ihit_first, ihit_last);
+        }
+        if (ihit_first < 0 || ihit_last < 0 || ihit_first == ihit_last) continue;
+
+        // check_file << "before conversion unit. ihit_first, ihit_last = " << ihit_first << ", " << ihit_last << std::endl;
+        
+        // convert the position from cm to mm
+        Vector3d first_hit_pos = Vector3d(hitcol[ihit_first].getPosition().x*10.0, hitcol[ihit_first].getPosition().y*10.0, hitcol[ihit_first].getPosition().z*10.0);
+        Vector3d last_hit_pos = Vector3d(hitcol[ihit_last].getPosition().x*10.0, hitcol[ihit_last].getPosition().y*10.0, hitcol[ihit_last].getPosition().z*10.0);
+        double len, p, cos;
+        m_pid_svc->getTrackPars(first_hit_pos, last_hit_pos, trk, 0, len, p, cos);
+
+        // check_file << "before dndx calc" << std::endl;
+
+        /// dN/dx reconstruction
+        if (m_method.value() == "Simple") { // Track level implementation
+            int particle_type = m_pid_svc->getParticleType(pdgid);
+            double bg = m_pid_svc->getBetaGamma(p_truth, particle_type);
+            double dndx_mean = m_pid_svc->getDndxMean(bg, cos);
+            double dndx_sigma = m_pid_svc->getDndxSigma(bg, cos, len);
+            double dndx_meas = m_pid_svc->getDndx(dndx_mean, dndx_sigma);
+            // std::cout << "pdgid: " << pdgid << " bg: " << bg << " dndx_mean: " << dndx_mean << " dndx_sigma: " << dndx_sigma << " dndx_meas: " << dndx_meas << std::endl;
+
+            Quantity q;
+            q.value = dndx_meas;
+            q.error = dndx_sigma;
+
+            //check_file << "before hypotheses loop" << std::endl;
+            std::array<Hypothesis, 5> hypotheses;
+            for (int pid = 0; pid < 5; pid++) {
+                bg = m_pid_svc->getBetaGamma(p_truth, pid);
+                dndx_mean = m_pid_svc->getDndxMean(bg, cos);
+                dndx_sigma = m_pid_svc->getDndxSigma(bg, cos, len);
+
+                Hypothesis h;
+                h.chi2 = m_pid_svc->getChi2(dndx_meas, dndx_mean, dndx_sigma);
+                h.expected = dndx_mean;
+                h.sigma = dndx_sigma;
+
+                hypotheses[pid] = h;
+            }
+            //check_file << "after hypotheses loop" << std::endl;
+
+            MutableRecDqdx dndx_track(q, particle_type, 0, hypotheses);
+            dndx_track.setTrack(trk);
+            outCol->push_back(dndx_track);
+        }
+        else if (m_method.value() == "Full") {
+            // Hit level implementation, loop over hits ...
+        }
+        else {
+            return StatusCode::FAILURE;
+        }
+    }
+
+
+    return StatusCode::SUCCESS;
+
+}
+
+StatusCode DCHDndxAlg::finalize() {
+
+    return StatusCode::SUCCESS;
+}
+
+void DCHDndxAlg::getFirstAndLastHitsByZ(const podio::RelationRange<edm4hep::TrackerHit>& hitcol, int& first, int& last) {
+    bool first_dc_hit = true;
+    double zmin, zmax;
+    for (size_t ihit = 0; ihit < hitcol.size(); ihit++) {
+        auto hit = hitcol[ihit];
+        double z = hit.getPosition()[2];
+        if (!isCDCHit(&hit)) {
+            continue;
+        }
+        if (first_dc_hit) {
+            zmin = z;
+            zmax = z;
+            first = ihit;
+            last = ihit;
+            first_dc_hit = false;
+        }
+        else {
+            if (z < zmin) {
+                zmin = z;
+                first = ihit;
+            }
+            if (z > zmax) {
+                zmax = z;
+                last = ihit;
+            }
+        }
+    }
+}
+
+void DCHDndxAlg::getFirstAndLastHitsByRadius(const podio::RelationRange<edm4hep::TrackerHit>& hitcol, int& first, int& last) {
+    bool first_dc_hit = true;
+    double rmin, rmax;
+    for (size_t ihit = 0; ihit < hitcol.size(); ihit++) {
+        auto hit = hitcol[ihit];
+        double x = hit.getPosition()[0];
+        double y = hit.getPosition()[1];
+        double r = sqrt(x*x + y*y);
+        if (!isCDCHit(&hit)) {
+            continue;
+        }
+        if (first_dc_hit) {
+            rmin = r;
+            rmax = r;
+            first = ihit;
+            last = ihit;
+            first_dc_hit = false;
+        }
+        else {
+            if (r < rmin) {
+                rmin = r;
+                first = ihit;
+            }
+            if (r > rmax) {
+                rmax = r;
+                last = ihit;
+            }
+        }
+    }
+}
+
+int DCHDndxAlg::getDetTypeID(unsigned long long cellID) const {
+    UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
+    encoder.setValue(cellID);
+    return encoder[lcio::ILDCellID0::subdet];
+}
+
+bool DCHDndxAlg::isCDCHit(edm4hep::TrackerHit* hit){
+    return m_geom_svc->lcdd()->constant<int>("DetID_DC")==
+        getDetTypeID(hit->getCellID());
+}
\ No newline at end of file

Reconstruction/ParticleID/src/DCHDndxAlg.h

+#ifndef DCHDndxAlg_h
+#define DCHDndxAlg_h 1
+
+#include "k4FWCore/DataHandle.h"
+
+#include "edm4hep/TrackerHit.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "edm4hep/RecDqdx.h"
+#include "edm4hep/RecDqdxCollection.h"
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+
+#include "SimplePIDSvc/ISimplePIDSvc.h"
+#include "DetInterface/IGeomSvc.h"
+
+/**
+ * @class DCHDndxAlg
+ * @brief Algorithm to calculate dN/dx for DCH
+ * @author Guang Zhao (zhaog@ihep.ac.cn)
+*/
+
+class DCHDndxAlg : public GaudiAlgorithm {
+public:
+    DCHDndxAlg(const std::string& name, ISvcLocator* svcLoc);
+
+    virtual StatusCode initialize();
+    virtual StatusCode execute();
+    virtual StatusCode finalize();
+
+protected:
+    Gaudi::Property<std::string> m_method{this, "Method", "Simple"};
+
+    DataHandle<edm4hep::TrackCollection> _trackCol{"SDTRecTrackCollection", Gaudi::DataHandle::Reader, this};
+    //DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _hitassCol{"DCHitAssociationCollection", Gaudi::DataHandle::Reader, this};
+    DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> _trkParAssCol{"SDTRecTrackCollectionParticleAssociation", Gaudi::DataHandle::Reader, this};
+    DataHandle<edm4hep::RecDqdxCollection> _dndxCol{"DndxTracks", Gaudi::DataHandle::Writer, this};
+
+private:
+    SmartIF<ISimplePIDSvc> m_pid_svc;
+    SmartIF<IGeomSvc> m_geom_svc;
+    void getFirstAndLastHitsByZ(const podio::RelationRange<edm4hep::TrackerHit>& hitcol, int& first, int& last);
+    void getFirstAndLastHitsByRadius(const podio::RelationRange<edm4hep::TrackerHit>& hitcol, int& first, int& last);
+    int getDetTypeID(unsigned long long cellID) const;
+    bool isCDCHit(edm4hep::TrackerHit* hit);
+};
+
+#endif

Reconstruction/ParticleID/src/TPCDndxAlg.cpp

+#include "TPCDndxAlg.h"
+#include "DataHelper/HelixClass.h"
+#include "DataHelper/Navigation.h"
+
+#include "edm4hep/Hypothesis.h"
+#include "edm4hep/MCParticle.h"
+#include "edm4hep/Quantity.h"
+#include "edm4hep/Vector3d.h"
+#include "lcio.h"
+
+#include <cmath>
+
+using namespace edm4hep;
+
+DECLARE_COMPONENT(TPCDndxAlg)
+
+TPCDndxAlg::TPCDndxAlg(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
+    // Input
+    declareProperty("ClupatraTracks", _trackCol, "handler of the input track collection");
+    declareProperty("ClupatraTracksParticleAssociation", _trkParAssCol, "handler of the input track particle association collection");
+
+    // Output
+    declareProperty("DndxTracks", _dndxCol, "handler of the collection of dN/dx tracks");
+}
+
+StatusCode TPCDndxAlg::initialize() {
+    info() << "Initilize DndxAlg ..." << endmsg;
+
+    if (m_method == "Simple") {
+        m_pid_svc = service("SimplePIDSvc");
+    }
+    else {
+        m_pid_svc = nullptr;
+    }
+
+    return GaudiAlgorithm::initialize();
+}
+
+StatusCode TPCDndxAlg::execute() {
+    info() << "Dndx reconstruction started" << endmsg;
+
+    const edm4hep::TrackCollection* trkcol = nullptr;
+    const edm4hep::MCRecoTrackParticleAssociationCollection* trkparasscol = nullptr;
+    try {
+        trkcol = _trackCol.get();
+        trkparasscol = _trkParAssCol.get();
+    }
+    catch(...) {
+        //
+    }
+
+    if (trkcol == nullptr || trkparasscol == nullptr) {
+        return StatusCode::SUCCESS;
+    }
+
+    RecDqdxCollection* outCol = _dndxCol.createAndPut();
+
+    // Navigation nav;
+    // nav.Initialize();
+    for (std::size_t i = 0; i < trkcol->size(); i++) {
+        Track trk(trkcol->at(i));
+        
+        /// MC truth information
+        int pdgid = -999;
+        double p_truth = -999.;
+
+        float max_weight = -999.;
+        int max_weight_idx = -1;
+        int ass_idx = 0;
+        for (auto ass : *trkparasscol) {
+            if (ass.getRec() == trk) {
+                float weight = ass.getWeight();
+                if (weight > max_weight) {
+                    max_weight = weight;
+                    max_weight_idx = ass_idx;
+                }
+            }
+            ass_idx++;
+        }
+        if (max_weight_idx < 0) continue;
+        pdgid = trkparasscol->at(max_weight_idx).getSim().getPDG();
+        double px = trkparasscol->at(max_weight_idx).getSim().getMomentum()[0];
+        double py = trkparasscol->at(max_weight_idx).getSim().getMomentum()[1];
+        double pz = trkparasscol->at(max_weight_idx).getSim().getMomentum()[2];
+        p_truth = sqrt(px*px + py*py + pz*pz);
+
+        /// Track parameters
+        podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
+        int nhits = hitcol.size();
+        if (nhits == 0) return StatusCode::SUCCESS;
+
+        int ihit_first = -1;
+        int ihit_last = -1;
+        TrackState trk_par;
+        for (auto it = trk.trackStates_begin(); it != trk.trackStates_end(); it++) {
+            if (it->location == TrackState::AtFirstHit) {
+                trk_par = *it;
+                break;
+            }
+        }
+        if (trk_par.tanLambda > 0.1) {
+            ihit_first = 0;
+            ihit_last = nhits - 1;
+        }
+        else {
+            getFirstAndLastHitsByRadius(hitcol, ihit_first, ihit_last);
+        }
+        if (ihit_first < 0 || ihit_last < 0 || ihit_first == ihit_last) continue;
+    
+        const Vector3d& first_hit_pos = hitcol[ihit_first].getPosition();
+        const Vector3d& last_hit_pos = hitcol[ihit_last].getPosition();
+        double len, p, cos;
+        m_pid_svc->getTrackPars(first_hit_pos, last_hit_pos, trk, TrackState::AtFirstHit, len, p, cos);
+
+        /// dN/dx reconstruction
+        if (m_method.value() == "Simple") { // Track level implementation
+            int particle_type = m_pid_svc->getParticleType(pdgid);
+            double bg = m_pid_svc->getBetaGamma(p_truth, particle_type);
+            double dndx_mean = m_pid_svc->getDndxMean(bg, cos);
+            double dndx_sigma = m_pid_svc->getDndxSigma(bg, cos, len);
+            double dndx_meas = m_pid_svc->getDndx(dndx_mean, dndx_sigma);
+
+            Quantity q;
+            q.value = dndx_meas;
+            q.error = dndx_sigma;
+
+            std::array<Hypothesis, 5> hypotheses;
+            for (int pid = 0; pid < 5; pid++) {
+                bg = m_pid_svc->getBetaGamma(p_truth, pid);
+                dndx_mean = m_pid_svc->getDndxMean(bg, cos);
+                dndx_sigma = m_pid_svc->getDndxSigma(bg, cos, len);
+
+                Hypothesis h;
+                h.chi2 = m_pid_svc->getChi2(dndx_meas, dndx_mean, dndx_sigma);
+                h.expected = dndx_mean;
+                h.sigma = dndx_sigma;
+
+                hypotheses[pid] = h;
+            }
+
+            MutableRecDqdx dndx_track(q, particle_type, 0, hypotheses);
+            dndx_track.setTrack(trk);
+            outCol->push_back(dndx_track);
+        }
+        else if (m_method.value() == "Full") {
+            // Hit level implementation, loop over hits ...
+        }
+        else {
+            return StatusCode::FAILURE;
+        }
+    }
+
+
+    return StatusCode::SUCCESS;
+
+}
+
+StatusCode TPCDndxAlg::finalize() {
+
+    return StatusCode::SUCCESS;
+}
+
+void TPCDndxAlg::getFirstAndLastHitsByRadius(const podio::RelationRange<edm4hep::TrackerHit>& hitcol, int& first, int& last) {
+    bool first_hit = true;
+    double rmin, rmax;
+    for (size_t ihit = 0; ihit < hitcol.size(); ihit++) {
+        auto hit = hitcol[ihit];
+        double x = hit.getPosition()[0];
+        double y = hit.getPosition()[1];
+        double r = sqrt(x*x + y*y);
+        // if (!isCDCHit(&hit)) {
+        //     continue;
+        // }
+        if (first_hit) {
+            rmin = r;
+            rmax = r;
+            first = ihit;
+            last = ihit;
+            first_hit = false;
+        }
+        else {
+            if (r < rmin) {
+                rmin = r;
+                first = ihit;
+            }
+            if (r > rmax) {
+                rmax = r;
+                last = ihit;
+            }
+        }
+    }
+}

Reconstruction/ParticleID/src/TPCDndxAlg.h

+#ifndef TPCDndxAlg_h
+#define TPCDndxAlg_h 1
+
+#include "k4FWCore/DataHandle.h"
+
+#include "edm4hep/TrackerHit.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "edm4hep/RecDqdx.h"
+#include "edm4hep/RecDqdxCollection.h"
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+
+#include "SimplePIDSvc/ISimplePIDSvc.h"
+
+/**
+ * @class TPCDndxAlg
+ * @brief Algorithm to calculate dN/dx for TPC
+ * @author Guang Zhao (zhaog@ihep.ac.cn)
+*/
+
+class TPCDndxAlg : public GaudiAlgorithm {
+public:
+    TPCDndxAlg(const std::string& name, ISvcLocator* svcLoc);
+
+    virtual StatusCode initialize();
+    virtual StatusCode execute();
+    virtual StatusCode finalize();
+
+protected:
+    Gaudi::Property<std::string> m_method{this, "Method", "Simple"};
+
+    DataHandle<edm4hep::TrackCollection> _trackCol{"ClupatraTracks", Gaudi::DataHandle::Reader, this};
+    DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> _trkParAssCol{"ClupatraTracksParticleAssociation", Gaudi::DataHandle::Reader, this};
+    DataHandle<edm4hep::RecDqdxCollection> _dndxCol{"DndxTracks", Gaudi::DataHandle::Writer, this};
+
+private:
+    SmartIF<ISimplePIDSvc> m_pid_svc;
+    void getFirstAndLastHitsByRadius(const podio::RelationRange<edm4hep::TrackerHit>& hitcol, int& first, int& last);
+};
+
+#endif

Service/CMakeLists.txt

 add_subdirectory(EventSeeder)
 add_subdirectory(GearSvc)
 add_subdirectory(TrackSystemSvc)
+add_subdirectory(SimplePIDSvc)

Service/SimplePIDSvc/CMakeLists.txt

+
+gaudi_add_header_only_library(SimplePIDSvc)
+
+gaudi_add_module(SimplePIDSvcPlugins
+                 SOURCES src/SimplePIDSvc.cpp
+                 LINK SimplePIDSvc
+                      GearSvc
+                      Gaudi::GaudiKernel
+                      DataHelperLib
+                      ${DD4hep_COMPONENT_LIBRARIES}
+                      ${ROOT_LIBRARIES}
+                      ${GEAR_LIBRARIES}
+)
+
+#target_include_directories(SimplePIDSvc PUBLIC
+#  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include
+#  $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
+install(TARGETS SimplePIDSvc SimplePIDSvcPlugins
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)
+

Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h

+#ifndef I_SimplePID_SVC_H
+#define I_SimplePID_SVC_H
+
+#include "GaudiKernel/IService.h"
+#include "edm4hep/Vector3d.h"
+#include "edm4hep/TrackerHit.h"
+#include "edm4hep/Track.h"
+#include "edm4hep/Vector3d.h"
+
+/**
+ * @class ISimplePIDSvc
+ * @brief Interface for Simple PID
+ * @author Guang Zhao (zhaog@ihep.ac.cn)
+*/
+
+class ISimplePIDSvc: virtual public IService {
+public:
+    DeclareInterfaceID(ISimplePIDSvc, 0, 1); // major/minor version
+    
+    virtual ~ISimplePIDSvc() = default;
+
+    virtual double getDndx(double mean, double sigma) = 0;
+    virtual double getDndxMean(double bg, double cos) = 0;
+    virtual double getDndxSigma(double bg, double cos, double len) = 0;
+    virtual double getChi2(double dndx_meas, double dndx_exp, double dndx_sigma) = 0;
+    virtual void getTrackPars(const edm4hep::Vector3d& hit1, const edm4hep::Vector3d& hit2, const edm4hep::Track& trk, int location,
+                              double& len, double& p, double& cos) = 0;
+    virtual double getBetaGamma(double p, int pid_type) = 0; // e, mu, pi, K, p: 0, 1, 2, 3, 4
+    virtual int getParticleType(int pdgid) = 0;
+};
+
+
+#endif

Service/SimplePIDSvc/src/SimplePIDSvc.cpp

+#include "SimplePIDSvc.h"
+#include "DataHelper/HelixClass.h"
+#include "GearSvc/IGearSvc.h"
+#include "gear/BField.h"
+#include "TFile.h"
+#include "TH2D.h"
+#include "TRandom.h"
+#include <iostream>
+#include <cmath>
+
+DECLARE_COMPONENT(SimplePIDSvc)
+
+using namespace edm4hep;
+
+SimplePIDSvc::SimplePIDSvc(const std::string& name, ISvcLocator* svc)
+    : base_class(name, svc)
+{
+    m_rootFile = nullptr;
+}
+
+SimplePIDSvc::~SimplePIDSvc()
+{
+    if (m_rootFile != nullptr) delete m_rootFile;
+}
+
+double SimplePIDSvc::getDndx(double mean, double sigma) {
+    return gRandom->Gaus(mean, sigma);
+}
+
+double SimplePIDSvc::getDndxMean(double bg, double cos)
+{
+    return interpolate(m_dndxMean, bg, cos);
+}
+
+double SimplePIDSvc::getDndxSigma(double bg, double cos, double len)
+{
+    return interpolate(m_dndxSigma, bg, cos)/sqrt(len*0.1); // len in mm, need to convert to cm
+}
+
+double SimplePIDSvc::getChi2(double dndx_meas, double dndx_exp, double dndx_sigma) {
+    double chi = (dndx_meas - dndx_exp)/dndx_sigma;
+    return chi*chi;
+}
+
+void SimplePIDSvc::getTrackPars(const edm4hep::Vector3d& first_hit_pos, const edm4hep::Vector3d& last_hit_pos, const edm4hep::Track& trk, int location,
+                                double& length, double& p, double& costheta) 
+{
+    TrackState trk_par;
+    for (auto it = trk.trackStates_begin(); it != trk.trackStates_end(); it++) {
+        if (it->location == location) {
+            trk_par = *it;
+            break;
+        }
+    }
+	auto _gear = service<IGearSvc>("GearSvc");
+	double bfield = _gear->getGearMgr()->getBField().at(gear::Vector3D(0.,0.0,0.)).z();
+    HelixClass hc = HelixClass();
+    hc.Initialize_Canonical(trk_par.phi, trk_par.D0, trk_par.Z0, trk_par.omega, trk_par.tanLambda, bfield);
+    double R = hc.getRadius();
+    double phi1 = atan2(first_hit_pos[1] - hc.getYC(), first_hit_pos[0] - hc.getXC());
+    double phi2 = atan2(last_hit_pos[1] - hc.getYC(), last_hit_pos[0] - hc.getXC());
+    double z1 = first_hit_pos[2];
+    double z2 = last_hit_pos[2];
+
+
+    // std::cout << "xc, yc = " << hc.getXC() << ", " << hc.getYC() << std::endl;
+    // std::cout << "r = " << R << std::endl;
+    // std::cout << "x1, y1 = " << first_hit_pos[0] << ", " << first_hit_pos[1] << std::endl;
+    // std::cout << "x2, y2 = " << last_hit_pos[0] << ", " << last_hit_pos[1] << std::endl;
+    // std::cout << "phi1, phi2 = " << phi1 << ", " << phi2 << std::endl;
+    // std::cout << "tanLambda = " << trk_par.tanLambda << std::endl;
+    // std::cout << "z1, z2 = " << z1 << ", " << z2 << std::endl;
+    // const float* p3 = hc.getMomentum();
+    // std::cout << "momentum = " << sqrt(p3[0]*p3[0] + p3[1]*p3[1] + p3[2]*p3[2]) << std::endl;
+    // std::cout << "charge = " << trk_par.omega << std::endl;
+
+    // if (fabs(z2 - z1) > fabs(2*M_PI*R/cos(atan(trk_par.tanLambda)))) { // length is larger than one cycle
+    if (fabs(trk_par.tanLambda) > 0.1) { // if the track is not vertical, use z
+        length = fabs((z1 - z2)/sin(atan(trk_par.tanLambda)));
+    }
+    else {
+        int outward_direction = int(fabs(z1) < fabs(z2));
+        int positive_charge = int(trk_par.omega > 0.);
+        if (outward_direction * positive_charge > 0) {
+            length = phi1 > phi2 ? (phi1 - phi2)*R/cos(atan(trk_par.tanLambda)) : (phi1 - phi2 + 2*M_PI)*R/cos(atan(trk_par.tanLambda));
+        }
+        else {
+            length = phi1 < phi2 ? (phi2 - phi1)*R/cos(atan(trk_par.tanLambda)) : (phi2 - phi1 + 2*M_PI)*R/cos(atan(trk_par.tanLambda));
+        }
+    }
+    // std::cout << "length (rphi): " << fabs((phi1 - phi2)*R/cos(atan(trk_par.tanLambda))) << std::endl;
+    // std::cout << "length (z): " << fabs((z1 - z2)/sin(atan(trk_par.tanLambda))) << std::endl << std::endl;
+    // std::cout << "length: " << length << std::endl;
+
+    double pt = hc.getPXY();
+    double pz = pt * hc.getTanLambda();
+    p = sqrt(pt*pt + pz*pz);
+    costheta = cos(M_PI/2 - atan(hc.getTanLambda()));
+}
+
+double SimplePIDSvc::getBetaGamma(double p, int pid_type) { // p: GeV/c
+    double mass;
+    switch (pid_type)
+    {
+    case 0:
+        mass = 0.511;
+        break;
+    case 1:
+        mass = 105.658;
+        break;
+    case 2:
+        mass = 139.570;
+        break;
+    case 3:
+        mass = 493.677;
+        break;
+    case 4:
+        mass = 938.272;
+        break;
+    
+    default:
+        mass = 105.658; // default is muon
+        break;
+    }
+    return p/mass*1000.;
+}
+
+int SimplePIDSvc::getParticleType(int pdgid) {
+    int type;
+    switch (abs(pdgid))
+    {
+    case 11:
+        type = 0;
+        break;
+    case 13:
+        type = 1;
+        break;
+    case 211:
+        type = 2;
+        break;
+    case 321:
+        type = 3;
+        break;
+    case 2212:
+        type = 4;
+        break;
+    
+    default:
+        type = -1;
+        break;
+    }
+    return type;
+}
+
+StatusCode SimplePIDSvc::initialize()
+{
+    m_rootFile = new TFile(m_parFile.value().c_str());
+    m_dndxMean = (TH2D*)m_rootFile->Get("dndx_mean");
+    m_dndxSigma = (TH2D*)m_rootFile->Get("dndx_sigma");
+    return StatusCode::SUCCESS;
+}
+
+StatusCode SimplePIDSvc::finalize()
+{
+    if (m_rootFile != nullptr) delete m_rootFile;
+    return StatusCode::SUCCESS;
+}
+
+double SimplePIDSvc::interpolate(TH2D* h, double x, double y) {
+    int nx = h->GetNbinsX();
+    int ny = h->GetNbinsY();
+    double xmax = h->GetXaxis()->GetBinCenter(nx);
+    double xmin = h->GetXaxis()->GetBinCenter(1);
+    double ymax = h->GetYaxis()->GetBinCenter(ny);
+    double ymin = h->GetYaxis()->GetBinCenter(1);
+    double dlogx = (log10(xmax) - log10(xmin)) * 0.01;
+    double dy = (ymax - ymin) * 0.01;
+
+    if (xmin <= x && x <= xmax && ymin <= y && y <= ymax) {
+        return h->Interpolate(x, y);
+    }
+
+    int xloc, yloc;
+    if (x > xmax) xloc = 1;
+    else if (x < xmin) xloc = -1;
+    else xloc = 0;
+    if (y > ymax) yloc = 1;
+    else if (y < ymin) yloc = -1;
+    else yloc = 0;
+
+    double x1, x2, y1, y2, z1, z2;
+    if (xloc == 1) {
+        x1 = pow(10, log10(xmax)-dlogx);
+        x2 = xmax;
+    }
+    else if (xloc == -1) {
+        x1 = xmin;
+        x2 = pow(10, log10(xmin) + dlogx);
+    }
+    if (yloc == 1) {
+        y1 = ymax - dy;
+        y2 = ymax;
+    }
+    else if (yloc == -1) {
+        y1 = ymin;
+        y2 = ymin + dy;
+    }
+
+    if (xloc != 0 && yloc == 0) {
+        z1 = h->Interpolate(x1, y);
+        z2 = h->Interpolate(x2, y);
+        return linear_interpolate(x1, x2, z1, z2, x);
+    }
+    else if (xloc == 0 && yloc != 0) {
+        z1 = h->Interpolate(x, y1);
+        z2 = h->Interpolate(x, y2);
+        return linear_interpolate(y1, y2, z1, z2, y);
+    }
+    else if (xloc != 0 && yloc != 0) {
+        double y_boundary = yloc > 0 ? ymax : ymin;
+        z1 = h->Interpolate(x1, y_boundary);
+        z2 = h->Interpolate(x2, y_boundary);
+        double z_tmp = linear_interpolate(x1, x2, z1, z2, x);
+
+        double x_boundary = xloc > 0 ? xmax : xmin;
+        z1 = h->Interpolate(x_boundary, y1);
+        z2 = h->Interpolate(x_boundary, y2);
+        double z_tmp2 = linear_interpolate(y1, y2, z1, z2, y);
+
+        return z_tmp + z_tmp2 - h->Interpolate(x_boundary, y_boundary);
+    }
+    else {
+        return 0.;
+    }
+}
+
+double SimplePIDSvc::linear_interpolate(double x1, double x2, double y1, double y2, double x) {
+    return (y2-y1)/(x2-x1)*(x-x1) + y1;
+}

Service/SimplePIDSvc/src/SimplePIDSvc.h

+#ifndef SIMPLEPID_SVC_H
+#define SIMPLEPID_SVC_H
+
+#include "SimplePIDSvc/ISimplePIDSvc.h"
+#include <GaudiKernel/Service.h>
+#include "DD4hep/Detector.h"
+#include "gear/GearMgr.h"
+#include "edm4hep/Vector3d.h"
+
+class TFile;
+class TH2D;
+
+/**
+ * @class SimplePIDSvc
+ * @brief Simple PID service
+ * @author Guang Zhao (zhaog@ihep.ac.cn)
+*/
+
+class SimplePIDSvc : public extends<Service, ISimplePIDSvc>
+{
+    public:
+        SimplePIDSvc(const std::string& name, ISvcLocator* svc);
+        virtual ~SimplePIDSvc();
+
+        double getDndx(double mean, double sigma) override;
+        double getDndxMean(double bg, double cos) override;
+        double getDndxSigma(double bg, double cos, double len) override;
+        double getChi2(double dndx_meas, double dndx_exp, double dndx_sigma) override;
+        void getTrackPars(const edm4hep::Vector3d& hit1, const edm4hep::Vector3d& hit2, const edm4hep::Track& trk, int location,
+                          double& len, double& p, double& cos) override;
+        double getBetaGamma(double p, int pid_type) override;
+        int getParticleType(int pdgid) override;
+
+        StatusCode initialize() override;
+        StatusCode finalize() override;
+
+    private:
+        Gaudi::Property<std::string> m_parFile{this, "ParFile", "dNdx_TPC.root"};
+        TFile* m_rootFile;
+        TH2D* m_dndxMean;
+        TH2D* m_dndxSigma;
+        double interpolate(TH2D* h, double bg, double cos);
+        double linear_interpolate(double x1, double x2, double y1, double y2, double x);
+};
+
+#endif