diff --git a/Reconstruction/RecPFACyber/CMakeLists.txt b/Reconstruction/RecPFACyber/CMakeLists.txt
index 5209dcef1784cfef3e0e2d3d8b95989c06c4a2f9..481f609d0dc3d6a55c0ecbc7d4db662506676675 100644
--- a/Reconstruction/RecPFACyber/CMakeLists.txt
+++ b/Reconstruction/RecPFACyber/CMakeLists.txt
@@ -23,6 +23,7 @@ gaudi_add_library(CrystalCaloRecLib
${LCIO_LIBRARIES}
${ROOT_LIBRARIES}
${TMVA_LIBRARIES}
+ EDM4CEPC::edm4cepc EDM4CEPC::edm4cepcDict
EDM4HEP::edm4hep EDM4HEP::edm4hepDict
${DD4hep_COMPONENT_LIBRARIES}
DetIdentifier
diff --git a/Reconstruction/RecPFACyber/include/Algorithm/AxisMergingAlg.h b/Reconstruction/RecPFACyber/include/Algorithm/AxisMergingAlg.h
index 50a5416c360cecb9897efbc36c629c532c0a7693..df153329979bedaf1a27b0692a874ec45116baf2 100644
--- a/Reconstruction/RecPFACyber/include/Algorithm/AxisMergingAlg.h
+++ b/Reconstruction/RecPFACyber/include/Algorithm/AxisMergingAlg.h
@@ -41,7 +41,18 @@ private:
std::vector<Cyber::CaloHalfCluster*> m_newAxisVCol;
static bool compLayer( const Cyber::CaloHalfCluster* sh1, const Cyber::CaloHalfCluster* sh2 )
- { return sh1->getBeginningDlayer() < sh2->getBeginningDlayer(); }
+ { if( sh1->getBeginningDlayer() != sh2->getBeginningDlayer() )
+ return sh1->getBeginningDlayer() < sh2->getBeginningDlayer();
+ else
+ return sh1->getEnergy() > sh2->getEnergy();
+ }
+ static bool compLayerPtr( const std::shared_ptr<Cyber::CaloHalfCluster> sh1, const std::shared_ptr<Cyber::CaloHalfCluster> sh2 )
+ { if( sh1->getBeginningDlayer() != sh2->getBeginningDlayer() )
+ return sh1->getBeginningDlayer() < sh2->getBeginningDlayer();
+ else
+ return sh1->getEnergy() > sh2->getEnergy();
+ }
+
};
diff --git a/Reconstruction/RecPFACyber/include/Algorithm/HcalClusteringAlg.h b/Reconstruction/RecPFACyber/include/Algorithm/HcalClusteringAlg.h
index bb915789ae0680790cceb4a3e607e91c5caf484f..41932c81c1205f16108509b960e09fb905084b66 100644
--- a/Reconstruction/RecPFACyber/include/Algorithm/HcalClusteringAlg.h
+++ b/Reconstruction/RecPFACyber/include/Algorithm/HcalClusteringAlg.h
@@ -25,7 +25,7 @@ public:
template<typename T1, typename T2> StatusCode Clustering(std::vector<T1*>& m_input, std::vector<std::shared_ptr<T2>>& m_output);
StatusCode LongiConeLinking( const std::map<int, std::vector<Cyber::CaloHit*> >& orderedHit, std::vector<std::shared_ptr<Cyber::Calo3DCluster> >& ClusterCol );
-
+ StatusCode MergeToCluster( std::vector<const Cyber::CaloHit*>& m_isohits, std::vector<std::shared_ptr<Cyber::Calo3DCluster>>& m_clusters );
private:
// std::vector<std::shared_ptr<Cyber::CaloHit>> m_hcalHits;
diff --git a/Reconstruction/RecPFACyber/include/CyberDataCol.h b/Reconstruction/RecPFACyber/include/CyberDataCol.h
index 5719eeb85ffe52086940cf42b31a93a8f33206b5..1d6854f6239eaa240bb5172984d17679f648b748 100644
--- a/Reconstruction/RecPFACyber/include/CyberDataCol.h
+++ b/Reconstruction/RecPFACyber/include/CyberDataCol.h
@@ -37,6 +37,11 @@
#include "edm4hep/Vertex.h"
#include "edm4hep/VertexCollection.h"
#include "edm4hep/ClusterCollection.h"
+#include "edm4cepc/RecTof.h"
+#include "edm4cepc/RecTofCollection.h"
+#include "edm4hep/RecDqdx.h"
+#include "edm4hep/RecDqdxCollection.h"
+#include "edm4hep/ParticleIDCollection.h"
#include "edm4hep/ReconstructedParticleCollection.h"
#include "edm4hep/MCRecoCaloAssociation.h"
#include "edm4hep/MCRecoTrackerAssociation.h"
@@ -77,5 +82,8 @@ public:
//Energy calibration service
SmartIF<ICrystalEcalSvc> EnergyCorrSvc;
+ //PID collections
+ edm4hep::RecTofCollection* tofCol = nullptr;
+ edm4hep::RecDqdxCollection* dNdxCol = nullptr;
};
#endif
diff --git a/Reconstruction/RecPFACyber/include/CyberPFAlg.h b/Reconstruction/RecPFACyber/include/CyberPFAlg.h
index c333c5a92588f3792f7df39078cdea3fbf680b35..6f1d195ee7fd38f79ca2e882d4d640249e881638 100644
--- a/Reconstruction/RecPFACyber/include/CyberPFAlg.h
+++ b/Reconstruction/RecPFACyber/include/CyberPFAlg.h
@@ -124,6 +124,8 @@ protected:
Gaudi::Property< std::string > name_MCParticleCol{ this, "MCParticleCollection", "MCParticle" };
Gaudi::Property< std::string > name_MCPTrkAssoCol{this, "MCRecoTrackParticleAssociationCollection", "MarlinTrkAssociation"};
Gaudi::Property< std::vector<std::string> > name_TrackCol{ this, "TrackCollections", {"MarlinTrkTracks"} };
+ Gaudi::Property< std::string > name_dNdxCol{this, "DndxCollection", "DndxTracks"};
+ Gaudi::Property< std::string > name_tofCol{this, "TOFCollection", "RecTofCollection"};
Gaudi::Property< std::vector<std::string> > name_EcalHits{ this, "ECalCaloHitCollections", {"ECALBarrel"} };
Gaudi::Property< std::vector<std::string> > name_EcalReadout{ this, "ECalReadOutNames", {"EcalBarrelCollection"} };
Gaudi::Property< std::vector<std::string> > name_EcalMCPAssociation{ this, "ECalMCPAssociationName", {"ECALBarrelParticleAssoCol"} };
@@ -143,6 +145,8 @@ protected:
//std::vector<CaloType*> r_HCalHitCols;
std::vector<CaloType*> r_CaloHitCols;
std::map<std::string, CaloParticleAssoType*> map_CaloMCPAssoCols;
+ DataHandle<edm4hep::RecTofCollection>* r_TofCol;
+ DataHandle<edm4hep::RecDqdxCollection>* r_dNdxCol;
//Global parameters.
Gaudi::Property<float> m_BField{this, "BField", 3., "Magnetic field"};
@@ -296,20 +300,21 @@ protected:
m_EcalClus_trk_d0, m_EcalClus_trk_z0, m_EcalClus_trk_phi, m_EcalClus_trk_tanL, m_EcalClus_trk_omega, m_EcalClus_trk_kappa,
m_EcalClus_truthMC_tag, m_EcalClus_truthMC_pid, m_EcalClus_truthMC_px, m_EcalClus_truthMC_py, m_EcalClus_truthMC_pz, m_EcalClus_truthMC_E,
m_EcalClus_truthMC_EPx, m_EcalClus_truthMC_EPy, m_EcalClus_truthMC_EPz, m_EcalClus_truthMC_weight;
- FloatVec m_HcalClus_x, m_HcalClus_y, m_HcalClus_z, m_HcalClus_E, m_HcalClus_nTrk, m_HcalClus_pTrk,
+ FloatVec m_HcalClus_x, m_HcalClus_y, m_HcalClus_z, m_HcalClus_E, m_HcalClus_nTrk, m_HcalClus_pTrk, m_HcalClus_nHit,
m_HcalClus_hit_x, m_HcalClus_hit_y, m_HcalClus_hit_z, m_HcalClus_hit_E, m_HcalClus_hit_tag,
m_HcalClus_truthMC_tag, m_HcalClus_truthMC_pid, m_HcalClus_truthMC_px, m_HcalClus_truthMC_py, m_HcalClus_truthMC_pz, m_HcalClus_truthMC_E,
m_HcalClus_truthMC_EPx, m_HcalClus_truthMC_EPy, m_HcalClus_truthMC_EPz, m_HcalClus_truthMC_weight;
- FloatVec m_SimpleHcalClus_x, m_SimpleHcalClus_y, m_SimpleHcalClus_z, m_SimpleHcalClus_E, m_SimpleHcalClus_nTrk, m_SimpleHcalClus_pTrk,
+ FloatVec m_SimpleHcalClus_x, m_SimpleHcalClus_y, m_SimpleHcalClus_z, m_SimpleHcalClus_E, m_SimpleHcalClus_nTrk, m_SimpleHcalClus_pTrk, m_SimpleHcalClus_nHit,
m_SimpleHcalClus_hit_x, m_SimpleHcalClus_hit_y, m_SimpleHcalClus_hit_z, m_SimpleHcalClus_hit_E, m_SimpleHcalClus_hit_tag,
m_SimpleHcalClus_truthMC_tag, m_SimpleHcalClus_truthMC_pid, m_SimpleHcalClus_truthMC_px, m_SimpleHcalClus_truthMC_py, m_SimpleHcalClus_truthMC_pz, m_SimpleHcalClus_truthMC_E,
m_SimpleHcalClus_truthMC_EPx, m_SimpleHcalClus_truthMC_EPy, m_SimpleHcalClus_truthMC_EPz, m_SimpleHcalClus_truthMC_weight;
TTree *t_Track;
int m_Ntrk;
+ FloatVec m_trk_p, m_trk_px, m_trk_py, m_trk_pz, m_trk_truthweight;
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_type, m_Nhit;
+ IntVec m_trkstate_tag, m_trkstate_location, m_type, m_Nhit, m_pid, m_pid_truth;
FloatVec m_trkstate_x_ECAL, m_trkstate_y_ECAL, m_trkstate_z_ECAL, m_trkstate_x_HCAL, m_trkstate_y_HCAL, m_trkstate_z_HCAL;
IntVec m_trkstate_tag_ECAL, m_trkstate_tag_HCAL;
diff --git a/Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h b/Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h
index 7dff502962f7c70fc70547047c8dbb32a4e76d24..ca4dcf96d6a660fa494ec345bb7ff24707ab7c0d 100644
--- a/Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h
+++ b/Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h
@@ -100,7 +100,12 @@ namespace Cyber {
TrackFitInEcal* track = new TrackFitInEcal();
static bool compLayer( const Cyber::Calo1DCluster* hit1, const Cyber::Calo1DCluster* hit2 )
- { return hit1->getDlayer() < hit2->getDlayer(); }
+ { if(hit1->getDlayer() != hit2->getDlayer())
+ return hit1->getDlayer() < hit2->getDlayer();
+ else
+ return hit1->getEnergy() > hit2->getEnergy();
+ }
+
};
diff --git a/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h b/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h
index 5c1904eb4fa09bb3aa48c3c0ec24687e356364b1..5d4eeb8dfca89e1b9c56e9e7f66ad16721db2e5a 100644
--- a/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h
+++ b/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h
@@ -83,8 +83,8 @@ namespace Cyber{
static int Nmodule;
static int Nstave;
static int Nlayer;
- static int NbarPhi_odd[14];
- static int NbarPhi_even[14];
+ static int NbarPhi_odd[9];
+ static int NbarPhi_even[9];
static int NbarZ;
static float barsize;
static float ecal_innerR;
diff --git a/Reconstruction/RecPFACyber/include/Objects/Track.h b/Reconstruction/RecPFACyber/include/Objects/Track.h
index 6af8328f6ea26185941bdeb3b5d8ce8d3913b1cc..01212021c96a4091d140afb58b47d06f70daeb17 100644
--- a/Reconstruction/RecPFACyber/include/Objects/Track.h
+++ b/Reconstruction/RecPFACyber/include/Objects/Track.h
@@ -44,8 +44,10 @@ namespace Cyber {
void addLinkedMCP( std::pair<edm4hep::MCParticle, float> _pair ) { MCParticleWeight.push_back(_pair); }
void setLinkedMCP( std::vector<std::pair<edm4hep::MCParticle, float>> _pairVec ) { MCParticleWeight.clear(); MCParticleWeight = _pairVec; }
+ void setPID(int _pid) {m_pid = _pid; }
void setType(int _type) { m_type=_type; }
int getType() const { return m_type; }
+ int getPID() const { return m_pid; }
private:
edm4hep::Track m_track;
@@ -56,6 +58,7 @@ namespace Cyber {
std::vector<Cyber::CaloHalfCluster*> m_halfClusterVCol;
int m_type;
+ int m_pid;
std::vector< std::pair<edm4hep::MCParticle, float> > MCParticleWeight;
};
diff --git a/Reconstruction/RecPFACyber/include/Tools/TrackCreator.h b/Reconstruction/RecPFACyber/include/Tools/TrackCreator.h
index b4f468418f7dcaa4d21cd77f8e802b07635b1c46..e95fedda058bee3e454cb73b35f666e8c8681e49 100644
--- a/Reconstruction/RecPFACyber/include/Tools/TrackCreator.h
+++ b/Reconstruction/RecPFACyber/include/Tools/TrackCreator.h
@@ -31,6 +31,13 @@ namespace Cyber{
StatusCode Reset(){};
private:
+ const std::map<int, int> PDGIDs = {
+ {0, -11},
+ {1, -13},
+ {2, 211},
+ {3, 321},
+ {4, 2212},
+ };
const Cyber::Settings settings;
Cyber::Algorithm* m_TrkExtraAlg;
Cyber::Settings m_TrkExtraSettings;
diff --git a/Reconstruction/RecPFACyber/script/rec.py b/Reconstruction/RecPFACyber/script/rec.py
index 20706f0849b9d20141ca31e11860373d96799862..e28309eef53944cba599e6755235837575b4ffab 100644
--- a/Reconstruction/RecPFACyber/script/rec.py
+++ b/Reconstruction/RecPFACyber/script/rec.py
@@ -60,7 +60,7 @@ CyberPFAlg.Seed = 1024
CyberPFAlg.BField = 3.
CyberPFAlg.Debug = 0
CyberPFAlg.SkipEvt = 0
-CyberPFAlg.WriteAna = 1
+CyberPFAlg.WriteAna = 0
CyberPFAlg.AnaFileName = "RecAnaTuple_TDR_o1_v01.root"
CyberPFAlg.UseMCPTrack = 0
CyberPFAlg.UseTruthMatchTrack = 0
diff --git a/Reconstruction/RecPFACyber/script/sim.py b/Reconstruction/RecPFACyber/script/sim.py
index 2b10df3feca785ce91b5a2ad27c45429ca74b61e..2c5f1fe1b3df8ccab594010d79eb943d5c7b8d06 100644
--- a/Reconstruction/RecPFACyber/script/sim.py
+++ b/Reconstruction/RecPFACyber/script/sim.py
@@ -84,6 +84,18 @@ detsimalg.RootDetElem = "WorldDetElemTool"
from Configurables import TimeProjectionChamberSensDetTool
tpc_sensdettool = TimeProjectionChamberSensDetTool("TimeProjectionChamberSensDetTool")
tpc_sensdettool.TypeOption = 1
+tpc_sensdettool.DoHeedSim = True
+dedxoption = "TrackHeedSimTool"
+tpc_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import TrackHeedSimTool
+dedx_simtool = TrackHeedSimTool("TrackHeedSimTool")
+dedx_simtool.detector = "TPC"
+dedx_simtool.only_primary = False#True
+dedx_simtool.use_max_step = False#True
+dedx_simtool.max_step = 1#mm
+dedx_simtool.save_mc = True
+#dedx_simtool.OutputLevel = DEBUG
from Configurables import CalorimeterSensDetTool
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/AxisMergingAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/AxisMergingAlg.cpp
index ae2a1cbd486f60c88ee9de31533eacb5b83ae006..64b601c06576723a23f646ab75b8ab63d66c5756 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/AxisMergingAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/AxisMergingAlg.cpp
@@ -75,6 +75,7 @@ StatusCode AxisMergingAlg::RunAlgorithm( CyberDataCol& m_datacol ){
if(m_newAxisUCol.size()<2){ //No need to merge, save into OutputAxis.
std::vector<const Cyber::CaloHalfCluster*> tmp_axisCol; tmp_axisCol.clear();
for(int ic=0; ic<m_newAxisUCol.size(); ic++) tmp_axisCol.push_back(m_newAxisUCol[ic]);
+ std::sort( tmp_axisCol.begin(), tmp_axisCol.end(), compLayer );
p_HalfClusterU->at(ih)->setHalfClusters(settings.map_stringPars["OutputAxisName"], tmp_axisCol);
continue;
}
@@ -207,7 +208,7 @@ cout<<endl;
for(int ic=0; ic<tmp_goodAxis.size(); ic++) tmp_axisCol.push_back(tmp_goodAxis[ic]);
// for(int ic=0; ic<m_newAxisUCol.size(); ic++) tmp_axisCol.push_back(m_newAxisUCol[ic]);
p_HalfClusterU->at(ih)->setHalfClusters(settings.map_stringPars["OutputAxisName"], tmp_axisCol);
-
+ p_HalfClusterU->at(ih)->sortBarShowersByLayer();
}
//Merge empty HalfCluster into the existing clusters.
@@ -221,6 +222,7 @@ cout<<endl;
}
}
}
+ std::sort( p_HalfClusterU->begin(), p_HalfClusterU->end(), compLayerPtr );
//cout << "yyy: Merge axis in V pHalfClusterV" << endl;
//Merge axis in HalfClusterV:
@@ -249,6 +251,7 @@ cout<<endl;
if(m_newAxisVCol.size()<2){
std::vector<const Cyber::CaloHalfCluster*> tmp_axisCol; tmp_axisCol.clear();
for(int ic=0; ic<m_newAxisVCol.size(); ic++) tmp_axisCol.push_back(m_newAxisVCol[ic]);
+ std::sort(tmp_axisCol.begin(), tmp_axisCol.end(), compLayer);
p_HalfClusterV->at(ih)->setHalfClusters(settings.map_stringPars["OutputAxisName"], tmp_axisCol);
continue;
}
@@ -315,9 +318,10 @@ cout<<endl;
//convert to constant object and save in HalfCluster:
std::vector<const Cyber::CaloHalfCluster*> tmp_axisCol; tmp_axisCol.clear();
for(int ic=0; ic<tmp_goodAxis.size(); ic++) tmp_axisCol.push_back(tmp_goodAxis[ic]);
+ std::sort(tmp_axisCol.begin(), tmp_axisCol.end(), compLayer);
// for(int ic=0; ic<m_newAxisVCol.size(); ic++) tmp_axisCol.push_back(m_newAxisVCol[ic]);
p_HalfClusterV->at(ih)->setHalfClusters(settings.map_stringPars["OutputAxisName"], tmp_axisCol);
-
+ p_HalfClusterV->at(ih)->sortBarShowersByLayer();
}
//Merge empty HalfClusters into the existing cluster.
@@ -330,6 +334,7 @@ cout<<endl;
m_datacol.map_HalfCluster["emptyHalfClusterV"].push_back(p_emptyHFClusV[ih]);
}
}
+ std::sort( p_HalfClusterV->begin(), p_HalfClusterV->end(), compLayerPtr );
return StatusCode::SUCCESS;
};
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/ConeClustering2DAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/ConeClustering2DAlg.cpp
index 38ea6252a1ed3d82071912c4054f3504c2a5da65..795bf7d42b7199a890b42593eafc821e54395d6b 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/ConeClustering2DAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/ConeClustering2DAlg.cpp
@@ -8,7 +8,7 @@ StatusCode ConeClustering2DAlg::ReadSettings(Settings& m_settings){
//Initialize parameters
if(settings.map_floatPars.find("th_beginLayer")==settings.map_floatPars.end()) settings.map_floatPars["th_beginLayer"] = 1;
- if(settings.map_floatPars.find("th_stopLayer")==settings.map_floatPars.end()) settings.map_floatPars["th_stopLayer"] = 15;
+ if(settings.map_floatPars.find("th_stopLayer")==settings.map_floatPars.end()) settings.map_floatPars["th_stopLayer"] = 10;
if(settings.map_floatPars.find("th_ConeTheta")==settings.map_floatPars.end()) settings.map_floatPars["th_ConeTheta"] = TMath::Pi()/4.;
if(settings.map_floatPars.find("th_ConeR")==settings.map_floatPars.end()) settings.map_floatPars["th_ConeR"] = 50;
if(settings.map_intPars.find("th_Nshowers")==settings.map_intPars.end()) settings.map_intPars["th_Nshowers"] = 4;
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/EnergySplittingAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/EnergySplittingAlg.cpp
index 3aced5c8e4441c999656103723a359e1736c971a..3e7248f5486f21f08c2733085aa4f39f276088c4 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/EnergySplittingAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/EnergySplittingAlg.cpp
@@ -227,6 +227,7 @@ cout<<"Empty half cluster energy "<<tmp_totE_U<<", "<<tmp_totE_V<<endl;
else m_axisVCol = p_HalfClusterV[ih]->getHalfClusterCol(settings.map_stringPars["ReadinAxisName"]);
m_1dShowerVCol.clear();
+ p_HalfClusterV[ih]->sortBarShowersByLayer();
std::vector<const Cyber::Calo1DCluster*> m_1dclusCol = p_HalfClusterV[ih]->getCluster();
//cout<<"1D Cluster size: "<<m_1dclusCol.size()<<", axis size: "<<m_axisVCol.size()<<endl;
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/GlobalClusteringAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/GlobalClusteringAlg.cpp
index 2eab8d9903cccebfd471d2aba728b8d126a0e308..ea4c2140cebccb54181eb377619287219bed8746 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/GlobalClusteringAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/GlobalClusteringAlg.cpp
@@ -38,16 +38,17 @@ StatusCode GlobalClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
//double totE_rest = 0.;
for(int ibar=0; ibar<m_bars.size(); ibar++)
{
- if(m_bars.at(ibar)->getEnergy()>settings.map_floatPars["unit_threshold"])
- {
- m_processbars.push_back(m_bars.at(ibar));
- //totE += m_bars.at(ibar)->getEnergy();
- }
- else
- {
- m_restbars.push_back(m_bars.at(ibar));
- //totE_rest += m_bars.at(ibar)->getEnergy();
- }
+ m_processbars.push_back(m_bars.at(ibar));
+ // if(m_bars.at(ibar)->getEnergy()>settings.map_floatPars["unit_threshold"])
+ // {
+ // m_processbars.push_back(m_bars.at(ibar));
+ // //totE += m_bars.at(ibar)->getEnergy();
+ // }
+ // else
+ // {
+ // m_restbars.push_back(m_bars.at(ibar));
+ // //totE_rest += m_bars.at(ibar)->getEnergy();
+ // }
}
//cout<<"Input bar after threshold: "<<m_processbars.size()<<", totE "<<totE<<", rest E "<<totE_rest<<endl;
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp
index 6bb6341851b33e377ce5c9ef1fd4da06e0f3e080..e7147853d4521b0c4cbf38486a26ce6cf511812b 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp
@@ -17,7 +17,10 @@ StatusCode HcalClusteringAlg::ReadSettings(Cyber::Settings& m_settings){
if(settings.map_floatPars.find("th_ConeR_l2")==settings.map_floatPars.end()) settings.map_floatPars["th_ConeR_l2"] = 120.;
if(settings.map_floatPars.find("th_ClusChi2")==settings.map_floatPars.end()) settings.map_floatPars["th_ClusChi2"] = 10e17;
if(settings.map_floatPars.find("fl_GoodClusLevel")==settings.map_floatPars.end()) settings.map_floatPars["fl_GoodClusLevel"] = 10;
-
+ if(settings.map_intPars.find("minNhit")==settings.map_intPars.end()) settings.map_intPars["minNhit"] = 3;
+ if(settings.map_intPars.find("minMergeLayer")==settings.map_intPars.end()) settings.map_intPars["minMergeLayer"] = 5;
+ if(settings.map_floatPars.find("maxMergeR")==settings.map_floatPars.end()) settings.map_floatPars["maxMergeR"] = 500;
+
return StatusCode::SUCCESS;
};
@@ -44,18 +47,53 @@ StatusCode HcalClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
m_hcalHits.push_back( m_datacol.map_CaloHit[settings.map_stringVecPars["InputHCALHits"][icl]][ih].get() );
}
+ //ordered hits by layer.
+ //TODO: this did not seperate barrel and endcap, need to do cone-clustering separately.
+ std::map<int, std::vector<Cyber::CaloHit*> > m_orderedHit;
+ m_orderedHit.clear();
+ for(int ih=0;ih<m_hcalHits.size();ih++)
+ m_orderedHit[m_hcalHits[ih]->getLayer()].push_back(m_hcalHits[ih]);
+
std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_clusterCol;
m_clusterCol.clear();
-// LongiConeLinking( m_orderedHit, m_clusterCol );
+ //LongiConeLinking( m_orderedHit, m_clusterCol );
Clustering(m_hcalHits, m_clusterCol);
+
+
// m_datacol.bk_Cluster3DCol.insert( m_datacol.bk_Cluster3DCol.end(), m_clusterCol.begin(), m_clusterCol.end() );
- // cout<<" Cluster size: "<<m_clusterCol.size()<<endl;
- // for(int ic=0; ic<m_clusterCol.size(); ic++)
- // {
- // cout<<" Cluster "<<ic<<":"<<m_clusterCol[ic]->getCaloHits().size()<<endl;
- // }
+ //Merge isolated hits.
+ //Do not merge hits at HCAL front 10 layers, they may need to be linked to ECAL clusters.
+ std::vector<const Cyber::CaloHit*> m_isoHits;
+ for(int ic=0; ic<m_clusterCol.size(); ic++){
+ if(m_clusterCol[ic]->getCaloHits().size()<=settings.map_intPars["minNhit"]){ //TODO: arbitrory criterion: cluster with <=3 hits are removed.
+ int minLayer = 999;
+ for(int ihit=0; ihit<m_clusterCol[ic]->getCaloHits().size(); ihit++){
+ if(m_clusterCol[ic]->getCaloHits()[ihit]->getLayer()<minLayer)
+ minLayer = m_clusterCol[ic]->getCaloHits()[ihit]->getLayer();
+ }
+
+ if(minLayer<=settings.map_intPars["minMergeLayer"]) continue;
+
+ for(int ihit=0; ihit<m_clusterCol[ic]->getCaloHits().size(); ihit++)
+ m_isoHits.push_back(m_clusterCol[ic]->getCaloHits()[ihit]);
+
+ m_clusterCol.erase(m_clusterCol.begin()+ic);
+ }
+ }
+
+ MergeToCluster(m_isoHits, m_clusterCol);
+
+ //cout<<" After HCAL clustering: Cluster size "<<m_clusterCol.size()<<endl;
+ //int Nisohits = 0;
+ //for(int ic=0; ic<m_clusterCol.size(); ic++)
+ //{
+ // cout<<" Cluster #"<<ic<<": hit size "<<m_clusterCol[ic]->getCaloHits().size()<<endl;
+ // if(m_clusterCol[ic]->getCaloHits().size()==1) Nisohits++;
+ //}
+ //cout<<" Isolated hit size: "<<Nisohits<<endl;
+
// m_datacol.map_CaloCluster[settings.map_stringPars["OutputCluster"]] = m_clusterCol;
m_datacol.map_CaloCluster[settings.map_stringPars["OutputHCALClusters"]]= m_clusterCol;
return StatusCode::SUCCESS;
@@ -174,4 +212,37 @@ template<typename T1, typename T2> StatusCode HcalClusteringAlg::Clustering(std:
//cout<<"how many neighbors: "<<number<<endl;
}
+StatusCode HcalClusteringAlg::MergeToCluster( std::vector<const Cyber::CaloHit*>& m_isohits, std::vector<std::shared_ptr<Cyber::Calo3DCluster>>& m_clusters ){
+
+ if(m_isohits.size()==0) return StatusCode::SUCCESS;
+
+ for(int ihit=0; ihit<m_isohits.size(); ihit++){
+ int index_closest = -1;
+ double minDistance = 1e6;
+ for(int icl=0; icl<m_clusters.size(); icl++){
+ for(auto jhit : m_clusters[icl]->getCaloHits()){
+ double tmp_R = sqrt( pow(jhit->getPosition().x()-m_isohits[ihit]->getPosition().x(),2) + pow(jhit->getPosition().y()-m_isohits[ihit]->getPosition().y(),2) + pow(jhit->getPosition().z()-m_isohits[ihit]->getPosition().z(),2) );
+ if( tmp_R>settings.map_floatPars["maxMergeR"] && tmp_R<minDistance)
+ { minDistance = tmp_R; index_closest = icl; }
+ }
+ }
+
+ if(index_closest>=0){
+ m_clusters[index_closest]->addHit( m_isohits[ihit] );
+ m_isohits.erase(m_isohits.begin()+ihit);
+ ihit--;
+ }
+ else{
+ std::shared_ptr<Cyber::Calo3DCluster> m_clus = std::make_shared<Cyber::Calo3DCluster>();
+ m_clus->addHit(m_isohits[ihit]);
+ m_clusters.push_back(m_clus);
+ m_isohits.erase(m_isohits.begin()+ihit);
+ ihit--;
+ }
+
+ }
+
+ return StatusCode::SUCCESS;
+}
+
#endif
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/HoughClusteringAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/HoughClusteringAlg.cpp
index 5c6aeb1ffd2502694d534aa42d429ede88cc06eb..e7de187a7ec74e77ece21b33ac81f11eb297f5d6 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/HoughClusteringAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/HoughClusteringAlg.cpp
@@ -62,7 +62,7 @@ StatusCode HoughClusteringAlg::ReadSettings(Cyber::Settings& m_settings){
// Algorithm parameter settings
if(settings.map_intPars.find("th_Layers")==settings.map_intPars.end())
- settings.map_intPars["th_Layers"] = 10;
+ settings.map_intPars["th_Layers"] = 7;
if(settings.map_intPars.find("th_peak")==settings.map_intPars.end())
settings.map_intPars["th_peak"] = 3;
if(settings.map_intPars.find("th_continueN")==settings.map_intPars.end())
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/LocalMaxFindingAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/LocalMaxFindingAlg.cpp
index 2adf2b98de56694701681156b2f981e160b6b0d7..502086633046420ea9de14acad6e82803b9dcda0 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/LocalMaxFindingAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/LocalMaxFindingAlg.cpp
@@ -7,7 +7,7 @@ using namespace Cyber;
StatusCode LocalMaxFindingAlg::ReadSettings(Cyber::Settings& m_settings){
settings = m_settings;
- if(settings.map_floatPars.find("Eth_localMax")==settings.map_floatPars.end()) settings.map_floatPars["Eth_localMax"] = 0.005;
+ if(settings.map_floatPars.find("Eth_localMax")==settings.map_floatPars.end()) settings.map_floatPars["Eth_localMax"] = 0.0075;
if(settings.map_floatPars.find("Eth_MaxWithNeigh")==settings.map_floatPars.end()) settings.map_floatPars["Eth_MaxWithNeigh"] = 0.;
if(settings.map_stringPars.find("OutputLocalMaxName")==settings.map_stringPars.end()) settings.map_stringPars["OutputLocalMaxName"] = "AllLocalMax";
return StatusCode::SUCCESS;
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/PFOCreatingAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/PFOCreatingAlg.cpp
index 8fe0dd3c105670f620e2156a432ea0c44818ea8b..6c07bd84831f818cead23e168d71e8d298f2d170 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/PFOCreatingAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/PFOCreatingAlg.cpp
@@ -64,6 +64,7 @@ StatusCode PFOCreatingAlg::RunAlgorithm( CyberDataCol& m_datacol ){
std::shared_ptr<Cyber::PFObject> tmp_pfo = std::make_shared<Cyber::PFObject>();
tmp_pfo->addTrack(ecal_cls_track[0]);
+ tmp_pfo->setPID(ecal_cls_track[0]->getPID());
tmp_pfo->addECALCluster(m_ecal_clusters[ie]);
for(int ic=0; ic<hcal_clus_candidate.size(); ic++){
tmp_pfo->addHCALCluster(hcal_clus_candidate[ic]);
@@ -92,6 +93,7 @@ StatusCode PFOCreatingAlg::RunAlgorithm( CyberDataCol& m_datacol ){
for(int it=0; it<m_tracks.size(); it++){
std::shared_ptr<Cyber::PFObject> tmp_pfo = std::make_shared<Cyber::PFObject>();
tmp_pfo->addTrack(m_tracks[it]);
+ tmp_pfo->setPID(m_tracks[it]->getPID());
m_pfobjects.push_back(tmp_pfo);
}
@@ -264,6 +266,7 @@ StatusCode PFOCreatingAlg::CreateLeftPFO(std::vector<Cyber::Track*>& _tracks,
for(int it=0; it<_tracks.size(); it++){
std::shared_ptr<Cyber::PFObject> tmp_pfo = std::make_shared<Cyber::PFObject>();
tmp_pfo->addTrack(_tracks[it]);
+ tmp_pfo->setPID(_tracks[it]->getPID());
_pfobjects.push_back(tmp_pfo);
}
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp
index 73e77e32dbe5ddc4a24a5d950d100cd32aa64a93..adc6ae4a360759982e509403f234e033668d85c9 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp
@@ -9,16 +9,18 @@ StatusCode PFOReclusteringAlg::ReadSettings(Settings& m_settings){
//Initialize parameters
if(settings.map_stringPars.find("ReadinPFOName")==settings.map_stringPars.end()) settings.map_stringPars["ReadinPFOName"] = "outputPFO";
if(settings.map_floatPars.find("ECALChargedCalib")==settings.map_floatPars.end()) settings.map_floatPars["ECALChargedCalib"] = 1.26;
+ if(settings.map_floatPars.find("ECALChargedEMCalib")==settings.map_floatPars.end()) settings.map_floatPars["ECALChargedEMCalib"] = 1.0;
if(settings.map_floatPars.find("HCALChargedCalib")==settings.map_floatPars.end()) settings.map_floatPars["HCALChargedCalib"] = 4.0;
if(settings.map_floatPars.find("ECALNeutralCalib")==settings.map_floatPars.end()) settings.map_floatPars["ECALNeutralCalib"] = 1.0;
if(settings.map_floatPars.find("HCALNeutralCalib")==settings.map_floatPars.end()) settings.map_floatPars["HCALNeutralCalib"] = 4.0;
+ if(settings.map_floatPars.find("ElecEoverP")==settings.map_floatPars.end()) settings.map_floatPars["ElecEoverP"] = 0.8;
if(settings.map_floatPars.find("EnergyRes")==settings.map_floatPars.end()) settings.map_floatPars["EnergyRes"] = 0.4;
if(settings.map_floatPars.find("SplitSigma")==settings.map_floatPars.end()) settings.map_floatPars["SplitSigma"] = 0.2;
if(settings.map_floatPars.find("NeutralMergeSigma")==settings.map_floatPars.end()) settings.map_floatPars["NeutralMergeSigma"] = 0.;
if(settings.map_floatPars.find("VirtualMergeSigma")==settings.map_floatPars.end()) settings.map_floatPars["VirtualMergeSigma"] = 0.4;
- if(settings.map_floatPars.find("MinAngleForNeuMerge")==settings.map_floatPars.end()) settings.map_floatPars["MinAngleForNeuMerge"] = 0.18;
- if(settings.map_floatPars.find("MinAngleForVirMerge")==settings.map_floatPars.end()) settings.map_floatPars["MinAngleForVirMerge"] = 0.20;
+ if(settings.map_floatPars.find("MinAngleForNeuMerge")==settings.map_floatPars.end()) settings.map_floatPars["MinAngleForNeuMerge"] = 0.20;
+ if(settings.map_floatPars.find("MinAngleForVirMerge")==settings.map_floatPars.end()) settings.map_floatPars["MinAngleForVirMerge"] = 0.24;
return StatusCode::SUCCESS;
@@ -477,8 +479,17 @@ StatusCode PFOReclusteringAlg::ReCluster_SplitFromChg( std::vector< std::shared_
if(m_chargedPFOs[ipfo]->getECALClusters().size()==0 && m_chargedPFOs[ipfo]->getHCALClusters().size()==0) continue;
double track_energy = m_chargedPFOs[ipfo]->getTrackMomentum();
- double ECAL_energy = settings.map_floatPars["ECALChargedCalib"]*m_chargedPFOs[ipfo]->getECALClusterEnergy();
+ double ECAL_energy = settings.map_floatPars["ECALChargedEMCalib"]*m_chargedPFOs[ipfo]->getECALClusterEnergy();
double HCAL_energy = settings.map_floatPars["HCALChargedCalib"]*m_chargedPFOs[ipfo]->getHCALClusterEnergy();
+
+//cout<<" Charged PFO #"<<ipfo<<": track P "<<track_energy<<", cluster E "<<ECAL_energy+HCAL_energy<<", PID "<<m_chargedPFOs[ipfo]->getPID()<<endl;
+
+ //FIXME: If PFO is identified as electron (E/p > 0.8): use ECALChargedEMCalib(1.0). Else: use ECALChargedCalib(1.26).
+ if( track_energy < 15. || ECAL_energy/track_energy < settings.map_floatPars["ElecEoverP"] )
+ //if( (track_energy<15. && abs(m_chargedPFOs[ipfo]->getPID())!=11) ||
+ // (track_energy>15. && ECAL_energy/track_energy < settings.map_floatPars["ElecEoverP"]) )
+ ECAL_energy = ECAL_energy*settings.map_floatPars["ECALChargedCalib"]/settings.map_floatPars["ECALChargedEMCalib"];
+
if(track_energy<0 || ECAL_energy<0 || HCAL_energy<0){
std::cout<<"ERROR: Charged PFO info break. Ptrk "<<track_energy<<", E_ecal "<<ECAL_energy<<", E_hcal "<<HCAL_energy<<endl;
continue;
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/TrackClusterConnectingAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/TrackClusterConnectingAlg.cpp
index 853fda17f6a3aba5756aa78c0f928775dd0f56ee..58c2ff6b1efdd409d7888b877329edd4f442d8fe 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/TrackClusterConnectingAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/TrackClusterConnectingAlg.cpp
@@ -18,9 +18,13 @@ StatusCode TrackClusterConnectingAlg::ReadSettings(Settings& m_settings){
if(settings.map_floatPars.find("th_ChFragDepth")==settings.map_floatPars.end()) settings.map_floatPars["th_ChFragDepth"] = 100.;
if(settings.map_floatPars.find("th_ChFragMinR")==settings.map_floatPars.end()) settings.map_floatPars["th_ChFragMinR"] = 200.;
if(settings.map_floatPars.find("th_HcalMatchingR")==settings.map_floatPars.end()) settings.map_floatPars["th_HcalMatchingR"] = 100.;
+ if(settings.map_floatPars.find("th_NeutralHcalMergeTheta")==settings.map_floatPars.end()) settings.map_floatPars["th_NeutralHcalMergeTheta"] = 0.3;
+ if(settings.map_floatPars.find("th_NeutralHcalMergeR")==settings.map_floatPars.end()) settings.map_floatPars["th_NeutralHcalMergeR"] = 300.;
if(settings.map_floatPars.find("th_MIPEnergy")==settings.map_floatPars.end()) settings.map_floatPars["th_MIPEnergy"] = 0.5;
if(settings.map_floatPars.find("th_AbsorbCone")==settings.map_floatPars.end()) settings.map_floatPars["th_AbsorbCone"] = 0.8;
+ if(settings.map_intPars.find("minNhit")==settings.map_intPars.end()) settings.map_intPars["minNhit"] = 5;
+
if(settings.map_stringPars.find("OutputMergedECALCluster")==settings.map_stringPars.end()) settings.map_stringPars["OutputMergedECALCluster"] = "TrkMergedECAL";
if(settings.map_stringPars.find("OutputCombPFO")==settings.map_stringPars.end()) settings.map_stringPars["OutputCombPFO"] = "outputPFO";
@@ -95,11 +99,12 @@ StatusCode TrackClusterConnectingAlg::RunAlgorithm( CyberDataCol& m_datacol ){
//cout<<"Print all HCAL cluster"<<endl;
//double totE_Hcal = 0;
//for(int i=0; i<m_HcalClusters.size(); i++){
-// printf(" HCAL Cluster #%d: En = %.6f, position (%.3f, %.3f, %.3f) \n", i,
+// printf(" HCAL Cluster #%d: En = %.6f, position (%.3f, %.3f, %.3f), address %p \n", i,
// m_HcalClusters[i]->getHitsE(),
// m_HcalClusters[i]->getHitCenter().x(),
// m_HcalClusters[i]->getHitCenter().y(),
-// m_HcalClusters[i]->getHitCenter().z() );
+// m_HcalClusters[i]->getHitCenter().z(),
+// m_HcalClusters[i] );
// totE_Hcal += m_HcalClusters[i]->getHitsE();
//}
//cout<<"Hcal cluster total E "<<totE_Hcal<<endl;
@@ -108,12 +113,15 @@ StatusCode TrackClusterConnectingAlg::RunAlgorithm( CyberDataCol& m_datacol ){
std::sort(m_PFObjects.begin(), m_PFObjects.end(), compTrkP);
HcalExtrapolatingMatch(m_HcalClusters, m_PFObjects);
-//cout<<" TrackClusterConnectingAlg: PFO size after HCAL matching: "<<m_PFObjects.size()<<endl;
-//for(int i=0; i<m_PFObjects.size(); i++){
-// cout<<" PFO #"<<i<<": track size "<<m_PFObjects[i]->getTracks().size()<<", leading P "<<m_PFObjects[i]->getTrackMomentum();
-// cout<<", ECAL cluster size "<<m_PFObjects[i]->getECALClusters().size()<<", totE "<<m_PFObjects[i]->getECALClusterEnergy();
-// cout<<", HCAL cluster size "<<m_PFObjects[i]->getHCALClusters().size()<<", totE "<<m_PFObjects[i]->getHCALClusterEnergy()<<endl;
-//}
+ //cout<<" TrackClusterConnectingAlg: PFO size after HCAL matching: "<<m_PFObjects.size()<<endl;
+ //for(int i=0; i<m_PFObjects.size(); i++){
+ // cout<<" PFO #"<<i<<": track size "<<m_PFObjects[i]->getTracks().size()<<", leading P "<<m_PFObjects[i]->getTrackMomentum();
+ // cout<<", ECAL cluster size "<<m_PFObjects[i]->getECALClusters().size()<<", totE "<<m_PFObjects[i]->getECALClusterEnergy();
+ // cout<<", HCAL cluster size "<<m_PFObjects[i]->getHCALClusters().size()<<", totE "<<m_PFObjects[i]->getHCALClusterEnergy()<<endl;
+ // for(int j=0; j<m_PFObjects[i]->getHCALClusters().size(); j++){
+ // printf(" HCAL #%d: energy %.6f, Nhit %d, address %p \n", j, m_PFObjects[i]->getHCALClusters()[j]->getHitsE(), m_PFObjects[i]->getHCALClusters()[j]->getCaloHits().size(), m_PFObjects[i]->getHCALClusters()[j]);
+ // }
+ //}
m_datacol.map_CaloCluster[ settings.map_stringPars["OutputMergedECALCluster"] ] = m_absorbedEcal;
m_datacol.map_PFObjects[settings.map_stringPars["OutputCombPFO"]] = m_PFObjects;
@@ -158,6 +166,7 @@ StatusCode TrackClusterConnectingAlg::PFOCreating( std::vector<const Cyber::Calo
std::vector<const Cyber::Track*> m_trkInClus = m_clusters[ic]->getAssociatedTracks();
if(m_trkInClus.size()!=0){
m_newPFO->addTrack( m_trkInClus[0] );
+ m_newPFO->setPID( m_trkInClus[0]->getPID() );
auto iter = find(m_leftTrks.begin(), m_leftTrks.end(), m_trkInClus[0]);
if( iter!=m_leftTrks.end() )
m_leftTrks.erase(iter);
@@ -170,6 +179,7 @@ StatusCode TrackClusterConnectingAlg::PFOCreating( std::vector<const Cyber::Calo
for(int itrk=0; itrk<m_leftTrks.size(); itrk++){
std::shared_ptr<Cyber::PFObject> m_newPFO = std::make_shared<Cyber::PFObject>();
m_newPFO->addTrack( m_leftTrks[itrk] );
+ m_newPFO->setPID( m_leftTrks[itrk]->getPID() );
m_PFOs.push_back(m_newPFO);
m_bkCol.map_PFObjects["bkPFO"].push_back(m_newPFO);
}
@@ -334,18 +344,38 @@ StatusCode TrackClusterConnectingAlg::EcalChFragAbsorption( std::vector<const Cy
StatusCode TrackClusterConnectingAlg::HcalExtrapolatingMatch(std::vector<const Cyber::Calo3DCluster*>& m_clusters, std::vector<std::shared_ptr<Cyber::PFObject>>& m_PFOs){
+ if(m_clusters.size()<=0){
+ std::cout<<"Error in TrackClusterConnectingAlg: No HCAL cluster. Skip this match "<<std::endl;
+ return StatusCode::SUCCESS;
+ }
+//cout<<" HcalExtrapolatingMatch: Readin cluster size "<<m_clusters.size()<<endl;
+//cout<<" Readin PFO: "<<m_PFOs.size()<<endl;
+//for(int i=0; i<m_PFOs.size(); i++){
+// cout<<" In PFO #"<<i<<": Ntrk "<<m_PFOs[i]->getTracks().size()<<", N Ecal "<<m_PFOs[i]->getECALClusters().size()<<", N Hcal "<<m_PFOs[i]->getHCALClusters().size()<<endl;
+//}
+
+
+/* std::vector<const Cyber::CaloHit*> m_isoHits;
+ for(int ic=0; ic<m_clusters.size(); ic++){
+ if(m_clusters[ic]->getCaloHits().size()<=3){
+ for(int ihit=0; ihit<m_clusters[ic]->getCaloHits().size(); ihit++)
+ m_isoHits.push_back(m_clusters[ic]->getCaloHits()[ihit]);
+
+ m_clusters.erase(m_clusters.begin()+ic);
+ }
+ }
+cout<<", real cluster size "<<m_clusters.size()<<", isolated hit size "<<m_isoHits.size()<<endl;
+*/
for(int ic=0; ic<m_clusters.size(); ic++){
std::vector<const Cyber::CaloHit*> hcal_hits = m_clusters[ic]->getCaloHits();
//cout<<"HCAL Cluster #"<<ic<<": Nhit "<<hcal_hits.size()<<", En "<<m_clusters[ic]->getHitsE()<<endl;
bool isInPfo = false;
- int index_selPfo = -1;
for(int ipfo=0; ipfo<m_PFOs.size(); ipfo++){
//Link HCAL cluster to charged PFO
if(m_PFOs[ipfo]->getTracks().size()!=0){
std::vector<TrackState> trk_points = m_PFOs[ipfo]->getTracks()[0]->getAllTrackStates();
-
bool is_candidate = false;
double minDistance = 99999;
for(int ihit=0; ihit<hcal_hits.size(); ihit++){
@@ -368,25 +398,56 @@ StatusCode TrackClusterConnectingAlg::HcalExtrapolatingMatch(std::vector<const C
//cout<<" Link cluster #"<<ic<<" to pfo #"<<ipfo<<endl;
m_PFOs[ipfo]->addHCALCluster( m_clusters[ic] );
isInPfo = true;
- index_selPfo = ipfo;
break;
}
-
}
//Link HCAL cluster to neutral PFO
+ else{
+ TVector3 pfo_pos(0., 0., 0.);
+ for(int jcl=0; jcl<m_PFOs[ipfo]->getECALClusters().size(); jcl++) pfo_pos += m_PFOs[ipfo]->getECALClusters()[jcl]->getShowerCenter() * m_PFOs[ipfo]->getECALClusters()[jcl]->getLongiE();
+ for(int jcl=0; jcl<m_PFOs[ipfo]->getHCALClusters().size(); jcl++) pfo_pos += m_PFOs[ipfo]->getHCALClusters()[jcl]->getHitCenter() * m_PFOs[ipfo]->getHCALClusters()[jcl]->getHitsE();
+ pfo_pos = pfo_pos*(1./(m_PFOs[ipfo]->getECALClusterEnergy() + m_PFOs[ipfo]->getHCALClusterEnergy()) );
+
+ TVector3 clus_pos = m_clusters[ic]->getHitCenter();
+ if( pfo_pos.Angle(clus_pos)<settings.map_floatPars["th_NeutralHcalMergeTheta"] && (clus_pos-pfo_pos).Mag()<settings.map_floatPars["th_NeutralHcalMergeR"] ){
+ m_PFOs[ipfo]->addHCALCluster( m_clusters[ic] );
+ isInPfo = true;
+ break;
+ }
+ }
}//end loop pfos
//if(isInPfo) cout<<" Merged into PFO: Ptrk = "<<m_PFOs[index_selPfo]->getTrackMomentum()<<endl;
//If HCAL cluster is not linked to any existing PFO: create a new one.
- if(!isInPfo){
+ if(!isInPfo && m_clusters[ic]->getCaloHits().size()>=settings.map_intPars["minNhit"]){
//cout<<" Create a new neutral PFO "<<endl;
std::shared_ptr<Cyber::PFObject> m_newPFO = std::make_shared<Cyber::PFObject>();
m_newPFO->addHCALCluster( m_clusters[ic] );
m_PFOs.push_back(m_newPFO);
m_bkCol.map_PFObjects["bkPFO"].push_back(m_newPFO);
}
+ //If only have few hits: merge to closest PFO.
+ else if(!isInPfo){
+
+ int index_closest = -1;
+ double minDistance = 1e6;
+ for(int ipfo=0; ipfo<m_PFOs.size(); ipfo++){
+ TVector3 pfo_pos(0., 0., 0.);
+ for(int jcl=0; jcl<m_PFOs[ipfo]->getECALClusters().size(); jcl++) pfo_pos += m_PFOs[ipfo]->getECALClusters()[jcl]->getShowerCenter() * m_PFOs[ipfo]->getECALClusters()[jcl]->getLongiE();
+ for(int jcl=0; jcl<m_PFOs[ipfo]->getHCALClusters().size(); jcl++) pfo_pos += m_PFOs[ipfo]->getHCALClusters()[jcl]->getHitCenter() * m_PFOs[ipfo]->getHCALClusters()[jcl]->getHitsE();
+ pfo_pos = pfo_pos*(1./(m_PFOs[ipfo]->getECALClusterEnergy() + m_PFOs[ipfo]->getHCALClusterEnergy()) );
+
+ double tmp_dis = (pfo_pos-m_clusters[ic]->getHitCenter()).Mag();
+ if(tmp_dis<minDistance){
+ minDistance = tmp_dis;
+ index_closest = ipfo;
+ }
+ }
+
+ if(index_closest>=0) m_PFOs[index_closest]->addHCALCluster( m_clusters[ic] );
+ }
}
return StatusCode::SUCCESS;
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/TrackExtrapolatingAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/TrackExtrapolatingAlg.cpp
index 2f5639a730950f5ae292ab00a64cadc1720903fe..c15e73e2a366649bada7702c93d8946e263d0eec 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/TrackExtrapolatingAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/TrackExtrapolatingAlg.cpp
@@ -22,7 +22,7 @@ StatusCode TrackExtrapolatingAlg::ReadSettings(Settings& m_settings){
if(settings.map_intPars.find("ECAL_Nlayers")==settings.map_intPars.end())
settings.map_intPars["ECAL_Nlayers"] = 28;
if(settings.map_floatPars.find("ECAL_layer_width")==settings.map_floatPars.end())
- settings.map_floatPars["ECAL_layer_width"] = 10;
+ settings.map_floatPars["ECAL_layer_width"] = Cyber::CaloUnit::barsize;
if(settings.map_floatPars.find("ECAL_half_length")==settings.map_floatPars.end())
settings.map_floatPars["ECAL_half_length"] = 2900;
if(settings.map_floatPars.find("ECAL_endcap_zmin")==settings.map_floatPars.end())
@@ -49,7 +49,7 @@ StatusCode TrackExtrapolatingAlg::ReadSettings(Settings& m_settings){
settings.map_floatPars["HCAL_endcap_zmax"] = 4575; // 3260 + (3455-2140)
if(settings.map_intPars.find("Nmodule")==settings.map_intPars.end())
- settings.map_intPars["Nmodule"] = 32;
+ settings.map_intPars["Nmodule"] = Cyber::CaloUnit::Nmodule;
if(settings.map_floatPars.find("B_field")==settings.map_floatPars.end())
settings.map_floatPars["B_field"] = 3.0;
@@ -339,7 +339,13 @@ vector<float> & barrel_delta_phi, vector<float> & endcap_delta_phi){
// Endcap
float z0 = CALO_trk_state.referencePoint.Z() + CALO_trk_state.Z0;
for(int iz=0; iz<layer_z.size(); iz++){
- float layer_delta_phi = ( layer_z[iz] - z0 ) / ( rho * CALO_trk_state.tanLambda );
+ float layer_delta_phi;
+ if (CALO_trk_state.tanLambda>0){
+ layer_delta_phi = ( layer_z[iz] - z0 ) / ( rho * CALO_trk_state.tanLambda );
+ }
+ else{
+ layer_delta_phi = ( -layer_z[iz] - z0 ) / ( rho * CALO_trk_state.tanLambda );
+ }
if(CALO_trk_state.Kappa < 0){
endcap_delta_phi.push_back(layer_delta_phi);
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/TruthClusterMergingAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/TruthClusterMergingAlg.cpp
index 1e1471892df6ca3382572e65b1777282d1e56aba..ae884ba2c967677ec85b7677a907344902ed3898 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/TruthClusterMergingAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/TruthClusterMergingAlg.cpp
@@ -153,8 +153,14 @@ StatusCode TruthClusterMergingAlg::RunAlgorithm( CyberDataCol& m_datacol ){
const Cyber::Track* p_HcalLeadingTrk = nullptr;
if(index>=0 && m_linkedHcalClus[index]->getAssociatedTracks().size()>0) p_HcalLeadingTrk = m_linkedHcalClus[index]->getAssociatedTracks()[0];
- if(!p_EcalLeadingTrk && p_HcalLeadingTrk) tmp_newpfo->addTrack(p_HcalLeadingTrk);
- else if(p_EcalLeadingTrk) tmp_newpfo->addTrack(p_EcalLeadingTrk);
+ if(!p_EcalLeadingTrk && p_HcalLeadingTrk){
+ tmp_newpfo->addTrack(p_HcalLeadingTrk);
+ tmp_newpfo->setPID(p_HcalLeadingTrk->getPID());
+ }
+ else if(p_EcalLeadingTrk){
+ tmp_newpfo->addTrack(p_EcalLeadingTrk);
+ tmp_newpfo->setPID(p_EcalLeadingTrk->getPID());
+ }
merged_CombClusterCol.push_back(tmp_newpfo);
m_datacol.map_PFObjects["bkPFO"].push_back(tmp_newpfo);
diff --git a/Reconstruction/RecPFACyber/src/CyberDataCol.cpp b/Reconstruction/RecPFACyber/src/CyberDataCol.cpp
index 98106f7ad23971b58918d4ec2b2efc43ac2e415c..f8386eab63e3b44fe11d4e9d354a064795cdaac3 100644
--- a/Reconstruction/RecPFACyber/src/CyberDataCol.cpp
+++ b/Reconstruction/RecPFACyber/src/CyberDataCol.cpp
@@ -31,7 +31,10 @@ StatusCode CyberDataCol::Clear(){
map_2DCluster.clear();
map_CaloCluster.clear();
map_PFObjects.clear();
-
+
+ //if(tofCol) delete tofCol;
+ //if(dNdxCol) delete dNdxCol;
+
return StatusCode::SUCCESS;
};
diff --git a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
index 354cb1289eb28f165e981051ff0b120e959ea522..9bdcc842c04f3119615139cfb5fc205d10483e1b 100644
--- a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
@@ -23,16 +23,16 @@ int Cyber::CaloUnit::System_Barrel = 20;
int Cyber::CaloUnit::System_Endcap = 29;
int Cyber::CaloUnit::Nmodule = 32;
int Cyber::CaloUnit::Nstave = 15;
-int Cyber::CaloUnit::Nlayer = 14;
-int Cyber::CaloUnit::NbarPhi_odd[14] = {39, 37, 37, 37, 37, 37, 35, 35, 35, 35, 33, 33, 33, 33};
-int Cyber::CaloUnit::NbarPhi_even[14] = {27, 29, 29, 31, 31, 33, 35, 35, 37, 37, 39, 41, 41, 43};
-int Cyber::CaloUnit::NbarZ = 36;
+int Cyber::CaloUnit::Nlayer = 9;
+int Cyber::CaloUnit::NbarPhi_odd[9] = {25, 25, 25, 23, 23, 23, 23, 23, 21};
+int Cyber::CaloUnit::NbarPhi_even[9] = {19, 19, 21, 23, 23, 25, 27, 27, 29};
+int Cyber::CaloUnit::NbarZ = 24;
//int Cyber::CaloUnit::over_module[28] = {13,15,16,18,19,21,22,24,25,26,28,29,30,32,33,35,36,38,39,41,42,43,45,46};
//int Cyber::CaloUnit::over_module_set = 2;
-float Cyber::CaloUnit::barsize = 10.; //mm
+float Cyber::CaloUnit::barsize = 15.2; //mm
float Cyber::CaloUnit::ecal_innerR = 1830; //mm
float Cyber::CaloUnit::ecal_endcap_deadarea = 8.5; //mm
-float Cyber::CaloUnit::ecal_endcap_barsize = 10.; //mm
+float Cyber::CaloUnit::ecal_endcap_barsize = 15.2; //mm
DECLARE_COMPONENT( CyberPFAlg )
@@ -110,8 +110,11 @@ StatusCode CyberPFAlg::initialize()
//---MC particle---
if(!name_MCParticleCol.empty()) r_MCParticleCol = new DataHandle<edm4hep::MCParticleCollection> (name_MCParticleCol, Gaudi::DataHandle::Reader, this);
- //---Tracks---
+ //---Tracks, dN/dx and TOF---
for(auto& _trk : name_TrackCol) if(!_trk.empty()) r_TrackCols.push_back( new TrackType(_trk, Gaudi::DataHandle::Reader, this) );
+ if(!name_dNdxCol.empty()) r_dNdxCol = new DataHandle<edm4hep::RecDqdxCollection> (name_dNdxCol, Gaudi::DataHandle::Reader, this);
+ if(!name_tofCol.empty()) r_TofCol = new DataHandle<edm4hep::RecTofCollection>(name_tofCol, Gaudi::DataHandle::Reader, this);
+
//---Calo Hits---
for(auto& _ecal : name_EcalHits){
@@ -615,6 +618,7 @@ StatusCode CyberPFAlg::initialize()
t_Cluster->Branch("HcalClus_y", &m_HcalClus_y);
t_Cluster->Branch("HcalClus_z", &m_HcalClus_z);
t_Cluster->Branch("HcalClus_E", &m_HcalClus_E);
+ t_Cluster->Branch("HcalClus_nHit", &m_HcalClus_nHit);
t_Cluster->Branch("HcalClus_nTrk", &m_HcalClus_nTrk);
t_Cluster->Branch("HcalClus_ptrk", &m_HcalClus_pTrk);
t_Cluster->Branch("HcalClus_hit_x", &m_HcalClus_hit_x);
@@ -636,6 +640,7 @@ StatusCode CyberPFAlg::initialize()
t_Cluster->Branch("SimpleHcalClus_y", &m_SimpleHcalClus_y);
t_Cluster->Branch("SimpleHcalClus_z", &m_SimpleHcalClus_z);
t_Cluster->Branch("SimpleHcalClus_E", &m_SimpleHcalClus_E);
+ t_Cluster->Branch("SimpleHcalClus_nHit", &m_SimpleHcalClus_nHit);
t_Cluster->Branch("SimpleHcalClus_nTrk", &m_SimpleHcalClus_nTrk);
t_Cluster->Branch("SimpleHcalClus_ptrk", &m_SimpleHcalClus_pTrk);
t_Cluster->Branch("SimpleHcalClus_hit_x", &m_SimpleHcalClus_hit_x);
@@ -658,6 +663,13 @@ StatusCode CyberPFAlg::initialize()
t_Track->Branch("m_Ntrk", &m_Ntrk);
t_Track->Branch("m_type", &m_type);
t_Track->Branch("m_Nhit", &m_Nhit);
+ t_Track->Branch("m_pid", &m_pid);
+ t_Track->Branch("m_pid_truth", &m_pid_truth);
+ t_Track->Branch("m_trk_px", &m_trk_px);
+ t_Track->Branch("m_trk_py", &m_trk_py);
+ t_Track->Branch("m_trk_pz", &m_trk_pz);
+ t_Track->Branch("m_trk_p", &m_trk_p);
+ t_Track->Branch("m_trk_truthweight", &m_trk_truthweight);
t_Track->Branch("m_trkstate_d0", &m_trkstate_d0);
t_Track->Branch("m_trkstate_z0", &m_trkstate_z0);
t_Track->Branch("m_trkstate_phi", &m_trkstate_phi);
@@ -722,6 +734,8 @@ StatusCode CyberPFAlg::execute()
CyberDataCol m_DataCol;
m_DataCol.Clear();
m_DataCol.EnergyCorrSvc = m_energycorsvc;
+ m_DataCol.tofCol = const_cast<edm4hep::RecTofCollection*>(r_TofCol->get());
+ m_DataCol.dNdxCol = const_cast<edm4hep::RecDqdxCollection*>(r_dNdxCol->get());
//Readin collections
@@ -1504,6 +1518,7 @@ StatusCode CyberPFAlg::execute()
m_HcalClus_pTrk.push_back(m_HcalClusterCol[icl]->getAssociatedTracks()[0]->getMomentum());
else
m_HcalClus_pTrk.push_back(-99);
+ m_HcalClus_nHit.push_back(m_HcalClusterCol[icl]->getCaloHits().size());
for(int ih=0; ih<m_HcalClusterCol[icl]->getCaloHits().size(); ih++){
m_HcalClus_hit_tag.push_back(icl);
@@ -1539,6 +1554,7 @@ StatusCode CyberPFAlg::execute()
m_SimpleHcalClus_pTrk.push_back(m_SimpleHcalClusterCol[icl]->getAssociatedTracks()[0]->getMomentum());
else
m_SimpleHcalClus_pTrk.push_back(-99);
+ m_SimpleHcalClus_nHit.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits().size());
for(int ih=0; ih<m_SimpleHcalClusterCol[icl]->getCaloHits().size(); ih++){
m_SimpleHcalClus_hit_tag.push_back(icl);
@@ -1574,6 +1590,13 @@ StatusCode CyberPFAlg::execute()
for(int itrk=0; itrk<m_Ntrk; itrk++){
m_type.push_back(m_trkCol[itrk]->getType());
m_Nhit.push_back(m_trkCol[itrk]->getTrackerHits());
+ m_pid.push_back(m_trkCol[itrk]->getPID());
+ m_pid_truth.push_back(m_trkCol[itrk]->getLeadingMCP().getPDG());
+ m_trk_truthweight.push_back( m_trkCol[itrk]->getLeadingMCPweight() );
+ m_trk_px.push_back( m_trkCol[itrk]->getP3().x() );
+ m_trk_py.push_back( m_trkCol[itrk]->getP3().y() );
+ m_trk_pz.push_back( m_trkCol[itrk]->getP3().z() );
+ m_trk_p.push_back( m_trkCol[itrk]->getMomentum() );
std::vector<TrackState> AllTrackStates = m_trkCol[itrk]->getAllTrackStates();
for(int istate=0; istate<AllTrackStates.size(); istate++){
m_trkstate_d0.push_back( AllTrackStates[istate].D0 );
@@ -2112,6 +2135,7 @@ void CyberPFAlg::ClearCluster(){
m_HcalClus_y.clear();
m_HcalClus_z.clear();
m_HcalClus_E.clear();
+ m_HcalClus_nHit.clear();
m_HcalClus_nTrk.clear();
m_HcalClus_pTrk.clear();
m_HcalClus_hit_x.clear();
@@ -2133,6 +2157,7 @@ void CyberPFAlg::ClearCluster(){
m_SimpleHcalClus_y.clear();
m_SimpleHcalClus_z.clear();
m_SimpleHcalClus_E.clear();
+ m_SimpleHcalClus_nHit.clear();
m_SimpleHcalClus_nTrk.clear();
m_SimpleHcalClus_pTrk.clear();
m_SimpleHcalClus_hit_x.clear();
@@ -2155,6 +2180,13 @@ void CyberPFAlg::ClearCluster(){
void CyberPFAlg::ClearTrack(){
m_type.clear();
m_Nhit.clear();
+ m_pid.clear();
+ m_pid_truth.clear();
+ m_trk_px.clear();
+ m_trk_py.clear();
+ m_trk_pz.clear();
+ m_trk_p.clear();
+ m_trk_truthweight.clear();
m_trkstate_d0.clear();
m_trkstate_z0.clear();
m_trkstate_phi.clear();
diff --git a/Reconstruction/RecPFACyber/src/Objects/Calo1DCluster.cc b/Reconstruction/RecPFACyber/src/Objects/Calo1DCluster.cc
index 46705b93cd3ad0744e5b25cbaa16fc67394ad329..708ab079eaa63bf3bc5f7b5c0f7346e2cfcf02af 100644
--- a/Reconstruction/RecPFACyber/src/Objects/Calo1DCluster.cc
+++ b/Reconstruction/RecPFACyber/src/Objects/Calo1DCluster.cc
@@ -174,8 +174,8 @@ namespace Cyber{
if( Bars[Bars.size()-1]->getSystem() == CaloUnit::System_Barrel ){
if( Bars[Bars.size()-1]->getSlayer()==0 ) edge = Bars[Bars.size()-1]->getBar() + Bars[Bars.size()-1]->getStave()*CaloUnit::NbarZ;
if( Bars[Bars.size()-1]->getSlayer()==1 ){
- if(Bars[Bars.size()-1]->getModule()%2==0) edge = Bars[Bars.size()-1]->getBar() + Bars[Bars.size()-1]->getModule()*(CaloUnit::NbarPhi_even[Bars[Bars.size()-1]->getDlayer()]);
- else edge = Bars[Bars.size()-1]->getBar() + Bars[Bars.size()-1]->getModule()*(CaloUnit::NbarPhi_odd[Bars[Bars.size()-1]->getDlayer()]);
+ if(Bars[Bars.size()-1]->getModule()%2==0) edge = Bars[Bars.size()-1]->getBar() + Bars[Bars.size()-1]->getModule()*(CaloUnit::NbarPhi_even[Bars[Bars.size()-1]->getDlayer()-1]);
+ else edge = Bars[Bars.size()-1]->getBar() + Bars[Bars.size()-1]->getModule()*(CaloUnit::NbarPhi_odd[Bars[Bars.size()-1]->getDlayer()-1]);
}
}
if( Bars[Bars.size()-1]->getSystem() == CaloUnit::System_Endcap ){
diff --git a/Reconstruction/RecPFACyber/src/Objects/CaloUnit.cc b/Reconstruction/RecPFACyber/src/Objects/CaloUnit.cc
index 41934d95c7e32977668b01618383bb7066312315..aedb0e6780e73be05839ac85cfd9e5545e976954 100644
--- a/Reconstruction/RecPFACyber/src/Objects/CaloUnit.cc
+++ b/Reconstruction/RecPFACyber/src/Objects/CaloUnit.cc
@@ -102,10 +102,10 @@ namespace Cyber{
}
else{
if(module%2==0){
- if(slayer==1 && bar==NbarPhi_even[dlayer-1]) isEdge = true;
+ if(slayer==1 && bar==NbarPhi_even[dlayer-1]-1) isEdge = true;
}
else{
- if(slayer==1 && bar==NbarPhi_odd[dlayer-1]) isEdge = true;
+ if(slayer==1 && bar==NbarPhi_odd[dlayer-1]-1) isEdge = true;
}
}
return isEdge;
diff --git a/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp b/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
index 060086543a038b0d1ca166cb8de564b4d93920ae..f105c6c74aa0f88dc28a84345a293167fde29827 100644
--- a/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
+++ b/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
@@ -71,6 +71,38 @@ namespace Cyber{
}
+ //Assign PID from dNdx and TOF
+ if( m_DataCol.dNdxCol && m_DataCol.tofCol ){
+ std::array<double, 5> chi2s; chi2s.fill(0);
+
+ for (auto dqdx : *m_DataCol.dNdxCol){
+ if (dqdx.getTrack() == const_TrkCol[itrk]){
+ for (int i = 0; i < 5; i++){
+ chi2s[i] = dqdx.getHypotheses(i).chi2;
+ }
+ }
+ }
+
+ for (auto tof : *m_DataCol.tofCol){
+ if (tof.getTrack() == const_TrkCol[itrk]){
+ double toft = tof.getTime();
+ std::array<float, 5> tofexpts = tof.getTimeExp();
+ double tofexpterr = tof.getSigma();
+ std::array<double, 5> tofchi2s;
+ for (int i = 0; i < 5; i++){
+ tofchi2s[i] = std::pow( (tofexpts[i] - toft) / tofexpterr, 2);
+ chi2s[i] += tofchi2s[i];
+ }
+ }
+ }
+
+ int minchi2idx = std::distance(chi2s.begin(), std::min_element(chi2s.begin(), chi2s.end()));
+ int pdgid = m_trk->getCharge() * PDGIDs.at(minchi2idx);
+
+ m_trk->setPID(pdgid);
+//cout<<" Readin track #"<<itrk<<": pid "<<pdgid<<", truth MC pid "<<m_trk->getLeadingMCP().getPDG()<<endl;
+ }
+
m_trkCol.push_back(m_trk);
}
}