diff --git a/Reconstruction/RecPFACyber/include/Algorithm/TrackMatchingAlg.h b/Reconstruction/RecPFACyber/include/Algorithm/TrackMatchingAlg.h
index f13345d624cec3e4eadace5b9d35e4b33caeaff6..542dca11a76b05256f1ec084ed299b39cf6de5b3 100644
--- a/Reconstruction/RecPFACyber/include/Algorithm/TrackMatchingAlg.h
+++ b/Reconstruction/RecPFACyber/include/Algorithm/TrackMatchingAlg.h
@@ -23,7 +23,7 @@ public:
StatusCode RunAlgorithm( CyberDataCol& m_datacol );
StatusCode ClearAlgorithm();
-
+ StatusCode GetGoodTracks(std::vector<std::shared_ptr<Cyber::Track>>& rawTracks, std::vector<std::shared_ptr<Cyber::Track>>& goodTracks);
StatusCode GetExtrpoECALPoints(const Cyber::Track* track, std::vector<TVector3>& extrapo_points);
StatusCode CreateTrackAxis(vector<TVector3>& extrapo_points, std::vector<const Cyber::Calo1DCluster*>& localMaxVCol,
Cyber::CaloHalfCluster* t_track_axis);
@@ -48,7 +48,8 @@ private:
// std::vector<const Cyber::CaloHalfCluster*> m_trackAxisVCol;
// std::vector<const Cyber::CaloHalfCluster*> m_trackAxisUCol;
-
-
+ static bool compTrkP( std::shared_ptr<Cyber::Track> trk1, std::shared_ptr<Cyber::Track> trk2 ){
+ return trk1->getMomentum() > trk2->getMomentum();
+ }
};
#endif
diff --git a/Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h b/Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h
index 8486c16e57d85fd3ff38c42863a2299610ba4132..bff03f0bbe30f94a8bd02263d844bbe5d7baa05b 100644
--- a/Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h
+++ b/Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h
@@ -27,7 +27,7 @@ namespace Cyber {
double getEnergy() const;
TVector3 getPos() const;
- TVector3 getAxis() const { return axis; }
+ TVector3 getAxis() const;
TVector3 getEnergyCenter() const;
std::vector<int> getEnergyCenterTower() const;
int getSlayer() const { return slayer; }
@@ -59,7 +59,7 @@ namespace Cyber {
bool isSubset(const CaloHalfCluster* clus) const;
double OverlapRatioE( const CaloHalfCluster* clus ) const;
- void fitAxis( std::string name );
+ void fitAxis( std::string name ) const;
void setType( int _type ) { type = _type; }
void sortBarShowersByLayer() { std::sort(m_1dclusters.begin(), m_1dclusters.end(), compLayer); }
void addUnit(const Calo1DCluster* _1dcluster);
@@ -85,9 +85,9 @@ namespace Cyber {
int type; // yyy: new definition: track: 10000, Hough: 100, cone: 1, merge: sum them
std::vector< std::vector<int> > towerID; //[module, part, stave]
int slayer;
- TVector3 axis;
- double trk_dr;
- double trk_dz;
+ mutable TVector3 axis = TVector3(99999., 99999., 99999.);
+ mutable double trk_dr;
+ mutable double trk_dz;
double Hough_alpha;
double Hough_rho;
double Hough_intercept;
diff --git a/Reconstruction/RecPFACyber/include/Objects/Track.h b/Reconstruction/RecPFACyber/include/Objects/Track.h
index 76937799e441ec38bae948f29e25e3f4599bc204..6af8328f6ea26185941bdeb3b5d8ce8d3913b1cc 100644
--- a/Reconstruction/RecPFACyber/include/Objects/Track.h
+++ b/Reconstruction/RecPFACyber/include/Objects/Track.h
@@ -16,6 +16,7 @@ namespace Cyber {
void Clear() { m_trackStates.clear(); m_halfClusterUCol.clear(); m_halfClusterVCol.clear(); }
void setOriginTrack(edm4hep::Track& _trk) { m_track = _trk; }
+ int getTrackerHits() const;
int trackStates_size(std::string name) const;
int trackStates_size() const ;
edm4hep::Track getOriginTrack() const { return m_track; }
@@ -26,12 +27,16 @@ namespace Cyber {
std::vector<Cyber::CaloHalfCluster*> getAssociatedHalfClustersV() const { return m_halfClusterVCol; }
std::vector< std::pair<edm4hep::MCParticle, float> > getLinkedMCP() const { return MCParticleWeight; }
edm4hep::MCParticle getLeadingMCP() const;
+ float getD0() const;
+ float getZ0() const;
float getLeadingMCPweight() const;
float getPt() const;
float getPz() const;
float getMomentum() const { return sqrt( getPt()*getPt() + getPz()*getPz() ); }
TVector3 getP3() const;
float getCharge() const;
+ TVector3 getStartPoint() const;
+ TVector3 getEndPoint() const;
void setTrackStates( std::string name, std::vector<TrackState>& _states ) { m_trackStates[name]=_states; }
void addAssociatedHalfClusterU( Cyber::CaloHalfCluster* _cl ) { m_halfClusterUCol.push_back(_cl); }
diff --git a/Reconstruction/RecPFACyber/include/Tools/TrackCreator.h b/Reconstruction/RecPFACyber/include/Tools/TrackCreator.h
index b1095099df66ccd4f3aed0c0471a4e767ccd19a0..16c48c81e8d312bf3b28a419fdf3ffa2c68083cb 100644
--- a/Reconstruction/RecPFACyber/include/Tools/TrackCreator.h
+++ b/Reconstruction/RecPFACyber/include/Tools/TrackCreator.h
@@ -22,7 +22,7 @@ namespace Cyber{
StatusCode CreateTracksFromMCParticle(CyberDataCol& m_DataCol,
DataHandle<edm4hep::MCParticleCollection>& r_MCParticleCol);
-
+ StatusCode SelectGoodTrack(std::vector<std::shared_ptr<Cyber::Track>>& trkCol);
StatusCode Reset(){};
private:
@@ -30,6 +30,15 @@ namespace Cyber{
Cyber::Algorithm* m_TrkExtraAlg;
Cyber::Settings m_TrkExtraSettings;
+ static bool compTrkIP( std::shared_ptr<Cyber::Track> trk1, std::shared_ptr<Cyber::Track> trk2 ){
+ float IP1 = sqrt(trk1->getD0()*trk1->getD0() + trk1->getZ0()*trk1->getZ0() );
+ float IP2 = sqrt(trk2->getD0()*trk2->getD0() + trk2->getZ0()*trk2->getZ0() );
+ return IP1<IP2;
+ }
+
+ static bool compTrkP( std::shared_ptr<Cyber::Track> trk1, std::shared_ptr<Cyber::Track> trk2 ){
+ return trk1->getMomentum() > trk2->getMomentum();
+ }
};
};
diff --git a/Reconstruction/RecPFACyber/script/digi.py b/Reconstruction/RecPFACyber/script/digi.py
index 8402b2ce64e70fcabab925c35b788c2f71c78e5b..a77368738eb3be0175758c4ceb16c7e8038386a2 100644
--- a/Reconstruction/RecPFACyber/script/digi.py
+++ b/Reconstruction/RecPFACyber/script/digi.py
@@ -73,6 +73,7 @@ EcalDigi.CalibrECAL = 1.
EcalDigi.AttenuationLength = 7e10
EcalDigi.TimeResolution = 0.5 #unit: ns
EcalDigi.ChargeThresholdFrac = 0.05
+EcalDigi.EcalMIP_Thre = 0.05
EcalDigi.Debug=1
EcalDigi.WriteNtuple = 1
EcalDigi.OutFileName = "Digi_ECAL.root"
@@ -95,7 +96,7 @@ HcalDigi.CalibrHCAL = 1.
HcalDigi.UseRealisticDigi = 1 # Flag to use digitization model.
HcalDigi.SiPMPixel = 57600 # 57600 for 6025PE (6*6 mm, 25 um pixel pitch)
HcalDigi.TileNonUniformity = 0.0
-HcalDigi.EffAttenLength = 20.
+HcalDigi.EffAttenLength = 1e7
HcalDigi.ADCError = 0.0
HcalDigi.MIPADCMean = 80.*30.0 # Light yield 80 pe/mip
HcalDigi.PeADCMean = 30.0
diff --git a/Reconstruction/RecPFACyber/script/rec.py b/Reconstruction/RecPFACyber/script/rec.py
index fd1de01db1f48229fb0ab16e7687788a73409c5b..6e273e209ab19d0ad227cfe8446cdb73b9c9c781 100644
--- a/Reconstruction/RecPFACyber/script/rec.py
+++ b/Reconstruction/RecPFACyber/script/rec.py
@@ -75,6 +75,7 @@ CyberPFAlg.HCalMCPAssociationName = ["HCALBarrelParticleAssoCol"]
##--- Output collections ---
CyberPFAlg.OutputPFO = "outputPFO";
+CyberPFAlg.RecoPFOCollection = "CyberPFO"
#----Algorithms----
CyberPFAlg.AlgList = ["GlobalClusteringAlg", #1
diff --git a/Reconstruction/RecPFACyber/script/sim.py b/Reconstruction/RecPFACyber/script/sim.py
index 057642463c141ac31f610c395246d243cdbb6575..5763949174aef40692cad33d7a255dade0bdae4f 100644
--- a/Reconstruction/RecPFACyber/script/sim.py
+++ b/Reconstruction/RecPFACyber/script/sim.py
@@ -1,5 +1,5 @@
#!/usr/bin/env python
-import os
+import os,sys
from Gaudi.Configuration import *
from Configurables import k4DataSvc
@@ -91,8 +91,8 @@ from Configurables import CalorimeterSensDetTool
from Configurables import DriftChamberSensDetTool
cal_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
cal_sensdettool.CalNamesMergeDisable = ["EcalBarrel", "CaloDetectorEndcap", "HcalBarrel", "HcalEndcaps"]
-cal_sensdettool.CalNamesApplyBirks = ["EcalBarrel", "CaloDetectorEndcap", "HcalBarrel","HcalEndcaps"]
-cal_sensdettool.CalNamesBirksConstants = [0.008415, 0.008415, 0.126, 0.126] # BGO and Glass scintillator
+#cal_sensdettool.CalNamesApplyBirks = ["EcalBarrel", "CaloDetectorEndcap", "HcalBarrel","HcalEndcaps"]
+#cal_sensdettool.CalNamesBirksConstants = [0.008415, 0.008415, 0.01, 0.01] # BGO and Glass scintillator
# output
from Configurables import PodioOutput
diff --git a/Reconstruction/RecPFACyber/script/tracking.py b/Reconstruction/RecPFACyber/script/tracking.py
index 69af563e0af67b1661b3891455edd8a5756dc3f4..a6d7538fa458435b2d21cf8007adb6d0add91e51 100644
--- a/Reconstruction/RecPFACyber/script/tracking.py
+++ b/Reconstruction/RecPFACyber/script/tracking.py
@@ -7,8 +7,8 @@ dsvc = k4DataSvc("EventDataSvc", input="Sim_TDR_o1_v01_E240_nnHgg.root")
from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
seed = [12340]
-# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
-rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4
rndmengine.SetSingleton = True
rndmengine.Seeds = seed
@@ -51,7 +51,8 @@ podioinput = PodioInput("PodioReader", collections=[
"VXDCollection",
"SITCollection",
"TPCCollection",
- "SETCollection",
+# "SETCollection",
+ "OTKBarrelCollection",
"FTDCollection"
])
@@ -59,23 +60,29 @@ podioinput = PodioInput("PodioReader", collections=[
vxdhitname = "VXDTrackerHits"
sithitname = "SITTrackerHits"
gashitname = "TPCTrackerHits"
-sethitname = "SETTrackerHits"
-setspname = "SETSpacePoints"
+sethitname = "OTKBarrelTrackerHits"
+setspname = "OTKBarrelSpacePoints"
ftdhitname = "FTDTrackerHits"
ftdspname = "FTDSpacePoints"
-from Configurables import PlanarDigiAlg
-digiVXD = PlanarDigiAlg("VXDDigi")
+from Configurables import SmearDigiTool
+vxdtool = SmearDigiTool("VXD")
+vxdtool.ResolutionU = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004]
+vxdtool.ResolutionV = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004]
+vxdtool.UsePlanarTag = True
+vxdtool.ParameterizeResolution = False
+vxdtool.ParametersU = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359]
+vxdtool.ParametersV = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359]
+#vxdtool.OutputLevel = DEBUG
+
+from Configurables import SiTrackerDigiAlg
+digiVXD = SiTrackerDigiAlg("VXDDigi")
digiVXD.SimTrackHitCollection = "VXDCollection"
digiVXD.TrackerHitCollection = vxdhitname
digiVXD.TrackerHitAssociationCollection = "VXDTrackerHitAssociation"
-digiVXD.ResolutionU = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004]
-digiVXD.ResolutionV = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004]
-digiVXD.UsePlanarTag = True
-digiVXD.ParameterizeResolution = False
-digiVXD.ParametersU = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359]
-digiVXD.ParametersV = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359]
+digiVXD.DigiTool = "SmearDigiTool/VXD"
#digiVXD.OutputLevel = DEBUG
+from Configurables import PlanarDigiAlg
digiSIT = PlanarDigiAlg("SITDigi")
digiSIT.IsStrip = False
digiSIT.SimTrackHitCollection = "SITCollection"
@@ -91,9 +98,9 @@ digiSIT.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.296
digiSET = PlanarDigiAlg("SETDigi")
digiSET.IsStrip = False
-digiSET.SimTrackHitCollection = "SETCollection"
+digiSET.SimTrackHitCollection = "OTKBarrelCollection"
digiSET.TrackerHitCollection = sethitname
-digiSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation"
+digiSET.TrackerHitAssociationCollection = "OTKBarrelTrackerHitAssociation"
digiSET.ResolutionU = [0.005]
digiSET.ResolutionV = [0.021]
digiSET.UsePlanarTag = True
@@ -120,6 +127,13 @@ digiTPC = TPCDigiAlg("TPCDigi")
digiTPC.TPCCollection = "TPCCollection"
digiTPC.TPCLowPtCollection = "TPCLowPtCollection"
digiTPC.TPCTrackerHitsCol = gashitname
+#default value, modify them according to future Garfield simulation results
+#digiTPC.PixelClustering = True
+#digiTPC.PointResolutionRPhi = 0.144
+#digiTPC.DiffusionCoeffRPhi = 0.0323
+#digiTPC.PointResolutionZ = 0.4
+#digiTPC.DiffusionCoeffZ = 0.23
+#digiTPC.N_eff = 30
#digiTPC.OutputLevel = DEBUG
# tracking
@@ -223,14 +237,14 @@ from Configurables import TrackParticleRelationAlg
tpr = TrackParticleRelationAlg("Track2Particle")
tpr.MCParticleCollection = "MCParticle"
tpr.TrackList = ["CompleteTracks", "ClupatraTracks"]
-tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
+tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
#tpr.OutputLevel = DEBUG
from Configurables import TrueMuonTagAlg
tmt = TrueMuonTagAlg("TrueMuonTag")
tmt.MCParticleCollection = "MCParticle"
tmt.TrackList = ["CompleteTracks"]
-tmt.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
+tmt.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
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
@@ -248,6 +262,6 @@ mgr = ApplicationMgr(
EvtSel = 'NONE',
EvtMax = 3,
ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc],
- HistogramPersistency = 'ROOT',
+ #HistogramPersistency = 'ROOT',
OutputLevel = INFO
)
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/HoughClusteringAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/HoughClusteringAlg.cpp
index db69d13a16eaf1032906fa0c27ad8b73b9023e56..dec300f17f22c17962151cb0d6f9f22dd6f877a0 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/HoughClusteringAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/HoughClusteringAlg.cpp
@@ -5,7 +5,7 @@
#include <algorithm>
#include <cmath>
#include <set>
-#include "TCanvas.h"
+
using namespace std;
StatusCode HoughClusteringAlg::ReadSettings(Cyber::Settings& m_settings){
@@ -104,7 +104,6 @@ StatusCode HoughClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
if(p_HalfClusterV.size()==0){ std::cout<<" HoughClusteringAlg: No HalfClusterV in present data collection! "<<std::endl; }
if(p_HalfClusterU.size()==0){ std::cout<<" HoughClusteringAlg: No HalfClusterU in present data collection! "<<std::endl; }
-//cout<<"Readin HalfCluster size: "<<p_HalfClusterV.size()<<", "<<p_HalfClusterU.size()<<endl;
//std::vector<const Cyber::CaloHalfCluster*> m_refHFClusVCol; m_refHFClusVCol.clear();
// Processing V(xy) plane
@@ -122,17 +121,13 @@ StatusCode HoughClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
continue;
}
-//cout<<" HoughClusteringAlg: Find Hough axis in HalfCluster "<<it<<". Local maximum size V = "<<m_localMaxVCol.size()<<endl;
- // cout<<" HoughClusteringAlg: Creating m_HoughObjectsV"<<endl;
std::vector<Cyber::HoughObject> m_HoughObjectsV; m_HoughObjectsV.clear();
for(int il=0; il<m_localMaxVCol.size(); il++){
Cyber::HoughObject m_obj(m_localMaxVCol[il], Cyber::CaloUnit::barsize, Cyber::CaloUnit::ecal_innerR);
m_HoughObjectsV.push_back(m_obj);
}
-//cout<<" HoughClusteringAlg: HoughObjectV size "<<m_HoughObjectsV.size()<<endl;
- // cout<<" HoughClusteringAlg: Hough transformation"<<endl;
HoughTransformation(m_HoughObjectsV);
// cout<<" HoughClusteringAlg: Creating hough_spaceV"<<endl;
@@ -141,16 +136,13 @@ StatusCode HoughClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
settings.map_floatPars["rho_low"], settings.map_floatPars["rho_high"],
settings.map_floatPars["bin_width_rho"], settings.map_intPars["Nbins_rho"]);
-//cout<<" HoughClusteringAlg: Filling hough_spaceV"<<endl;
FillHoughSpace(m_HoughObjectsV, hough_spaceV);
-//cout<<" HoughClusteringAlg: Finding clusters from Hough space"<<endl;
-
//Create output HoughClusters
m_longiClusVCol.clear();
ClusterFinding(m_HoughObjectsV, hough_spaceV, m_longiClusVCol );
+
CleanClusters(m_longiClusVCol);
-//cout << " HoughClusteringAlg: final output m_longiClusVCol.size() = " << m_longiClusVCol.size() << endl;
m_datacol.map_HalfCluster["bkHalfCluster"].insert( m_datacol.map_HalfCluster["bkHalfCluster"].end(), m_longiClusVCol.begin(), m_longiClusVCol.end() );
std::vector<const Cyber::CaloHalfCluster*> m_constHoughCluster; m_constHoughCluster.clear();
@@ -168,7 +160,6 @@ StatusCode HoughClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
for(int ic=0; ic<m_longiClusVCol.size(); ic++)
m_constHoughCluster.push_back(m_longiClusVCol[ic].get());
- //m_refHFClusVCol.insert(m_refHFClusVCol.end(), m_constHoughCluster.begin(), m_constHoughCluster.end());
p_HalfClusterV[it]->setLocalMax("HoughLocalMax", m_houghMax);
p_HalfClusterV[it]->setLocalMax(settings.map_stringPars["LeftLocalMaxName"], left_localMaxVCol);
@@ -177,7 +168,6 @@ StatusCode HoughClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
left_localMaxVCol.clear();
} // end of V plane
-//cout<<"Finish Hough in V. Reference HFClusterV size "<<m_refHFClusVCol.size()<<endl;
// Processing U(r-phi) plane
@@ -191,21 +181,16 @@ StatusCode HoughClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
}
if(m_localMaxUCol.size()<settings.map_intPars["th_peak"]){
- //std::cout << " yyy: m_localMaxUCol.size()<th_peak, continue" << std::endl;
continue;
}
-//cout<<" HoughClusteringAlg: Find Hough axis in HalfCluster "<<it<<". Local maximum size U = "<<m_localMaxUCol.size()<<endl;
- // cout<<" HoughClusteringAlg: Creating m_HoughObjectsU"<<endl;
std::vector<Cyber::HoughObject> m_HoughObjectsU; m_HoughObjectsU.clear();
for(int il=0; il<m_localMaxUCol.size(); il++){
Cyber::HoughObject m_obj(m_localMaxUCol[il], Cyber::CaloUnit::barsize, Cyber::CaloUnit::ecal_innerR);
m_HoughObjectsU.push_back(m_obj);
}
-//cout<<" HoughClusteringAlg: HoughObjectU size "<<m_HoughObjectsU.size()<<endl;
- // cout<<" HoughClusteringAlg: Hough transformation"<<endl;
HoughTransformation(m_HoughObjectsU);
// cout<<" HoughClusteringAlg: Creating hough_spaceU"<<endl;
@@ -213,16 +198,13 @@ StatusCode HoughClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
settings.map_floatPars["bin_width_alphaU"], settings.map_intPars["Nbins_alphaU"],
settings.map_floatPars["rho_low"], settings.map_floatPars["rho_high"],
settings.map_floatPars["bin_width_rho"], settings.map_intPars["Nbins_rho"]);
-
-//cout<<" HoughClusteringAlg: Filling hough_spaceU"<<endl;
FillHoughSpace(m_HoughObjectsU, hough_spaceU);
-//cout<<" HoughClusteringAlg: Finding clusters from Hough space"<<endl;
//Create output HoughClusters
m_longiClusUCol.clear();
ClusterFinding(m_HoughObjectsU, hough_spaceU, m_longiClusUCol );
CleanClusters(m_longiClusUCol);
-//cout << " HoughClusteringAlg: final output m_longiClusUCol.size() = " << m_longiClusUCol.size() << endl;
+
m_datacol.map_HalfCluster["bkHalfCluster"].insert( m_datacol.map_HalfCluster["bkHalfCluster"].end(), m_longiClusUCol.begin(), m_longiClusUCol.end() );
std::vector<const Cyber::CaloHalfCluster*> m_constHoughCluster; m_constHoughCluster.clear();
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp
index acc9165da38f61223c3e9a4f84610054da499633..f7d59919e200afab5a4136f32e904963ed6f4cd7 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp
@@ -14,8 +14,8 @@ StatusCode PFOReclusteringAlg::ReadSettings(Settings& m_settings){
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.;
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.6;
- if(settings.map_floatPars.find("MinAngleForNeuMerge")==settings.map_floatPars.end()) settings.map_floatPars["MinAngleForNeuMerge"] = 0.18;
+ if(settings.map_floatPars.find("VirtualMergeSigma")==settings.map_floatPars.end()) settings.map_floatPars["VirtualMergeSigma"] = 0.5;
+ if(settings.map_floatPars.find("MinAngleForNeuMerge")==settings.map_floatPars.end()) settings.map_floatPars["MinAngleForNeuMerge"] = 0.12;
if(settings.map_floatPars.find("MinAngleForVirMerge")==settings.map_floatPars.end()) settings.map_floatPars["MinAngleForVirMerge"] = 0.12;
@@ -66,7 +66,9 @@ StatusCode PFOReclusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
//If P_trk < E_cluster, create a virtual neutral PFO.
ReCluster_SplitFromChg(m_chargedPFOs, m_neutralPFOs);
+
/*
+ cout<<" PFO after ReCluster_SplitFromChg: "<<p_PFObjects->size()<<", charged "<<m_chargedPFOs.size()<<", neutral "<<m_neutralPFOs.size()<<endl;
totE_Ecal = 0;
totE_Hcal = 0;
cout<<"After split from Ch: charged "<<m_chargedPFOs.size()<<", neutral "<<m_neutralPFOs.size()<<", total "<<p_PFObjects->size()<<endl;
@@ -96,9 +98,10 @@ StatusCode PFOReclusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
m_datacol.map_CaloHit["bkHit"].insert( m_datacol.map_CaloHit["bkHit"].end(), m_bkCol.map_CaloHit["bkHit"].begin(), m_bkCol.map_CaloHit["bkHit"].end() );
m_datacol.map_CaloCluster["bk3DCluster"].insert( m_datacol.map_CaloCluster["bk3DCluster"].end(), m_bkCol.map_CaloCluster["bk3DCluster"].begin(), m_bkCol.map_CaloCluster["bk3DCluster"].end() );
m_datacol.map_PFObjects["bkPFO"].insert( m_datacol.map_PFObjects["bkPFO"].end(), m_bkCol.map_PFObjects["bkPFO"].begin(), m_bkCol.map_PFObjects["bkPFO"].end() );
+
/*
- totE_Ecal = 0;
- totE_Hcal = 0;
+ double totE_Ecal = 0;
+ double totE_Hcal = 0;
cout<<"After merge all virtual to Ch: charged "<<m_chargedPFOs.size()<<", neutral "<<m_neutralPFOs.size()<<", total "<<p_PFObjects->size()<<endl;
for(int i=0; i<m_neutralPFOs.size(); i++){
cout<<" PFO #"<<i<<": track size "<<m_neutralPFOs[i]->getTracks().size()<<", leading P "<<m_neutralPFOs[i]->getTrackMomentum();
@@ -365,10 +368,11 @@ StatusCode PFOReclusteringAlg::ReCluster_MergeToChg(std::vector< std::shared_ptr
double tmp_delta_E = delta_energy + settings.map_floatPars["HCALCalib"]*all_neutral_HCAL_clus[clus_index]->getHitsE();
// cout<<" If include this cluster: new deltaE "<<tmp_delta_E<<", merge = "<<(tmp_delta_E > sigmaE * settings.map_floatPars["VirtualMergeSigma"])<<endl;
+
if(tmp_delta_E > sigmaE * settings.map_floatPars["VirtualMergeSigma"]){
double absorbed_energy = sigmaE*settings.map_floatPars["VirtualMergeSigma"] - delta_energy;
delta_energy = delta_energy + absorbed_energy;
-
+// cout<<" Absorberd energy: "<<absorbed_energy<<", current delta energy "<<delta_energy<<endl;
//Create a new virtual neutral cluster with energy = absorbed_energy.
@@ -396,7 +400,6 @@ StatusCode PFOReclusteringAlg::ReCluster_MergeToChg(std::vector< std::shared_ptr
if(tmp_HCAL_clus[0]==all_neutral_HCAL_clus[clus_index]){
std::shared_ptr<Cyber::CaloHit> m_newhit = all_neutral_HCAL_clus[clus_index]->getCaloHits()[0]->Clone();
m_newhit->setEnergy(all_neutral_HCAL_clus[clus_index]->getHitsE() - absorbed_energy/settings.map_floatPars["HCALCalib"] );
-
std::shared_ptr<Cyber::Calo3DCluster> m_newclus = std::make_shared<Cyber::Calo3DCluster>();
m_newclus->addHit(m_newhit.get());
m_newclus->setType(-1);
@@ -405,7 +408,7 @@ StatusCode PFOReclusteringAlg::ReCluster_MergeToChg(std::vector< std::shared_ptr
m_bkCol.map_CaloCluster["bk3DCluster"].push_back(m_newclus);
std::vector<const Calo3DCluster*> tmp_clusters; tmp_clusters.clear();
- tmp_clusters.push_back(m_clus.get());
+ tmp_clusters.push_back(m_newclus.get());
m_neutralPFOs[ip]->setHCALCluster( tmp_clusters );
is_found = true;
diff --git a/Reconstruction/RecPFACyber/src/Algorithm/TrackClusterConnectingAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/TrackClusterConnectingAlg.cpp
index fa79d1c3acce328f0653a2e2f8ef02e5f9cad3cb..5d3ff5807bdc6ce9f73fb70930b19d1501d2f393 100644
--- a/Reconstruction/RecPFACyber/src/Algorithm/TrackClusterConnectingAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/Algorithm/TrackClusterConnectingAlg.cpp
@@ -13,7 +13,7 @@ StatusCode TrackClusterConnectingAlg::ReadSettings(Settings& m_settings){
if(settings.map_floatPars.find("th_ChFragEn")==settings.map_floatPars.end()) settings.map_floatPars["th_ChFragEn"] = 2.;
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_HcalMatchingaR")==settings.map_floatPars.end()) settings.map_floatPars["th_HcalMatchingR"] = 100.;
+ if(settings.map_floatPars.find("th_HcalMatchingR")==settings.map_floatPars.end()) settings.map_floatPars["th_HcalMatchingR"] = 150.;
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;
@@ -44,6 +44,9 @@ StatusCode TrackClusterConnectingAlg::Initialize( CyberDataCol& m_datacol ){
}
//cout<<"Readin Track size: "<<m_tracks.size()<<", ECAL cluster size: "<<m_EcalClusters.size()<<", HCAL cluster size "<<m_HcalClusters.size()<<endl;
+//cout<<"Print track"<<endl;
+//for(int i=0; i<m_tracks.size(); i++)
+// cout<<"Track #"<<i<<": P = "<<m_tracks[i]->getMomentum()<<", Pt = "<<m_tracks[i]->getPt()<<endl;
//cout<<"Print all ECAL cluster "<<endl;
//for(int ic=0; ic<m_EcalClusters.size(); ic++){
// cout<<" ECAL Cluster #"<<ic<<": En = "<<m_EcalClusters[ic]->getLongiE()<<", track size "<<m_EcalClusters[ic]->getAssociatedTracks().size();
@@ -64,6 +67,12 @@ StatusCode TrackClusterConnectingAlg::RunAlgorithm( CyberDataCol& m_datacol ){
m_absorbedEcal.clear();
EcalChFragAbsorption(m_EcalClusters, m_tracks, m_absorbedEcal);
//cout<<" TrackClusterConnectingAlg: After ECAL charged fragment absorption: cluster size "<<m_absorbedEcal.size()<<endl;
+//cout<<"Print merged ECAL cluster "<<endl;
+//for(int ic=0; ic<m_absorbedEcal.size(); ic++){
+// cout<<" ECAL Cluster #"<<ic<<": En = "<<m_absorbedEcal[ic]->getLongiE()<<", track size "<<m_absorbedEcal[ic]->getAssociatedTracks().size();
+// if(m_absorbedEcal[ic]->getAssociatedTracks().size()>0) cout<<", Leading track P = "<<m_absorbedEcal[ic]->getAssociatedTracks()[0]->getMomentum()<<endl;
+// else cout<<endl;
+//}
//2. Create PFObject with ECAL cluster and track
std::vector<const Cyber::Calo3DCluster*> tmp_constClus;
@@ -77,6 +86,18 @@ StatusCode TrackClusterConnectingAlg::RunAlgorithm( CyberDataCol& m_datacol ){
// cout<<", HCAL cluster size "<<m_PFObjects[i]->getHCALClusters().size()<<", totE "<<m_PFObjects[i]->getHCALClusterEnergy()<<endl;
//}
+//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,
+// m_HcalClusters[i]->getHitsE(),
+// m_HcalClusters[i]->getHitCenter().x(),
+// m_HcalClusters[i]->getHitCenter().y(),
+// m_HcalClusters[i]->getHitCenter().z() );
+// totE_Hcal += m_HcalClusters[i]->getHitsE();
+//}
+//cout<<"Hcal cluster total E "<<totE_Hcal<<endl;
+
//3. Add HCAL clusters into the PFObject.
std::sort(m_PFObjects.begin(), m_PFObjects.end(), compTrkP);
HcalExtrapolatingMatch(m_HcalClusters, m_PFObjects);
diff --git a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
index 06158fed9d5503a566b1eb1f801e14b20bcd39dd..70185d2c88b8102f2a1f95efc8b3509bb11ed055 100644
--- a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
+++ b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp
@@ -52,6 +52,15 @@ StatusCode CyberPFAlg::initialize()
m_pTrackCreatorSettings.map_stringVecPars["trackCollections"] = name_TrackCol.value();
m_pTrackCreatorSettings.map_floatPars["BField"] = m_BField;
+ m_pTrackCreatorSettings.map_floatPars["TrkEndZCut"] = 200.;
+ m_pTrackCreatorSettings.map_floatPars["TrkEndRCutMin"] = 700.;
+ m_pTrackCreatorSettings.map_floatPars["TrkEndRCutMax"] = 1600.;
+ m_pTrackCreatorSettings.map_floatPars["TrkStartRCutMin"] = 635.;
+ m_pTrackCreatorSettings.map_floatPars["TrkStartRCutMax"] = 640.;
+ m_pTrackCreatorSettings.map_floatPars["TrkLengthCut"] = 130.;
+ m_pTrackCreatorSettings.map_floatPars["BrokenTrkMinP"] = 0.5;
+ m_pTrackCreatorSettings.map_floatPars["BrokenTrkDeltaPCut"] = 0.15;
+ m_pTrackCreatorSettings.map_floatPars["BrokenTrkDistance"] = 10.;
std::vector<std::string> name_CaloHits = name_EcalHits;
std::vector<std::string> name_CaloReadout = name_EcalReadout;
@@ -659,7 +668,7 @@ StatusCode CyberPFAlg::execute()
if(_nEvt==0) std::cout<<"CyberPFAlg::execute Start"<<std::endl;
std::cout<<"Processing event: "<<_nEvt<<std::endl;
- if(_nEvt<m_Nskip){ _nEvt++; return GaudiAlgorithm::initialize(); }
+ if(_nEvt<m_Nskip){ _nEvt++; return StatusCode::SUCCESS; }
//InitializeForNewEvent();
CyberDataCol m_DataCol;
@@ -669,10 +678,8 @@ StatusCode CyberPFAlg::execute()
//Readin collections
m_pMCParticleCreator->CreateMCParticle( m_DataCol, *r_MCParticleCol );
-
if(m_useTruthTrk) m_pTrackCreator->CreateTracksFromMCParticle(m_DataCol, *r_MCParticleCol);
else m_pTrackCreator->CreateTracks( m_DataCol, r_TrackCols, r_MCPTrkAssoCol );
-
m_pCaloHitsCreator->CreateCaloHits( m_DataCol, r_CaloHitCols, map_readout_decoder, map_CaloMCPAssoCols );
//Perform PFA algorithm
@@ -688,912 +695,915 @@ StatusCode CyberPFAlg::execute()
// yyy_endrec = clock(); // 閲嶅缓缁撴潫鐨勬椂闂�
-cout<<"Write tuples"<<endl;
- //---------------------Write Ana tuples-------------------------
- // MC particles
- ClearMCParticle();
- std::vector<edm4hep::MCParticle> m_MCPCol = m_DataCol.collectionMap_MC[name_MCParticleCol.value()];
- for(int imc=0; imc<m_MCPCol.size(); imc++){
- m_mcPdgid.push_back( m_MCPCol[imc].getPDG() );
- m_mcStatus.push_back( m_MCPCol[imc].getGeneratorStatus() );
- m_mcPx.push_back( m_MCPCol[imc].getMomentum()[0] );
- m_mcPy.push_back( m_MCPCol[imc].getMomentum()[1] );
- m_mcPz.push_back( m_MCPCol[imc].getMomentum()[2] );
- m_mcEn.push_back( m_MCPCol[imc].getEnergy() );
- m_mcMass.push_back( m_MCPCol[imc].getMass() );
- m_mcCharge.push_back( m_MCPCol[imc].getCharge() );
- m_mcEPx.push_back( m_MCPCol[imc].getEndpoint()[0] );
- m_mcEPy.push_back( m_MCPCol[imc].getEndpoint()[1] );
- m_mcEPz.push_back( m_MCPCol[imc].getEndpoint()[2] );
- //double tmp_phi = std::atan2(m_MCPCol[imc].getMomentum()[1], m_MCPCol[imc].getMomentum()[0])* 180.0 / M_PI;
- //if (tmp_phi < 0) tmp_phi += 360.0;
- //double tmp_theta = std::atan2(m_MCPCol[imc].getMomentum()[2], sqrt(m_MCPCol[imc].getMomentum()[1]*m_MCPCol[imc].getMomentum()[1]+m_MCPCol[imc].getMomentum()[0]*m_MCPCol[imc].getMomentum()[0]))* 180.0 / M_PI + 90;
- //cout<<"MCParticle: "<<imc<<" PDG: "<<m_MCPCol[imc].getPDG()<<" Theta: "<<tmp_theta<<" Phi: "<<tmp_phi<<endl;
-
- double EnDep_ecal = GetParticleDepEnergy(m_MCPCol[imc], m_DataCol.map_BarCol["BarCol"]);
- double EnDep_hcal = GetParticleDepEnergy(m_MCPCol[imc], m_DataCol.map_CaloHit["HCALBarrel"]);
- m_depEn_ecal.push_back(EnDep_ecal);
- m_depEn_hcal.push_back(EnDep_hcal);
- }
- t_MCParticle->Fill();
-
-
- //Save Raw bars information
- ClearBar();
- m_totE_EcalSim = 0.;
- for(int ibar=0;ibar<m_DataCol.map_BarCol["BarCol"].size();ibar++){
- auto p_hitbar = m_DataCol.map_BarCol["BarCol"][ibar].get();
- m_simBar_x.push_back(p_hitbar->getPosition().x());
- m_simBar_y.push_back(p_hitbar->getPosition().y());
- m_simBar_z.push_back(p_hitbar->getPosition().z());
- m_simBar_Q1.push_back(p_hitbar->getQ1());
- m_simBar_Q2.push_back(p_hitbar->getQ2());
- m_simBar_T1.push_back(p_hitbar->getT1());
- m_simBar_T2.push_back(p_hitbar->getT2());
- m_simBar_module.push_back(p_hitbar->getModule());
- m_simBar_dlayer.push_back(p_hitbar->getDlayer());
- m_simBar_stave.push_back(p_hitbar->getStave());
- m_simBar_slayer.push_back(p_hitbar->getSlayer());
- m_simBar_bar.push_back(p_hitbar->getBar());
- m_totE_EcalSim += (p_hitbar->getQ1()+p_hitbar->getQ2())/2.;
-
- auto truthMap = p_hitbar->getLinkedMCP();
- for(auto iter: truthMap){
- m_simBar_truthMC_tag.push_back(ibar);
- m_simBar_truthMC_pid.push_back(iter.first.getPDG());
- m_simBar_truthMC_px.push_back(iter.first.getMomentum().x);
- m_simBar_truthMC_py.push_back(iter.first.getMomentum().y);
- m_simBar_truthMC_pz.push_back(iter.first.getMomentum().z);
- m_simBar_truthMC_E.push_back(iter.first.getEnergy());
- m_simBar_truthMC_EPx.push_back(iter.first.getEndpoint().x);
- m_simBar_truthMC_EPy.push_back(iter.first.getEndpoint().y);
- m_simBar_truthMC_EPz.push_back(iter.first.getEndpoint().z);
- m_simBar_truthMC_weight.push_back(iter.second);
+ if(m_WriteAna){
+ cout<<"Write tuples"<<endl;
+ //---------------------Write Ana tuples-------------------------
+ // MC particles
+ ClearMCParticle();
+ std::vector<edm4hep::MCParticle> m_MCPCol = m_DataCol.collectionMap_MC[name_MCParticleCol.value()];
+ for(int imc=0; imc<m_MCPCol.size(); imc++){
+ m_mcPdgid.push_back( m_MCPCol[imc].getPDG() );
+ m_mcStatus.push_back( m_MCPCol[imc].getGeneratorStatus() );
+ m_mcPx.push_back( m_MCPCol[imc].getMomentum()[0] );
+ m_mcPy.push_back( m_MCPCol[imc].getMomentum()[1] );
+ m_mcPz.push_back( m_MCPCol[imc].getMomentum()[2] );
+ m_mcEn.push_back( m_MCPCol[imc].getEnergy() );
+ m_mcMass.push_back( m_MCPCol[imc].getMass() );
+ m_mcCharge.push_back( m_MCPCol[imc].getCharge() );
+ m_mcEPx.push_back( m_MCPCol[imc].getEndpoint()[0] );
+ m_mcEPy.push_back( m_MCPCol[imc].getEndpoint()[1] );
+ m_mcEPz.push_back( m_MCPCol[imc].getEndpoint()[2] );
+ //double tmp_phi = std::atan2(m_MCPCol[imc].getMomentum()[1], m_MCPCol[imc].getMomentum()[0])* 180.0 / M_PI;
+ //if (tmp_phi < 0) tmp_phi += 360.0;
+ //double tmp_theta = std::atan2(m_MCPCol[imc].getMomentum()[2], sqrt(m_MCPCol[imc].getMomentum()[1]*m_MCPCol[imc].getMomentum()[1]+m_MCPCol[imc].getMomentum()[0]*m_MCPCol[imc].getMomentum()[0]))* 180.0 / M_PI + 90;
+ //cout<<"MCParticle: "<<imc<<" PDG: "<<m_MCPCol[imc].getPDG()<<" Theta: "<<tmp_theta<<" Phi: "<<tmp_phi<<endl;
+
+ double EnDep_ecal = GetParticleDepEnergy(m_MCPCol[imc], m_DataCol.map_BarCol["BarCol"]);
+ double EnDep_hcal = GetParticleDepEnergy(m_MCPCol[imc], m_DataCol.map_CaloHit["HCALBarrel"]);
+ m_depEn_ecal.push_back(EnDep_ecal);
+ m_depEn_hcal.push_back(EnDep_hcal);
}
- }
-
- std::vector<Cyber::CaloHit*> m_hcalHitsCol; m_hcalHitsCol.clear();
- for(int ih=0; ih<m_DataCol.map_CaloHit["HCALBarrel"].size(); ih++)
- m_hcalHitsCol.push_back( m_DataCol.map_CaloHit["HCALBarrel"][ih].get() );
-
- m_totE_HcalSim = 0.;
- for(int ihit=0; ihit<m_hcalHitsCol.size(); ihit++){
- m_HcalHit_x.push_back( m_hcalHitsCol[ihit]->getPosition().x() );
- m_HcalHit_y.push_back( m_hcalHitsCol[ihit]->getPosition().y() );
- m_HcalHit_z.push_back( m_hcalHitsCol[ihit]->getPosition().z() );
- m_HcalHit_E.push_back( m_hcalHitsCol[ihit]->getEnergy() );
- m_HcalHit_layer.push_back( m_hcalHitsCol[ihit]->getLayer() );
- m_totE_HcalSim += m_hcalHitsCol[ihit]->getEnergy();
-
- auto truthMap = m_hcalHitsCol[ihit]->getLinkedMCP();
- for(auto iter: truthMap){
- m_HcalHit_truthMC_tag.push_back(ihit);
- m_HcalHit_truthMC_pid.push_back(iter.first.getPDG());
- m_HcalHit_truthMC_px.push_back(iter.first.getMomentum().x);
- m_HcalHit_truthMC_py.push_back(iter.first.getMomentum().y);
- m_HcalHit_truthMC_pz.push_back(iter.first.getMomentum().z);
- m_HcalHit_truthMC_E.push_back(iter.first.getEnergy());
- m_HcalHit_truthMC_EPx.push_back(iter.first.getEndpoint().x);
- m_HcalHit_truthMC_EPy.push_back(iter.first.getEndpoint().y);
- m_HcalHit_truthMC_EPz.push_back(iter.first.getEndpoint().z);
- m_HcalHit_truthMC_weight.push_back(iter.second);
- }
-
- }
- t_SimBar->Fill();
-
- //Save localMax
- ClearLocalMax();
- std::vector<Cyber::CaloHalfCluster*> m_halfclusters; m_halfclusters.clear();
- for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColU"].size(); i++)
- m_halfclusters.push_back( m_DataCol.map_HalfCluster["HalfClusterColU"][i].get() );
-
- std::vector<const Calo1DCluster*> m_local_max; m_local_max.clear();
- for(int ic=0;ic<m_halfclusters.size(); ic++){
- std::vector<const Calo1DCluster*> tmp_shower = m_halfclusters[ic]->getLocalMaxCol("AllLocalMax");
- m_local_max.insert(m_local_max.end(), tmp_shower.begin(), tmp_shower.end());
- }
- for(int il=0; il<m_local_max.size(); il++){
- m_localMaxU_tag.push_back( il );
- m_localMaxU_x.push_back( m_local_max[il]->getPos().x() );
- m_localMaxU_y.push_back( m_local_max[il]->getPos().y() );
- m_localMaxU_z.push_back( m_local_max[il]->getPos().z() );
- m_localMaxU_E.push_back( m_local_max[il]->getEnergy() );
-
- auto truthMap = m_local_max[il]->getLinkedMCP();
- for(auto iter: truthMap){
- m_localMaxU_mc_tag.push_back(il);
- m_localMaxU_mc_pdg.push_back(iter.first.getPDG());
- m_localMaxU_mc_px.push_back(iter.first.getMomentum().x);
- m_localMaxU_mc_py.push_back(iter.first.getMomentum().y);
- m_localMaxU_mc_pz.push_back(iter.first.getMomentum().z);
- m_localMaxU_mc_weight.push_back(iter.second);
+ t_MCParticle->Fill();
+
+
+ //Save Raw bars information
+ ClearBar();
+ m_totE_EcalSim = 0.;
+ for(int ibar=0;ibar<m_DataCol.map_BarCol["BarCol"].size();ibar++){
+ auto p_hitbar = m_DataCol.map_BarCol["BarCol"][ibar].get();
+ m_simBar_x.push_back(p_hitbar->getPosition().x());
+ m_simBar_y.push_back(p_hitbar->getPosition().y());
+ m_simBar_z.push_back(p_hitbar->getPosition().z());
+ m_simBar_Q1.push_back(p_hitbar->getQ1());
+ m_simBar_Q2.push_back(p_hitbar->getQ2());
+ m_simBar_T1.push_back(p_hitbar->getT1());
+ m_simBar_T2.push_back(p_hitbar->getT2());
+ m_simBar_module.push_back(p_hitbar->getModule());
+ m_simBar_dlayer.push_back(p_hitbar->getDlayer());
+ m_simBar_stave.push_back(p_hitbar->getStave());
+ m_simBar_slayer.push_back(p_hitbar->getSlayer());
+ m_simBar_bar.push_back(p_hitbar->getBar());
+ m_totE_EcalSim += (p_hitbar->getQ1()+p_hitbar->getQ2())/2.;
+
+ auto truthMap = p_hitbar->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_simBar_truthMC_tag.push_back(ibar);
+ m_simBar_truthMC_pid.push_back(iter.first.getPDG());
+ m_simBar_truthMC_px.push_back(iter.first.getMomentum().x);
+ m_simBar_truthMC_py.push_back(iter.first.getMomentum().y);
+ m_simBar_truthMC_pz.push_back(iter.first.getMomentum().z);
+ m_simBar_truthMC_E.push_back(iter.first.getEnergy());
+ m_simBar_truthMC_EPx.push_back(iter.first.getEndpoint().x);
+ m_simBar_truthMC_EPy.push_back(iter.first.getEndpoint().y);
+ m_simBar_truthMC_EPz.push_back(iter.first.getEndpoint().z);
+ m_simBar_truthMC_weight.push_back(iter.second);
+ }
}
- }
- m_halfclusters.clear();
- m_local_max.clear();
- for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColV"].size(); i++)
- m_halfclusters.push_back( m_DataCol.map_HalfCluster["HalfClusterColV"][i].get() );
-
- for(int ic=0;ic<m_halfclusters.size(); ic++){
- std::vector<const Calo1DCluster*> tmp_shower = m_halfclusters[ic]->getLocalMaxCol("AllLocalMax");
- m_local_max.insert(m_local_max.end(), tmp_shower.begin(), tmp_shower.end());
- }
- for(int il=0; il<m_local_max.size(); il++){
- m_localMaxV_tag.push_back( il );
- m_localMaxV_x.push_back( m_local_max[il]->getPos().x() );
- m_localMaxV_y.push_back( m_local_max[il]->getPos().y() );
- m_localMaxV_z.push_back( m_local_max[il]->getPos().z() );
- m_localMaxV_E.push_back( m_local_max[il]->getEnergy() );
-
- auto truthMap = m_local_max[il]->getLinkedMCP();
- for(auto iter: truthMap){
- m_localMaxV_mc_tag.push_back(il);
- m_localMaxV_mc_pdg.push_back(iter.first.getPDG());
- m_localMaxV_mc_px.push_back(iter.first.getMomentum().x);
- m_localMaxV_mc_py.push_back(iter.first.getMomentum().y);
- m_localMaxV_mc_pz.push_back(iter.first.getMomentum().z);
- m_localMaxV_mc_weight.push_back(iter.second);
+
+ std::vector<Cyber::CaloHit*> m_hcalHitsCol; m_hcalHitsCol.clear();
+ for(int ih=0; ih<m_DataCol.map_CaloHit["HCALBarrel"].size(); ih++)
+ m_hcalHitsCol.push_back( m_DataCol.map_CaloHit["HCALBarrel"][ih].get() );
+
+ m_totE_HcalSim = 0.;
+ for(int ihit=0; ihit<m_hcalHitsCol.size(); ihit++){
+ m_HcalHit_x.push_back( m_hcalHitsCol[ihit]->getPosition().x() );
+ m_HcalHit_y.push_back( m_hcalHitsCol[ihit]->getPosition().y() );
+ m_HcalHit_z.push_back( m_hcalHitsCol[ihit]->getPosition().z() );
+ m_HcalHit_E.push_back( m_hcalHitsCol[ihit]->getEnergy() );
+ m_HcalHit_layer.push_back( m_hcalHitsCol[ihit]->getLayer() );
+ m_totE_HcalSim += m_hcalHitsCol[ihit]->getEnergy();
+
+ auto truthMap = m_hcalHitsCol[ihit]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_HcalHit_truthMC_tag.push_back(ihit);
+ m_HcalHit_truthMC_pid.push_back(iter.first.getPDG());
+ m_HcalHit_truthMC_px.push_back(iter.first.getMomentum().x);
+ m_HcalHit_truthMC_py.push_back(iter.first.getMomentum().y);
+ m_HcalHit_truthMC_pz.push_back(iter.first.getMomentum().z);
+ m_HcalHit_truthMC_E.push_back(iter.first.getEnergy());
+ m_HcalHit_truthMC_EPx.push_back(iter.first.getEndpoint().x);
+ m_HcalHit_truthMC_EPy.push_back(iter.first.getEndpoint().y);
+ m_HcalHit_truthMC_EPz.push_back(iter.first.getEndpoint().z);
+ m_HcalHit_truthMC_weight.push_back(iter.second);
+ }
+
}
- }
- t_LocalMax->Fill();
-
- //Save 1DCluster
- ClearLayer();
- m_halfclusters.clear();
- for(int i=0; i<m_DataCol.map_HalfCluster["ESHalfClusterU"].size(); i++)
- m_halfclusters.push_back( m_DataCol.map_HalfCluster["ESHalfClusterU"][i].get() );
-
- m_local_max.clear();
- for(int ic=0;ic<m_halfclusters.size(); ic++){
- //std::vector<const Calo1DCluster*> tmp_shower = m_halfclusters[ic]->getLocalMaxCol("AllLocalMax");
- std::vector<const CaloHalfCluster*> m_axis = m_halfclusters[ic]->getHalfClusterCol("MergedAxis");
- if(m_axis.size()>0){
- std::vector<const Calo1DCluster*> tmp_shower = m_axis[0]->getCluster();
+ t_SimBar->Fill();
+
+ //Save localMax
+ ClearLocalMax();
+ std::vector<Cyber::CaloHalfCluster*> m_halfclusters; m_halfclusters.clear();
+ for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColU"].size(); i++)
+ m_halfclusters.push_back( m_DataCol.map_HalfCluster["HalfClusterColU"][i].get() );
+
+ std::vector<const Calo1DCluster*> m_local_max; m_local_max.clear();
+ for(int ic=0;ic<m_halfclusters.size(); ic++){
+ std::vector<const Calo1DCluster*> tmp_shower = m_halfclusters[ic]->getLocalMaxCol("AllLocalMax");
m_local_max.insert(m_local_max.end(), tmp_shower.begin(), tmp_shower.end());
}
- }
- for(int il=0; il<m_local_max.size(); il++){
- m_barShowerU_tag.push_back( il );
- m_barShowerU_x.push_back( m_local_max[il]->getPos().x() );
- m_barShowerU_y.push_back( m_local_max[il]->getPos().y() );
- m_barShowerU_z.push_back( m_local_max[il]->getPos().z() );
- m_barShowerU_E.push_back( m_local_max[il]->getEnergy() );
-
- auto truthMap = m_local_max[il]->getLinkedMCP();
- for(auto iter: truthMap){
- m_barShowerU_mc_tag.push_back(il);
- m_barShowerU_mc_pdg.push_back(iter.first.getPDG());
- m_barShowerU_mc_px.push_back(iter.first.getMomentum().x);
- m_barShowerU_mc_py.push_back(iter.first.getMomentum().y);
- m_barShowerU_mc_pz.push_back(iter.first.getMomentum().z);
- m_barShowerU_mc_weight.push_back(iter.second);
+ for(int il=0; il<m_local_max.size(); il++){
+ m_localMaxU_tag.push_back( il );
+ m_localMaxU_x.push_back( m_local_max[il]->getPos().x() );
+ m_localMaxU_y.push_back( m_local_max[il]->getPos().y() );
+ m_localMaxU_z.push_back( m_local_max[il]->getPos().z() );
+ m_localMaxU_E.push_back( m_local_max[il]->getEnergy() );
+
+ auto truthMap = m_local_max[il]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_localMaxU_mc_tag.push_back(il);
+ m_localMaxU_mc_pdg.push_back(iter.first.getPDG());
+ m_localMaxU_mc_px.push_back(iter.first.getMomentum().x);
+ m_localMaxU_mc_py.push_back(iter.first.getMomentum().y);
+ m_localMaxU_mc_pz.push_back(iter.first.getMomentum().z);
+ m_localMaxU_mc_weight.push_back(iter.second);
+ }
}
- }
- m_halfclusters.clear();
- m_local_max.clear();
- for(int i=0; i<m_DataCol.map_HalfCluster["ESHalfClusterV"].size(); i++)
- m_halfclusters.push_back( m_DataCol.map_HalfCluster["ESHalfClusterV"][i].get() );
-
- for(int ic=0;ic<m_halfclusters.size(); ic++){
- //std::vector<const Calo1DCluster*> tmp_shower = m_halfclusters[ic]->getLocalMaxCol("AllLocalMax");
- std::vector<const CaloHalfCluster*> m_axis = m_halfclusters[ic]->getHalfClusterCol("MergedAxis");
- if(m_axis.size()>0){
- std::vector<const Calo1DCluster*> tmp_shower = m_axis[0]->getCluster();
+ m_halfclusters.clear();
+ m_local_max.clear();
+ for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColV"].size(); i++)
+ m_halfclusters.push_back( m_DataCol.map_HalfCluster["HalfClusterColV"][i].get() );
+
+ for(int ic=0;ic<m_halfclusters.size(); ic++){
+ std::vector<const Calo1DCluster*> tmp_shower = m_halfclusters[ic]->getLocalMaxCol("AllLocalMax");
m_local_max.insert(m_local_max.end(), tmp_shower.begin(), tmp_shower.end());
}
- }
- for(int il=0; il<m_local_max.size(); il++){
- m_barShowerV_tag.push_back( il );
- m_barShowerV_x.push_back( m_local_max[il]->getPos().x() );
- m_barShowerV_y.push_back( m_local_max[il]->getPos().y() );
- m_barShowerV_z.push_back( m_local_max[il]->getPos().z() );
- m_barShowerV_E.push_back( m_local_max[il]->getEnergy() );
-
- auto truthMap = m_local_max[il]->getLinkedMCP();
- for(auto iter: truthMap){
- m_barShowerV_mc_tag.push_back(il);
- m_barShowerV_mc_pdg.push_back(iter.first.getPDG());
- m_barShowerV_mc_px.push_back(iter.first.getMomentum().x);
- m_barShowerV_mc_py.push_back(iter.first.getMomentum().y);
- m_barShowerV_mc_pz.push_back(iter.first.getMomentum().z);
- m_barShowerV_mc_weight.push_back(iter.second);
- }
- }
- t_Layers->Fill();
-
- std::vector<const Cyber::CaloHalfCluster*> m_halfclusterV; m_halfclusterV.clear();
- std::vector<const Cyber::CaloHalfCluster*> m_halfclusterU; m_halfclusterU.clear();
- for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColU"].size(); i++){
- m_halfclusterU.push_back( m_DataCol.map_HalfCluster["HalfClusterColU"][i].get() );
- }
- for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColV"].size(); i++){
- m_halfclusterV.push_back( m_DataCol.map_HalfCluster["HalfClusterColV"][i].get() );
- }
- // Hough
- ClearHough();
- int houghU_index=0;
- int houghV_index=0;
- for(int i=0; i<m_halfclusterU.size(); i++){ // loop half cluster U
- std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisU = m_halfclusterU[i]->getHalfClusterCol("HoughAxis");
- for(int ita=0; ita<m_mergedaxisU.size(); ita++){ // loop axis U
- // General information of the axis
- m_houghU_tag.push_back(houghU_index);
- m_houghU_type.push_back(m_mergedaxisU[ita]->getType());
- m_houghU_x.push_back(m_mergedaxisU[ita]->getPos().x());
- m_houghU_y.push_back(m_mergedaxisU[ita]->getPos().y());
- m_houghU_z.push_back(m_mergedaxisU[ita]->getPos().z());
- m_houghU_E.push_back(m_mergedaxisU[ita]->getEnergy());
-
- // MC truth information of the Axis
- auto truthMap = m_mergedaxisU[ita]->getLinkedMCP();
+ for(int il=0; il<m_local_max.size(); il++){
+ m_localMaxV_tag.push_back( il );
+ m_localMaxV_x.push_back( m_local_max[il]->getPos().x() );
+ m_localMaxV_y.push_back( m_local_max[il]->getPos().y() );
+ m_localMaxV_z.push_back( m_local_max[il]->getPos().z() );
+ m_localMaxV_E.push_back( m_local_max[il]->getEnergy() );
+
+ auto truthMap = m_local_max[il]->getLinkedMCP();
for(auto iter: truthMap){
- m_houghU_truth_tag.push_back(houghU_index);
- m_houghU_truth_MC_px.push_back(iter.first.getMomentum().x);
- m_houghU_truth_MC_py.push_back(iter.first.getMomentum().y);
- m_houghU_truth_MC_pz.push_back(iter.first.getMomentum().z);
- m_houghU_truth_MC_E.push_back(iter.first.getEnergy());
- m_houghU_truth_MC_weight.push_back(iter.second);
+ m_localMaxV_mc_tag.push_back(il);
+ m_localMaxV_mc_pdg.push_back(iter.first.getPDG());
+ m_localMaxV_mc_px.push_back(iter.first.getMomentum().x);
+ m_localMaxV_mc_py.push_back(iter.first.getMomentum().y);
+ m_localMaxV_mc_pz.push_back(iter.first.getMomentum().z);
+ m_localMaxV_mc_weight.push_back(iter.second);
}
-
- // Hits on axis
- for(int ilm=0; ilm<m_mergedaxisU[ita]->getCluster().size(); ilm++){ // loop local max
- m_houghU_hit_tag.push_back(houghU_index);
- m_houghU_hit_x.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().x() );
- m_houghU_hit_y.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().y() );
- m_houghU_hit_z.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().z() );
- m_houghU_hit_E.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getEnergy() );
- }
-
- houghU_index++;
}
- }
- for(int i=0; i<m_halfclusterV.size(); i++){ // loop half cluster V
- std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisV = m_halfclusterV[i]->getHalfClusterCol("HoughAxis");
- for(int ita=0; ita<m_mergedaxisV.size(); ita++){ // loop axis V
- // General information of the axis
- m_houghV_tag.push_back(houghV_index);
- m_houghV_type.push_back(m_mergedaxisV[ita]->getType());
- m_houghV_x.push_back(m_mergedaxisV[ita]->getPos().x());
- m_houghV_y.push_back(m_mergedaxisV[ita]->getPos().y());
- m_houghV_z.push_back(m_mergedaxisV[ita]->getPos().z());
- m_houghV_E.push_back(m_mergedaxisV[ita]->getEnergy());
- m_houghV_alpha.push_back(m_mergedaxisV[ita]->getHoughAlpha());
- m_houghV_rho.push_back(m_mergedaxisV[ita]->getHoughRho());
-
- // MC truth information of the Axis
- auto truthMap = m_mergedaxisV[ita]->getLinkedMCP();
- for(auto iter: truthMap){
- m_houghV_truth_tag.push_back(houghV_index);
- m_houghV_truth_MC_px.push_back(iter.first.getMomentum().x);
- m_houghV_truth_MC_py.push_back(iter.first.getMomentum().y);
- m_houghV_truth_MC_pz.push_back(iter.first.getMomentum().z);
- m_houghV_truth_MC_E.push_back(iter.first.getEnergy());
- m_houghV_truth_MC_weight.push_back(iter.second);
+ t_LocalMax->Fill();
+
+ //Save 1DCluster
+ ClearLayer();
+ m_halfclusters.clear();
+ for(int i=0; i<m_DataCol.map_HalfCluster["ESHalfClusterU"].size(); i++)
+ m_halfclusters.push_back( m_DataCol.map_HalfCluster["ESHalfClusterU"][i].get() );
+
+ m_local_max.clear();
+ for(int ic=0;ic<m_halfclusters.size(); ic++){
+ //std::vector<const Calo1DCluster*> tmp_shower = m_halfclusters[ic]->getLocalMaxCol("AllLocalMax");
+ std::vector<const CaloHalfCluster*> m_axis = m_halfclusters[ic]->getHalfClusterCol("MergedAxis");
+ if(m_axis.size()>0){
+ std::vector<const Calo1DCluster*> tmp_shower = m_axis[0]->getCluster();
+ m_local_max.insert(m_local_max.end(), tmp_shower.begin(), tmp_shower.end());
}
-
- // Hits on axis
- for(int ilm=0; ilm<m_mergedaxisV[ita]->getCluster().size(); ilm++){
- m_houghV_hit_tag.push_back(houghV_index);
- m_houghV_hit_x.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().x());
- m_houghV_hit_y.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().y());
- m_houghV_hit_z.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().z());
- m_houghV_hit_E.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getEnergy());
- }
-
- houghV_index++;
}
- }
- t_Hough->Fill();
- // Cone
- ClearCone();
- int coneU_index=0;
- int coneV_index=0;
- for(int i=0; i<m_halfclusterU.size(); i++){ // loop half cluster U
- std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisU = m_halfclusterU[i]->getHalfClusterCol("ConeAxis");
- for(int ita=0; ita<m_mergedaxisU.size(); ita++){ // loop axis U
- // General information of the axis
- m_coneU_tag.push_back(coneU_index);
- m_coneU_type.push_back(m_mergedaxisU[ita]->getType());
- m_coneU_x.push_back(m_mergedaxisU[ita]->getPos().x());
- m_coneU_y.push_back(m_mergedaxisU[ita]->getPos().y());
- m_coneU_z.push_back(m_mergedaxisU[ita]->getPos().z());
- m_coneU_E.push_back(m_mergedaxisU[ita]->getEnergy());
-
- // MC truth information of the Axis
- auto truthMap = m_mergedaxisU[ita]->getLinkedMCP();
+ for(int il=0; il<m_local_max.size(); il++){
+ m_barShowerU_tag.push_back( il );
+ m_barShowerU_x.push_back( m_local_max[il]->getPos().x() );
+ m_barShowerU_y.push_back( m_local_max[il]->getPos().y() );
+ m_barShowerU_z.push_back( m_local_max[il]->getPos().z() );
+ m_barShowerU_E.push_back( m_local_max[il]->getEnergy() );
+
+ auto truthMap = m_local_max[il]->getLinkedMCP();
for(auto iter: truthMap){
- m_coneU_truth_tag.push_back(coneU_index);
- m_coneU_truth_MC_px.push_back(iter.first.getMomentum().x);
- m_coneU_truth_MC_py.push_back(iter.first.getMomentum().y);
- m_coneU_truth_MC_pz.push_back(iter.first.getMomentum().z);
- m_coneU_truth_MC_E.push_back(iter.first.getEnergy());
- m_coneU_truth_MC_weight.push_back(iter.second);
+ m_barShowerU_mc_tag.push_back(il);
+ m_barShowerU_mc_pdg.push_back(iter.first.getPDG());
+ m_barShowerU_mc_px.push_back(iter.first.getMomentum().x);
+ m_barShowerU_mc_py.push_back(iter.first.getMomentum().y);
+ m_barShowerU_mc_pz.push_back(iter.first.getMomentum().z);
+ m_barShowerU_mc_weight.push_back(iter.second);
}
-
- // Hits on axis
- for(int ilm=0; ilm<m_mergedaxisU[ita]->getCluster().size(); ilm++){ // loop local max
- m_coneU_hit_tag.push_back(coneU_index);
- m_coneU_hit_x.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().x() );
- m_coneU_hit_y.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().y() );
- m_coneU_hit_z.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().z() );
- m_coneU_hit_E.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getEnergy() );
- }
-
- coneU_index++;
}
- }
- for(int i=0; i<m_halfclusterV.size(); i++){ // loop half cluster V
- std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisV = m_halfclusterV[i]->getHalfClusterCol("ConeAxis");
- for(int ita=0; ita<m_mergedaxisV.size(); ita++){ // loop axis V
- // General information of the axis
- m_coneV_tag.push_back(coneV_index);
- m_coneV_type.push_back(m_mergedaxisV[ita]->getType());
- m_coneV_x.push_back(m_mergedaxisV[ita]->getPos().x());
- m_coneV_y.push_back(m_mergedaxisV[ita]->getPos().y());
- m_coneV_z.push_back(m_mergedaxisV[ita]->getPos().z());
- m_coneV_E.push_back(m_mergedaxisV[ita]->getEnergy());
-
- // MC truth information of the Axis
- auto truthMap = m_mergedaxisV[ita]->getLinkedMCP();
- for(auto iter: truthMap){
- m_coneV_truth_tag.push_back(coneV_index);
- m_coneV_truth_MC_px.push_back(iter.first.getMomentum().x);
- m_coneV_truth_MC_py.push_back(iter.first.getMomentum().y);
- m_coneV_truth_MC_pz.push_back(iter.first.getMomentum().z);
- m_coneV_truth_MC_E.push_back(iter.first.getEnergy());
- m_coneV_truth_MC_weight.push_back(iter.second);
+ m_halfclusters.clear();
+ m_local_max.clear();
+ for(int i=0; i<m_DataCol.map_HalfCluster["ESHalfClusterV"].size(); i++)
+ m_halfclusters.push_back( m_DataCol.map_HalfCluster["ESHalfClusterV"][i].get() );
+
+ for(int ic=0;ic<m_halfclusters.size(); ic++){
+ //std::vector<const Calo1DCluster*> tmp_shower = m_halfclusters[ic]->getLocalMaxCol("AllLocalMax");
+ std::vector<const CaloHalfCluster*> m_axis = m_halfclusters[ic]->getHalfClusterCol("MergedAxis");
+ if(m_axis.size()>0){
+ std::vector<const Calo1DCluster*> tmp_shower = m_axis[0]->getCluster();
+ m_local_max.insert(m_local_max.end(), tmp_shower.begin(), tmp_shower.end());
}
-
- // Hits on axis
- for(int ilm=0; ilm<m_mergedaxisV[ita]->getCluster().size(); ilm++){
- m_coneV_hit_tag.push_back(coneV_index);
- m_coneV_hit_x.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().x());
- m_coneV_hit_y.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().y());
- m_coneV_hit_z.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().z());
- m_coneV_hit_E.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getEnergy());
- }
-
- coneV_index++;
}
- }
- t_Cone->Fill();
- // Track axis
- ClearTrackAxis();
- int trackU_index=0;
- int trackV_index=0;
- for(int i=0; i<m_halfclusterU.size(); i++){ // loop half cluster U
- std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisU = m_halfclusterU[i]->getHalfClusterCol("TrackAxis");
- for(int ita=0; ita<m_mergedaxisU.size(); ita++){ // loop axis U
- // General information of the axis
- m_trackU_tag.push_back(trackU_index);
- m_trackU_type.push_back(m_mergedaxisU[ita]->getType());
- m_trackU_x.push_back(m_mergedaxisU[ita]->getPos().x());
- m_trackU_y.push_back(m_mergedaxisU[ita]->getPos().y());
- m_trackU_z.push_back(m_mergedaxisU[ita]->getPos().z());
- m_trackU_E.push_back(m_mergedaxisU[ita]->getEnergy());
-
- // MC truth information of the Axis
- auto truthMap = m_mergedaxisU[ita]->getLinkedMCP();
+ for(int il=0; il<m_local_max.size(); il++){
+ m_barShowerV_tag.push_back( il );
+ m_barShowerV_x.push_back( m_local_max[il]->getPos().x() );
+ m_barShowerV_y.push_back( m_local_max[il]->getPos().y() );
+ m_barShowerV_z.push_back( m_local_max[il]->getPos().z() );
+ m_barShowerV_E.push_back( m_local_max[il]->getEnergy() );
+
+ auto truthMap = m_local_max[il]->getLinkedMCP();
for(auto iter: truthMap){
- m_trackU_truth_tag.push_back(trackU_index);
- m_trackU_truth_MC_px.push_back(iter.first.getMomentum().x);
- m_trackU_truth_MC_py.push_back(iter.first.getMomentum().y);
- m_trackU_truth_MC_pz.push_back(iter.first.getMomentum().z);
- m_trackU_truth_MC_E.push_back(iter.first.getEnergy());
- m_trackU_truth_MC_weight.push_back(iter.second);
+ m_barShowerV_mc_tag.push_back(il);
+ m_barShowerV_mc_pdg.push_back(iter.first.getPDG());
+ m_barShowerV_mc_px.push_back(iter.first.getMomentum().x);
+ m_barShowerV_mc_py.push_back(iter.first.getMomentum().y);
+ m_barShowerV_mc_pz.push_back(iter.first.getMomentum().z);
+ m_barShowerV_mc_weight.push_back(iter.second);
}
-
- // Hits on axis
- for(int ilm=0; ilm<m_mergedaxisU[ita]->getCluster().size(); ilm++){ // loop local max
- m_trackU_hit_tag.push_back(trackU_index);
- m_trackU_hit_x.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().x() );
- m_trackU_hit_y.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().y() );
- m_trackU_hit_z.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().z() );
- m_trackU_hit_E.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getEnergy() );
- }
-
- trackU_index++;
}
- }
- for(int i=0; i<m_halfclusterV.size(); i++){ // loop half cluster V
- std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisV = m_halfclusterV[i]->getHalfClusterCol("TrackAxis");
- for(int ita=0; ita<m_mergedaxisV.size(); ita++){ // loop axis V
- // General information of the axis
- m_trackV_tag.push_back(trackV_index);
- m_trackV_type.push_back(m_mergedaxisV[ita]->getType());
- m_trackV_x.push_back(m_mergedaxisV[ita]->getPos().x());
- m_trackV_y.push_back(m_mergedaxisV[ita]->getPos().y());
- m_trackV_z.push_back(m_mergedaxisV[ita]->getPos().z());
- m_trackV_E.push_back(m_mergedaxisV[ita]->getEnergy());
-
- // MC truth information of the Axis
- auto truthMap = m_mergedaxisV[ita]->getLinkedMCP();
- for(auto iter: truthMap){
- m_trackV_truth_tag.push_back(trackV_index);
- m_trackV_truth_MC_px.push_back(iter.first.getMomentum().x);
- m_trackV_truth_MC_py.push_back(iter.first.getMomentum().y);
- m_trackV_truth_MC_pz.push_back(iter.first.getMomentum().z);
- m_trackV_truth_MC_E.push_back(iter.first.getEnergy());
- m_trackV_truth_MC_weight.push_back(iter.second);
- }
-
- // Hits on axis
- for(int ilm=0; ilm<m_mergedaxisV[ita]->getCluster().size(); ilm++){
- m_trackV_hit_tag.push_back(trackV_index);
- m_trackV_hit_x.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().x());
- m_trackV_hit_y.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().y());
- m_trackV_hit_z.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().z());
- m_trackV_hit_E.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getEnergy());
- }
-
- trackV_index++;
+ t_Layers->Fill();
+
+ std::vector<const Cyber::CaloHalfCluster*> m_halfclusterV; m_halfclusterV.clear();
+ std::vector<const Cyber::CaloHalfCluster*> m_halfclusterU; m_halfclusterU.clear();
+ for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColU"].size(); i++){
+ m_halfclusterU.push_back( m_DataCol.map_HalfCluster["HalfClusterColU"][i].get() );
}
- }
- t_TrackAxis->Fill();
- //Axis
- ClearAxis();
- int axisU_index=0;
- int axisV_index=0;
- for(int i=0; i<m_halfclusterU.size(); i++){ // loop half cluster U
- std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisU = m_halfclusterU[i]->getHalfClusterCol("MergedAxis");
- for(int ita=0; ita<m_mergedaxisU.size(); ita++){ // loop axis U
- // General information of the axis
- m_axisU_tag.push_back(axisU_index);
- m_axisU_type.push_back(m_mergedaxisU[ita]->getType());
- m_axisU_x.push_back(m_mergedaxisU[ita]->getPos().x());
- m_axisU_y.push_back(m_mergedaxisU[ita]->getPos().y());
- m_axisU_z.push_back(m_mergedaxisU[ita]->getPos().z());
- m_axisU_E.push_back(m_mergedaxisU[ita]->getEnergy());
-
- // MC truth information of the Axis
- auto truthMap = m_mergedaxisU[ita]->getLinkedMCP();
- for(auto iter: truthMap){
- m_axisU_truth_tag.push_back(axisU_index);
- m_axisU_truth_MC_px.push_back(iter.first.getMomentum().x);
- m_axisU_truth_MC_py.push_back(iter.first.getMomentum().y);
- m_axisU_truth_MC_pz.push_back(iter.first.getMomentum().z);
- m_axisU_truth_MC_E.push_back(iter.first.getEnergy());
- m_axisU_truth_MC_weight.push_back(iter.second);
- }
- // Hits on axis
- for(int ilm=0; ilm<m_mergedaxisU[ita]->getCluster().size(); ilm++){ // loop local max
- m_axisU_hit_tag.push_back(axisU_index);
- m_axisU_hit_x.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().x() );
- m_axisU_hit_y.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().y() );
- m_axisU_hit_z.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().z() );
- m_axisU_hit_E.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getEnergy() );
- }
-
- axisU_index++;
+ for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColV"].size(); i++){
+ m_halfclusterV.push_back( m_DataCol.map_HalfCluster["HalfClusterColV"][i].get() );
}
- }
- for(int i=0; i<m_halfclusterV.size(); i++){ // loop half cluster V
- std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisV = m_halfclusterV[i]->getHalfClusterCol("MergedAxis");
- for(int ita=0; ita<m_mergedaxisV.size(); ita++){ // loop axis V
- // General information of the axis
- m_axisV_tag.push_back(axisV_index);
- m_axisV_type.push_back(m_mergedaxisV[ita]->getType());
- m_axisV_x.push_back(m_mergedaxisV[ita]->getPos().x());
- m_axisV_y.push_back(m_mergedaxisV[ita]->getPos().y());
- m_axisV_z.push_back(m_mergedaxisV[ita]->getPos().z());
- m_axisV_E.push_back(m_mergedaxisV[ita]->getEnergy());
-
- // MC truth information of the Axis
- auto truthMap = m_mergedaxisV[ita]->getLinkedMCP();
- for(auto iter: truthMap){
- m_axisV_truth_tag.push_back(axisV_index);
- m_axisV_truth_MC_px.push_back(iter.first.getMomentum().x);
- m_axisV_truth_MC_py.push_back(iter.first.getMomentum().y);
- m_axisV_truth_MC_pz.push_back(iter.first.getMomentum().z);
- m_axisV_truth_MC_E.push_back(iter.first.getEnergy());
- m_axisV_truth_MC_weight.push_back(iter.second);
- }
- // Hits on axis
- for(int ilm=0; ilm<m_mergedaxisV[ita]->getCluster().size(); ilm++){
- m_axisV_hit_tag.push_back(axisV_index);
- m_axisV_hit_x.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().x());
- m_axisV_hit_y.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().y());
- m_axisV_hit_z.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().z());
- m_axisV_hit_E.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getEnergy());
+ // Hough
+ ClearHough();
+ int houghU_index=0;
+ int houghV_index=0;
+ for(int i=0; i<m_halfclusterU.size(); i++){ // loop half cluster U
+ std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisU = m_halfclusterU[i]->getHalfClusterCol("HoughAxis");
+ for(int ita=0; ita<m_mergedaxisU.size(); ita++){ // loop axis U
+ // General information of the axis
+ m_houghU_tag.push_back(houghU_index);
+ m_houghU_type.push_back(m_mergedaxisU[ita]->getType());
+ m_houghU_x.push_back(m_mergedaxisU[ita]->getPos().x());
+ m_houghU_y.push_back(m_mergedaxisU[ita]->getPos().y());
+ m_houghU_z.push_back(m_mergedaxisU[ita]->getPos().z());
+ m_houghU_E.push_back(m_mergedaxisU[ita]->getEnergy());
+
+ // MC truth information of the Axis
+ auto truthMap = m_mergedaxisU[ita]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_houghU_truth_tag.push_back(houghU_index);
+ m_houghU_truth_MC_px.push_back(iter.first.getMomentum().x);
+ m_houghU_truth_MC_py.push_back(iter.first.getMomentum().y);
+ m_houghU_truth_MC_pz.push_back(iter.first.getMomentum().z);
+ m_houghU_truth_MC_E.push_back(iter.first.getEnergy());
+ m_houghU_truth_MC_weight.push_back(iter.second);
+ }
+
+ // Hits on axis
+ for(int ilm=0; ilm<m_mergedaxisU[ita]->getCluster().size(); ilm++){ // loop local max
+ m_houghU_hit_tag.push_back(houghU_index);
+ m_houghU_hit_x.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().x() );
+ m_houghU_hit_y.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().y() );
+ m_houghU_hit_z.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().z() );
+ m_houghU_hit_E.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getEnergy() );
+ }
+
+ houghU_index++;
}
-
- axisV_index++;
}
- }
- m_halfclusterU.clear();
- m_halfclusterV.clear();
- for(int i=0; i<m_DataCol.map_HalfCluster["emptyHalfClusterU"].size(); i++){
- m_halfclusterU.push_back( m_DataCol.map_HalfCluster["emptyHalfClusterU"][i].get() );
- }
- for(int i=0; i<m_DataCol.map_HalfCluster["emptyHalfClusterV"].size(); i++){
- m_halfclusterV.push_back( m_DataCol.map_HalfCluster["emptyHalfClusterV"][i].get() );
- }
- for(int i=0; i<m_halfclusterU.size(); i++){
- m_emptyAxisU_tag.push_back(m_halfclusterU[i]->getType());
- m_emptyAxisU_x.push_back(m_halfclusterU[i]->getPos().x());
- m_emptyAxisU_y.push_back(m_halfclusterU[i]->getPos().y());
- m_emptyAxisU_z.push_back(m_halfclusterU[i]->getPos().z());
- m_emptyAxisU_E.push_back(m_halfclusterU[i]->getEnergy());
- }
- for(int i=0; i<m_halfclusterV.size(); i++){
- m_emptyAxisV_tag.push_back(m_halfclusterV[i]->getType());
- m_emptyAxisV_x.push_back(m_halfclusterV[i]->getPos().x());
- m_emptyAxisV_y.push_back(m_halfclusterV[i]->getPos().y());
- m_emptyAxisV_z.push_back(m_halfclusterV[i]->getPos().z());
- m_emptyAxisV_E.push_back(m_halfclusterV[i]->getEnergy());
- }
-
- t_Axis->Fill();
-
-
- //Half cluster
- ClearHalfCluster();
- m_halfclusterV.clear();
- m_halfclusterU.clear();
- m_totE_HFClusV = 0;
- m_totE_HFClusU = 0;
- //for(int i=0; i<m_DataCol.map_HalfCluster["ESHalfClusterU"].size(); i++){
- // m_halfclusterU.push_back( m_DataCol.map_HalfCluster["ESHalfClusterU"][i]->getHalfClusterCol("MergedAxis")[0] );
- //}
- //for(int i=0; i<m_DataCol.map_HalfCluster["ESHalfClusterV"].size(); i++){
- // m_halfclusterV.push_back( m_DataCol.map_HalfCluster["ESHalfClusterV"][i]->getHalfClusterCol("MergedAxis")[0] );
- //}
- for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColU"].size(); i++){
- m_halfclusterU.push_back( m_DataCol.map_HalfCluster["HalfClusterColU"][i].get() );
- }
- for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColV"].size(); i++){
- m_halfclusterV.push_back( m_DataCol.map_HalfCluster["HalfClusterColV"][i].get() );
- }
- for(int i=0; i<m_halfclusterV.size(); i++){
- m_HalfClusterV_x.push_back(m_halfclusterV[i]->getPos().x());
- m_HalfClusterV_y.push_back(m_halfclusterV[i]->getPos().y());
- m_HalfClusterV_z.push_back(m_halfclusterV[i]->getPos().z());
- m_HalfClusterV_E.push_back(m_halfclusterV[i]->getEnergy());
- m_HalfClusterV_tag.push_back(i);
- m_HalfClusterV_type.push_back(m_halfclusterV[i]->getType());
- m_totE_HFClusV += m_halfclusterV[i]->getEnergy();
-
- // MC truth information of the HFCluster
- auto truthMap = m_halfclusterV[i]->getLinkedMCP();
- for(auto iter: truthMap){
- m_HalfClusterV_truth_tag.push_back(i);
- m_HalfClusterV_truthMC_px.push_back(iter.first.getMomentum().x);
- m_HalfClusterV_truthMC_py.push_back(iter.first.getMomentum().y);
- m_HalfClusterV_truthMC_pz.push_back(iter.first.getMomentum().z);
- m_HalfClusterV_truthMC_E.push_back(iter.first.getEnergy());
- m_HalfClusterV_truthMC_weight.push_back(iter.second);
- }
-
- // 1DClusters (hits)
- for(int ilm=0; ilm<m_halfclusterV[i]->getCluster().size(); ilm++){
- m_HalfClusterV_hit_tag.push_back(i);
- m_HalfClusterV_hit_x.push_back( m_halfclusterV[i]->getCluster()[ilm]->getPos().x() );
- m_HalfClusterV_hit_y.push_back( m_halfclusterV[i]->getCluster()[ilm]->getPos().y() );
- m_HalfClusterV_hit_z.push_back( m_halfclusterV[i]->getCluster()[ilm]->getPos().z() );
- m_HalfClusterV_hit_E.push_back( m_halfclusterV[i]->getCluster()[ilm]->getEnergy() );
+ for(int i=0; i<m_halfclusterV.size(); i++){ // loop half cluster V
+ std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisV = m_halfclusterV[i]->getHalfClusterCol("HoughAxis");
+ for(int ita=0; ita<m_mergedaxisV.size(); ita++){ // loop axis V
+ // General information of the axis
+ m_houghV_tag.push_back(houghV_index);
+ m_houghV_type.push_back(m_mergedaxisV[ita]->getType());
+ m_houghV_x.push_back(m_mergedaxisV[ita]->getPos().x());
+ m_houghV_y.push_back(m_mergedaxisV[ita]->getPos().y());
+ m_houghV_z.push_back(m_mergedaxisV[ita]->getPos().z());
+ m_houghV_E.push_back(m_mergedaxisV[ita]->getEnergy());
+ m_houghV_alpha.push_back(m_mergedaxisV[ita]->getHoughAlpha());
+ m_houghV_rho.push_back(m_mergedaxisV[ita]->getHoughRho());
+
+ // MC truth information of the Axis
+ auto truthMap = m_mergedaxisV[ita]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_houghV_truth_tag.push_back(houghV_index);
+ m_houghV_truth_MC_px.push_back(iter.first.getMomentum().x);
+ m_houghV_truth_MC_py.push_back(iter.first.getMomentum().y);
+ m_houghV_truth_MC_pz.push_back(iter.first.getMomentum().z);
+ m_houghV_truth_MC_E.push_back(iter.first.getEnergy());
+ m_houghV_truth_MC_weight.push_back(iter.second);
+ }
+
+ // Hits on axis
+ for(int ilm=0; ilm<m_mergedaxisV[ita]->getCluster().size(); ilm++){
+ m_houghV_hit_tag.push_back(houghV_index);
+ m_houghV_hit_x.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().x());
+ m_houghV_hit_y.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().y());
+ m_houghV_hit_z.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().z());
+ m_houghV_hit_E.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getEnergy());
+ }
+
+ houghV_index++;
}
- }
- for(int i=0; i<m_halfclusterU.size(); i++){
- m_HalfClusterU_x.push_back(m_halfclusterU[i]->getPos().x());
- m_HalfClusterU_y.push_back(m_halfclusterU[i]->getPos().y());
- m_HalfClusterU_z.push_back(m_halfclusterU[i]->getPos().z());
- m_HalfClusterU_E.push_back(m_halfclusterU[i]->getEnergy());
- m_HalfClusterU_tag.push_back(i);
- m_HalfClusterU_type.push_back(m_halfclusterU[i]->getType());
- m_totE_HFClusU += m_halfclusterU[i]->getEnergy();
-
- // MC truth information of the HFCluster
- auto truthMap = m_halfclusterU[i]->getLinkedMCP();
- for(auto iter: truthMap){
- m_HalfClusterU_truth_tag.push_back(i);
- m_HalfClusterU_truthMC_px.push_back(iter.first.getMomentum().x);
- m_HalfClusterU_truthMC_py.push_back(iter.first.getMomentum().y);
- m_HalfClusterU_truthMC_pz.push_back(iter.first.getMomentum().z);
- m_HalfClusterU_truthMC_E.push_back(iter.first.getEnergy());
- m_HalfClusterU_truthMC_weight.push_back(iter.second);
+ }
+ t_Hough->Fill();
+ // Cone
+ ClearCone();
+ int coneU_index=0;
+ int coneV_index=0;
+ for(int i=0; i<m_halfclusterU.size(); i++){ // loop half cluster U
+ std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisU = m_halfclusterU[i]->getHalfClusterCol("ConeAxis");
+ for(int ita=0; ita<m_mergedaxisU.size(); ita++){ // loop axis U
+ // General information of the axis
+ m_coneU_tag.push_back(coneU_index);
+ m_coneU_type.push_back(m_mergedaxisU[ita]->getType());
+ m_coneU_x.push_back(m_mergedaxisU[ita]->getPos().x());
+ m_coneU_y.push_back(m_mergedaxisU[ita]->getPos().y());
+ m_coneU_z.push_back(m_mergedaxisU[ita]->getPos().z());
+ m_coneU_E.push_back(m_mergedaxisU[ita]->getEnergy());
+
+ // MC truth information of the Axis
+ auto truthMap = m_mergedaxisU[ita]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_coneU_truth_tag.push_back(coneU_index);
+ m_coneU_truth_MC_px.push_back(iter.first.getMomentum().x);
+ m_coneU_truth_MC_py.push_back(iter.first.getMomentum().y);
+ m_coneU_truth_MC_pz.push_back(iter.first.getMomentum().z);
+ m_coneU_truth_MC_E.push_back(iter.first.getEnergy());
+ m_coneU_truth_MC_weight.push_back(iter.second);
+ }
+
+ // Hits on axis
+ for(int ilm=0; ilm<m_mergedaxisU[ita]->getCluster().size(); ilm++){ // loop local max
+ m_coneU_hit_tag.push_back(coneU_index);
+ m_coneU_hit_x.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().x() );
+ m_coneU_hit_y.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().y() );
+ m_coneU_hit_z.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().z() );
+ m_coneU_hit_E.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getEnergy() );
+ }
+
+ coneU_index++;
}
-
- // 1DClusters (hits)
- for(int ilm=0; ilm<m_halfclusterU[i]->getCluster().size(); ilm++){
- m_HalfClusterU_hit_tag.push_back(i);
- m_HalfClusterU_hit_x.push_back( m_halfclusterU[i]->getCluster()[ilm]->getPos().x() );
- m_HalfClusterU_hit_y.push_back( m_halfclusterU[i]->getCluster()[ilm]->getPos().y() );
- m_HalfClusterU_hit_z.push_back( m_halfclusterU[i]->getCluster()[ilm]->getPos().z() );
- m_HalfClusterU_hit_E.push_back( m_halfclusterU[i]->getCluster()[ilm]->getEnergy() );
+ }
+ for(int i=0; i<m_halfclusterV.size(); i++){ // loop half cluster V
+ std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisV = m_halfclusterV[i]->getHalfClusterCol("ConeAxis");
+ for(int ita=0; ita<m_mergedaxisV.size(); ita++){ // loop axis V
+ // General information of the axis
+ m_coneV_tag.push_back(coneV_index);
+ m_coneV_type.push_back(m_mergedaxisV[ita]->getType());
+ m_coneV_x.push_back(m_mergedaxisV[ita]->getPos().x());
+ m_coneV_y.push_back(m_mergedaxisV[ita]->getPos().y());
+ m_coneV_z.push_back(m_mergedaxisV[ita]->getPos().z());
+ m_coneV_E.push_back(m_mergedaxisV[ita]->getEnergy());
+
+ // MC truth information of the Axis
+ auto truthMap = m_mergedaxisV[ita]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_coneV_truth_tag.push_back(coneV_index);
+ m_coneV_truth_MC_px.push_back(iter.first.getMomentum().x);
+ m_coneV_truth_MC_py.push_back(iter.first.getMomentum().y);
+ m_coneV_truth_MC_pz.push_back(iter.first.getMomentum().z);
+ m_coneV_truth_MC_E.push_back(iter.first.getEnergy());
+ m_coneV_truth_MC_weight.push_back(iter.second);
+ }
+
+ // Hits on axis
+ for(int ilm=0; ilm<m_mergedaxisV[ita]->getCluster().size(); ilm++){
+ m_coneV_hit_tag.push_back(coneV_index);
+ m_coneV_hit_x.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().x());
+ m_coneV_hit_y.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().y());
+ m_coneV_hit_z.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().z());
+ m_coneV_hit_E.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getEnergy());
+ }
+
+ coneV_index++;
}
- }
- t_HalfCluster->Fill();
-
-
- //Tower
- ClearTower();
- std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_tower = m_DataCol.map_CaloCluster["ESTower"];
- for(int it=0; it<m_tower.size(); it++){
- ClearTower();
- towerID[0] = m_tower[it]->getTowerID()[0][0];
- towerID[1] = m_tower[it]->getTowerID()[0][1];
- towerID[2] = m_tower[it]->getTowerID()[0][2];
-
- std::vector<const CaloHalfCluster*> m_HFClusU = m_tower[it]->getHalfClusterUCol("ESHalfClusterU");
- std::vector<const CaloHalfCluster*> m_HFClusV = m_tower[it]->getHalfClusterVCol("ESHalfClusterV");
-
- m_NclusU = m_HFClusU.size();
- m_NclusV = m_HFClusV.size();
- m_totEn = m_tower[it]->getEnergy();
- m_totEn_U = 0.;
- m_totEn_V = 0.;
- for(int ic=0; ic<m_NclusU; ic++){
- m_HalfClusterU_tag.push_back(it);
- m_HalfClusterU_x.push_back(m_HFClusU[ic]->getPos().x());
- m_HalfClusterU_y.push_back(m_HFClusU[ic]->getPos().y());
- m_HalfClusterU_z.push_back(m_HFClusU[ic]->getPos().z());
- m_HalfClusterU_E.push_back(m_HFClusU[ic]->getEnergy());
- m_HalfClusterU_type.push_back(m_HFClusU[ic]->getType());
- m_HalfClusterU_nTrk.push_back(m_HFClusU[ic]->getAssociatedTracks().size());
- m_totEn_U += m_HFClusU[ic]->getEnergy();
}
-
- for(int ic=0; ic<m_NclusV; ic++){
- m_HalfClusterV_tag.push_back(it);
- m_HalfClusterV_x.push_back(m_HFClusV[ic]->getPos().x());
- m_HalfClusterV_y.push_back(m_HFClusV[ic]->getPos().y());
- m_HalfClusterV_z.push_back(m_HFClusV[ic]->getPos().z());
- m_HalfClusterV_E.push_back(m_HFClusV[ic]->getEnergy());
- m_HalfClusterV_type.push_back(m_HFClusV[ic]->getType());
- m_HalfClusterV_nTrk.push_back(m_HFClusV[ic]->getAssociatedTracks().size());
- m_totEn_V += m_HFClusV[ic]->getEnergy();
+ t_Cone->Fill();
+ // Track axis
+ ClearTrackAxis();
+ int trackU_index=0;
+ int trackV_index=0;
+ for(int i=0; i<m_halfclusterU.size(); i++){ // loop half cluster U
+ std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisU = m_halfclusterU[i]->getHalfClusterCol("TrackAxis");
+ for(int ita=0; ita<m_mergedaxisU.size(); ita++){ // loop axis U
+ // General information of the axis
+ m_trackU_tag.push_back(trackU_index);
+ m_trackU_type.push_back(m_mergedaxisU[ita]->getType());
+ m_trackU_x.push_back(m_mergedaxisU[ita]->getPos().x());
+ m_trackU_y.push_back(m_mergedaxisU[ita]->getPos().y());
+ m_trackU_z.push_back(m_mergedaxisU[ita]->getPos().z());
+ m_trackU_E.push_back(m_mergedaxisU[ita]->getEnergy());
+
+ // MC truth information of the Axis
+ auto truthMap = m_mergedaxisU[ita]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_trackU_truth_tag.push_back(trackU_index);
+ m_trackU_truth_MC_px.push_back(iter.first.getMomentum().x);
+ m_trackU_truth_MC_py.push_back(iter.first.getMomentum().y);
+ m_trackU_truth_MC_pz.push_back(iter.first.getMomentum().z);
+ m_trackU_truth_MC_E.push_back(iter.first.getEnergy());
+ m_trackU_truth_MC_weight.push_back(iter.second);
+ }
+
+ // Hits on axis
+ for(int ilm=0; ilm<m_mergedaxisU[ita]->getCluster().size(); ilm++){ // loop local max
+ m_trackU_hit_tag.push_back(trackU_index);
+ m_trackU_hit_x.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().x() );
+ m_trackU_hit_y.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().y() );
+ m_trackU_hit_z.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().z() );
+ m_trackU_hit_E.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getEnergy() );
+ }
+
+ trackU_index++;
+ }
}
-
- t_Tower->Fill();
- }
-
-cout<<" Write 3D cluster"<<endl;
- //3D cluster
- ClearCluster();
- std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_EcalClusterCol = m_DataCol.map_CaloCluster["TrkMergedECAL"];
- std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_HcalClusterCol = m_DataCol.map_CaloCluster["HCALCluster"];
- std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_SimpleHcalClusterCol = m_DataCol.map_CaloCluster["SimpleHCALCluster"];
- m_totE_Ecal = 0.;
- m_totE_Hcal = 0.;
- m_Nclus_Ecal = m_EcalClusterCol.size();
- m_Nclus_Hcal = m_HcalClusterCol.size();
- for(int icl=0; icl<m_EcalClusterCol.size(); icl++){
- m_EcalClus_x.push_back(m_EcalClusterCol[icl]->getShowerCenter().x());
- m_EcalClus_y.push_back(m_EcalClusterCol[icl]->getShowerCenter().y());
- m_EcalClus_z.push_back(m_EcalClusterCol[icl]->getShowerCenter().z());
- m_EcalClus_E.push_back(m_EcalClusterCol[icl]->getLongiE());
- m_EcalClus_nTrk.push_back(m_EcalClusterCol[icl]->getAssociatedTracks().size());
-
- double tmp_phi = std::atan2(m_EcalClusterCol[icl]->getShowerCenter().y(), m_EcalClusterCol[icl]->getShowerCenter().x())* 180.0 / M_PI;
- if (tmp_phi < 0) tmp_phi += 360.0;
- double tmp_theta = std::atan2(m_EcalClusterCol[icl]->getShowerCenter().z(), m_EcalClusterCol[icl]->getShowerCenter().Perp())* 180.0 / M_PI + 90;
- //cout<<" Theta: "<<tmp_theta<<" Phi: "<<tmp_phi<<endl;
- m_EcalClus_Escale.push_back(m_energycorsvc->energyCorrection(m_EcalClusterCol[icl]->getLongiE(), tmp_phi, tmp_theta));
-
-
- if(m_EcalClusterCol[icl]->getAssociatedTracks().size()==1){
- const Track* trk = m_EcalClusterCol[icl]->getAssociatedTracks()[0];
- m_EcalClus_pTrk.push_back(trk->getMomentum());
-
- std::vector<TrackState> AllTrackStates = trk->getAllTrackStates();
- for(int istate=0; istate<AllTrackStates.size(); istate++){
- m_EcalClus_trk_tag.push_back(icl);
- m_EcalClus_trk_d0.push_back(AllTrackStates[istate].D0);
- m_EcalClus_trk_z0.push_back(AllTrackStates[istate].Z0);
- m_EcalClus_trk_phi.push_back(AllTrackStates[istate].phi0);
- m_EcalClus_trk_tanL.push_back( AllTrackStates[istate].tanLambda );
- m_EcalClus_trk_kappa.push_back( AllTrackStates[istate].Kappa);
- m_EcalClus_trk_omega.push_back( AllTrackStates[istate].Omega );
- m_EcalClus_trk_location.push_back( AllTrackStates[istate].location );
+ for(int i=0; i<m_halfclusterV.size(); i++){ // loop half cluster V
+ std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisV = m_halfclusterV[i]->getHalfClusterCol("TrackAxis");
+ for(int ita=0; ita<m_mergedaxisV.size(); ita++){ // loop axis V
+ // General information of the axis
+ m_trackV_tag.push_back(trackV_index);
+ m_trackV_type.push_back(m_mergedaxisV[ita]->getType());
+ m_trackV_x.push_back(m_mergedaxisV[ita]->getPos().x());
+ m_trackV_y.push_back(m_mergedaxisV[ita]->getPos().y());
+ m_trackV_z.push_back(m_mergedaxisV[ita]->getPos().z());
+ m_trackV_E.push_back(m_mergedaxisV[ita]->getEnergy());
+
+ // MC truth information of the Axis
+ auto truthMap = m_mergedaxisV[ita]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_trackV_truth_tag.push_back(trackV_index);
+ m_trackV_truth_MC_px.push_back(iter.first.getMomentum().x);
+ m_trackV_truth_MC_py.push_back(iter.first.getMomentum().y);
+ m_trackV_truth_MC_pz.push_back(iter.first.getMomentum().z);
+ m_trackV_truth_MC_E.push_back(iter.first.getEnergy());
+ m_trackV_truth_MC_weight.push_back(iter.second);
+ }
+
+ // Hits on axis
+ for(int ilm=0; ilm<m_mergedaxisV[ita]->getCluster().size(); ilm++){
+ m_trackV_hit_tag.push_back(trackV_index);
+ m_trackV_hit_x.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().x());
+ m_trackV_hit_y.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().y());
+ m_trackV_hit_z.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().z());
+ m_trackV_hit_E.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getEnergy());
+ }
+
+ trackV_index++;
}
-
}
- else
- m_EcalClus_pTrk.push_back(-99);
-
- m_EcalClus_typeU.push_back(m_EcalClusterCol[icl]->getHalfClusterUCol("LinkedLongiCluster")[0]->getType());
- m_EcalClus_typeV.push_back(m_EcalClusterCol[icl]->getHalfClusterVCol("LinkedLongiCluster")[0]->getType());
- for(int ii=0; ii<m_EcalClusterCol[icl]->getHalfClusterUCol("LinkedLongiCluster").size(); ii++){
- for(int ihit=0; ihit<m_EcalClusterCol[icl]->getHalfClusterUCol("LinkedLongiCluster")[ii]->getBars().size(); ihit++){
- auto shower = m_EcalClusterCol[icl]->getHalfClusterUCol("LinkedLongiCluster")[ii]->getBars()[ihit];
- m_EcalClus_hitU_tag.push_back(icl);
- m_EcalClus_hitU_x.push_back(shower->getPosition().x());
- m_EcalClus_hitU_y.push_back(shower->getPosition().y());
- m_EcalClus_hitU_z.push_back(shower->getPosition().z());
- m_EcalClus_hitU_E.push_back(shower->getEnergy());
+ t_TrackAxis->Fill();
+ //Axis
+ ClearAxis();
+ int axisU_index=0;
+ int axisV_index=0;
+ for(int i=0; i<m_halfclusterU.size(); i++){ // loop half cluster U
+ std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisU = m_halfclusterU[i]->getHalfClusterCol("MergedAxis");
+ for(int ita=0; ita<m_mergedaxisU.size(); ita++){ // loop axis U
+ // General information of the axis
+ m_axisU_tag.push_back(axisU_index);
+ m_axisU_type.push_back(m_mergedaxisU[ita]->getType());
+ m_axisU_x.push_back(m_mergedaxisU[ita]->getPos().x());
+ m_axisU_y.push_back(m_mergedaxisU[ita]->getPos().y());
+ m_axisU_z.push_back(m_mergedaxisU[ita]->getPos().z());
+ m_axisU_E.push_back(m_mergedaxisU[ita]->getEnergy());
+
+ // MC truth information of the Axis
+ auto truthMap = m_mergedaxisU[ita]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_axisU_truth_tag.push_back(axisU_index);
+ m_axisU_truth_MC_px.push_back(iter.first.getMomentum().x);
+ m_axisU_truth_MC_py.push_back(iter.first.getMomentum().y);
+ m_axisU_truth_MC_pz.push_back(iter.first.getMomentum().z);
+ m_axisU_truth_MC_E.push_back(iter.first.getEnergy());
+ m_axisU_truth_MC_weight.push_back(iter.second);
+ }
+ // Hits on axis
+ for(int ilm=0; ilm<m_mergedaxisU[ita]->getCluster().size(); ilm++){ // loop local max
+ m_axisU_hit_tag.push_back(axisU_index);
+ m_axisU_hit_x.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().x() );
+ m_axisU_hit_y.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().y() );
+ m_axisU_hit_z.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getPos().z() );
+ m_axisU_hit_E.push_back( m_mergedaxisU[ita]->getCluster()[ilm]->getEnergy() );
+ }
+
+ axisU_index++;
}
}
- for(int ii=0; ii<m_EcalClusterCol[icl]->getHalfClusterVCol("LinkedLongiCluster").size(); ii++){
- for(int ihit=0; ihit<m_EcalClusterCol[icl]->getHalfClusterVCol("LinkedLongiCluster")[ii]->getBars().size(); ihit++){
- auto shower = m_EcalClusterCol[icl]->getHalfClusterVCol("LinkedLongiCluster")[ii]->getBars()[ihit];
- m_EcalClus_hitV_tag.push_back(icl);
- m_EcalClus_hitV_x.push_back(shower->getPosition().x());
- m_EcalClus_hitV_y.push_back(shower->getPosition().y());
- m_EcalClus_hitV_z.push_back(shower->getPosition().z());
- m_EcalClus_hitV_E.push_back(shower->getEnergy());
+ for(int i=0; i<m_halfclusterV.size(); i++){ // loop half cluster V
+ std::vector<const Cyber::CaloHalfCluster*> m_mergedaxisV = m_halfclusterV[i]->getHalfClusterCol("MergedAxis");
+ for(int ita=0; ita<m_mergedaxisV.size(); ita++){ // loop axis V
+ // General information of the axis
+ m_axisV_tag.push_back(axisV_index);
+ m_axisV_type.push_back(m_mergedaxisV[ita]->getType());
+ m_axisV_x.push_back(m_mergedaxisV[ita]->getPos().x());
+ m_axisV_y.push_back(m_mergedaxisV[ita]->getPos().y());
+ m_axisV_z.push_back(m_mergedaxisV[ita]->getPos().z());
+ m_axisV_E.push_back(m_mergedaxisV[ita]->getEnergy());
+
+ // MC truth information of the Axis
+ auto truthMap = m_mergedaxisV[ita]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_axisV_truth_tag.push_back(axisV_index);
+ m_axisV_truth_MC_px.push_back(iter.first.getMomentum().x);
+ m_axisV_truth_MC_py.push_back(iter.first.getMomentum().y);
+ m_axisV_truth_MC_pz.push_back(iter.first.getMomentum().z);
+ m_axisV_truth_MC_E.push_back(iter.first.getEnergy());
+ m_axisV_truth_MC_weight.push_back(iter.second);
+ }
+ // Hits on axis
+ for(int ilm=0; ilm<m_mergedaxisV[ita]->getCluster().size(); ilm++){
+ m_axisV_hit_tag.push_back(axisV_index);
+ m_axisV_hit_x.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().x());
+ m_axisV_hit_y.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().y());
+ m_axisV_hit_z.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getPos().z());
+ m_axisV_hit_E.push_back(m_mergedaxisV[ita]->getCluster()[ilm]->getEnergy());
+ }
+
+ axisV_index++;
}
}
-
- m_totE_Ecal += m_EcalClusterCol[icl]->getLongiE();
- auto truthMap = m_EcalClusterCol[icl]->getLinkedMCP();
- for(auto iter: truthMap){
- m_EcalClus_truthMC_tag.push_back(icl);
- m_EcalClus_truthMC_pid.push_back(iter.first.getPDG() );
- m_EcalClus_truthMC_px.push_back(iter.first.getMomentum().x);
- m_EcalClus_truthMC_py.push_back(iter.first.getMomentum().y);
- m_EcalClus_truthMC_pz.push_back(iter.first.getMomentum().z);
- m_EcalClus_truthMC_E.push_back(iter.first.getEnergy());
- m_EcalClus_truthMC_EPx.push_back(iter.first.getEndpoint().x);
- m_EcalClus_truthMC_EPy.push_back(iter.first.getEndpoint().y);
- m_EcalClus_truthMC_EPz.push_back(iter.first.getEndpoint().z);
- m_EcalClus_truthMC_weight.push_back(iter.second);
+ m_halfclusterU.clear();
+ m_halfclusterV.clear();
+ for(int i=0; i<m_DataCol.map_HalfCluster["emptyHalfClusterU"].size(); i++){
+ m_halfclusterU.push_back( m_DataCol.map_HalfCluster["emptyHalfClusterU"][i].get() );
}
- }
-
- for(int icl=0; icl<m_HcalClusterCol.size(); icl++){
- m_HcalClus_x.push_back(m_HcalClusterCol[icl]->getHitCenter().x());
- m_HcalClus_y.push_back(m_HcalClusterCol[icl]->getHitCenter().y());
- m_HcalClus_z.push_back(m_HcalClusterCol[icl]->getHitCenter().z());
- m_HcalClus_E.push_back(m_HcalClusterCol[icl]->getHitsE());
- m_HcalClus_nTrk.push_back(m_HcalClusterCol[icl]->getAssociatedTracks().size());
- if(m_HcalClusterCol[icl]->getAssociatedTracks().size()==1)
- m_HcalClus_pTrk.push_back(m_HcalClusterCol[icl]->getAssociatedTracks()[0]->getMomentum());
- else
- m_HcalClus_pTrk.push_back(-99);
-
- for(int ih=0; ih<m_HcalClusterCol[icl]->getCaloHits().size(); ih++){
- m_HcalClus_hit_tag.push_back(icl);
- m_HcalClus_hit_x.push_back(m_HcalClusterCol[icl]->getCaloHits()[ih]->getPosition().x());
- m_HcalClus_hit_y.push_back(m_HcalClusterCol[icl]->getCaloHits()[ih]->getPosition().y());
- m_HcalClus_hit_z.push_back(m_HcalClusterCol[icl]->getCaloHits()[ih]->getPosition().z());
- m_HcalClus_hit_E.push_back(m_HcalClusterCol[icl]->getCaloHits()[ih]->getEnergy());
+ for(int i=0; i<m_DataCol.map_HalfCluster["emptyHalfClusterV"].size(); i++){
+ m_halfclusterV.push_back( m_DataCol.map_HalfCluster["emptyHalfClusterV"][i].get() );
}
-
- m_totE_Hcal += m_HcalClusterCol[icl]->getHitsE();
- auto truthMap = m_HcalClusterCol[icl]->getLinkedMCP();
- for(auto iter: truthMap){
- m_HcalClus_truthMC_tag.push_back(icl);
- m_HcalClus_truthMC_pid.push_back(iter.first.getPDG() );
- m_HcalClus_truthMC_px.push_back(iter.first.getMomentum().x);
- m_HcalClus_truthMC_py.push_back(iter.first.getMomentum().y);
- m_HcalClus_truthMC_pz.push_back(iter.first.getMomentum().z);
- m_HcalClus_truthMC_E.push_back(iter.first.getEnergy());
- m_HcalClus_truthMC_EPx.push_back(iter.first.getEndpoint().x);
- m_HcalClus_truthMC_EPy.push_back(iter.first.getEndpoint().y);
- m_HcalClus_truthMC_EPz.push_back(iter.first.getEndpoint().z);
- m_HcalClus_truthMC_weight.push_back(iter.second);
+ for(int i=0; i<m_halfclusterU.size(); i++){
+ m_emptyAxisU_tag.push_back(m_halfclusterU[i]->getType());
+ m_emptyAxisU_x.push_back(m_halfclusterU[i]->getPos().x());
+ m_emptyAxisU_y.push_back(m_halfclusterU[i]->getPos().y());
+ m_emptyAxisU_z.push_back(m_halfclusterU[i]->getPos().z());
+ m_emptyAxisU_E.push_back(m_halfclusterU[i]->getEnergy());
}
- }
-
- for(int icl=0; icl<m_SimpleHcalClusterCol.size(); icl++){
- m_SimpleHcalClus_x.push_back(m_SimpleHcalClusterCol[icl]->getHitCenter().x());
- m_SimpleHcalClus_y.push_back(m_SimpleHcalClusterCol[icl]->getHitCenter().y());
- m_SimpleHcalClus_z.push_back(m_SimpleHcalClusterCol[icl]->getHitCenter().z());
- m_SimpleHcalClus_E.push_back(m_SimpleHcalClusterCol[icl]->getHitsE());
- m_SimpleHcalClus_nTrk.push_back(m_SimpleHcalClusterCol[icl]->getAssociatedTracks().size());
- if(m_SimpleHcalClusterCol[icl]->getAssociatedTracks().size()==1)
- m_SimpleHcalClus_pTrk.push_back(m_SimpleHcalClusterCol[icl]->getAssociatedTracks()[0]->getMomentum());
- else
- m_SimpleHcalClus_pTrk.push_back(-99);
-
- for(int ih=0; ih<m_SimpleHcalClusterCol[icl]->getCaloHits().size(); ih++){
- m_SimpleHcalClus_hit_tag.push_back(icl);
- m_SimpleHcalClus_hit_x.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits()[ih]->getPosition().x());
- m_SimpleHcalClus_hit_y.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits()[ih]->getPosition().y());
- m_SimpleHcalClus_hit_z.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits()[ih]->getPosition().z());
- m_SimpleHcalClus_hit_E.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits()[ih]->getEnergy());
+ for(int i=0; i<m_halfclusterV.size(); i++){
+ m_emptyAxisV_tag.push_back(m_halfclusterV[i]->getType());
+ m_emptyAxisV_x.push_back(m_halfclusterV[i]->getPos().x());
+ m_emptyAxisV_y.push_back(m_halfclusterV[i]->getPos().y());
+ m_emptyAxisV_z.push_back(m_halfclusterV[i]->getPos().z());
+ m_emptyAxisV_E.push_back(m_halfclusterV[i]->getEnergy());
}
-
- auto truthMap = m_SimpleHcalClusterCol[icl]->getLinkedMCP();
- for(auto iter: truthMap){
- m_SimpleHcalClus_truthMC_tag.push_back(icl);
- m_SimpleHcalClus_truthMC_pid.push_back(iter.first.getPDG() );
- m_SimpleHcalClus_truthMC_px.push_back(iter.first.getMomentum().x);
- m_SimpleHcalClus_truthMC_py.push_back(iter.first.getMomentum().y);
- m_SimpleHcalClus_truthMC_pz.push_back(iter.first.getMomentum().z);
- m_SimpleHcalClus_truthMC_E.push_back(iter.first.getEnergy());
- m_SimpleHcalClus_truthMC_EPx.push_back(iter.first.getEndpoint().x);
- m_SimpleHcalClus_truthMC_EPy.push_back(iter.first.getEndpoint().y);
- m_SimpleHcalClus_truthMC_EPz.push_back(iter.first.getEndpoint().z);
- m_SimpleHcalClus_truthMC_weight.push_back(iter.second);
+
+ t_Axis->Fill();
+
+
+ //Half cluster
+ ClearHalfCluster();
+ m_halfclusterV.clear();
+ m_halfclusterU.clear();
+ m_totE_HFClusV = 0;
+ m_totE_HFClusU = 0;
+ //for(int i=0; i<m_DataCol.map_HalfCluster["ESHalfClusterU"].size(); i++){
+ // m_halfclusterU.push_back( m_DataCol.map_HalfCluster["ESHalfClusterU"][i]->getHalfClusterCol("MergedAxis")[0] );
+ //}
+ //for(int i=0; i<m_DataCol.map_HalfCluster["ESHalfClusterV"].size(); i++){
+ // m_halfclusterV.push_back( m_DataCol.map_HalfCluster["ESHalfClusterV"][i]->getHalfClusterCol("MergedAxis")[0] );
+ //}
+ for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColU"].size(); i++){
+ m_halfclusterU.push_back( m_DataCol.map_HalfCluster["HalfClusterColU"][i].get() );
}
- }
- t_Cluster->Fill();
-
- // Save Track info
- ClearTrack();
- std::vector<Cyber::Track*> m_trkCol;
- for(int it=0; it<m_DataCol.TrackCol.size(); it++)
- m_trkCol.push_back( m_DataCol.TrackCol[it].get() );
-
- m_Ntrk = m_trkCol.size();
- for(int itrk=0; itrk<m_Ntrk; itrk++){
- m_type.push_back(m_trkCol[itrk]->getType());
- std::vector<TrackState> AllTrackStates = m_trkCol[itrk]->getAllTrackStates();
- for(int istate=0; istate<AllTrackStates.size(); istate++){
- m_trkstate_d0.push_back( AllTrackStates[istate].D0 );
- m_trkstate_z0.push_back( AllTrackStates[istate].Z0 );
- m_trkstate_phi.push_back( AllTrackStates[istate].phi0 );
- m_trkstate_tanL.push_back( AllTrackStates[istate].tanLambda );
- m_trkstate_kappa.push_back( AllTrackStates[istate].Kappa);
- m_trkstate_omega.push_back( AllTrackStates[istate].Omega );
- m_trkstate_refx.push_back( AllTrackStates[istate].referencePoint.X() );
- m_trkstate_refy.push_back( AllTrackStates[istate].referencePoint.Y() );
- m_trkstate_refz.push_back( AllTrackStates[istate].referencePoint.Z() );
- m_trkstate_location.push_back( AllTrackStates[istate].location );
- m_trkstate_tag.push_back(itrk);
+ for(int i=0; i<m_DataCol.map_HalfCluster["HalfClusterColV"].size(); i++){
+ m_halfclusterV.push_back( m_DataCol.map_HalfCluster["HalfClusterColV"][i].get() );
}
- std::vector<TrackState> EcalTrackStates = m_trkCol[itrk]->getTrackStates("Ecal");
- for(int istate=0; istate<EcalTrackStates.size(); istate++){
- m_trkstate_x_ECAL.push_back(EcalTrackStates[istate].referencePoint.X());
- m_trkstate_y_ECAL.push_back(EcalTrackStates[istate].referencePoint.Y());
- m_trkstate_z_ECAL.push_back(EcalTrackStates[istate].referencePoint.Z());
- m_trkstate_tag_ECAL.push_back(itrk);
+ for(int i=0; i<m_halfclusterV.size(); i++){
+ m_HalfClusterV_x.push_back(m_halfclusterV[i]->getPos().x());
+ m_HalfClusterV_y.push_back(m_halfclusterV[i]->getPos().y());
+ m_HalfClusterV_z.push_back(m_halfclusterV[i]->getPos().z());
+ m_HalfClusterV_E.push_back(m_halfclusterV[i]->getEnergy());
+ m_HalfClusterV_tag.push_back(i);
+ m_HalfClusterV_type.push_back(m_halfclusterV[i]->getType());
+ m_totE_HFClusV += m_halfclusterV[i]->getEnergy();
+
+ // MC truth information of the HFCluster
+ auto truthMap = m_halfclusterV[i]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_HalfClusterV_truth_tag.push_back(i);
+ m_HalfClusterV_truthMC_px.push_back(iter.first.getMomentum().x);
+ m_HalfClusterV_truthMC_py.push_back(iter.first.getMomentum().y);
+ m_HalfClusterV_truthMC_pz.push_back(iter.first.getMomentum().z);
+ m_HalfClusterV_truthMC_E.push_back(iter.first.getEnergy());
+ m_HalfClusterV_truthMC_weight.push_back(iter.second);
+ }
+
+ // 1DClusters (hits)
+ for(int ilm=0; ilm<m_halfclusterV[i]->getCluster().size(); ilm++){
+ m_HalfClusterV_hit_tag.push_back(i);
+ m_HalfClusterV_hit_x.push_back( m_halfclusterV[i]->getCluster()[ilm]->getPos().x() );
+ m_HalfClusterV_hit_y.push_back( m_halfclusterV[i]->getCluster()[ilm]->getPos().y() );
+ m_HalfClusterV_hit_z.push_back( m_halfclusterV[i]->getCluster()[ilm]->getPos().z() );
+ m_HalfClusterV_hit_E.push_back( m_halfclusterV[i]->getCluster()[ilm]->getEnergy() );
+ }
}
- std::vector<TrackState> HcalTrackStates = m_trkCol[itrk]->getTrackStates("Hcal");
- for(int istate=0; istate<HcalTrackStates.size(); istate++){
- m_trkstate_x_ECAL.push_back(HcalTrackStates[istate].referencePoint.X());
- m_trkstate_y_ECAL.push_back(HcalTrackStates[istate].referencePoint.Y());
- m_trkstate_z_ECAL.push_back(HcalTrackStates[istate].referencePoint.Z());
- m_trkstate_tag_HCAL.push_back(itrk);
+ for(int i=0; i<m_halfclusterU.size(); i++){
+ m_HalfClusterU_x.push_back(m_halfclusterU[i]->getPos().x());
+ m_HalfClusterU_y.push_back(m_halfclusterU[i]->getPos().y());
+ m_HalfClusterU_z.push_back(m_halfclusterU[i]->getPos().z());
+ m_HalfClusterU_E.push_back(m_halfclusterU[i]->getEnergy());
+ m_HalfClusterU_tag.push_back(i);
+ m_HalfClusterU_type.push_back(m_halfclusterU[i]->getType());
+ m_totE_HFClusU += m_halfclusterU[i]->getEnergy();
+
+ // MC truth information of the HFCluster
+ auto truthMap = m_halfclusterU[i]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_HalfClusterU_truth_tag.push_back(i);
+ m_HalfClusterU_truthMC_px.push_back(iter.first.getMomentum().x);
+ m_HalfClusterU_truthMC_py.push_back(iter.first.getMomentum().y);
+ m_HalfClusterU_truthMC_pz.push_back(iter.first.getMomentum().z);
+ m_HalfClusterU_truthMC_E.push_back(iter.first.getEnergy());
+ m_HalfClusterU_truthMC_weight.push_back(iter.second);
+ }
+
+ // 1DClusters (hits)
+ for(int ilm=0; ilm<m_halfclusterU[i]->getCluster().size(); ilm++){
+ m_HalfClusterU_hit_tag.push_back(i);
+ m_HalfClusterU_hit_x.push_back( m_halfclusterU[i]->getCluster()[ilm]->getPos().x() );
+ m_HalfClusterU_hit_y.push_back( m_halfclusterU[i]->getCluster()[ilm]->getPos().y() );
+ m_HalfClusterU_hit_z.push_back( m_halfclusterU[i]->getCluster()[ilm]->getPos().z() );
+ m_HalfClusterU_hit_E.push_back( m_halfclusterU[i]->getCluster()[ilm]->getEnergy() );
+ }
}
- }
- t_Track->Fill();
-
- // yyy: pfo
- ClearPFO();
- std::vector<Cyber::PFObject*> m_pfobjects; m_pfobjects.clear();
- for(int ip=0; ip<m_DataCol.map_PFObjects["outputPFO"].size(); ip++)
- m_pfobjects.push_back(m_DataCol.map_PFObjects["outputPFO"][ip].get());
-
-
- for(int ip=0; ip<m_pfobjects.size(); ip++){
- std::vector<const Track*> t_tracks = m_pfobjects[ip]->getTracks();
- std::vector<const Calo3DCluster*> t_ecal_clusters = m_pfobjects[ip]->getECALClusters();
- std::vector<const Calo3DCluster*> t_hcal_clusters = m_pfobjects[ip]->getHCALClusters();
-
- pfo_tag.push_back(ip);
- pfo_n_track.push_back(t_tracks.size());
- pfo_n_ecal_clus.push_back(t_ecal_clusters.size());
- pfo_n_hcal_clus.push_back(t_hcal_clusters.size());
-
- for(int it=0; it<t_tracks.size(); it++){
- std::vector<TrackState> AllTrackStates = t_tracks[it]->getAllTrackStates();
- for(int istate=0; istate<AllTrackStates.size(); istate++){
- pfo_trk_tag.push_back(ip);
- pfo_trk_d0.push_back( AllTrackStates[istate].D0 );
- pfo_trk_z0.push_back( AllTrackStates[istate].Z0 );
- pfo_trk_phi.push_back( AllTrackStates[istate].phi0 );
- pfo_trk_tanL.push_back( AllTrackStates[istate].tanLambda );
- pfo_trk_kappa.push_back( AllTrackStates[istate].Kappa);
- pfo_trk_omega.push_back( AllTrackStates[istate].Omega );
- pfo_trk_location.push_back( AllTrackStates[istate].location );
+ t_HalfCluster->Fill();
+
+
+ //Tower
+ ClearTower();
+ std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_tower = m_DataCol.map_CaloCluster["ESTower"];
+ for(int it=0; it<m_tower.size(); it++){
+ ClearTower();
+ towerID[0] = m_tower[it]->getTowerID()[0][0];
+ towerID[1] = m_tower[it]->getTowerID()[0][1];
+ towerID[2] = m_tower[it]->getTowerID()[0][2];
+
+ std::vector<const CaloHalfCluster*> m_HFClusU = m_tower[it]->getHalfClusterUCol("ESHalfClusterU");
+ std::vector<const CaloHalfCluster*> m_HFClusV = m_tower[it]->getHalfClusterVCol("ESHalfClusterV");
+
+ m_NclusU = m_HFClusU.size();
+ m_NclusV = m_HFClusV.size();
+ m_totEn = m_tower[it]->getEnergy();
+ m_totEn_U = 0.;
+ m_totEn_V = 0.;
+ for(int ic=0; ic<m_NclusU; ic++){
+ m_HalfClusterU_tag.push_back(it);
+ m_HalfClusterU_x.push_back(m_HFClusU[ic]->getPos().x());
+ m_HalfClusterU_y.push_back(m_HFClusU[ic]->getPos().y());
+ m_HalfClusterU_z.push_back(m_HFClusU[ic]->getPos().z());
+ m_HalfClusterU_E.push_back(m_HFClusU[ic]->getEnergy());
+ m_HalfClusterU_type.push_back(m_HFClusU[ic]->getType());
+ m_HalfClusterU_nTrk.push_back(m_HFClusU[ic]->getAssociatedTracks().size());
+ m_totEn_U += m_HFClusU[ic]->getEnergy();
+ }
+
+ for(int ic=0; ic<m_NclusV; ic++){
+ m_HalfClusterV_tag.push_back(it);
+ m_HalfClusterV_x.push_back(m_HFClusV[ic]->getPos().x());
+ m_HalfClusterV_y.push_back(m_HFClusV[ic]->getPos().y());
+ m_HalfClusterV_z.push_back(m_HFClusV[ic]->getPos().z());
+ m_HalfClusterV_E.push_back(m_HFClusV[ic]->getEnergy());
+ m_HalfClusterV_type.push_back(m_HFClusV[ic]->getType());
+ m_HalfClusterV_nTrk.push_back(m_HFClusV[ic]->getAssociatedTracks().size());
+ m_totEn_V += m_HFClusV[ic]->getEnergy();
}
+
+ t_Tower->Fill();
}
- for(int ie=0; ie<t_ecal_clusters.size(); ie++){
- pfo_ecal_tag.push_back(ip);
- pfo_ecal_clus_x.push_back(t_ecal_clusters[ie]->getShowerCenter().x());
- pfo_ecal_clus_y.push_back(t_ecal_clusters[ie]->getShowerCenter().y());
- pfo_ecal_clus_z.push_back(t_ecal_clusters[ie]->getShowerCenter().z());
- pfo_ecal_clus_E.push_back(t_ecal_clusters[ie]->getLongiE());
-
- double tmp_phi = std::atan2(t_ecal_clusters[ie]->getShowerCenter().y(), t_ecal_clusters[ie]->getShowerCenter().x())* 180.0 / M_PI;
+ cout<<" Write 3D cluster"<<endl;
+ //3D cluster
+ ClearCluster();
+ std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_EcalClusterCol = m_DataCol.map_CaloCluster["TrkMergedECAL"];
+ std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_HcalClusterCol = m_DataCol.map_CaloCluster["HCALCluster"];
+ std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_SimpleHcalClusterCol = m_DataCol.map_CaloCluster["SimpleHCALCluster"];
+ m_totE_Ecal = 0.;
+ m_totE_Hcal = 0.;
+ m_Nclus_Ecal = m_EcalClusterCol.size();
+ m_Nclus_Hcal = m_HcalClusterCol.size();
+ for(int icl=0; icl<m_EcalClusterCol.size(); icl++){
+ m_EcalClus_x.push_back(m_EcalClusterCol[icl]->getShowerCenter().x());
+ m_EcalClus_y.push_back(m_EcalClusterCol[icl]->getShowerCenter().y());
+ m_EcalClus_z.push_back(m_EcalClusterCol[icl]->getShowerCenter().z());
+ m_EcalClus_E.push_back(m_EcalClusterCol[icl]->getLongiE());
+ m_EcalClus_nTrk.push_back(m_EcalClusterCol[icl]->getAssociatedTracks().size());
+
+ double tmp_phi = std::atan2(m_EcalClusterCol[icl]->getShowerCenter().y(), m_EcalClusterCol[icl]->getShowerCenter().x())* 180.0 / M_PI;
if (tmp_phi < 0) tmp_phi += 360.0;
- double tmp_theta = std::atan2(t_ecal_clusters[ie]->getShowerCenter().z(), t_ecal_clusters[ie]->getShowerCenter().Perp())* 180.0 / M_PI + 90;
- pfo_ecal_clus_Escale.push_back(m_energycorsvc->energyCorrection(t_ecal_clusters[ie]->getLongiE(), tmp_phi, tmp_theta));
-
+ double tmp_theta = std::atan2(m_EcalClusterCol[icl]->getShowerCenter().z(), m_EcalClusterCol[icl]->getShowerCenter().Perp())* 180.0 / M_PI + 90;
+ //cout<<" Theta: "<<tmp_theta<<" Phi: "<<tmp_phi<<endl;
+ m_EcalClus_Escale.push_back(m_energycorsvc->energyCorrection(m_EcalClusterCol[icl]->getLongiE(), tmp_phi, tmp_theta));
+
+
+ if(m_EcalClusterCol[icl]->getAssociatedTracks().size()==1){
+ const Track* trk = m_EcalClusterCol[icl]->getAssociatedTracks()[0];
+ m_EcalClus_pTrk.push_back(trk->getMomentum());
+
+ std::vector<TrackState> AllTrackStates = trk->getAllTrackStates();
+ for(int istate=0; istate<AllTrackStates.size(); istate++){
+ m_EcalClus_trk_tag.push_back(icl);
+ m_EcalClus_trk_d0.push_back(AllTrackStates[istate].D0);
+ m_EcalClus_trk_z0.push_back(AllTrackStates[istate].Z0);
+ m_EcalClus_trk_phi.push_back(AllTrackStates[istate].phi0);
+ m_EcalClus_trk_tanL.push_back( AllTrackStates[istate].tanLambda );
+ m_EcalClus_trk_kappa.push_back( AllTrackStates[istate].Kappa);
+ m_EcalClus_trk_omega.push_back( AllTrackStates[istate].Omega );
+ m_EcalClus_trk_location.push_back( AllTrackStates[istate].location );
+ }
+
+ }
+ else
+ m_EcalClus_pTrk.push_back(-99);
+
+ m_EcalClus_typeU.push_back(m_EcalClusterCol[icl]->getHalfClusterUCol("LinkedLongiCluster")[0]->getType());
+ m_EcalClus_typeV.push_back(m_EcalClusterCol[icl]->getHalfClusterVCol("LinkedLongiCluster")[0]->getType());
+ for(int ii=0; ii<m_EcalClusterCol[icl]->getHalfClusterUCol("LinkedLongiCluster").size(); ii++){
+ for(int ihit=0; ihit<m_EcalClusterCol[icl]->getHalfClusterUCol("LinkedLongiCluster")[ii]->getBars().size(); ihit++){
+ auto shower = m_EcalClusterCol[icl]->getHalfClusterUCol("LinkedLongiCluster")[ii]->getBars()[ihit];
+ m_EcalClus_hitU_tag.push_back(icl);
+ m_EcalClus_hitU_x.push_back(shower->getPosition().x());
+ m_EcalClus_hitU_y.push_back(shower->getPosition().y());
+ m_EcalClus_hitU_z.push_back(shower->getPosition().z());
+ m_EcalClus_hitU_E.push_back(shower->getEnergy());
+ }
+ }
+ for(int ii=0; ii<m_EcalClusterCol[icl]->getHalfClusterVCol("LinkedLongiCluster").size(); ii++){
+ for(int ihit=0; ihit<m_EcalClusterCol[icl]->getHalfClusterVCol("LinkedLongiCluster")[ii]->getBars().size(); ihit++){
+ auto shower = m_EcalClusterCol[icl]->getHalfClusterVCol("LinkedLongiCluster")[ii]->getBars()[ihit];
+ m_EcalClus_hitV_tag.push_back(icl);
+ m_EcalClus_hitV_x.push_back(shower->getPosition().x());
+ m_EcalClus_hitV_y.push_back(shower->getPosition().y());
+ m_EcalClus_hitV_z.push_back(shower->getPosition().z());
+ m_EcalClus_hitV_E.push_back(shower->getEnergy());
+ }
+ }
+
+ m_totE_Ecal += m_EcalClusterCol[icl]->getLongiE();
+ auto truthMap = m_EcalClusterCol[icl]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_EcalClus_truthMC_tag.push_back(icl);
+ m_EcalClus_truthMC_pid.push_back(iter.first.getPDG() );
+ m_EcalClus_truthMC_px.push_back(iter.first.getMomentum().x);
+ m_EcalClus_truthMC_py.push_back(iter.first.getMomentum().y);
+ m_EcalClus_truthMC_pz.push_back(iter.first.getMomentum().z);
+ m_EcalClus_truthMC_E.push_back(iter.first.getEnergy());
+ m_EcalClus_truthMC_EPx.push_back(iter.first.getEndpoint().x);
+ m_EcalClus_truthMC_EPy.push_back(iter.first.getEndpoint().y);
+ m_EcalClus_truthMC_EPz.push_back(iter.first.getEndpoint().z);
+ m_EcalClus_truthMC_weight.push_back(iter.second);
+ }
}
- for(int ih=0; ih<t_hcal_clusters.size(); ih++){
- pfo_hcal_tag.push_back(ip);
- pfo_hcal_clus_x.push_back(t_hcal_clusters[ih]->getHitCenter().x());
- pfo_hcal_clus_y.push_back(t_hcal_clusters[ih]->getHitCenter().y());
- pfo_hcal_clus_z.push_back(t_hcal_clusters[ih]->getHitCenter().z());
- pfo_hcal_clus_E.push_back(t_hcal_clusters[ih]->getHitsE());
+
+ for(int icl=0; icl<m_HcalClusterCol.size(); icl++){
+ m_HcalClus_x.push_back(m_HcalClusterCol[icl]->getHitCenter().x());
+ m_HcalClus_y.push_back(m_HcalClusterCol[icl]->getHitCenter().y());
+ m_HcalClus_z.push_back(m_HcalClusterCol[icl]->getHitCenter().z());
+ m_HcalClus_E.push_back(m_HcalClusterCol[icl]->getHitsE());
+ m_HcalClus_nTrk.push_back(m_HcalClusterCol[icl]->getAssociatedTracks().size());
+ if(m_HcalClusterCol[icl]->getAssociatedTracks().size()==1)
+ m_HcalClus_pTrk.push_back(m_HcalClusterCol[icl]->getAssociatedTracks()[0]->getMomentum());
+ else
+ m_HcalClus_pTrk.push_back(-99);
+
+ for(int ih=0; ih<m_HcalClusterCol[icl]->getCaloHits().size(); ih++){
+ m_HcalClus_hit_tag.push_back(icl);
+ m_HcalClus_hit_x.push_back(m_HcalClusterCol[icl]->getCaloHits()[ih]->getPosition().x());
+ m_HcalClus_hit_y.push_back(m_HcalClusterCol[icl]->getCaloHits()[ih]->getPosition().y());
+ m_HcalClus_hit_z.push_back(m_HcalClusterCol[icl]->getCaloHits()[ih]->getPosition().z());
+ m_HcalClus_hit_E.push_back(m_HcalClusterCol[icl]->getCaloHits()[ih]->getEnergy());
+ }
+
+ m_totE_Hcal += m_HcalClusterCol[icl]->getHitsE();
+ auto truthMap = m_HcalClusterCol[icl]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_HcalClus_truthMC_tag.push_back(icl);
+ m_HcalClus_truthMC_pid.push_back(iter.first.getPDG() );
+ m_HcalClus_truthMC_px.push_back(iter.first.getMomentum().x);
+ m_HcalClus_truthMC_py.push_back(iter.first.getMomentum().y);
+ m_HcalClus_truthMC_pz.push_back(iter.first.getMomentum().z);
+ m_HcalClus_truthMC_E.push_back(iter.first.getEnergy());
+ m_HcalClus_truthMC_EPx.push_back(iter.first.getEndpoint().x);
+ m_HcalClus_truthMC_EPy.push_back(iter.first.getEndpoint().y);
+ m_HcalClus_truthMC_EPz.push_back(iter.first.getEndpoint().z);
+ m_HcalClus_truthMC_weight.push_back(iter.second);
+ }
}
+
+ for(int icl=0; icl<m_SimpleHcalClusterCol.size(); icl++){
+ m_SimpleHcalClus_x.push_back(m_SimpleHcalClusterCol[icl]->getHitCenter().x());
+ m_SimpleHcalClus_y.push_back(m_SimpleHcalClusterCol[icl]->getHitCenter().y());
+ m_SimpleHcalClus_z.push_back(m_SimpleHcalClusterCol[icl]->getHitCenter().z());
+ m_SimpleHcalClus_E.push_back(m_SimpleHcalClusterCol[icl]->getHitsE());
+ m_SimpleHcalClus_nTrk.push_back(m_SimpleHcalClusterCol[icl]->getAssociatedTracks().size());
+ if(m_SimpleHcalClusterCol[icl]->getAssociatedTracks().size()==1)
+ m_SimpleHcalClus_pTrk.push_back(m_SimpleHcalClusterCol[icl]->getAssociatedTracks()[0]->getMomentum());
+ else
+ m_SimpleHcalClus_pTrk.push_back(-99);
+
+ for(int ih=0; ih<m_SimpleHcalClusterCol[icl]->getCaloHits().size(); ih++){
+ m_SimpleHcalClus_hit_tag.push_back(icl);
+ m_SimpleHcalClus_hit_x.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits()[ih]->getPosition().x());
+ m_SimpleHcalClus_hit_y.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits()[ih]->getPosition().y());
+ m_SimpleHcalClus_hit_z.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits()[ih]->getPosition().z());
+ m_SimpleHcalClus_hit_E.push_back(m_SimpleHcalClusterCol[icl]->getCaloHits()[ih]->getEnergy());
+ }
+
+ auto truthMap = m_SimpleHcalClusterCol[icl]->getLinkedMCP();
+ for(auto iter: truthMap){
+ m_SimpleHcalClus_truthMC_tag.push_back(icl);
+ m_SimpleHcalClus_truthMC_pid.push_back(iter.first.getPDG() );
+ m_SimpleHcalClus_truthMC_px.push_back(iter.first.getMomentum().x);
+ m_SimpleHcalClus_truthMC_py.push_back(iter.first.getMomentum().y);
+ m_SimpleHcalClus_truthMC_pz.push_back(iter.first.getMomentum().z);
+ m_SimpleHcalClus_truthMC_E.push_back(iter.first.getEnergy());
+ m_SimpleHcalClus_truthMC_EPx.push_back(iter.first.getEndpoint().x);
+ m_SimpleHcalClus_truthMC_EPy.push_back(iter.first.getEndpoint().y);
+ m_SimpleHcalClus_truthMC_EPz.push_back(iter.first.getEndpoint().z);
+ m_SimpleHcalClus_truthMC_weight.push_back(iter.second);
+ }
+ }
+ t_Cluster->Fill();
+
+ // Save Track info
+ ClearTrack();
+ std::vector<Cyber::Track*> m_trkCol;
+ for(int it=0; it<m_DataCol.TrackCol.size(); it++)
+ m_trkCol.push_back( m_DataCol.TrackCol[it].get() );
+
+ m_Ntrk = m_trkCol.size();
+ for(int itrk=0; itrk<m_Ntrk; itrk++){
+ m_type.push_back(m_trkCol[itrk]->getType());
+ std::vector<TrackState> AllTrackStates = m_trkCol[itrk]->getAllTrackStates();
+ for(int istate=0; istate<AllTrackStates.size(); istate++){
+ m_trkstate_d0.push_back( AllTrackStates[istate].D0 );
+ m_trkstate_z0.push_back( AllTrackStates[istate].Z0 );
+ m_trkstate_phi.push_back( AllTrackStates[istate].phi0 );
+ m_trkstate_tanL.push_back( AllTrackStates[istate].tanLambda );
+ m_trkstate_kappa.push_back( AllTrackStates[istate].Kappa);
+ m_trkstate_omega.push_back( AllTrackStates[istate].Omega );
+ m_trkstate_refx.push_back( AllTrackStates[istate].referencePoint.X() );
+ m_trkstate_refy.push_back( AllTrackStates[istate].referencePoint.Y() );
+ m_trkstate_refz.push_back( AllTrackStates[istate].referencePoint.Z() );
+ m_trkstate_location.push_back( AllTrackStates[istate].location );
+ m_trkstate_tag.push_back(itrk);
+ }
+ std::vector<TrackState> EcalTrackStates = m_trkCol[itrk]->getTrackStates("Ecal");
+ for(int istate=0; istate<EcalTrackStates.size(); istate++){
+ m_trkstate_x_ECAL.push_back(EcalTrackStates[istate].referencePoint.X());
+ m_trkstate_y_ECAL.push_back(EcalTrackStates[istate].referencePoint.Y());
+ m_trkstate_z_ECAL.push_back(EcalTrackStates[istate].referencePoint.Z());
+ m_trkstate_tag_ECAL.push_back(itrk);
+ }
+ std::vector<TrackState> HcalTrackStates = m_trkCol[itrk]->getTrackStates("Hcal");
+ for(int istate=0; istate<HcalTrackStates.size(); istate++){
+ m_trkstate_x_ECAL.push_back(HcalTrackStates[istate].referencePoint.X());
+ m_trkstate_y_ECAL.push_back(HcalTrackStates[istate].referencePoint.Y());
+ m_trkstate_z_ECAL.push_back(HcalTrackStates[istate].referencePoint.Z());
+ m_trkstate_tag_HCAL.push_back(itrk);
+ }
+ }
+ t_Track->Fill();
+
+ // yyy: pfo
+ ClearPFO();
+ std::vector<Cyber::PFObject*> m_pfobjects; m_pfobjects.clear();
+ for(int ip=0; ip<m_DataCol.map_PFObjects["outputPFO"].size(); ip++)
+ m_pfobjects.push_back(m_DataCol.map_PFObjects["outputPFO"][ip].get());
+
+
+ for(int ip=0; ip<m_pfobjects.size(); ip++){
+ std::vector<const Track*> t_tracks = m_pfobjects[ip]->getTracks();
+ std::vector<const Calo3DCluster*> t_ecal_clusters = m_pfobjects[ip]->getECALClusters();
+ std::vector<const Calo3DCluster*> t_hcal_clusters = m_pfobjects[ip]->getHCALClusters();
+
+ pfo_tag.push_back(ip);
+ pfo_n_track.push_back(t_tracks.size());
+ pfo_n_ecal_clus.push_back(t_ecal_clusters.size());
+ pfo_n_hcal_clus.push_back(t_hcal_clusters.size());
+
+ for(int it=0; it<t_tracks.size(); it++){
+ std::vector<TrackState> AllTrackStates = t_tracks[it]->getAllTrackStates();
+ for(int istate=0; istate<AllTrackStates.size(); istate++){
+ pfo_trk_tag.push_back(ip);
+ pfo_trk_d0.push_back( AllTrackStates[istate].D0 );
+ pfo_trk_z0.push_back( AllTrackStates[istate].Z0 );
+ pfo_trk_phi.push_back( AllTrackStates[istate].phi0 );
+ pfo_trk_tanL.push_back( AllTrackStates[istate].tanLambda );
+ pfo_trk_kappa.push_back( AllTrackStates[istate].Kappa);
+ pfo_trk_omega.push_back( AllTrackStates[istate].Omega );
+ pfo_trk_location.push_back( AllTrackStates[istate].location );
+ }
+ }
+
+ for(int ie=0; ie<t_ecal_clusters.size(); ie++){
+ pfo_ecal_tag.push_back(ip);
+ pfo_ecal_clus_x.push_back(t_ecal_clusters[ie]->getShowerCenter().x());
+ pfo_ecal_clus_y.push_back(t_ecal_clusters[ie]->getShowerCenter().y());
+ pfo_ecal_clus_z.push_back(t_ecal_clusters[ie]->getShowerCenter().z());
+ pfo_ecal_clus_E.push_back(t_ecal_clusters[ie]->getLongiE());
+
+ double tmp_phi = std::atan2(t_ecal_clusters[ie]->getShowerCenter().y(), t_ecal_clusters[ie]->getShowerCenter().x())* 180.0 / M_PI;
+ if (tmp_phi < 0) tmp_phi += 360.0;
+ double tmp_theta = std::atan2(t_ecal_clusters[ie]->getShowerCenter().z(), t_ecal_clusters[ie]->getShowerCenter().Perp())* 180.0 / M_PI + 90;
+ pfo_ecal_clus_Escale.push_back(m_energycorsvc->energyCorrection(t_ecal_clusters[ie]->getLongiE(), tmp_phi, tmp_theta));
+
+ }
+ for(int ih=0; ih<t_hcal_clusters.size(); ih++){
+ pfo_hcal_tag.push_back(ip);
+ pfo_hcal_clus_x.push_back(t_hcal_clusters[ih]->getHitCenter().x());
+ pfo_hcal_clus_y.push_back(t_hcal_clusters[ih]->getHitCenter().y());
+ pfo_hcal_clus_z.push_back(t_hcal_clusters[ih]->getHitCenter().z());
+ pfo_hcal_clus_E.push_back(t_hcal_clusters[ih]->getHitsE());
+ }
+ }
+ t_PFO->Fill();
+
}
- t_PFO->Fill();
//Clean Events
//system("/cefs/higgs/songwz/winter22/CEPCSW/workarea/memory/memory_test.sh before_clean");
@@ -1608,22 +1618,24 @@ cout<<" Write 3D cluster"<<endl;
StatusCode CyberPFAlg::finalize()
{
- m_wfile->cd();
- t_MCParticle->Write();
- t_SimBar->Write();
- t_LocalMax->Write();
- t_Layers->Write();
- t_Hough->Write();
- t_Cone->Write();
- t_TrackAxis->Write();
- t_Axis->Write();
- t_HalfCluster->Write();
- t_Tower->Write();
- t_Cluster->Write();
- t_Track->Write();
- t_PFO->Write();
- m_wfile->Close();
- delete m_wfile, t_MCParticle, t_SimBar, t_LocalMax, t_Layers, t_Hough, t_Cone, t_TrackAxis, t_Axis, t_HalfCluster, t_Tower, t_Cluster, t_Track, t_PFO;
+ if(m_WriteAna){
+ m_wfile->cd();
+ t_MCParticle->Write();
+ t_SimBar->Write();
+ t_LocalMax->Write();
+ t_Layers->Write();
+ t_Hough->Write();
+ t_Cone->Write();
+ t_TrackAxis->Write();
+ t_Axis->Write();
+ t_HalfCluster->Write();
+ t_Tower->Write();
+ t_Cluster->Write();
+ t_Track->Write();
+ t_PFO->Write();
+ m_wfile->Close();
+ delete m_wfile, t_MCParticle, t_SimBar, t_LocalMax, t_Layers, t_Hough, t_Cone, t_TrackAxis, t_Axis, t_HalfCluster, t_Tower, t_Cluster, t_Track, t_PFO;
+ }
delete m_pMCParticleCreator;
delete m_pTrackCreator;
diff --git a/Reconstruction/RecPFACyber/src/Objects/CaloHalfCluster.cc b/Reconstruction/RecPFACyber/src/Objects/CaloHalfCluster.cc
index 70555e6e1b0cb18e024189a692e6f02325796aec..68889cd38df04cf3f961f77ba48d177fab5c8f32 100644
--- a/Reconstruction/RecPFACyber/src/Objects/CaloHalfCluster.cc
+++ b/Reconstruction/RecPFACyber/src/Objects/CaloHalfCluster.cc
@@ -80,7 +80,6 @@ namespace Cyber{
for(int ii=0; ii<id.size(); ii++)
if( find(towerID.begin(), towerID.end(), id[ii])==towerID.end() ) towerID.push_back(id[ii]);
- fitAxis("");
}
}
@@ -90,7 +89,6 @@ namespace Cyber{
auto iter = find( m_1dclusters.begin(), m_1dclusters.end(), _1dcluster);
if( iter != m_1dclusters.end() ){
m_1dclusters.erase(iter);
- fitAxis("");
}
}
@@ -129,6 +127,11 @@ namespace Cyber{
return pos;
}
+ TVector3 CaloHalfCluster::getAxis() const{
+ if (axis.Mag() > 10.) fitAxis("");
+ return axis;
+ }
+
TVector3 CaloHalfCluster::getEnergyCenter() const{
TVector3 pos = getPos();
double maxEn = -99;
@@ -267,7 +270,7 @@ namespace Cyber{
}
- void CaloHalfCluster::fitAxis( std::string name ){
+ void CaloHalfCluster::fitAxis( std::string name ) const{
std::vector<const Cyber::Calo1DCluster*> barShowerCol; barShowerCol.clear();
if(!name.empty() && map_localMax.find(name)!=map_localMax.end() ) barShowerCol = map_localMax.at(name);
else barShowerCol = m_1dclusters;
@@ -337,7 +340,6 @@ namespace Cyber{
}
}
- fitAxis("");
}
diff --git a/Reconstruction/RecPFACyber/src/Objects/Track.cc b/Reconstruction/RecPFACyber/src/Objects/Track.cc
index bb083dd024f9a9b54a9c463420a2e99c6ff0b9fd..204f3e3e85f1e054725df899ab2c4cedd85990b7 100644
--- a/Reconstruction/RecPFACyber/src/Objects/Track.cc
+++ b/Reconstruction/RecPFACyber/src/Objects/Track.cc
@@ -9,6 +9,11 @@ namespace Cyber{
const double Track::B = 3.;
+ int Track::getTrackerHits() const {
+ if(m_track.isAvailable()) return m_track.trackerHits_size();
+ else return 0;
+ }
+
int Track::trackStates_size(std::string name) const{
std::vector<TrackState> emptyCol; emptyCol.clear();
if(m_trackStates.find(name)!=m_trackStates.end()) emptyCol = m_trackStates.at(name);
@@ -33,12 +38,37 @@ namespace Cyber{
return emptyCol;
}
+ float Track::getD0() const{
+ 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)
+ d0 = it.D0;
+ }
+ }
+
+ return d0;
+ }
+
+ float Track::getZ0() const{
+ 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)
+ z0 = it.Z0;
+ }
+ }
+
+ return z0;
+ }
+
+
float Track::getPt() const{
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)
- pt = 1./it.Kappa;
+ pt = 1./fabs(it.Kappa);
}
}
@@ -49,8 +79,8 @@ namespace Cyber{
std::vector<TrackState> trkStates = getTrackStates("Input");
float pz = -99.;
for(auto it: trkStates){
- if(it.location==4 || it.location==1){ //Calorimeter(for real track) or IP (for truth track)
- pz = it.tanLambda/it.Kappa;
+ if( (m_type!=0 && it.location==4) || (m_type==0 && it.location==1)){ //Calorimeter(for real track) or IP (for truth track)
+ pz = it.tanLambda/fabs(it.Kappa);
}
}
@@ -63,10 +93,10 @@ namespace Cyber{
float pt = -99.;
float pz = -99.;
for(auto it: trkStates){
- if(it.location==4 || it.location==1){ //Calorimeter(for real track) or IP (for truth track)
- pt = 1./it.Kappa;
+ if((m_type!=0 && it.location==4) || (m_type==0 && it.location==1)){ //Calorimeter(for real track) or IP (for truth track)
+ pt = 1./fabs(it.Kappa);
phi = it.phi0;
- pz = it.tanLambda/it.Kappa;
+ pz = it.tanLambda/fabs(it.Kappa);
}
}
@@ -78,7 +108,7 @@ namespace Cyber{
std::vector<TrackState> trkStates = getTrackStates("Input");
float omega = -99.;
for(auto it: trkStates){
- if(it.location==4 || it.location==1){ //Calorimeter(for real track) or IP (for truth track)
+ if((m_type!=0 && it.location==4) || (m_type==0 && it.location==1)){ //Calorimeter(for real track) or IP (for truth track)
omega = it.Omega;
}
}
@@ -86,6 +116,24 @@ namespace Cyber{
return omega/fabs(omega);
}
+ TVector3 Track::getStartPoint() const{
+ std::vector<TrackState> trkStates = getTrackStates("Input");
+ TVector3 startpoint (0.,0.,0.);
+ for(auto it: trkStates)
+ if(it.location==2) startpoint = it.referencePoint;
+
+ return startpoint;
+ }
+
+ TVector3 Track::getEndPoint() const{
+ std::vector<TrackState> trkStates = getTrackStates("Input");
+ TVector3 endpoint (0.,0.,0.);
+ for(auto it: trkStates)
+ if(it.location==3) endpoint = it.referencePoint;
+
+ return endpoint;
+ }
+
edm4hep::MCParticle Track::getLeadingMCP() const{
float maxWeight = -1.;
edm4hep::MCParticle mcp;
diff --git a/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp b/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp
index b36d83257e7f3c1d7589c8afef0150fb18a9fdb0..225315215f41fa3850fb43984b515160bdcf477a 100644
--- a/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp
+++ b/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp
@@ -33,7 +33,20 @@ namespace Cyber{
//PFO
std::vector<std::shared_ptr<Cyber::PFObject>> p_pfos = m_DataCol.map_PFObjects[settings.map_stringPars.at("OutputPFO")];
+/*
std::cout<<" Input PFO size: "<<p_pfos.size()<<std::endl;
+ double totE_Ecal = 0;
+ double totE_Hcal = 0;
+ for(int i=0; i<p_pfos.size(); i++){
+ cout<<" PFO #"<<i<<": track size "<<p_pfos[i]->getTracks().size()<<", leading P "<<p_pfos[i]->getTrackMomentum();
+ cout<<", ECAL cluster size "<<p_pfos[i]->getECALClusters().size()<<", totE "<<p_pfos[i]->getECALClusterEnergy();
+ cout<<", HCAL cluster size "<<p_pfos[i]->getHCALClusters().size()<<", totE "<<p_pfos[i]->getHCALClusterEnergy()<<endl;
+ totE_Ecal += p_pfos[i]->getECALClusterEnergy();
+ totE_Hcal += p_pfos[i]->getHCALClusterEnergy();
+ }
+ cout<<"-----Neutral cluster Ecal total energy: "<<totE_Ecal<<", Hcal total energy: "<<totE_Hcal<<endl;
+*/
+
for(int ip=0; ip<p_pfos.size(); ip++){
auto m_pfo = m_pfocol->create();
std::vector<const Track*> vec_trks = p_pfos[ip]->getTracks();
@@ -118,7 +131,7 @@ std::cout<<" Input PFO size: "<<p_pfos.size()<<std::endl;
m_clus.addToHits(_hit);
}
- double tmp_clusE = p_clus->getEnergy()*settings.map_floatPars.at("HCALCalib");
+ double tmp_clusE = p_clus->getHitsE()*settings.map_floatPars.at("HCALCalib");
m_clus.setEnergy( tmp_clusE );
edm4hep::Vector3f pos( p_clus->getHitCenter().x(), p_clus->getHitCenter().y(), p_clus->getHitCenter().z() );
m_clus.setPosition( pos );
@@ -172,9 +185,26 @@ std::cout<<" Input PFO size: "<<p_pfos.size()<<std::endl;
}
}
+/*
+double totE = 0;
+for(int i=0; i<m_pfocol->size(); i++){
+ auto m_pfo = m_pfocol->at(i);
+ if(m_pfo.getCharge()!=0) continue;
+ cout<<" PFO #"<<i<<": track size "<<m_pfo.tracks_size()<<", cluster size "<<m_pfo.clusters_size()<<", energy "<<m_pfo.getEnergy()<<endl;
+ totE += m_pfo.getEnergy();
+}
+cout<<"-----Neutral cluster total energy: "<<totE<<endl;
+totE = 0;
+for(int i=0; i<m_pfocol->size(); i++){
+ auto m_pfo = m_pfocol->at(i);
+ if(m_pfo.getCharge()==0) continue;
+ cout<<" PFO #"<<i<<": track size "<<m_pfo.tracks_size()<<", cluster size "<<m_pfo.clusters_size()<<", energy "<<m_pfo.getEnergy()<<endl;
+ totE += m_pfo.getEnergy();
+}
+cout<<"-----Charged cluster Ecal total energy: "<<totE<<endl;
std::cout<<" Created PFO size: "<<m_pfocol->size()<<std::endl;
-
+*/
return StatusCode::SUCCESS;
}
diff --git a/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp b/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
index aca3e85ca297d8cefee9d587bdbe3d986a00adef..7b09a9547c65498cf6dd2ad3768e014e5c870920 100644
--- a/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
+++ b/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp
@@ -24,6 +24,7 @@ namespace Cyber{
std::vector<edm4hep::Track> m_TrkCol; m_TrkCol.clear();
for(unsigned int icol=0; icol<const_TrkCol->size(); icol++){
edm4hep::Track m_trk = const_TrkCol->at(icol);
+ if(!m_trk.isAvailable() || m_trk.getType()==0) continue;
m_TrkCol.push_back(m_trk);
}
@@ -38,6 +39,7 @@ namespace Cyber{
for(auto iter : m_DataCol.collectionMap_Track){
auto const_TrkCol = iter.second;
for(int itrk=0; itrk<const_TrkCol.size(); itrk++){
+
//Cyber::Track* m_trk = new Cyber::Track();
std::shared_ptr<Cyber::Track> m_trk = std::make_shared<Cyber::Track>();
std::vector<Cyber::TrackState> m_trkstates;
@@ -49,7 +51,7 @@ namespace Cyber{
m_trkst.phi0 = const_TrkCol[itrk].getTrackStates(its).phi;
m_trkst.tanLambda = const_TrkCol[itrk].getTrackStates(its).tanLambda;
m_trkst.Omega = const_TrkCol[itrk].getTrackStates(its).omega;
- m_trkst.Kappa = m_trkst.Omega*1000./(0.3*settings.map_floatPars.at("BField"));
+ m_trkst.Kappa = m_trkst.Omega*1000./(0.299792458*settings.map_floatPars.at("BField"));
m_trkst.location = const_TrkCol[itrk].getTrackStates(its).location;
m_trkst.referencePoint.SetXYZ( const_TrkCol[itrk].getTrackStates(its).referencePoint[0],
const_TrkCol[itrk].getTrackStates(its).referencePoint[1],
@@ -68,9 +70,12 @@ namespace Cyber{
}
}
+
m_trkCol.push_back(m_trk);
}
}
+
+ //SelectGoodTrack(m_trkCol);
m_DataCol.TrackCol = m_trkCol;
@@ -150,5 +155,58 @@ namespace Cyber{
return StatusCode::SUCCESS;
}
+
+ StatusCode TrackCreator::SelectGoodTrack(std::vector<std::shared_ptr<Cyber::Track>>& trkCol){
+ if(trkCol.size()<2) return StatusCode::SUCCESS;
+
+ for(int itrk=0; itrk<trkCol.size(); itrk++){
+ //Endpoint cut
+ if( fabs( trkCol[itrk]->getEndPoint().z() )<settings.map_floatPars.at("TrkEndZCut") &&
+ trkCol[itrk]->getEndPoint().Perp()>settings.map_floatPars.at("TrkEndRCutMin") &&
+ trkCol[itrk]->getEndPoint().Perp()<settings.map_floatPars.at("TrkEndRCutMax") ){
+
+ trkCol.erase(trkCol.begin() + itrk);
+ itrk--;
+ continue;
+ }
+ //Start point cut
+ if( trkCol[itrk]->getStartPoint().Perp()>settings.map_floatPars.at("TrkStartRCutMin") &&
+ trkCol[itrk]->getStartPoint().Perp()<settings.map_floatPars.at("TrkStartRCutMax") ){
+
+ trkCol.erase(trkCol.begin() + itrk);
+ itrk--;
+ continue;
+ }
+ //Track length cut
+ if( (trkCol[itrk]->getEndPoint().Perp()-trkCol[itrk]->getStartPoint().Perp() )<settings.map_floatPars.at("TrkLengthCut") ){
+
+ trkCol.erase(trkCol.begin() + itrk);
+ itrk--;
+ continue;
+ }
+ }
+
+ if(trkCol.size()<2) return StatusCode::SUCCESS;
+
+ //Remove the broken tracks
+ std::sort(trkCol.begin(), trkCol.end(), compTrkIP);
+ for(int itrk=0; itrk<trkCol.size()-1; itrk++){
+ for(int jtrk=itrk+1; jtrk<trkCol.size(); jtrk++){
+ double deltaP = (trkCol[itrk]->getP3() - trkCol[jtrk]->getP3()).Mag();
+ if( trkCol[jtrk]->getP3().Perp() < settings.map_floatPars.at("BrokenTrkMinP") &&
+ ( deltaP/max(trkCol[itrk]->getMomentum(), trkCol[jtrk]->getMomentum()) < settings.map_floatPars.at("BrokenTrkDeltaPCut") ||
+ (trkCol[itrk]->getEndPoint()-trkCol[jtrk]->getStartPoint()).Mag() < settings.map_floatPars.at("BrokenTrkDistance") ) ){
+ trkCol.erase(trkCol.begin() + jtrk);
+ jtrk--;
+ if(itrk>jtrk+1) itrk--;
+ }
+ }
+ }
+
+ std::sort(trkCol.begin(), trkCol.end(), compTrkP);
+
+ return StatusCode::SUCCESS;
+ }
+
};
#endif