From d5ec7e5c7b6738e43b03831b799cd31313d4f1eb Mon Sep 17 00:00:00 2001
From: "guofangyi@ihep.ac.cn" <guofangyi@ihep.ac.cn>
Date: Wed, 16 Oct 2024 09:46:55 +0000
Subject: [PATCH] Fix track reference points in PFO and remove redundant
collections
---
Analysis/ReadDigi/src/ReadDigiAlg.cpp | 446 ++++++++++++------
Analysis/ReadDigi/src/ReadDigiAlg.h | 35 +-
Reconstruction/RecPFACyber/script/digi.py | 25 +-
Reconstruction/RecPFACyber/script/rec.py | 14 +-
Reconstruction/RecPFACyber/script/sim.py | 24 +-
Reconstruction/RecPFACyber/script/tracking.py | 27 +-
.../RecPFACyber/src/Objects/Track.cc | 12 +-
7 files changed, 395 insertions(+), 188 deletions(-)
diff --git a/Analysis/ReadDigi/src/ReadDigiAlg.cpp b/Analysis/ReadDigi/src/ReadDigiAlg.cpp
index 186a5dd7..c2cdb4ed 100644
--- a/Analysis/ReadDigi/src/ReadDigiAlg.cpp
+++ b/Analysis/ReadDigi/src/ReadDigiAlg.cpp
@@ -12,21 +12,24 @@ ReadDigiAlg::ReadDigiAlg(const std::string& name, ISvcLocator* svcLoc)
declareProperty("MCParticle", m_MCParticleCol, "MCParticle collection (input)");
//Tracker hit
- declareProperty("SITRawHits", m_SITRawColHdl, "SIT Raw Hit Collection of SpacePoints");
- declareProperty("SETRawHits", m_SETRawColHdl, "SET Raw Hit Collection of SpacePoints");
- declareProperty("FTDRawHits", m_FTDRawColHdl, "FTD Raw Hit Collection of SpacePoints");
- declareProperty("VTXTrackerHits", m_VTXTrackerHitColHdl, "VTX Hit Collection");
- declareProperty("SITTrackerHits", m_SITTrackerHitColHdl, "SIT Hit Collection");
- declareProperty("SETTrackerHits", m_SETTrackerHitColHdl, "SET Hit Collection");
- declareProperty("TPCTrackerHits", m_TPCTrackerHitColHdl, "TPC Hit Collection");
- declareProperty("FTDSpacePoints", m_FTDSpacePointColHdl, "FTD FTDSpacePoint Collection");
- declareProperty("FTDPixelTrackerHits", m_FTDPixelTrackerHitColHdl, "handler of FTD Pixel Hit Collection");
+ //declareProperty("SITRawHits", m_SITRawColHdl, "SIT Raw Hit Collection of SpacePoints");
+ //declareProperty("SETRawHits", m_SETRawColHdl, "SET Raw Hit Collection of SpacePoints");
+ //declareProperty("FTDRawHits", m_FTDRawColHdl, "FTD Raw Hit Collection of SpacePoints");
+ //declareProperty("VTXTrackerHits", m_VTXTrackerHitColHdl, "VTX Hit Collection");
+ //declareProperty("SITTrackerHits", m_SITTrackerHitColHdl, "SIT Hit Collection");
+ //declareProperty("SETTrackerHits", m_SETTrackerHitColHdl, "SET Hit Collection");
+ //declareProperty("TPCTrackerHits", m_TPCTrackerHitColHdl, "TPC Hit Collection");
+ //declareProperty("FTDSpacePoints", m_FTDSpacePointColHdl, "FTD FTDSpacePoint Collection");
+ //declareProperty("FTDPixelTrackerHits", m_FTDPixelTrackerHitColHdl, "handler of FTD Pixel Hit Collection");
//Track
declareProperty("FullTracks", m_fullTrk, "Full Track Collection");
declareProperty("TPCTracks", m_TPCTrk, "TPC Track Collection");
declareProperty("SiTracks", m_SiTrk, "Si Track Collection");
+ declareProperty("TPCTracksAssociation", m_TPCTrkAssoHdl, "TPC Track - MCParticle Collection");
+ declareProperty("FullTracksAssociation", m_fullTrkAssoHdl, "Full Track - MCParticle Collection");
+
//declareProperty("EcalBarrelCollection", m_ECalBarrelHitCol, "ECal Barrel");
//declareProperty("EcalEndcapsCollection", m_ECalEndcapHitCol, "ECal Endcap");
//declareProperty("HcalBarrelCollection", m_HCalBarrelHitCol, "HCal Barrel");
@@ -40,7 +43,9 @@ StatusCode ReadDigiAlg::initialize()
std::string s_outfile = _filename;
m_wfile = new TFile(s_outfile.c_str(), "recreate");
m_mctree = new TTree("MCParticle", "MCParticle");
- m_trktree = new TTree("RecoTrk", "RecoTrk");
+ m_sitrktree = new TTree("RecoSiTrk", "RecoSiTrk");
+ m_tpctrktree = new TTree("RecoTPCTrk", "RecoTPCTrk");
+ m_fulltrktree = new TTree("RecoFullTrk", "RecoFullTrk");
m_mctree->Branch("N_MCP", &N_MCP);
m_mctree->Branch("MCP_px", &MCP_px);
@@ -53,28 +58,103 @@ StatusCode ReadDigiAlg::initialize()
m_mctree->Branch("MCP_pdgid", &MCP_pdgid);
m_mctree->Branch("MCP_gStatus", &MCP_gStatus);
- m_trktree->Branch("N_SiTrk", &N_SiTrk);
- m_trktree->Branch("N_TPCTrk", &N_TPCTrk);
- m_trktree->Branch("N_fullTrk", &N_fullTrk);
- m_trktree->Branch("N_VTXhit", &N_VTXhit);
- m_trktree->Branch("N_SIThit", &N_SIThit);
- m_trktree->Branch("N_TPChit", &N_TPChit);
- m_trktree->Branch("N_SEThit", &N_SEThit);
- m_trktree->Branch("N_FTDhit", &N_FTDhit);
- m_trktree->Branch("N_SITrawhit", &N_SITrawhit);
- m_trktree->Branch("N_SETrawhit", &N_SETrawhit);
- m_trktree->Branch("trk_px", &m_trk_px);
- m_trktree->Branch("trk_py", &m_trk_py);
- m_trktree->Branch("trk_pz", &m_trk_pz);
- m_trktree->Branch("trk_p", &m_trk_p);
- m_trktree->Branch("trk_startpoint_x", &m_trk_startpoint_x);
- m_trktree->Branch("trk_startpoint_y", &m_trk_startpoint_y);
- m_trktree->Branch("trk_startpoint_z", &m_trk_startpoint_z);
- m_trktree->Branch("trk_endpoint_x", &m_trk_endpoint_x);
- m_trktree->Branch("trk_endpoint_y", &m_trk_endpoint_y);
- m_trktree->Branch("trk_endpoint_z", &m_trk_endpoint_z);
- m_trktree->Branch("trk_type", &m_trk_type);
- m_trktree->Branch("trk_Nhit", &m_trk_Nhit);
+ //m_trktree->Branch("N_VTXhit", &N_VTXhit);
+ //m_trktree->Branch("N_SIThit", &N_SIThit);
+ //m_trktree->Branch("N_TPChit", &N_TPChit);
+ //m_trktree->Branch("N_SEThit", &N_SEThit);
+ //m_trktree->Branch("N_FTDhit", &N_FTDhit);
+ //m_trktree->Branch("N_SITrawhit", &N_SITrawhit);
+ //m_trktree->Branch("N_SETrawhit", &N_SETrawhit);
+ m_sitrktree->Branch("N_SiTrk", &N_SiTrk);
+ m_sitrktree->Branch("N_TPCTrk", &N_TPCTrk);
+ m_sitrktree->Branch("N_fullTrk", &N_fullTrk);
+ m_sitrktree->Branch("trk_px", &m_trk_px);
+ m_sitrktree->Branch("trk_py", &m_trk_py);
+ m_sitrktree->Branch("trk_pz", &m_trk_pz);
+ m_sitrktree->Branch("trk_p", &m_trk_p);
+ m_sitrktree->Branch("trk_Nhit", &m_trk_Nhit);
+ m_sitrktree->Branch("trkstate_d0", &m_trkstate_d0 );
+ m_sitrktree->Branch("trkstate_z0", &m_trkstate_z0 );
+ m_sitrktree->Branch("trkstate_phi", &m_trkstate_phi );
+ m_sitrktree->Branch("trkstate_tanL", &m_trkstate_tanL );
+ m_sitrktree->Branch("trkstate_omega", &m_trkstate_omega );
+ m_sitrktree->Branch("trkstate_kappa", &m_trkstate_kappa );
+ m_sitrktree->Branch("trkstate_refx", &m_trkstate_refx );
+ m_sitrktree->Branch("trkstate_refy", &m_trkstate_refy );
+ m_sitrktree->Branch("trkstate_refz", &m_trkstate_refz );
+ m_sitrktree->Branch("trkstate_tag", &m_trkstate_tag );
+ m_sitrktree->Branch("trkstate_location", &m_trkstate_location );
+ m_sitrktree->Branch("truthMC_tag", &m_truthMC_tag);
+ m_sitrktree->Branch("truthMC_pid", &m_truthMC_pid);
+ m_sitrktree->Branch("truthMC_px", &m_truthMC_px);
+ m_sitrktree->Branch("truthMC_py", &m_truthMC_py);
+ m_sitrktree->Branch("truthMC_pz", &m_truthMC_pz);
+ m_sitrktree->Branch("truthMC_E", &m_truthMC_E);
+ m_sitrktree->Branch("truthMC_EPx", &m_truthMC_EPx);
+ m_sitrktree->Branch("truthMC_EPy", &m_truthMC_EPy);
+ m_sitrktree->Branch("truthMC_EPz", &m_truthMC_EPz);
+ m_sitrktree->Branch("truthMC_weight", &m_truthMC_weight);
+
+ m_tpctrktree->Branch("N_SiTrk", &N_SiTrk);
+ m_tpctrktree->Branch("N_TPCTrk", &N_TPCTrk);
+ m_tpctrktree->Branch("N_fullTrk", &N_fullTrk);
+ m_tpctrktree->Branch("trk_px", &m_trk_px);
+ m_tpctrktree->Branch("trk_py", &m_trk_py);
+ m_tpctrktree->Branch("trk_pz", &m_trk_pz);
+ m_tpctrktree->Branch("trk_p", &m_trk_p);
+ m_tpctrktree->Branch("trk_Nhit", &m_trk_Nhit);
+ m_tpctrktree->Branch("trkstate_d0", &m_trkstate_d0 );
+ m_tpctrktree->Branch("trkstate_z0", &m_trkstate_z0 );
+ m_tpctrktree->Branch("trkstate_phi", &m_trkstate_phi );
+ m_tpctrktree->Branch("trkstate_tanL", &m_trkstate_tanL );
+ m_tpctrktree->Branch("trkstate_omega", &m_trkstate_omega );
+ m_tpctrktree->Branch("trkstate_kappa", &m_trkstate_kappa );
+ m_tpctrktree->Branch("trkstate_refx", &m_trkstate_refx );
+ m_tpctrktree->Branch("trkstate_refy", &m_trkstate_refy );
+ m_tpctrktree->Branch("trkstate_refz", &m_trkstate_refz );
+ m_tpctrktree->Branch("trkstate_tag", &m_trkstate_tag );
+ m_tpctrktree->Branch("trkstate_location", &m_trkstate_location );
+ m_tpctrktree->Branch("truthMC_tag", &m_truthMC_tag);
+ m_tpctrktree->Branch("truthMC_pid", &m_truthMC_pid);
+ m_tpctrktree->Branch("truthMC_px", &m_truthMC_px);
+ m_tpctrktree->Branch("truthMC_py", &m_truthMC_py);
+ m_tpctrktree->Branch("truthMC_pz", &m_truthMC_pz);
+ m_tpctrktree->Branch("truthMC_E", &m_truthMC_E);
+ m_tpctrktree->Branch("truthMC_EPx", &m_truthMC_EPx);
+ m_tpctrktree->Branch("truthMC_EPy", &m_truthMC_EPy);
+ m_tpctrktree->Branch("truthMC_EPz", &m_truthMC_EPz);
+ m_tpctrktree->Branch("truthMC_weight", &m_truthMC_weight);
+
+ m_fulltrktree->Branch("N_SiTrk", &N_SiTrk);
+ m_fulltrktree->Branch("N_TPCTrk", &N_TPCTrk);
+ m_fulltrktree->Branch("N_fullTrk", &N_fullTrk);
+ m_fulltrktree->Branch("trk_px", &m_trk_px);
+ m_fulltrktree->Branch("trk_py", &m_trk_py);
+ m_fulltrktree->Branch("trk_pz", &m_trk_pz);
+ m_fulltrktree->Branch("trk_p", &m_trk_p);
+ m_fulltrktree->Branch("trk_Nhit", &m_trk_Nhit);
+ m_fulltrktree->Branch("trk_type", &m_trk_type);
+ m_fulltrktree->Branch("trkstate_d0", &m_trkstate_d0 );
+ m_fulltrktree->Branch("trkstate_z0", &m_trkstate_z0 );
+ m_fulltrktree->Branch("trkstate_phi", &m_trkstate_phi );
+ m_fulltrktree->Branch("trkstate_tanL", &m_trkstate_tanL );
+ m_fulltrktree->Branch("trkstate_omega", &m_trkstate_omega );
+ m_fulltrktree->Branch("trkstate_kappa", &m_trkstate_kappa );
+ m_fulltrktree->Branch("trkstate_refx", &m_trkstate_refx );
+ m_fulltrktree->Branch("trkstate_refy", &m_trkstate_refy );
+ m_fulltrktree->Branch("trkstate_refz", &m_trkstate_refz );
+ m_fulltrktree->Branch("trkstate_tag", &m_trkstate_tag );
+ m_fulltrktree->Branch("trkstate_location", &m_trkstate_location );
+ m_fulltrktree->Branch("truthMC_tag", &m_truthMC_tag);
+ m_fulltrktree->Branch("truthMC_pid", &m_truthMC_pid);
+ m_fulltrktree->Branch("truthMC_px", &m_truthMC_px);
+ m_fulltrktree->Branch("truthMC_py", &m_truthMC_py);
+ m_fulltrktree->Branch("truthMC_pz", &m_truthMC_pz);
+ m_fulltrktree->Branch("truthMC_E", &m_truthMC_E);
+ m_fulltrktree->Branch("truthMC_EPx", &m_truthMC_EPx);
+ m_fulltrktree->Branch("truthMC_EPy", &m_truthMC_EPy);
+ m_fulltrktree->Branch("truthMC_EPz", &m_truthMC_EPz);
+ m_fulltrktree->Branch("truthMC_weight", &m_truthMC_weight);
return GaudiAlgorithm::initialize();
}
@@ -106,56 +186,58 @@ StatusCode ReadDigiAlg::execute()
}catch(GaudiException &e){
debug()<<"MC Particle is not available "<<endmsg;
}
+ m_mctree->Fill();
- try{
- if(m_VTXTrackerHitColHdl.get())
- N_VTXhit = m_VTXTrackerHitColHdl.get()->size();
- }catch(GaudiException &e){
- debug()<<"VTX hit is not available "<<endmsg;
- }
-
- try{
- if(m_SITTrackerHitColHdl.get())
- N_SIThit = m_SITTrackerHitColHdl.get()->size();
- }catch(GaudiException &e){
- debug()<<"SIT hit is not available "<<endmsg;
- }
-
- try{
- if(m_TPCTrackerHitColHdl.get())
- N_TPChit = m_TPCTrackerHitColHdl.get()->size();
- }catch(GaudiException &e){
- debug()<<"TPC hit is not available "<<endmsg;
- }
-
- try{
- if(m_SETTrackerHitColHdl.get())
- N_SEThit = m_SETTrackerHitColHdl.get()->size();
- }catch(GaudiException &e){
- debug()<<"SET hit is not available "<<endmsg;
- }
-
- try{
- if(m_FTDSpacePointColHdl.get())
- N_FTDhit = m_FTDSpacePointColHdl.get()->size();
- }catch(GaudiException &e){
- debug()<<"FDT space point is not available "<<endmsg;
- }
-
- try{
- if(m_SITRawColHdl.get())
- N_SITrawhit = m_SITRawColHdl.get()->size();
- }catch(GaudiException &e){
- debug()<<"SIT raw hit is not available "<<endmsg;
- }
-
- try{
- if(m_SETRawColHdl.get())
- N_SETrawhit = m_SETRawColHdl.get()->size();
- }catch(GaudiException &e){
- debug()<<"SET raw hit is not available "<<endmsg;
- }
+ //try{
+ //if(m_VTXTrackerHitColHdl.get())
+ // N_VTXhit = m_VTXTrackerHitColHdl.get()->size();
+ //}catch(GaudiException &e){
+ // debug()<<"VTX hit is not available "<<endmsg;
+ //}
+
+ //try{
+ //if(m_SITTrackerHitColHdl.get())
+ // N_SIThit = m_SITTrackerHitColHdl.get()->size();
+ //}catch(GaudiException &e){
+ // debug()<<"SIT hit is not available "<<endmsg;
+ //}
+
+ //try{
+ //if(m_TPCTrackerHitColHdl.get())
+ // N_TPChit = m_TPCTrackerHitColHdl.get()->size();
+ //}catch(GaudiException &e){
+ // debug()<<"TPC hit is not available "<<endmsg;
+ //}
+
+ //try{
+ //if(m_SETTrackerHitColHdl.get())
+ // N_SEThit = m_SETTrackerHitColHdl.get()->size();
+ //}catch(GaudiException &e){
+ // debug()<<"SET hit is not available "<<endmsg;
+ //}
+
+ //try{
+ //if(m_FTDSpacePointColHdl.get())
+ // N_FTDhit = m_FTDSpacePointColHdl.get()->size();
+ //}catch(GaudiException &e){
+ // debug()<<"FDT space point is not available "<<endmsg;
+ //}
+
+ //try{
+ //if(m_SITRawColHdl.get())
+ // N_SITrawhit = m_SITRawColHdl.get()->size();
+ //}catch(GaudiException &e){
+ // debug()<<"SIT raw hit is not available "<<endmsg;
+ //}
+
+ //try{
+ //if(m_SETRawColHdl.get())
+ // N_SETrawhit = m_SETRawColHdl.get()->size();
+ //}catch(GaudiException &e){
+ // debug()<<"SET raw hit is not available "<<endmsg;
+ //}
+ Clear();
try{
auto const_SiTrkCol = m_SiTrk.get();
if(const_SiTrkCol){
@@ -164,6 +246,23 @@ StatusCode ReadDigiAlg::execute()
auto m_trk = const_SiTrkCol->at(i);
if(m_trk.trackStates_size()==0) continue;
if(m_trk.trackerHits_size()==0) continue;
+
+ for(int istat=0; istat<m_trk.trackStates_size(); istat++){
+ edm4hep::TrackState m_trkstate = m_trk.getTrackStates(istat);
+
+ m_trkstate_d0.push_back( m_trkstate.D0 );
+ m_trkstate_z0.push_back( m_trkstate.Z0 );
+ m_trkstate_phi.push_back( m_trkstate.phi );
+ m_trkstate_tanL.push_back( m_trkstate.tanLambda );
+ //m_trkstate_kappa.push_back( m_trkstate.Kappa);
+ m_trkstate_omega.push_back( m_trkstate.omega );
+ m_trkstate_refx.push_back( m_trkstate.referencePoint.x );
+ m_trkstate_refy.push_back( m_trkstate.referencePoint.y );
+ m_trkstate_refz.push_back( m_trkstate.referencePoint.z );
+ m_trkstate_location.push_back( m_trkstate.location );
+ m_trkstate_tag.push_back(i);
+ }
+
edm4hep::TrackState m_trkstate = m_trk.getTrackStates(0);
HelixClassD * TrkInit_Helix = new HelixClassD();
@@ -171,21 +270,11 @@ StatusCode ReadDigiAlg::execute()
TVector3 TrkMom(TrkInit_Helix->getMomentum()[0],TrkInit_Helix->getMomentum()[1],TrkInit_Helix->getMomentum()[2]);
int NTrkHit = m_trk.trackerHits_size();
- TVector3 StartPointPos = TVector3((m_trk.getTrackerHits(0)).getPosition().x,(m_trk.getTrackerHits(0)).getPosition().y,(m_trk.getTrackerHits(0)).getPosition().z);
- TVector3 EndPointPos = TVector3((m_trk.getTrackerHits(NTrkHit - 1)).getPosition().x,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().y,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().z);
-
m_trk_Nhit.push_back(NTrkHit);
m_trk_px.push_back(TrkMom.x());
m_trk_py.push_back(TrkMom.y());
m_trk_pz.push_back(TrkMom.z());
m_trk_p.push_back(TrkMom.Mag());
- m_trk_endpoint_x.push_back(EndPointPos.x());
- m_trk_endpoint_y.push_back(EndPointPos.y());
- m_trk_endpoint_z.push_back(EndPointPos.z());
- m_trk_startpoint_x.push_back(StartPointPos.x());
- m_trk_startpoint_y.push_back(StartPointPos.y());
- m_trk_startpoint_z.push_back(StartPointPos.z());
- m_trk_type.push_back(1);
delete TrkInit_Helix;
}
@@ -193,37 +282,65 @@ StatusCode ReadDigiAlg::execute()
}catch(GaudiException &e){
debug()<<"Si track is not available "<<endmsg;
}
+ m_sitrktree->Fill();
+ Clear();
try{
auto const_TPCTrkCol = m_TPCTrk.get();
+ auto const_TPCTrkLinkCol = m_TPCTrkAssoHdl.get();
if(const_TPCTrkCol){
N_TPCTrk = const_TPCTrkCol->size();
for(int i=0; i<N_TPCTrk; i++){
auto m_trk = const_TPCTrkCol->at(i);
if(m_trk.trackStates_size()==0) continue;
if(m_trk.trackerHits_size()==0) continue;
+ int NTrkHit = m_trk.trackerHits_size();
+
+ for(int istat=0; istat<m_trk.trackStates_size(); istat++){
+ edm4hep::TrackState m_trkstate = m_trk.getTrackStates(istat);
+
+ m_trkstate_d0.push_back( m_trkstate.D0 );
+ m_trkstate_z0.push_back( m_trkstate.Z0 );
+ m_trkstate_phi.push_back( m_trkstate.phi );
+ m_trkstate_tanL.push_back( m_trkstate.tanLambda );
+ //m_trkstate_kappa.push_back( m_trkstate.Kappa);
+ m_trkstate_omega.push_back( m_trkstate.omega );
+ m_trkstate_refx.push_back( m_trkstate.referencePoint.x );
+ m_trkstate_refy.push_back( m_trkstate.referencePoint.y );
+ m_trkstate_refz.push_back( m_trkstate.referencePoint.z );
+ m_trkstate_location.push_back( m_trkstate.location );
+ m_trkstate_tag.push_back(i);
+ }
+
+ if(const_TPCTrkLinkCol){
+ for(auto iter: *const_TPCTrkLinkCol){
+ if(iter.getRec()==m_trk){
+ edm4hep::MCParticle m_mcp = iter.getSim();
+ m_truthMC_pid.push_back(m_mcp.getPDG());
+ m_truthMC_px.push_back(m_mcp.getMomentum().x);
+ m_truthMC_py.push_back(m_mcp.getMomentum().y);
+ m_truthMC_pz.push_back(m_mcp.getMomentum().z);
+ m_truthMC_E.push_back(m_mcp.getEnergy());
+ m_truthMC_EPx.push_back(m_mcp.getEndpoint().x);
+ m_truthMC_EPy.push_back(m_mcp.getEndpoint().y);
+ m_truthMC_EPz.push_back(m_mcp.getEndpoint().z);
+ m_truthMC_weight.push_back(iter.getWeight()/NTrkHit);
+ m_truthMC_tag.push_back(i);
+ }
+ }
+ }
+
edm4hep::TrackState m_trkstate = m_trk.getTrackStates(0);
HelixClassD * TrkInit_Helix = new HelixClassD();
TrkInit_Helix->Initialize_Canonical(m_trkstate.phi, m_trkstate.D0, m_trkstate.Z0, m_trkstate.omega, m_trkstate.tanLambda, _Bfield);
TVector3 TrkMom(TrkInit_Helix->getMomentum()[0],TrkInit_Helix->getMomentum()[1],TrkInit_Helix->getMomentum()[2]);
- int NTrkHit = m_trk.trackerHits_size();
- TVector3 StartPointPos = TVector3((m_trk.getTrackerHits(0)).getPosition().x,(m_trk.getTrackerHits(0)).getPosition().y,(m_trk.getTrackerHits(0)).getPosition().z);
- TVector3 EndPointPos = TVector3((m_trk.getTrackerHits(NTrkHit - 1)).getPosition().x,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().y,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().z);
-
m_trk_Nhit.push_back(NTrkHit);
m_trk_px.push_back(TrkMom.x());
m_trk_py.push_back(TrkMom.y());
m_trk_pz.push_back(TrkMom.z());
m_trk_p.push_back(TrkMom.Mag());
- m_trk_endpoint_x.push_back(EndPointPos.x());
- m_trk_endpoint_y.push_back(EndPointPos.y());
- m_trk_endpoint_z.push_back(EndPointPos.z());
- m_trk_startpoint_x.push_back(StartPointPos.x());
- m_trk_startpoint_y.push_back(StartPointPos.y());
- m_trk_startpoint_z.push_back(StartPointPos.z());
- m_trk_type.push_back(2);
delete TrkInit_Helix;
}
@@ -231,45 +348,93 @@ StatusCode ReadDigiAlg::execute()
}catch(GaudiException &e){
debug()<<"TPC track is not available "<<endmsg;
}
+ m_tpctrktree->Fill();
+ Clear();
try{
auto const_fullTrkCol = m_fullTrk.get();
+ auto const_fullTrkLinkCol = m_fullTrkAssoHdl.get();
if(const_fullTrkCol){
N_fullTrk = const_fullTrkCol->size();
for(int i=0; i<N_fullTrk; i++){
auto m_trk = const_fullTrkCol->at(i);
if(m_trk.trackStates_size()==0) continue;
if(m_trk.trackerHits_size()==0) continue;
- edm4hep::TrackState m_trkstate = m_trk.getTrackStates(0);
+ int NTrkHit = m_trk.trackerHits_size();
- HelixClassD * TrkInit_Helix = new HelixClassD();
- TrkInit_Helix->Initialize_Canonical(m_trkstate.phi, m_trkstate.D0, m_trkstate.Z0, m_trkstate.omega, m_trkstate.tanLambda, _Bfield);
- TVector3 TrkMom(TrkInit_Helix->getMomentum()[0],TrkInit_Helix->getMomentum()[1],TrkInit_Helix->getMomentum()[2]);
+ for(int istat=0; istat<m_trk.trackStates_size(); istat++){
+ edm4hep::TrackState m_trkstate = m_trk.getTrackStates(istat);
+
+ m_trkstate_d0.push_back( m_trkstate.D0 );
+ m_trkstate_z0.push_back( m_trkstate.Z0 );
+ m_trkstate_phi.push_back( m_trkstate.phi );
+ m_trkstate_tanL.push_back( m_trkstate.tanLambda );
+ //m_trkstate_kappa.push_back( m_trkstate.Kappa);
+ m_trkstate_omega.push_back( m_trkstate.omega );
+ m_trkstate_refx.push_back( m_trkstate.referencePoint.x );
+ m_trkstate_refy.push_back( m_trkstate.referencePoint.y );
+ m_trkstate_refz.push_back( m_trkstate.referencePoint.z );
+ m_trkstate_location.push_back( m_trkstate.location );
+ m_trkstate_tag.push_back(i);
+ }
+
+ int NSihit = 0;
+ int NTPChit = 0;
+ int Nelse = 0;
+ for(int ihit=0; ihit<NTrkHit; ihit++){
+ auto m_hit = m_trk.getTrackerHits(ihit);
+ if(m_hit.getType()==8) NSihit++;
+ else if(m_hit.getType()==0) NTPChit++;
+ else Nelse++;
+ }
+ //cout<<"In track #"<<i<<": Nhit "<<NTrkHit<<", Si Hit "<<NSihit<<", TPC hit "<<NTPChit<<", else "<<Nelse<<endl;
+
+ if(const_fullTrkLinkCol){
+ for(auto iter: *const_fullTrkLinkCol){
+ if(iter.getRec()==m_trk){
+ edm4hep::MCParticle m_mcp = iter.getSim();
+ m_truthMC_pid.push_back(m_mcp.getPDG());
+ m_truthMC_px.push_back(m_mcp.getMomentum().x);
+ m_truthMC_py.push_back(m_mcp.getMomentum().y);
+ m_truthMC_pz.push_back(m_mcp.getMomentum().z);
+ m_truthMC_E.push_back(m_mcp.getEnergy());
+ m_truthMC_EPx.push_back(m_mcp.getEndpoint().x);
+ m_truthMC_EPy.push_back(m_mcp.getEndpoint().y);
+ m_truthMC_EPz.push_back(m_mcp.getEndpoint().z);
+ m_truthMC_weight.push_back(iter.getWeight()/NTrkHit);
+ m_truthMC_tag.push_back(i);
+ }
+ }
+ }
- int NTrkHit = m_trk.trackerHits_size();
- TVector3 StartPointPos = TVector3((m_trk.getTrackerHits(0)).getPosition().x,(m_trk.getTrackerHits(0)).getPosition().y,(m_trk.getTrackerHits(0)).getPosition().z);
- TVector3 EndPointPos = TVector3((m_trk.getTrackerHits(NTrkHit - 1)).getPosition().x,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().y,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().z);
+ edm4hep::TrackState m_trkstate = m_trk.getTrackStates(0);
+ if(m_trkstate.location!=1){
+ std::cout<<"ERROR: first track state is not IP! Will use this for track momentum "<<std::endl;
+ }
+
+ double m_pt = (0.299792458*_Bfield)/fabs(m_trkstate.omega*1000.);
+ double m_phi = m_trkstate.phi;
+ double m_pz = m_trkstate.tanLambda * m_pt;
+ TVector3 TrkMom(m_pt*cos(m_phi), m_pt*sin(m_phi), m_pz);
+ //HelixClassD * TrkInit_Helix = new HelixClassD();
+ //TrkInit_Helix->Initialize_Canonical(m_trkstate.phi, m_trkstate.D0, m_trkstate.Z0, m_trkstate.omega, m_trkstate.tanLambda, _Bfield);
+ //TVector3 TrkMom(TrkInit_Helix->getMomentum()[0],TrkInit_Helix->getMomentum()[1],TrkInit_Helix->getMomentum()[2]);
m_trk_Nhit.push_back(NTrkHit);
+ m_trk_type.push_back(m_trk.getType());
m_trk_px.push_back(TrkMom.x());
m_trk_py.push_back(TrkMom.y());
m_trk_pz.push_back(TrkMom.z());
m_trk_p.push_back(TrkMom.Mag());
- m_trk_endpoint_x.push_back(EndPointPos.x());
- m_trk_endpoint_y.push_back(EndPointPos.y());
- m_trk_endpoint_z.push_back(EndPointPos.z());
- m_trk_startpoint_x.push_back(StartPointPos.x());
- m_trk_startpoint_y.push_back(StartPointPos.y());
- m_trk_startpoint_z.push_back(StartPointPos.z());
- m_trk_type.push_back(3);
- delete TrkInit_Helix;
+ //delete TrkInit_Helix;
}
}
}catch(GaudiException &e){
debug()<<"Full track is not available "<<endmsg;
}
+ m_fulltrktree->Fill();
/* const edm4hep::SimCalorimeterHitCollection* ECalBHitCol = m_ECalBarrelHitCol.get();
Nhit_EcalB = ECalBHitCol->size();
@@ -330,8 +495,6 @@ StatusCode ReadDigiAlg::execute()
}
*/
- m_mctree->Fill();
- m_trktree->Fill();
_nEvt++;
return StatusCode::SUCCESS;
}
@@ -341,7 +504,9 @@ StatusCode ReadDigiAlg::finalize()
debug() << "begin finalize ReadDigiAlg" << endmsg;
m_wfile->cd();
m_mctree->Write();
- m_trktree->Write();
+ m_sitrktree->Write();
+ m_tpctrktree->Write();
+ m_fulltrktree->Write();
m_wfile->Close();
return GaudiAlgorithm::finalize();
@@ -364,25 +529,40 @@ StatusCode ReadDigiAlg::Clear()
N_SiTrk = -99;
N_TPCTrk = -99;
N_fullTrk = -99;
- N_VTXhit = -99;
- N_SIThit = -99;
- N_TPChit = -99;
- N_SEThit = -99;
- N_FTDhit = -99;
- N_SITrawhit = -99;
- N_SETrawhit = -99;
+ //N_VTXhit = -99;
+ //N_SIThit = -99;
+ //N_TPChit = -99;
+ //N_SEThit = -99;
+ //N_FTDhit = -99;
+ //N_SITrawhit = -99;
+ //N_SETrawhit = -99;
m_trk_px.clear();
m_trk_py.clear();
m_trk_pz.clear();
m_trk_p.clear();
- m_trk_endpoint_x.clear();
- m_trk_endpoint_y.clear();
- m_trk_endpoint_z.clear();
- m_trk_startpoint_x.clear();
- m_trk_startpoint_y.clear();
- m_trk_startpoint_z.clear();
- m_trk_type.clear();
m_trk_Nhit.clear();
+ m_trk_type.clear();
+ m_trkstate_d0.clear();
+ m_trkstate_z0.clear();
+ m_trkstate_phi.clear();
+ m_trkstate_tanL.clear();
+ m_trkstate_omega.clear();
+ m_trkstate_kappa.clear();
+ m_trkstate_refx.clear();
+ m_trkstate_refy.clear();
+ m_trkstate_refz.clear();
+ m_trkstate_tag.clear();
+ m_trkstate_location.clear();
+ m_truthMC_tag.clear();
+ m_truthMC_pid.clear();
+ m_truthMC_px.clear();
+ m_truthMC_py.clear();
+ m_truthMC_pz.clear();
+ m_truthMC_E.clear();
+ m_truthMC_EPx.clear();
+ m_truthMC_EPy.clear();
+ m_truthMC_EPz.clear();
+ m_truthMC_weight.clear();
return StatusCode::SUCCESS;
}
diff --git a/Analysis/ReadDigi/src/ReadDigiAlg.h b/Analysis/ReadDigi/src/ReadDigiAlg.h
index abfcb759..ea79a581 100644
--- a/Analysis/ReadDigi/src/ReadDigiAlg.h
+++ b/Analysis/ReadDigi/src/ReadDigiAlg.h
@@ -9,6 +9,7 @@
#include "edm4hep/TrackerHitCollection.h"
#include "edm4hep/TrackCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
#include "HelixClassD.hh"
#include "TFile.h"
@@ -32,20 +33,23 @@ public :
private :
DataHandle<edm4hep::MCParticleCollection> m_MCParticleCol{"MCParticle", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_SITRawColHdl{"SITTrackerHits", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_SETRawColHdl{"SETTrackerHits", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_FTDRawColHdl{"FTDStripTrackerHits", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_SITRawColHdl{"SITTrackerHits", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_SETRawColHdl{"SETTrackerHits", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_FTDRawColHdl{"FTDStripTrackerHits", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_VTXTrackerHitColHdl{"VTXTrackerHits", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_SITTrackerHitColHdl{"SITSpacePoints", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_TPCTrackerHitColHdl{"TPCTrackerHits", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_SETTrackerHitColHdl{"SETSpacePoints", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_FTDSpacePointColHdl{"FTDSpacePoints", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_FTDPixelTrackerHitColHdl{"FTDPixelTrackerHits", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_VTXTrackerHitColHdl{"VTXTrackerHits", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_SITTrackerHitColHdl{"SITSpacePoints", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_TPCTrackerHitColHdl{"TPCTrackerHits", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_SETTrackerHitColHdl{"SETSpacePoints", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_FTDSpacePointColHdl{"FTDSpacePoints", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::TrackerHitCollection> m_FTDPixelTrackerHitColHdl{"FTDPixelTrackerHits", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::TrackCollection> m_SiTrk{"SiTracks", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::TrackCollection> m_TPCTrk{"ClupatraTracks", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackCollection> m_fullTrk{"MarlinTrkTracks", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackCollection> m_fullTrk{"CompleteTracks", Gaudi::DataHandle::Reader, this};
+
+ DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> m_TPCTrkAssoHdl{"ClupatraTracksParticleAssociation", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> m_fullTrkAssoHdl{"CompleteTracksParticleAssociation", Gaudi::DataHandle::Reader, this};
//DataHandle<edm4hep::SimCalorimeterHitCollection> m_ECalBarrelHitCol{"EcalBarrelCollection", Gaudi::DataHandle::Reader, this};
//DataHandle<edm4hep::SimCalorimeterHitCollection> m_ECalEndcapHitCol{"EcalEndcapCollection", Gaudi::DataHandle::Reader, this};
@@ -62,7 +66,7 @@ private :
int _nEvt;
TFile *m_wfile;
- TTree *m_mctree, *m_trktree;
+ TTree *m_mctree, *m_sitrktree, *m_tpctrktree, *m_fulltrktree;
//MCParticle
int N_MCP;
@@ -71,10 +75,13 @@ private :
//Tracker
int N_SiTrk, N_TPCTrk, N_fullTrk;
- int N_VTXhit, N_SIThit, N_TPChit, N_SEThit, N_FTDhit, N_SITrawhit, N_SETrawhit;
+ //int N_VTXhit, N_SIThit, N_TPChit, N_SEThit, N_FTDhit, N_SITrawhit, N_SETrawhit;
FloatVec m_trk_px, m_trk_py, m_trk_pz, m_trk_p;
- FloatVec m_trk_endpoint_x, m_trk_endpoint_y, m_trk_endpoint_z, m_trk_startpoint_x, m_trk_startpoint_y, m_trk_startpoint_z;
- IntVec m_trk_type, m_trk_Nhit;
+ FloatVec m_trkstate_d0, m_trkstate_z0, m_trkstate_phi, m_trkstate_tanL, m_trkstate_omega, m_trkstate_kappa;
+ FloatVec m_trkstate_refx, m_trkstate_refy, m_trkstate_refz;
+ IntVec m_trkstate_tag, m_trkstate_location, m_trk_Nhit, m_trk_type;
+ FloatVec m_truthMC_tag, m_truthMC_pid, m_truthMC_px, m_truthMC_py, m_truthMC_pz, m_truthMC_E;
+ FloatVec m_truthMC_EPx, m_truthMC_EPy, m_truthMC_EPz, m_truthMC_weight;
//ECal
diff --git a/Reconstruction/RecPFACyber/script/digi.py b/Reconstruction/RecPFACyber/script/digi.py
index a7736873..c3ea459f 100644
--- a/Reconstruction/RecPFACyber/script/digi.py
+++ b/Reconstruction/RecPFACyber/script/digi.py
@@ -3,7 +3,7 @@ import os
from Gaudi.Configuration import *
from Configurables import k4DataSvc
-dsvc = k4DataSvc("EventDataSvc", input="Tracking_TDR_o1_v01_E240_nnHgg.root")
+dsvc = k4DataSvc("EventDataSvc", input="Sim_TDR_o1_v01.root")
from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
seed = [12340]
@@ -52,8 +52,6 @@ podioinput = PodioInput("PodioReader", collections=[
"EcalBarrelContributionCollection",
"HcalBarrelCollection",
"HcalBarrelContributionCollection",
- "CompleteTracks",
- "CompleteTracksParticleAssociation"
])
########## Digitalization ################
@@ -75,7 +73,7 @@ EcalDigi.TimeResolution = 0.5 #unit: ns
EcalDigi.ChargeThresholdFrac = 0.05
EcalDigi.EcalMIP_Thre = 0.05
EcalDigi.Debug=1
-EcalDigi.WriteNtuple = 1
+EcalDigi.WriteNtuple = 0
EcalDigi.OutFileName = "Digi_ECAL.root"
#########################################
@@ -108,22 +106,33 @@ HcalDigi.ADCBaselineSigmaLG = 0.
HcalDigi.ADCHLRatio = 1
HcalDigi.ADCSwitch = 1e7
HcalDigi.ADCLimit = 1e7
-HcalDigi.WriteNtuple = 1
+HcalDigi.WriteNtuple = 0
HcalDigi.OutFileName = "Digi_HCAL.root"
# output
from Configurables import PodioOutput
out = PodioOutput("outputalg")
-out.filename = "CaloDigi_TDR_o1_v01_E240_nnHgg.root"
-out.outputCommands = ["keep *"]
+out.filename = "CaloDigi_TDR_o1_v01.root"
+out.outputCommands = ["drop *",
+ "keep MCParticle",
+ "keep VXDCollection",
+ "keep SITCollection",
+ "keep TPCCollection",
+ "keep OTKBarrelCollection",
+ "keep FTDCollection",
+ "keep ECALBarrel",
+ "keep HCALBarrel",
+ "keep ECALBarrelParticleAssoCol",
+ "keep HCALBarrelParticleAssoCol" ]
+
# ApplicationMgr
from Configurables import ApplicationMgr
mgr = ApplicationMgr(
TopAlg = [podioinput, EcalDigi, HcalDigi, out],
EvtSel = 'NONE',
- EvtMax = 3,
+ EvtMax = 10,
ExtSvc = [dsvc, rndmengine, rndmgensvc, geosvc],
HistogramPersistency = 'ROOT',
OutputLevel = ERROR
diff --git a/Reconstruction/RecPFACyber/script/rec.py b/Reconstruction/RecPFACyber/script/rec.py
index 6e273e20..e844f945 100644
--- a/Reconstruction/RecPFACyber/script/rec.py
+++ b/Reconstruction/RecPFACyber/script/rec.py
@@ -20,7 +20,7 @@ geomsvc.compact = geometry_path
########### k4DataSvc ####################
from Configurables import k4DataSvc
-podioevent = k4DataSvc("EventDataSvc", input="CaloDigi_TDR_o1_v01_E240_nnHgg.root")
+podioevent = k4DataSvc("EventDataSvc", input="Tracking_TDR_o1_v01.root")
##########################################
########## CEPCSWData #################
@@ -36,12 +36,8 @@ crystalecalcorr.CorrectionFile = os.path.join(cepcswdatatop, "CEPCSWData/offline
########## Podio Input ###################
from Configurables import PodioInput
inp = PodioInput("InputReader")
-inp.collections = [ "EcalBarrelCollection",
- "EcalBarrelContributionCollection",
- "ECALBarrel",
+inp.collections = [ "ECALBarrel",
"ECALBarrelParticleAssoCol",
- "HcalBarrelCollection",
- "HcalBarrelContributionCollection",
"HCALBarrel",
"HCALBarrelParticleAssoCol",
"MCParticle",
@@ -58,7 +54,7 @@ CyberPFAlg.BField = 3.
CyberPFAlg.Debug = 0
CyberPFAlg.SkipEvt = 0
CyberPFAlg.WriteAna = 1
-CyberPFAlg.AnaFileName = "RecAnaTuple_TDR_o1_v01_E240_nnHgg.root"
+CyberPFAlg.AnaFileName = "RecAnaTuple_TDR_o1_v01_mu.root"
CyberPFAlg.UseTruthTrack = 0
CyberPFAlg.EcalGlobalCalib = 1.05
CyberPFAlg.HcalGlobalCalib = 4.5
@@ -133,7 +129,7 @@ CyberPFAlg.AlgParValues = [ ["BarCol","Cluster1DCol","HalfClusterCol"],#1
##############################################################################
from Configurables import PodioOutput
out = PodioOutput("outputalg")
-out.filename = "Rec_TDR_o1_v01_E240_nnHgg.root"
+out.filename = "Rec_TDR_o1_v01.root"
out.outputCommands = ["keep *"]
@@ -143,7 +139,7 @@ from Configurables import ApplicationMgr
ApplicationMgr(
TopAlg=[inp, CyberPFAlg, out ],
EvtSel="NONE",
- EvtMax=3,
+ EvtMax=10,
ExtSvc=[podioevent, geomsvc],
#OutputLevel=DEBUG
)
diff --git a/Reconstruction/RecPFACyber/script/sim.py b/Reconstruction/RecPFACyber/script/sim.py
index 57639491..a8b14c52 100644
--- a/Reconstruction/RecPFACyber/script/sim.py
+++ b/Reconstruction/RecPFACyber/script/sim.py
@@ -42,16 +42,16 @@ from Configurables import HepMCRdr
from Configurables import GenPrinter
#gun = GtGunTool("GtGunTool")
-#gun.PositionXs = [0]
-#gun.PositionYs = [0]
-#gun.PositionZs = [0]
-#gun.Particles = ["pi-"]
-#gun.EnergyMins = [10]
-#gun.EnergyMaxs = [10]
-#gun.ThetaMins = [60]
-#gun.ThetaMaxs = [120]
-#gun.PhiMins = [0]
-#gun.PhiMaxs = [360]
+#gun.PositionXs = [0, 0]
+#gun.PositionYs = [0, 0]
+#gun.PositionZs = [0, 0]
+#gun.Particles = ["mu-", "mu+"]
+#gun.EnergyMins = [5, 10]
+#gun.EnergyMaxs = [5, 10]
+#gun.ThetaMins = [60, 60]
+#gun.ThetaMaxs = [120, 120]
+#gun.PhiMins = [0, 0]
+#gun.PhiMaxs = [360, 360]
#genprinter = GenPrinter("GenPrinter")
stdheprdr = StdHepRdr("StdHepRdr")
@@ -97,7 +97,7 @@ cal_sensdettool.CalNamesMergeDisable = ["EcalBarrel", "CaloDetectorEndcap", "Hca
# output
from Configurables import PodioOutput
out = PodioOutput("outputalg")
-out.filename = "Sim_TDR_o1_v01_E240_nnHgg.root"
+out.filename = "Sim_TDR_o1_v01.root"
out.outputCommands = ["keep *"]
# ApplicationMgr
@@ -105,7 +105,7 @@ from Configurables import ApplicationMgr
mgr = ApplicationMgr(
TopAlg = [genalg, detsimalg, out],
EvtSel = 'NONE',
- EvtMax = 3,
+ EvtMax = 10,
ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc],
HistogramPersistency = 'ROOT',
OutputLevel = INFO
diff --git a/Reconstruction/RecPFACyber/script/tracking.py b/Reconstruction/RecPFACyber/script/tracking.py
index a6d7538f..7e07f554 100644
--- a/Reconstruction/RecPFACyber/script/tracking.py
+++ b/Reconstruction/RecPFACyber/script/tracking.py
@@ -3,7 +3,7 @@ import os
from Gaudi.Configuration import *
from Configurables import k4DataSvc
-dsvc = k4DataSvc("EventDataSvc", input="Sim_TDR_o1_v01_E240_nnHgg.root")
+dsvc = k4DataSvc("EventDataSvc", input="CaloDigi_TDR_o1_v01.root")
from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
seed = [12340]
@@ -51,7 +51,6 @@ podioinput = PodioInput("PodioReader", collections=[
"VXDCollection",
"SITCollection",
"TPCCollection",
-# "SETCollection",
"OTKBarrelCollection",
"FTDCollection"
])
@@ -249,18 +248,34 @@ tmt.MuonTagEfficiency = 0.95 # muon true tag efficiency, default is 1.0 (100%)
tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta)
#tmt.OutputLevel = DEBUG
+from Configurables import ReadDigiAlg
+readtrk = ReadDigiAlg("ReadDigiAlg")
+readtrk.SiTracks = "SubsetTracks"
+readtrk.TPCTracks = "ClupatraTracks"
+readtrk.FullTracks = "CompleteTracks"
+readtrk.TPCTracksAssociation = "ClupatraTracksParticleAssociation"
+readtrk.FullTracksAssociation = "CompleteTracksParticleAssociation"
+readtrk.OutFileName = "TrackAnaTuple_mu.root"
+
# output
from Configurables import PodioOutput
out = PodioOutput("outputalg")
-out.filename = "Tracking_TDR_o1_v01_E240_nnHgg.root"
-out.outputCommands = ["keep *"]
+out.filename = "Tracking_TDR_o1_v01.root"
+out.outputCommands = ["drop *",
+ "keep ECALBarrel",
+ "keep ECALBarrelParticleAssoCol",
+ "keep HCALBarrel",
+ "keep HCALBarrelParticleAssoCol",
+ "keep MCParticle",
+ "keep CompleteTracks",
+ "keep CompleteTracksParticleAssociation" ]
# ApplicationMgr
from Configurables import ApplicationMgr
mgr = ApplicationMgr(
- TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tmt, out],
+ TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tmt, readtrk, out],
EvtSel = 'NONE',
- EvtMax = 3,
+ EvtMax = 10,
ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc],
#HistogramPersistency = 'ROOT',
OutputLevel = INFO
diff --git a/Reconstruction/RecPFACyber/src/Objects/Track.cc b/Reconstruction/RecPFACyber/src/Objects/Track.cc
index 204f3e3e..1d67ef0b 100644
--- a/Reconstruction/RecPFACyber/src/Objects/Track.cc
+++ b/Reconstruction/RecPFACyber/src/Objects/Track.cc
@@ -42,7 +42,7 @@ namespace Cyber{
std::vector<TrackState> trkStates = getTrackStates("Input");
float d0 = -99.;
for(auto it: trkStates){
- if(it.location==4 || it.location==1){ //Calorimeter(for real track) or IP (for truth track)
+ if(it.location==Cyber::TrackState::AtIP){
d0 = it.D0;
}
}
@@ -54,7 +54,7 @@ namespace Cyber{
std::vector<TrackState> trkStates = getTrackStates("Input");
float z0 = -99.;
for(auto it: trkStates){
- if(it.location==4 || it.location==1){ //Calorimeter(for real track) or IP (for truth track)
+ if(it.location==Cyber::TrackState::AtIP){
z0 = it.Z0;
}
}
@@ -67,7 +67,7 @@ namespace Cyber{
std::vector<TrackState> trkStates = getTrackStates("Input");
float pt = -99.;
for(auto it: trkStates){
- if(it.location==4 || it.location==1){ //Calorimeter(for real track) or IP (for truth track)
+ if(it.location==Cyber::TrackState::AtIP){
pt = 1./fabs(it.Kappa);
}
}
@@ -79,7 +79,7 @@ namespace Cyber{
std::vector<TrackState> trkStates = getTrackStates("Input");
float pz = -99.;
for(auto it: trkStates){
- if( (m_type!=0 && it.location==4) || (m_type==0 && it.location==1)){ //Calorimeter(for real track) or IP (for truth track)
+ if( it.location==Cyber::TrackState::AtIP ){
pz = it.tanLambda/fabs(it.Kappa);
}
}
@@ -93,7 +93,7 @@ namespace Cyber{
float pt = -99.;
float pz = -99.;
for(auto it: trkStates){
- if((m_type!=0 && it.location==4) || (m_type==0 && it.location==1)){ //Calorimeter(for real track) or IP (for truth track)
+ if(it.location==Cyber::TrackState::AtIP){
pt = 1./fabs(it.Kappa);
phi = it.phi0;
pz = it.tanLambda/fabs(it.Kappa);
@@ -108,7 +108,7 @@ namespace Cyber{
std::vector<TrackState> trkStates = getTrackStates("Input");
float omega = -99.;
for(auto it: trkStates){
- if((m_type!=0 && it.location==4) || (m_type==0 && it.location==1)){ //Calorimeter(for real track) or IP (for truth track)
+ if(it.location==Cyber::TrackState::AtIP){
omega = it.Omega;
}
}
--
GitLab