From 9bddeb0ed00a162dbffc78c7c471d75311404102 Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Mon, 27 May 2024 21:12:53 +0800
Subject: [PATCH 01/10] Implementation of dN/dx for gaseous detectors
---
Reconstruction/CMakeLists.txt | 1 +
Reconstruction/ParticleID/CMakeLists.txt | 28 +++
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 182 ++++++++++++++
Reconstruction/ParticleID/src/DCHDndxAlg.h | 40 ++++
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 147 ++++++++++++
Reconstruction/ParticleID/src/TPCDndxAlg.h | 36 +++
Service/CMakeLists.txt | 1 +
Service/SimplePIDSvc/CMakeLists.txt | 24 ++
.../include/SimplePIDSvc/ISimplePIDSvc.h | 26 ++
Service/SimplePIDSvc/src/SimplePIDSvc.cpp | 224 ++++++++++++++++++
Service/SimplePIDSvc/src/SimplePIDSvc.h | 39 +++
11 files changed, 748 insertions(+)
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/Reconstruction/CMakeLists.txt b/Reconstruction/CMakeLists.txt
index 6bb53e06..793d8828 100644
--- a/Reconstruction/CMakeLists.txt
+++ b/Reconstruction/CMakeLists.txt
@@ -5,3 +5,4 @@ add_subdirectory(SiliconTracking)
add_subdirectory(Tracking)
add_subdirectory(RecGenfitAlg)
add_subdirectory(RecAssociationMaker)
+add_subdirectory(ParticleID)
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..7b66d17e
--- /dev/null
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -0,0 +1,182 @@
+#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>
+
+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("DCHitAssociationCollection", _hitassCol, "handler of the input hit 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;
+
+ const edm4hep::TrackCollection* trkcol = nullptr;
+ const edm4hep::MCRecoTrackerAssociationCollection* hitasscol = nullptr;
+ try {
+ trkcol = _trackCol.get();
+ hitasscol = _hitassCol.get();
+ }
+ catch(...) {
+ //
+ }
+
+ if (trkcol == nullptr || hitasscol == 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.;
+ bool found_mc = false;
+
+ podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
+ for (auto hit : hitcol) {
+ for (auto ha : *hitasscol) {
+ if (ha.getRec() == hit) {
+ pdgid = ha.getSim().getMCParticle().getPDG();
+ double px = ha.getSim().getMCParticle().getMomentum()[0];
+ double py = ha.getSim().getMCParticle().getMomentum()[1];
+ double pz = ha.getSim().getMCParticle().getMomentum()[2];
+ p_truth = sqrt(px*px + py*py + pz*pz);
+ // std::cout << "pdgid: " << pdgid << " p_truth: " << p_truth << std::endl;
+ found_mc = true;
+ break;
+ }
+ }
+ if (found_mc) break;
+ }
+ if (!found_mc) continue;
+
+ /// Track parameters
+ if (hitcol.size() == 0) return StatusCode::SUCCESS;
+ bool first_dc_hit = true;
+ int ihit_first, ihit_last;
+ double zmin, zmax;
+ for (size_t ihit = 0; ihit < hitcol.size(); ihit++) {
+ auto hit = hitcol[ihit];
+ if (!isCDCHit(&hit)) {
+ continue;
+ }
+ if (first_dc_hit) {
+ zmin = hit.getPosition()[2];
+ zmax = hit.getPosition()[2];
+ first_dc_hit = false;
+ }
+ else {
+ if (hit.getPosition()[2] < zmin) {
+ zmin = hit.getPosition()[2];
+ ihit_first = ihit;
+ }
+ if (hit.getPosition()[2] > zmax) {
+ zmax = hit.getPosition()[2];
+ ihit_last = ihit;
+ }
+ }
+ }
+
+ TrackerHit first_hit = hitcol[ihit_first];
+ TrackerHit last_hit = hitcol[ihit_last];
+ double len, p, cos;
+ m_pid_svc->getTrackPars(first_hit, last_hit, trk, 0, len, p, cos);
+ len = len * 10.; // convert from cm to mm
+
+ /// 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;
+
+ 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;
+ }
+
+ Quantity q;
+ q.value = dndx_meas;
+ q.error = dndx_sigma;
+
+ MutableRecDqdx dndx_track(q, pdgid, 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;
+}
+
+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..56840ca8
--- /dev/null
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.h
@@ -0,0 +1,40 @@
+#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/RecDqdx.h"
+#include "edm4hep/RecDqdxCollection.h"
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+
+#include "SimplePIDSvc/ISimplePIDSvc.h"
+#include "DetInterface/IGeomSvc.h"
+
+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::RecDqdxCollection> _dndxCol{"DndxTracks", Gaudi::DataHandle::Writer, this};
+
+private:
+ SmartIF<ISimplePIDSvc> m_pid_svc;
+ SmartIF<IGeomSvc> m_geom_svc;
+ 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..801bd52b
--- /dev/null
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -0,0 +1,147 @@
+#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("TPCTrackerHitAss", _hitassCol, "handler of the input hit 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::MCRecoTrackerAssociationCollection* hitasscol = nullptr;
+ try {
+ trkcol = _trackCol.get();
+ hitasscol = _hitassCol.get();
+ }
+ catch(...) {
+ //
+ }
+
+ if (trkcol == nullptr || hitasscol == 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.;
+
+ podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
+ for (auto hit : hitcol) {
+ std::vector<SimTrackerHit> simhits = nav.GetRelatedTrackerHit(hit, hitasscol);
+ if (simhits.size() > 0) {
+ pdgid = simhits[0].getMCParticle().getPDG();
+ double px = simhits[0].getMCParticle().getMomentum()[0];
+ double py = simhits[0].getMCParticle().getMomentum()[1];
+ double pz = simhits[0].getMCParticle().getMomentum()[2];
+ p_truth = sqrt(px*px + py*py + pz*pz);
+ break;
+ }
+
+ //std::cout << hit.getCellID() << std::endl;
+ //double x = hit.getPosition()[0];
+ //double y = hit.getPosition()[1];
+ //double z = hit.getPosition()[2];
+ //double r = sqrt(x*x + y*y);
+ //std::cout << "r, z = " << r << ", " << z << std::endl;
+ //std::cout << "# of sim hits = " << simhits.size() << std::endl;
+ //for (auto s : simhits) {
+ // std::cout << s.getMCParticle().getPDG() << std::endl;
+ //}
+ }
+ if (p_truth < 0) continue;
+
+ /// Track parameters
+ int nhits = hitcol.size();
+ if (nhits == 0) return StatusCode::SUCCESS;
+ TrackerHit first_hit = hitcol[0];
+ TrackerHit last_hit = hitcol[nhits - 1];
+ double len, p, cos;
+ m_pid_svc->getTrackPars(first_hit, last_hit, 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);
+
+ 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;
+ }
+
+ Quantity q;
+ q.value = dndx_meas;
+ q.error = dndx_sigma;
+
+ MutableRecDqdx dndx_track(q, pdgid, 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;
+}
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.h b/Reconstruction/ParticleID/src/TPCDndxAlg.h
new file mode 100644
index 00000000..225ce161
--- /dev/null
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.h
@@ -0,0 +1,36 @@
+#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/RecDqdx.h"
+#include "edm4hep/RecDqdxCollection.h"
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+
+#include "SimplePIDSvc/ISimplePIDSvc.h"
+
+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::MCRecoTrackerAssociationCollection> _hitassCol{"TPCTrackerHitAss", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::RecDqdxCollection> _dndxCol{"DndxTracks", Gaudi::DataHandle::Writer, this};
+
+private:
+ SmartIF<ISimplePIDSvc> m_pid_svc;
+};
+
+#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..dcf1a80b
--- /dev/null
+++ b/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
@@ -0,0 +1,26 @@
+#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"
+
+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::TrackerHit& hit1, const edm4hep::TrackerHit& 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..e8890d80
--- /dev/null
+++ b/Service/SimplePIDSvc/src/SimplePIDSvc.cpp
@@ -0,0 +1,224 @@
+#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::TrackerHit& hit1, const edm4hep::TrackerHit& hit2, const edm4hep::Track& trk, int location,
+ double& length, double& p, double& costheta)
+{
+ Vector3d first_hit_pos = hit1.getPosition();
+ Vector3d last_hit_pos = hit2.getPosition();
+ 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); // WARNING: hard-coded magnetic field
+ 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];
+
+ // if (fabs(z2 - z1) > fabs(2*M_PI*R/cos(atan(trk_par.tanLambda)))) { // length is larger than one cycle
+ // length = fabs((z1 - z2)/sin(atan(trk_par.tanLambda)));
+ // //std::cout << "The trajectory longer than a cycle: " << z1 << ", " << z2 << std::endl;
+ // }
+ // else { // otherwise use rphi, because rphi is more accurate
+ // length = fabs((phi1 - phi2)*R/cos(atan(trk_par.tanLambda)));
+ // }
+ // std::cout << "phi1: " << phi1 << ", phi2: " << phi2 << ", z1: " << z1 << ", z2: " << z2 << std::endl;
+ // std::cout << "length (rphi): " << length << std::endl;
+ // std::cout << "length (z): " << fabs((z1 - z2)/sin(atan(trk_par.tanLambda))) << std::endl;
+
+ length = fabs((z1 - z2)/sin(atan(trk_par.tanLambda)));
+
+ 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..12438bcb
--- /dev/null
+++ b/Service/SimplePIDSvc/src/SimplePIDSvc.h
@@ -0,0 +1,39 @@
+#ifndef SIMPLEPID_SVC_H
+#define SIMPLEPID_SVC_H
+
+#include "SimplePIDSvc/ISimplePIDSvc.h"
+#include <GaudiKernel/Service.h>
+#include "DD4hep/Detector.h"
+#include "gear/GearMgr.h"
+
+class TFile;
+class TH2D;
+
+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::TrackerHit& hit1, const edm4hep::TrackerHit& 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", "/cefs/higgs/zhaog/dndx_files/dNdx_DCH.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
From e3f5e2cead332231df33f7d6a3268d833089e892 Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Tue, 28 May 2024 22:03:03 +0800
Subject: [PATCH 02/10] updata mc track association
---
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 38 +++++++++++----
Reconstruction/ParticleID/src/DCHDndxAlg.h | 4 +-
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 51 +++++++++-----------
Reconstruction/ParticleID/src/TPCDndxAlg.h | 3 +-
4 files changed, 57 insertions(+), 39 deletions(-)
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
index 7b66d17e..ea0b5d54 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -20,7 +20,7 @@ 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("DCHitAssociationCollection", _hitassCol, "handler of the input hit association collection");
+ declareProperty("SDTRecTrackCollectionParticleAssociation", _trkParAssCol, "handler of the input track particle association collection");
// Output
declareProperty("DndxTracks", _dndxCol, "handler of the collection of dN/dx tracks");
@@ -45,16 +45,16 @@ StatusCode DCHDndxAlg::execute() {
info() << "Dndx reconstruction started" << endmsg;
const edm4hep::TrackCollection* trkcol = nullptr;
- const edm4hep::MCRecoTrackerAssociationCollection* hitasscol = nullptr;
+ const edm4hep::MCRecoTrackParticleAssociationCollection* trkparasscol = nullptr;
try {
trkcol = _trackCol.get();
- hitasscol = _hitassCol.get();
+ trkparasscol = _trkParAssCol.get();
}
catch(...) {
//
}
- if (trkcol == nullptr || hitasscol == nullptr) {
+ if (trkcol == nullptr || trkparasscol == nullptr) {
return StatusCode::SUCCESS;
}
@@ -68,9 +68,28 @@ StatusCode DCHDndxAlg::execute() {
/// MC truth information
int pdgid = -999;
double p_truth = -999.;
- bool found_mc = false;
-
- podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
+
+ 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);
+
+ /*bool found_mc = false;
for (auto hit : hitcol) {
for (auto ha : *hitasscol) {
if (ha.getRec() == hit) {
@@ -86,9 +105,10 @@ StatusCode DCHDndxAlg::execute() {
}
if (found_mc) break;
}
- if (!found_mc) continue;
+ if (!found_mc) continue;*/
/// Track parameters
+ podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
if (hitcol.size() == 0) return StatusCode::SUCCESS;
bool first_dc_hit = true;
int ihit_first, ihit_last;
@@ -148,7 +168,7 @@ StatusCode DCHDndxAlg::execute() {
q.value = dndx_meas;
q.error = dndx_sigma;
- MutableRecDqdx dndx_track(q, pdgid, 0, hypotheses);
+ MutableRecDqdx dndx_track(q, particle_type, 0, hypotheses);
dndx_track.setTrack(trk);
outCol->push_back(dndx_track);
}
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.h b/Reconstruction/ParticleID/src/DCHDndxAlg.h
index 56840ca8..bd6f01d8 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.h
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.h
@@ -7,6 +7,7 @@
#include "edm4hep/TrackCollection.h"
#include "edm4hep/TrackerHitCollection.h"
#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
#include "edm4hep/RecDqdx.h"
#include "edm4hep/RecDqdxCollection.h"
@@ -27,7 +28,8 @@ 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::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:
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
index 801bd52b..42db714f 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -17,7 +17,7 @@ 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("TPCTrackerHitAss", _hitassCol, "handler of the input hit association collection");
+ declareProperty("ClupatraTracksParticleAssociation", _trkParAssCol, "handler of the input track particle association collection");
// Output
declareProperty("DndxTracks", _dndxCol, "handler of the collection of dN/dx tracks");
@@ -40,22 +40,21 @@ StatusCode TPCDndxAlg::execute() {
info() << "Dndx reconstruction started" << endmsg;
const edm4hep::TrackCollection* trkcol = nullptr;
- const edm4hep::MCRecoTrackerAssociationCollection* hitasscol = nullptr;
+ const edm4hep::MCRecoTrackParticleAssociationCollection* trkparasscol = nullptr;
try {
trkcol = _trackCol.get();
- hitasscol = _hitassCol.get();
+ trkparasscol = _trkParAssCol.get();
}
catch(...) {
//
}
- if (trkcol == nullptr || hitasscol == nullptr) {
+ 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++) {
@@ -65,32 +64,28 @@ StatusCode TPCDndxAlg::execute() {
int pdgid = -999;
double p_truth = -999.;
- podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
- for (auto hit : hitcol) {
- std::vector<SimTrackerHit> simhits = nav.GetRelatedTrackerHit(hit, hitasscol);
- if (simhits.size() > 0) {
- pdgid = simhits[0].getMCParticle().getPDG();
- double px = simhits[0].getMCParticle().getMomentum()[0];
- double py = simhits[0].getMCParticle().getMomentum()[1];
- double pz = simhits[0].getMCParticle().getMomentum()[2];
- p_truth = sqrt(px*px + py*py + pz*pz);
- break;
+ 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;
+ }
}
-
- //std::cout << hit.getCellID() << std::endl;
- //double x = hit.getPosition()[0];
- //double y = hit.getPosition()[1];
- //double z = hit.getPosition()[2];
- //double r = sqrt(x*x + y*y);
- //std::cout << "r, z = " << r << ", " << z << std::endl;
- //std::cout << "# of sim hits = " << simhits.size() << std::endl;
- //for (auto s : simhits) {
- // std::cout << s.getMCParticle().getPDG() << std::endl;
- //}
+ ass_idx++;
}
- if (p_truth < 0) continue;
+ 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;
TrackerHit first_hit = hitcol[0];
@@ -124,7 +119,7 @@ StatusCode TPCDndxAlg::execute() {
q.value = dndx_meas;
q.error = dndx_sigma;
- MutableRecDqdx dndx_track(q, pdgid, 0, hypotheses);
+ MutableRecDqdx dndx_track(q, particle_type, 0, hypotheses);
dndx_track.setTrack(trk);
outCol->push_back(dndx_track);
}
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.h b/Reconstruction/ParticleID/src/TPCDndxAlg.h
index 225ce161..76c95a84 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.h
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.h
@@ -7,6 +7,7 @@
#include "edm4hep/TrackCollection.h"
#include "edm4hep/TrackerHitCollection.h"
#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
#include "edm4hep/RecDqdx.h"
#include "edm4hep/RecDqdxCollection.h"
@@ -26,7 +27,7 @@ protected:
Gaudi::Property<std::string> m_method{this, "Method", "Simple"};
DataHandle<edm4hep::TrackCollection> _trackCol{"ClupatraTracks", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _hitassCol{"TPCTrackerHitAss", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> _trkParAssCol{"ClupatraTracksParticleAssociation", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::RecDqdxCollection> _dndxCol{"DndxTracks", Gaudi::DataHandle::Writer, this};
private:
--
GitLab
From 9110abbddd6f2fdf8cff05661efc26f45759bb8c Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Wed, 29 May 2024 09:44:27 +0800
Subject: [PATCH 03/10] Minor updates on comments and default parameters
---
Reconstruction/ParticleID/src/DCHDndxAlg.h | 6 ++++++
Reconstruction/ParticleID/src/TPCDndxAlg.h | 6 ++++++
Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h | 6 ++++++
Service/SimplePIDSvc/src/SimplePIDSvc.h | 8 +++++++-
4 files changed, 25 insertions(+), 1 deletion(-)
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.h b/Reconstruction/ParticleID/src/DCHDndxAlg.h
index bd6f01d8..8301662f 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.h
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.h
@@ -16,6 +16,12 @@
#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);
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.h b/Reconstruction/ParticleID/src/TPCDndxAlg.h
index 76c95a84..9c508e0b 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.h
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.h
@@ -15,6 +15,12 @@
#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);
diff --git a/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h b/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
index dcf1a80b..5e048a03 100644
--- a/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
+++ b/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
@@ -6,6 +6,12 @@
#include "edm4hep/TrackerHit.h"
#include "edm4hep/Track.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
diff --git a/Service/SimplePIDSvc/src/SimplePIDSvc.h b/Service/SimplePIDSvc/src/SimplePIDSvc.h
index 12438bcb..aec3da99 100644
--- a/Service/SimplePIDSvc/src/SimplePIDSvc.h
+++ b/Service/SimplePIDSvc/src/SimplePIDSvc.h
@@ -9,6 +9,12 @@
class TFile;
class TH2D;
+/**
+ * @class SimplePIDSvc
+ * @brief Simple PID service
+ * @author Guang Zhao (zhaog@ihep.ac.cn)
+*/
+
class SimplePIDSvc : public extends<Service, ISimplePIDSvc>
{
public:
@@ -28,7 +34,7 @@ class SimplePIDSvc : public extends<Service, ISimplePIDSvc>
StatusCode finalize() override;
private:
- Gaudi::Property<std::string> m_parFile{this, "ParFile", "/cefs/higgs/zhaog/dndx_files/dNdx_DCH.root"};
+ Gaudi::Property<std::string> m_parFile{this, "ParFile", "dNdx_TPC.root"};
TFile* m_rootFile;
TH2D* m_dndxMean;
TH2D* m_dndxSigma;
--
GitLab
From 00de4c5adcf3ca888837e1a08d8d0f38970ef91f Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Wed, 29 May 2024 16:06:31 +0800
Subject: [PATCH 04/10] update the script and bug fix
---
Detector/DetCRD/scripts/TDR_o1_v01/tracking.py | 15 ++++++++++++---
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 8 ++++----
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 8 ++++----
3 files changed, 20 insertions(+), 11 deletions(-)
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/ParticleID/src/DCHDndxAlg.cpp b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
index ea0b5d54..afd8cadf 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -150,6 +150,10 @@ StatusCode DCHDndxAlg::execute() {
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;
+
std::array<Hypothesis, 5> hypotheses;
for (int pid = 0; pid < 5; pid++) {
bg = m_pid_svc->getBetaGamma(p_truth, pid);
@@ -164,10 +168,6 @@ StatusCode DCHDndxAlg::execute() {
hypotheses[pid] = h;
}
- Quantity q;
- q.value = dndx_meas;
- q.error = dndx_sigma;
-
MutableRecDqdx dndx_track(q, particle_type, 0, hypotheses);
dndx_track.setTrack(trk);
outCol->push_back(dndx_track);
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
index 42db714f..454c4754 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -101,6 +101,10 @@ StatusCode TPCDndxAlg::execute() {
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);
@@ -115,10 +119,6 @@ StatusCode TPCDndxAlg::execute() {
hypotheses[pid] = h;
}
- Quantity q;
- q.value = dndx_meas;
- q.error = dndx_sigma;
-
MutableRecDqdx dndx_track(q, particle_type, 0, hypotheses);
dndx_track.setTrack(trk);
outCol->push_back(dndx_track);
--
GitLab
From be0e3fcbf526f60fc0a9b98456182f23af96607d Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Thu, 30 May 2024 00:30:45 +0800
Subject: [PATCH 05/10] Update track length calculation; fix DCH unit
conversion
---
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 9 ++--
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 4 +-
.../include/SimplePIDSvc/ISimplePIDSvc.h | 3 +-
Service/SimplePIDSvc/src/SimplePIDSvc.cpp | 45 ++++++++++++-------
Service/SimplePIDSvc/src/SimplePIDSvc.h | 3 +-
5 files changed, 41 insertions(+), 23 deletions(-)
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
index afd8cadf..f7b26c3b 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -135,11 +135,12 @@ StatusCode DCHDndxAlg::execute() {
}
}
- TrackerHit first_hit = hitcol[ihit_first];
- TrackerHit last_hit = hitcol[ihit_last];
+ // 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, last_hit, trk, 0, len, p, cos);
- len = len * 10.; // convert from cm to mm
+ m_pid_svc->getTrackPars(first_hit_pos, last_hit_pos, trk, 0, len, p, cos);
+ //len = len * 10.; // convert from cm to mm
/// dN/dx reconstruction
if (m_method.value() == "Simple") { // Track level implementation
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
index 454c4754..3d5c38fb 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -88,8 +88,8 @@ StatusCode TPCDndxAlg::execute() {
podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
int nhits = hitcol.size();
if (nhits == 0) return StatusCode::SUCCESS;
- TrackerHit first_hit = hitcol[0];
- TrackerHit last_hit = hitcol[nhits - 1];
+ const Vector3d& first_hit = hitcol[0].getPosition();
+ const Vector3d& last_hit = hitcol[nhits - 1].getPosition();
double len, p, cos;
m_pid_svc->getTrackPars(first_hit, last_hit, trk, TrackState::AtFirstHit, len, p, cos);
diff --git a/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h b/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
index 5e048a03..5a543229 100644
--- a/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
+++ b/Service/SimplePIDSvc/include/SimplePIDSvc/ISimplePIDSvc.h
@@ -5,6 +5,7 @@
#include "edm4hep/Vector3d.h"
#include "edm4hep/TrackerHit.h"
#include "edm4hep/Track.h"
+#include "edm4hep/Vector3d.h"
/**
* @class ISimplePIDSvc
@@ -22,7 +23,7 @@ public:
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::TrackerHit& hit1, const edm4hep::TrackerHit& hit2, const edm4hep::Track& trk, int location,
+ 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;
diff --git a/Service/SimplePIDSvc/src/SimplePIDSvc.cpp b/Service/SimplePIDSvc/src/SimplePIDSvc.cpp
index e8890d80..5ccece5b 100644
--- a/Service/SimplePIDSvc/src/SimplePIDSvc.cpp
+++ b/Service/SimplePIDSvc/src/SimplePIDSvc.cpp
@@ -42,11 +42,9 @@ double SimplePIDSvc::getChi2(double dndx_meas, double dndx_exp, double dndx_sigm
return chi*chi;
}
-void SimplePIDSvc::getTrackPars(const edm4hep::TrackerHit& hit1, const edm4hep::TrackerHit& hit2, const edm4hep::Track& trk, int location,
+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)
{
- Vector3d first_hit_pos = hit1.getPosition();
- Vector3d last_hit_pos = hit2.getPosition();
TrackState trk_par;
for (auto it = trk.trackStates_begin(); it != trk.trackStates_end(); it++) {
if (it->location == location) {
@@ -57,25 +55,42 @@ void SimplePIDSvc::getTrackPars(const edm4hep::TrackerHit& hit1, const edm4hep::
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); // WARNING: hard-coded magnetic field
+ 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
- // length = fabs((z1 - z2)/sin(atan(trk_par.tanLambda)));
- // //std::cout << "The trajectory longer than a cycle: " << z1 << ", " << z2 << std::endl;
- // }
- // else { // otherwise use rphi, because rphi is more accurate
- // length = fabs((phi1 - phi2)*R/cos(atan(trk_par.tanLambda)));
- // }
- // std::cout << "phi1: " << phi1 << ", phi2: " << phi2 << ", z1: " << z1 << ", z2: " << z2 << std::endl;
- // std::cout << "length (rphi): " << length << std::endl;
- // std::cout << "length (z): " << fabs((z1 - z2)/sin(atan(trk_par.tanLambda))) << std::endl;
-
- length = fabs((z1 - z2)/sin(atan(trk_par.tanLambda)));
+ 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();
diff --git a/Service/SimplePIDSvc/src/SimplePIDSvc.h b/Service/SimplePIDSvc/src/SimplePIDSvc.h
index aec3da99..661accb7 100644
--- a/Service/SimplePIDSvc/src/SimplePIDSvc.h
+++ b/Service/SimplePIDSvc/src/SimplePIDSvc.h
@@ -5,6 +5,7 @@
#include <GaudiKernel/Service.h>
#include "DD4hep/Detector.h"
#include "gear/GearMgr.h"
+#include "edm4hep/Vector3d.h"
class TFile;
class TH2D;
@@ -25,7 +26,7 @@ class SimplePIDSvc : public extends<Service, ISimplePIDSvc>
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::TrackerHit& hit1, const edm4hep::TrackerHit& hit2, const edm4hep::Track& trk, int location,
+ 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;
--
GitLab
From adf9646e02470cca186b30d57e3d72d08b5035c8 Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Thu, 30 May 2024 21:12:46 +0800
Subject: [PATCH 06/10] update track length calculation and bug fixes
---
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 111 ++++++++++++-------
Reconstruction/ParticleID/src/DCHDndxAlg.h | 2 +
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 47 +++++++-
Reconstruction/ParticleID/src/TPCDndxAlg.h | 1 +
4 files changed, 118 insertions(+), 43 deletions(-)
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
index f7b26c3b..29110f87 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -12,6 +12,7 @@
#include "lcio.h"
#include <cmath>
+#include <fstream>
using namespace edm4hep;
@@ -43,6 +44,7 @@ StatusCode DCHDndxAlg::initialize() {
StatusCode DCHDndxAlg::execute() {
info() << "Dndx reconstruction started" << endmsg;
+ static std::ofstream check_file("log.txt");
const edm4hep::TrackCollection* trkcol = nullptr;
const edm4hep::MCRecoTrackParticleAssociationCollection* trkparasscol = nullptr;
@@ -63,6 +65,7 @@ StatusCode DCHDndxAlg::execute() {
//Navigation nav;
//nav.Initialize();
for (std::size_t i = 0; i < trkcol->size(); i++) {
+ //check_file << "Track " << i << std::endl;
Track trk(trkcol->at(i));
/// MC truth information
@@ -89,50 +92,16 @@ StatusCode DCHDndxAlg::execute() {
double pz = trkparasscol->at(max_weight_idx).getSim().getMomentum()[2];
p_truth = sqrt(px*px + py*py + pz*pz);
- /*bool found_mc = false;
- for (auto hit : hitcol) {
- for (auto ha : *hitasscol) {
- if (ha.getRec() == hit) {
- pdgid = ha.getSim().getMCParticle().getPDG();
- double px = ha.getSim().getMCParticle().getMomentum()[0];
- double py = ha.getSim().getMCParticle().getMomentum()[1];
- double pz = ha.getSim().getMCParticle().getMomentum()[2];
- p_truth = sqrt(px*px + py*py + pz*pz);
- // std::cout << "pdgid: " << pdgid << " p_truth: " << p_truth << std::endl;
- found_mc = true;
- break;
- }
- }
- if (found_mc) break;
- }
- if (!found_mc) continue;*/
-
/// Track parameters
podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
if (hitcol.size() == 0) return StatusCode::SUCCESS;
- bool first_dc_hit = true;
+
int ihit_first, ihit_last;
- double zmin, zmax;
- for (size_t ihit = 0; ihit < hitcol.size(); ihit++) {
- auto hit = hitcol[ihit];
- if (!isCDCHit(&hit)) {
- continue;
- }
- if (first_dc_hit) {
- zmin = hit.getPosition()[2];
- zmax = hit.getPosition()[2];
- first_dc_hit = false;
- }
- else {
- if (hit.getPosition()[2] < zmin) {
- zmin = hit.getPosition()[2];
- ihit_first = ihit;
- }
- if (hit.getPosition()[2] > zmax) {
- zmax = hit.getPosition()[2];
- ihit_last = ihit;
- }
- }
+ if (trk.getTrackStates(0).tanLambda > 0.1) {
+ getFirstAndLastHitsByZ(hitcol, ihit_first, ihit_last);
+ }
+ else {
+ getFirstAndLastHitsByRadius(hitcol, ihit_first, ihit_last);
}
// convert the position from cm to mm
@@ -155,6 +124,7 @@ StatusCode DCHDndxAlg::execute() {
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);
@@ -168,6 +138,7 @@ StatusCode DCHDndxAlg::execute() {
hypotheses[pid] = h;
}
+ //check_file << "after hypotheses loop" << std::endl;
MutableRecDqdx dndx_track(q, particle_type, 0, hypotheses);
dndx_track.setTrack(trk);
@@ -191,6 +162,66 @@ 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);
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.h b/Reconstruction/ParticleID/src/DCHDndxAlg.h
index 8301662f..3601d755 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.h
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.h
@@ -41,6 +41,8 @@ protected:
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);
};
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
index 3d5c38fb..64e7a088 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -88,10 +88,20 @@ StatusCode TPCDndxAlg::execute() {
podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
int nhits = hitcol.size();
if (nhits == 0) return StatusCode::SUCCESS;
- const Vector3d& first_hit = hitcol[0].getPosition();
- const Vector3d& last_hit = hitcol[nhits - 1].getPosition();
+
+ int ihit_first, ihit_last;
+ if (trk.getTrackStates(0).tanLambda > 0.1) {
+ ihit_first = 0;
+ ihit_last = nhits - 1;
+ }
+ else {
+ getFirstAndLastHitsByRadius(hitcol, ihit_first, ihit_last);
+ }
+
+ 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, last_hit, trk, TrackState::AtFirstHit, 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
@@ -140,3 +150,34 @@ 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
index 9c508e0b..c8cd699e 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.h
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.h
@@ -38,6 +38,7 @@ protected:
private:
SmartIF<ISimplePIDSvc> m_pid_svc;
+ void getFirstAndLastHitsByRadius(const podio::RelationRange<edm4hep::TrackerHit>& hitcol, int& first, int& last);
};
#endif
--
GitLab
From 4259ef8fd0cee0441ff3538a10d843687fccb431 Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Fri, 31 May 2024 15:46:30 +0800
Subject: [PATCH 07/10] bug fixes
---
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 17 +++++++++++++++--
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 4 +++-
2 files changed, 18 insertions(+), 3 deletions(-)
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
index 29110f87..b5703b12 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -45,6 +45,9 @@ StatusCode DCHDndxAlg::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;
+ //if (ievent <= 262) return StatusCode::SUCCESS;
const edm4hep::TrackCollection* trkcol = nullptr;
const edm4hep::MCRecoTrackParticleAssociationCollection* trkparasscol = nullptr;
@@ -62,10 +65,11 @@ StatusCode DCHDndxAlg::execute() {
RecDqdxCollection* outCol = _dndxCol.createAndPut();
+ check_file << "before track loop" << std::endl;
//Navigation nav;
//nav.Initialize();
for (std::size_t i = 0; i < trkcol->size(); i++) {
- //check_file << "Track " << i << std::endl;
+ check_file << " Track " << i << std::endl;
Track trk(trkcol->at(i));
/// MC truth information
@@ -92,18 +96,25 @@ StatusCode DCHDndxAlg::execute() {
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();
+ check_file << "hitcol size = " << hitcol.size() << std::endl;
if (hitcol.size() == 0) return StatusCode::SUCCESS;
- int ihit_first, ihit_last;
+ check_file << "before get track state" << std::endl;
+ int ihit_first = -1;
+ int ihit_last = -1;
if (trk.getTrackStates(0).tanLambda > 0.1) {
getFirstAndLastHitsByZ(hitcol, ihit_first, ihit_last);
}
else {
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);
@@ -111,6 +122,8 @@ StatusCode DCHDndxAlg::execute() {
m_pid_svc->getTrackPars(first_hit_pos, last_hit_pos, trk, 0, len, p, cos);
//len = len * 10.; // convert from cm to mm
+ 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);
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
index 64e7a088..ca43923d 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -89,7 +89,8 @@ StatusCode TPCDndxAlg::execute() {
int nhits = hitcol.size();
if (nhits == 0) return StatusCode::SUCCESS;
- int ihit_first, ihit_last;
+ int ihit_first = -1;
+ int ihit_last = -1;
if (trk.getTrackStates(0).tanLambda > 0.1) {
ihit_first = 0;
ihit_last = nhits - 1;
@@ -97,6 +98,7 @@ StatusCode TPCDndxAlg::execute() {
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();
--
GitLab
From 3c187cd61c9055a1dd9f81b71d5e8bac732042d8 Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Fri, 31 May 2024 16:43:21 +0800
Subject: [PATCH 08/10] bug fix
---
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 11 ++++++++++-
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 9 ++++++++-
2 files changed, 18 insertions(+), 2 deletions(-)
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
index b5703b12..0b1660cd 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -106,11 +106,20 @@ StatusCode DCHDndxAlg::execute() {
check_file << "before get track state" << std::endl;
int ihit_first = -1;
int ihit_last = -1;
- if (trk.getTrackStates(0).tanLambda > 0.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);
+ check_file << "after get first and last by r" << std::endl;
}
if (ihit_first < 0 || ihit_last < 0 || ihit_first == ihit_last) continue;
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
index ca43923d..2bd63ab5 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -91,7 +91,14 @@ StatusCode TPCDndxAlg::execute() {
int ihit_first = -1;
int ihit_last = -1;
- if (trk.getTrackStates(0).tanLambda > 0.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;
}
--
GitLab
From 77b33c493554239f483727ef9f2c285e801c3e4d Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Fri, 31 May 2024 23:41:20 +0800
Subject: [PATCH 09/10] remove comments
---
Reconstruction/ParticleID/src/DCHDndxAlg.cpp | 28 +++++++++-----------
1 file changed, 12 insertions(+), 16 deletions(-)
diff --git a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
index 0b1660cd..ed9cf4af 100644
--- a/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/DCHDndxAlg.cpp
@@ -44,10 +44,9 @@ StatusCode DCHDndxAlg::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;
- //if (ievent <= 262) return StatusCode::SUCCESS;
+ // 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;
@@ -65,11 +64,10 @@ StatusCode DCHDndxAlg::execute() {
RecDqdxCollection* outCol = _dndxCol.createAndPut();
- check_file << "before track loop" << std::endl;
- //Navigation nav;
- //nav.Initialize();
+ // check_file << "before track loop" << std::endl;
+
for (std::size_t i = 0; i < trkcol->size(); i++) {
- check_file << " Track " << i << std::endl;
+ // check_file << " Track " << i << std::endl;
Track trk(trkcol->at(i));
/// MC truth information
@@ -96,14 +94,13 @@ StatusCode DCHDndxAlg::execute() {
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;
+ // check_file << "before track par" << std::endl;
/// Track parameters
podio::RelationRange<edm4hep::TrackerHit> hitcol = trk.getTrackerHits();
- check_file << "hitcol size = " << hitcol.size() << std::endl;
if (hitcol.size() == 0) return StatusCode::SUCCESS;
- check_file << "before get track state" << std::endl;
+ // check_file << "before get track state" << std::endl;
int ihit_first = -1;
int ihit_last = -1;
TrackState trk_par;
@@ -117,21 +114,20 @@ StatusCode DCHDndxAlg::execute() {
getFirstAndLastHitsByZ(hitcol, ihit_first, ihit_last);
}
else {
- check_file << "before get first and last by r" << std::endl;
+ // check_file << "before get first and last by r" << std::endl;
getFirstAndLastHitsByRadius(hitcol, ihit_first, ihit_last);
- check_file << "after get first and last by r" << std::endl;
}
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;
+ // 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);
- //len = len * 10.; // convert from cm to mm
- check_file << "before dndx calc" << std::endl;
+ // check_file << "before dndx calc" << std::endl;
/// dN/dx reconstruction
if (m_method.value() == "Simple") { // Track level implementation
--
GitLab
From 61c07b089ae5a5e88bbfce439ecee64e5f7ad263 Mon Sep 17 00:00:00 2001
From: Guang Zhao <zhaog@ihep.ac.cn>
Date: Sun, 2 Jun 2024 10:57:32 +0800
Subject: [PATCH 10/10] remove uneccesary functions
---
Reconstruction/ParticleID/src/TPCDndxAlg.cpp | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
index 2bd63ab5..0af4dbc2 100644
--- a/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
+++ b/Reconstruction/ParticleID/src/TPCDndxAlg.cpp
@@ -55,8 +55,8 @@ StatusCode TPCDndxAlg::execute() {
RecDqdxCollection* outCol = _dndxCol.createAndPut();
- Navigation nav;
- nav.Initialize();
+ // Navigation nav;
+ // nav.Initialize();
for (std::size_t i = 0; i < trkcol->size(); i++) {
Track trk(trkcol->at(i));
--
GitLab