From a3b474b19f71e71dcfd95b94f3c4c7a9bd39ad69 Mon Sep 17 00:00:00 2001
From: "zhaog@ihep.ac.cn" <zhaog@ihep.ac.cn>
Date: Thu, 13 Jun 2024 06:43:57 +0000
Subject: [PATCH] dN/dx for TPC and DCH
---
.../DetCRD/scripts/TDR_o1_v01/tracking.py | 15 +-
Reconstruction/CMakeLists.txt | 1 +
Reconstruction/ParticleID/CMakeLists.txt | 28 ++
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 252 ++++++++++++++++++
Reconstruction/ParticleID/src/DCHDndxAlg.h | 50 ++++
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 192 +++++++++++++
Reconstruction/ParticleID/src/TPCDndxAlg.h | 44 +++
Service/CMakeLists.txt | 1 +
Service/SimplePIDSvc/CMakeLists.txt | 24 ++
.../include/SimplePIDSvc/ISimplePIDSvc.h | 33 +++
Service/SimplePIDSvc/src/SimplePIDSvc.cpp | 239 +++++++++++++++++
Service/SimplePIDSvc/src/SimplePIDSvc.h | 46 ++++
12 files changed, 922 insertions(+), 3 deletions(-)
create mode 100644 Reconstruction/ParticleID/CMakeLists.txt
create mode 100644 Reconstruction/ParticleID/src/DCHDndxAlg.cpp
create mode 100644 Reconstruction/ParticleID/src/DCHDndxAlg.h
create mode 100644 Reconstruction/ParticleID/src/TPCDndxAlg.cpp
create mode 100644 Reconstruction/ParticleID/src/TPCDndxAlg.h
create mode 100644 Service/SimplePIDSvc/CMakeLists.txt
create mode 100644 Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
create mode 100644 Service/SimplePIDSvc/src/SimplePIDSvc.cpp
create mode 100644 Service/SimplePIDSvc/src/SimplePIDSvc.h
diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
index 6c7a4e89..5ce28166 100644
--- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
+++ b/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
)
diff --git a/Reconstruction/CMakeLists.txt b/Reconstruction/CMakeLists.txt
index 414e1160..b13381ae 100644
--- a/Reconstruction/CMakeLists.txt
+++ b/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)
diff --git a/Reconstruction/ParticleID/CMakeLists.txt b/Reconstruction/ParticleID/CMakeLists.txt
new file mode 100644
index 00000000..6530b72f
--- /dev/null
+++ b/Reconstruction/ParticleID/CMakeLists.txt
@@ -0,0 +1,28 @@
+# 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)
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
new file mode 100644
index 00000000..ed9cf4af
--- /dev/null
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -0,0 +1,252 @@
+#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
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.h b/Reconstruction/ParticleID/src/DCHDndxAlg.h
new file mode 100644
index 00000000..3601d755
--- /dev/null
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.h
@@ -0,0 +1,50 @@
+#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
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
new file mode 100644
index 00000000..0af4dbc2
--- /dev/null
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -0,0 +1,192 @@
+#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;
+ }
+ }
+ }
+}
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.h b/Reconstruction/ParticleID/src/TPCDndxAlg.h
new file mode 100644
index 00000000..c8cd699e
--- /dev/null
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.h
@@ -0,0 +1,44 @@
+#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
diff --git a/Service/CMakeLists.txt b/Service/CMakeLists.txt
index 24ff600b..14c79538 100644
--- a/Service/CMakeLists.txt
+++ b/Service/CMakeLists.txt
@@ -1,3 +1,4 @@
add_subdirectory(EventSeeder)
add_subdirectory(GearSvc)
add_subdirectory(TrackSystemSvc)
+add_subdirectory(SimplePIDSvc)
diff --git a/Service/SimplePIDSvc/CMakeLists.txt b/Service/SimplePIDSvc/CMakeLists.txt
new file mode 100644
index 00000000..6e8d14b9
--- /dev/null
+++ b/Service/SimplePIDSvc/CMakeLists.txt
@@ -0,0 +1,24 @@
+
+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)
+
diff --git a/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h b/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
new file mode 100644
index 00000000..5a543229
--- /dev/null
+++ b/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
@@ -0,0 +1,33 @@
+#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
diff --git a/Service/SimplePIDSvc/src/SimplePIDSvc.cpp b/Service/SimplePIDSvc/src/SimplePIDSvc.cpp
new file mode 100644
index 00000000..5ccece5b
--- /dev/null
+++ b/Service/SimplePIDSvc/src/SimplePIDSvc.cpp
@@ -0,0 +1,239 @@
+#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;
+}
diff --git a/Service/SimplePIDSvc/src/SimplePIDSvc.h b/Service/SimplePIDSvc/src/SimplePIDSvc.h
new file mode 100644
index 00000000..661accb7
--- /dev/null
+++ b/Service/SimplePIDSvc/src/SimplePIDSvc.h
@@ -0,0 +1,46 @@
+#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
--
GitLab