Reconstruction/RecPFACyber/include/Objects/Calo2DCluster.h

@@ -46,14 +46,14 @@ namespace Cyber {
     void setShowerUCol(std::vector<const Calo1DCluster*> _sh) { barShowerUCol=_sh; }
     void setShowerVCol(std::vector<const Calo1DCluster*> _sh) { barShowerVCol=_sh; }
     void addUnit(const Calo1DCluster* _1dcluster);
-    void addTowerID(int _m, int _p, int _s) { std::vector<int> id(3); id[0] = _m; id[1] = _p; id[2] = _s; towerID.push_back(id); }
+    void addTowerID(int _sys, int _m,  int _s, int _p) { std::vector<int> id(4); id[0] = _sys; id[1] = _m; id[2] = _s; id[3] = _p; towerID.push_back(id); }
     void addTowerID(std::vector<int> id) { towerID.push_back(id); }
     void setTowerID(std::vector<int> id) { towerID.clear(); towerID.push_back(id); }
     void setPos( TVector3 _vec ) { pos = _vec; }
     void setPos( double _x, double _y, double _z ) { pos.SetXYZ(_x, _y, _z); }
 
   private:
-    std::vector< std::vector<int> > towerID; //[module, stave]
+    std::vector< std::vector<int> > towerID; //[system, module, stave, part]
 
     TVector3 pos = TVector3(0.,0.,0.);
     std::vector<const CaloHit*> hits; 

Reconstruction/RecPFACyber/include/Objects/Calo3DCluster.h

@@ -58,6 +58,7 @@ namespace Cyber {
     int getType() const { return type; }
     double getEnergyScale() const { return Escale; }
 
+    void setAxis(TVector3 _axis ) { axis = _axis; }
     void setCaloHits( std::vector<const Cyber::CaloHit*> _hits ) { hits = _hits; }
     void setCaloHitsFrom2DCluster(); 
     void setTowers(std::vector<const Calo3DCluster*> _t) { m_towers = _t; }

Reconstruction/RecPFACyber/include/Objects/CaloHalfCluster.h

@@ -72,7 +72,7 @@ namespace Cyber {
     void setHoughPars(double _a, double _r) { Hough_alpha=_a; Hough_rho=_r; }
     void setIntercept(double _in) { Hough_intercept=_in; }
     void mergeHalfCluster( const CaloHalfCluster* clus );
-    void addTowerID(int _m, int _p, int _s) { std::vector<int> id(3); id[0] = _m; id[1] = _p; id[2] = _s; towerID.push_back(id); }
+    void addTowerID(int _sys, int _m,  int _s, int _p) { std::vector<int> id(4); id[0] = _sys; id[1] = _m; id[2] = _s; id[3] = _p; towerID.push_back(id); }
     void addTowerID(std::vector<int> id) { towerID.push_back(id); }
     void setTowerID(std::vector<int> id) { towerID.clear(); towerID.push_back(id); }
     void addAssociatedTrack(const Cyber::Track* _track){ m_TrackCol.push_back(_track); }
@@ -83,7 +83,7 @@ namespace Cyber {
 
   private:
     int type; // yyy: new definition: track: 10000, Hough: 100, cone: 1, merge: sum them
-    std::vector< std::vector<int> > towerID; //[module, part, stave]
+    std::vector< std::vector<int> > towerID; //[system, module, stave, part]
     int slayer;
     mutable TVector3 axis = TVector3(99999., 99999., 99999.);
     mutable double trk_dr;

Reconstruction/RecPFACyber/include/Objects/CaloHit.h

@@ -14,7 +14,7 @@ namespace Cyber{
     CaloHit() {};
     ~CaloHit() { Clear(); };
 
-    void Clear() { cellID=0; position.SetXYZ(0.,0.,0.); energy=-1; module=-1; layer=-1; ParentShower=nullptr; }
+    void Clear() { cellID=0; position.SetXYZ(0.,0.,0.); energy=-1; system = -1; module=-1; layer=-1; ParentShower=nullptr; }
     std::shared_ptr<CaloHit> Clone() const;
 
     void setOriginHit( edm4hep::CalorimeterHit& _hit ) { m_hit = _hit; }
@@ -22,15 +22,17 @@ namespace Cyber{
     double   getEnergy() const { return energy; } 
     int getLayer() const {return layer;}
     int getModule() const {return module;}
+    int getSystem() const {return system;}
     std::vector< std::pair<edm4hep::MCParticle, float> > getLinkedMCP() const { return MCParticleWeight; }
     edm4hep::MCParticle getLeadingMCP() const;
     float getLeadingMCPweight() const;
     edm4hep::CalorimeterHit getOriginHit() const { return m_hit; }
 
     void setcellID(unsigned long long _id) { cellID=_id; }
-    void setcellID(int _m, int _l) { module=_m; layer=_l; }
+    void setcellID(int _s, int _m, int _l) { system=_s; module=_m; layer=_l; }
     void setEnergy(double _en) { energy=_en; }
     void setPosition( TVector3 _vec ) { position=_vec; }
+    void setSystem(int _s) { system = _s; }
     void setModule(int _m ) { module = _m; }
     void setLayer(int _l) { layer = _l; }
     void setParentShower( Cyber::Calo2DCluster* _p ) { ParentShower=_p; }
@@ -38,6 +40,7 @@ namespace Cyber{
     void setLinkedMCP( std::vector<std::pair<edm4hep::MCParticle, float>> _pairVec ) { MCParticleWeight.clear(); MCParticleWeight = _pairVec; }
 
   private: 
+    int system; 
     int module; 
     int layer; 
     unsigned long long cellID; 

Reconstruction/RecPFACyber/include/Objects/CaloUnit.h

@@ -10,8 +10,8 @@ namespace Cyber{
   class CaloUnit{
 
   public:
-    CaloUnit(unsigned long long _cellID, int _system, int _module, int _slayer, int _dlayer, int _stave, int _bar, TVector3 _pos, double _Q1, double _Q2, double _T1, double _T2)
-    : cellID(_cellID), system(_system), module(_module), stave(_stave), dlayer(_dlayer), slayer(_slayer), bar(_bar), position(_pos), Q1(_Q1), Q2(_Q2), T1(_T1), T2(_T2) {}; 
+    CaloUnit(unsigned long long _cellID, int _system, int _module, int _slayer, int _dlayer, int _stave, int _part, int _bar, TVector3 _pos, double _Q1, double _Q2, double _T1, double _T2)
+    : cellID(_cellID), system(_system), module(_module), stave(_stave), part(_part), dlayer(_dlayer), slayer(_slayer), bar(_bar), position(_pos), Q1(_Q1), Q2(_Q2), T1(_T1), T2(_T2) {}; 
     CaloUnit() {};
     ~CaloUnit() { Clear(); };
     void Clear() { position.SetXYZ(0.,0.,0.); Q1=-1; Q2=-1; T1=-99; T2=-99; }
@@ -19,14 +19,20 @@ namespace Cyber{
     inline bool operator < (const CaloUnit &x) const {
       if(x.cellID==cellID) return false; 
 
-      if(slayer==0) return ( (stave<x.stave) || (stave==x.stave && bar<x.bar) );
-      else{
-        if( module==Nmodule && x.module==0 ) return true;
-        else if( module==0 && x.module==Nmodule ) return false; 
+      if(system==System_Barrel){
+        if(slayer==0) return ( (stave<x.stave) || (stave==x.stave && bar<x.bar) );
         else{
-          return ( (module<x.module) || (module==x.module && bar<x.bar) );
+          if( module==Nmodule && x.module==0 ) return true;
+          else if( module==0 && x.module==Nmodule ) return false; 
+          else{
+            return ( (module<x.module) || (module==x.module && bar<x.bar) );
+          }
         }
       }
+      else{
+        if(slayer==0) return (position.x() < x.getPosition().x());
+        else return (position.y() < x.getPosition().y());
+      }
     }
 
     inline bool operator == (const CaloUnit &x) const{
@@ -36,6 +42,7 @@ namespace Cyber{
     int getSystem() const { return system; }
     int getModule() const { return module; }
     int getStave()  const { return stave;  }
+    int getPart()   const { return part;   }
     int getDlayer() const { return dlayer; }
     int getSlayer() const { return slayer; }
     int getBar()    const { return bar;    }
@@ -44,6 +51,7 @@ namespace Cyber{
     double getT1()  const { return T1;     }
     double getT2()  const { return T2;     }
     double getBarLength() const {return barLength; }
+    int getNBarInLayer() const { return nBarInLayer; }
 
     TVector3 getPosition() const { return position; }
     double getEnergy() const { return (Q1+Q2); }
@@ -67,8 +75,11 @@ namespace Cyber{
     void setQ(double _q1, double _q2) { Q1=_q1; Q2=_q2; }
     void setT(double _t1, double _t2) { T1=_t1; T2=_t2; }
     void setBarLength(double _barLength) { barLength=_barLength; }
+    void setNBarInLayer(int _nBarInLayer) { nBarInLayer=_nBarInLayer; }
     std::shared_ptr<CaloUnit> Clone() const;
 
+    static int System_Barrel;
+    static int System_Endcap;
     static int Nmodule;
     static int Nstave;
     static int Nlayer;
@@ -77,6 +88,8 @@ namespace Cyber{
     static int NbarZ;
     static float barsize;
     static float ecal_innerR;
+    static float ecal_endcap_deadarea;
+    static float ecal_endcap_barsize;
 
   private:
 		unsigned long long cellID;
@@ -96,6 +109,7 @@ namespace Cyber{
     double barLength;
     std::vector< std::pair<edm4hep::MCParticle, float> > MCParticleWeight; 
 
+    int nBarInLayer;
   };
   
 }

Reconstruction/RecPFACyber/include/Objects/HoughObject.h

@@ -11,7 +11,7 @@ namespace Cyber {
   class HoughObject{
   public:
     HoughObject() {};
-    HoughObject( const Cyber::Calo1DCluster* _localmax, double _cellSize, double _ecal_inner_radius, double _phi=0.);
+    HoughObject( const Cyber::Calo1DCluster* _localmax, double _cellSize, double _ecal_inner_radius);
     ~HoughObject() { };
 
 
@@ -34,12 +34,13 @@ namespace Cyber {
 
     //int getModule() const { return (m_local_max->getTowerID())[0][0]; }
     int getSlayer() const { return m_local_max->getSlayer(); }
+    int getSystem() const { return m_local_max->getBars()[0]->getSystem(); }
     double getE() const { return m_local_max->getEnergy(); }
     double getCellSize() const { return m_cell_size; }
     const Cyber::Calo1DCluster* getLocalMax() const { return m_local_max; }
 
     void setCellSize(double _cs) { m_cell_size=_cs; }
-    void setCenterPoint(double& _ecal_inner_radius, double _phi=0.);
+    void setCenterPoint(double& _ecal_inner_radius);
     void setHoughLine(TF1& line1, TF1& line2);
 
 

Reconstruction/RecPFACyber/include/Tools/CaloHitsCreator.h

@@ -21,7 +21,8 @@ namespace Cyber{
                                std::vector<DataHandle<edm4hep::CalorimeterHitCollection>*>& r_CaloHitCols, 
                                std::map<std::string, dd4hep::DDSegmentation::BitFieldCoder*>& map_decoder, 
                                std::map<std::string, DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>*>& map_CaloParticleAssoCol,
-                               const dd4hep::VolumeManager& m_volumeManager ); 
+                               const dd4hep::VolumeManager& m_volumeManager,
+                               std::map<std::tuple<int, int, int, int, int>, int>& barNumberMapEndcapMap ); 
 
     //StatusCode CreateMCParticleCaloHitsAsso( std::vector<DataHandle<edm4hep::CalorimeterHitCollection>*>& r_CaloHitCols, 
     //                                         DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>* r_MCParticleRecoCaloCol );

Reconstruction/RecPFACyber/script/ana.py

@@ -14,13 +14,14 @@ cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest"
 ########## Podio Input ###################
 from Configurables import PodioInput
 inp = PodioInput("InputReader")
-inp.collections = [ "CyberPFO", "MCParticle" ]
+inp.collections = [ "CyberPFO", "CyberPFOPID", "MCParticle" ]
 ##########################################
 
 
 
 from Configurables import GenMatch
 genmatch = GenMatch("GenMatch")
+genmatch.InputPFOs = "CyberPFOPID"
 genmatch.nJets = 2
 genmatch.R = 0.6
 genmatch.OutputFile = "Jets_TDR_o1_v01.root"

Reconstruction/RecPFACyber/script/digi.py

@@ -86,16 +86,18 @@ EcalDigi.EcalSiPMPDE = 0.25              # NDL-EQR06, PDE 0.25
 EcalDigi.EcalSiPMDCR = 0                 # NDL-EQR06, dark count rate 2500000 [Hz]
 EcalDigi.EcalTimeInterval = 0.           # Time interval 0.000002 [s]. DCR*TimeInterval = dark count noise
 EcalDigi.EcalSiPMCT = 0.                 # SiPM crosstalk Probability 12%
-EcalDigi.EcalSiPMGainMean = 50           # 50 [ADC/p.e.]
+EcalDigi.EcalSiPMGainMean = 5            # 5 [ADC/p.e.]
 EcalDigi.EcalSiPMGainSigma = 0.08        # 0.08
 #EcalDigi.EcalSiPMNoiseSigma = 0          # 0
 # ADC
 EcalDigi.ADC = 8192                      # 13-bit, 8192
 EcalDigi.ADCSwitch = 8000                # 8000
 EcalDigi.Pedestal = 50                   # Pedestal 50 ADC
-EcalDigi.GainRatio_12 = 50               # Gain ratio 50
-EcalDigi.GainRatio_23 = 60               # Gain ratio 60
-#EcalDigi.EcalNoiseADCSigma = 4           # 4
+EcalDigi.GainRatio_12 = 30               # Gain ratio 30
+EcalDigi.GainRatio_23 = 10               # Gain ratio 10
+EcalDigi.EcalASICNoiseSigma = 4
+EcalDigi.EcalFEENoiseSigma = 5
+EcalDigi.ADCNonLinearity = 0             # ADC non-linearity 0
 # temperature control
 EcalDigi.UseCryTemp = 0
 EcalDigi.UseCryTempCor = 0

Reconstruction/RecPFACyber/script/rec.py

@@ -39,12 +39,12 @@ inp = PodioInput("InputReader")
 inp.collections = [ 
                     "ECALBarrel",
                     "ECALBarrelParticleAssoCol",
-#                    "ECALEndcaps",
-#                    "ECALEndcapsParticleAssoCol",
+                    "ECALEndcaps",
+                    "ECALEndcapsParticleAssoCol",
                     "HCALBarrel",
                     "HCALBarrelParticleAssoCol",
-#                    "HCALEndcaps",
-#                    "HCALEndcapsParticleAssoCol",
+                    "HCALEndcaps",
+                    "HCALEndcapsParticleAssoCol",
                     "MCParticle", 
                     "CompleteTracks", 
                     "CompleteTracksParticleAssociation",
@@ -75,18 +75,12 @@ CyberPFAlg.HcalNeutralCalib = 4.0
 CyberPFAlg.MCParticleCollection = "MCParticle"
 CyberPFAlg.TrackCollections = ["CompleteTracks"]
 CyberPFAlg.MCRecoTrackParticleAssociationCollection = "CompleteTracksParticleAssociation"
-#CyberPFAlg.ECalCaloHitCollections = ["ECALBarrel","ECALEndcaps"]
-#CyberPFAlg.ECalReadOutNames = ["EcalBarrelCollection","EcalEndcapsCollection"]
-#CyberPFAlg.ECalMCPAssociationName = ["ECALBarrelParticleAssoCol", "ECALEndcapsParticleAssoCol"]
-#CyberPFAlg.HCalCaloHitCollections = ["HCALBarrel", "HCALEndcaps"]
-#CyberPFAlg.HCalReadOutNames = ["HcalBarrelCollection", "HcalEndcapsCollection"]
-#CyberPFAlg.HCalMCPAssociationName = ["HCALBarrelParticleAssoCol", "HCALEndcapsParticleAssoCol"]
-CyberPFAlg.ECalCaloHitCollections = ["ECALBarrel"]
-CyberPFAlg.ECalReadOutNames = ["EcalBarrelCollection"]
-CyberPFAlg.ECalMCPAssociationName = ["ECALBarrelParticleAssoCol"]
-CyberPFAlg.HCalCaloHitCollections = ["HCALBarrel"]
-CyberPFAlg.HCalReadOutNames = ["HcalBarrelCollection"]
-CyberPFAlg.HCalMCPAssociationName = ["HCALBarrelParticleAssoCol"]
+CyberPFAlg.ECalCaloHitCollections = ["ECALBarrel","ECALEndcaps"]
+CyberPFAlg.ECalReadOutNames = ["EcalBarrelCollection","EcalEndcapsCollection"]
+CyberPFAlg.ECalMCPAssociationName = ["ECALBarrelParticleAssoCol", "ECALEndcapsParticleAssoCol"]
+CyberPFAlg.HCalCaloHitCollections = ["HCALBarrel", "HCALEndcaps"]
+CyberPFAlg.HCalReadOutNames = ["HcalBarrelCollection", "HcalEndcapsCollection"]
+CyberPFAlg.HCalMCPAssociationName = ["HCALBarrelParticleAssoCol", "HCALEndcapsParticleAssoCol"]
 
 ##--- Output collections ---
 CyberPFAlg.OutputPFO = "outputPFO";
@@ -113,8 +107,8 @@ CyberPFAlg.AlgParNames = [ ["InputECALBars","OutputECAL1DClusters","OutputECALHa
                                  ["OutputAxisName"], #6
                                  ["ReadinAxisName", "OutputClusName", "OutputTowerName"],  #9
                                  ["ReadinHFClusterName", "ReadinTowerName","OutputClusterName"], #11
-                                 ["InputHCALHits", "OutputHCALClusters"], #12
-                                 ["DoECALClustering","DoHCALClustering","InputHCALHits","OutputHCALClusters"], #15
+                                 ["OutputHCALClusters"], #12
+                                 ["DoECALClustering","DoHCALClustering","OutputHCALClusters"], #15
                                  ["ReadinECALClusterName", "ReadinHCALClusterName", "OutputCombPFO"],  #16
                                  ["ECALChargedCalib", "HCALChargedCalib", "ECALNeutralCalib", "HCALNeutralCalib"] ]#17
 CyberPFAlg.AlgParTypes = [ ["string","string","string"],#1
@@ -125,8 +119,8 @@ CyberPFAlg.AlgParTypes = [ ["string","string","string"],#1
                                  ["string"], #6
                                  ["string","string","string"],  #9
                                  ["string","string","string"], #11
-                                 ["string", "string"], #12
-                                 ["bool","bool","string","string"], #15
+                                 ["string"], #12
+                                 ["bool","bool","string"], #15
                                  ["string","string","string"],  #16
                                  ["double","double", "double","double"] ]#17
 CyberPFAlg.AlgParValues = [ ["BarCol","Cluster1DCol","HalfClusterCol"],#1
@@ -137,8 +131,8 @@ CyberPFAlg.AlgParValues = [ ["BarCol","Cluster1DCol","HalfClusterCol"],#1
                                   ["MergedAxis"], #6
                                   ["MergedAxis","ESHalfCluster","ESTower"],  #9
                                   ["ESHalfCluster","ESTower","EcalCluster"], #11
-                                  ["HCALBarrel", "SimpleHCALCluster"], #12
-                                  ["0","1","HCALBarrel","HCALCluster"], #15
+                                  ["SimpleHCALCluster"], #12
+                                  ["0","1","HCALCluster"], #15
                                   ["EcalCluster", "SimpleHCALCluster", "outputPFO"],  #16
                                   ["1.26","4.", "1.", "4."]  ]#17
 

Reconstruction/RecPFACyber/script/tracking.py

@@ -6,7 +6,7 @@ from Configurables import k4DataSvc
 dsvc = k4DataSvc("EventDataSvc", input="CaloDigi_TDR_o1_v01.root")
 
 from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
-seed = [1024]
+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.SetSingleton = True
@@ -42,7 +42,7 @@ tracksystemsvc = TrackSystemSvc("TrackSystemSvc")
 from Configurables import SimplePIDSvc
 pidsvc = SimplePIDSvc("SimplePIDSvc")
 cepcswdatatop = "/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest"
-pidsvc.ParFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Service/SimplePIDSvc/data/dNdx_TPC.root")
+pidsvc.ParFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Service/SimplePIDSvc/data/tdr25.1.0/dNdx_TPC.root")
 
 from Configurables import PodioInput
 podioinput = PodioInput("PodioReader", collections=[
@@ -51,13 +51,19 @@ podioinput = PodioInput("PodioReader", collections=[
     "VXDCollection",
     "SITCollection",
     "TPCCollection",
+#    "SETCollection",
     "OTKBarrelCollection",
     "FTDCollection",
     "MuonBarrelCollection",
     "MuonEndcapCollection"
     ])
 
-# digitization
+
+##################
+# Digitization
+##################
+
+## Config ##
 vxdhitname  = "VXDTrackerHits"
 sithitname  = "SITTrackerHits"
 gashitname  = "TPCTrackerHits"
@@ -65,17 +71,14 @@ sethitname  = "OTKBarrelTrackerHits"
 setspname   = "OTKBarrelSpacePoints"
 ftdhitname  = "FTDTrackerHits"
 ftdspname   = "FTDSpacePoints"
-from Configurables import SmearDigiTool
+from Configurables import SmearDigiTool,SiTrackerDigiAlg
+
+## VXD ##
 vxdtool = SmearDigiTool("VXD")
 vxdtool.ResolutionU = [0.005]
 vxdtool.ResolutionV = [0.005]
-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
@@ -83,46 +86,46 @@ digiVXD.TrackerHitAssociationCollection = "VXDTrackerHitAssociation"
 digiVXD.DigiTool = "SmearDigiTool/VXD"
 #digiVXD.OutputLevel = DEBUG
 
-from Configurables import PlanarDigiAlg
-digiSIT = PlanarDigiAlg("SITDigi")
-digiSIT.IsStrip = False
+## SIT ##
+sittool = SmearDigiTool("SIT")
+sittool.ResolutionU = [0.0098]
+sittool.ResolutionV = [0.0433]
+#sittool.OutputLevel = DEBUG
+
+digiSIT = SiTrackerDigiAlg("SITDigi")
 digiSIT.SimTrackHitCollection = "SITCollection"
 digiSIT.TrackerHitCollection = sithitname
 digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation"
-digiSIT.ResolutionU = [0.0098]
-digiSIT.ResolutionV = [0.0433]
-digiSIT.UsePlanarTag = True
-digiSIT.ParameterizeResolution = False
-digiSIT.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452]
-digiSIT.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01]
+digiSIT.DigiTool = "SmearDigiTool/SIT"
 #digiSIT.OutputLevel = DEBUG
 
-digiSET = PlanarDigiAlg("SETDigi")
-digiSET.IsStrip = False
-digiSET.SimTrackHitCollection = "OTKBarrelCollection"
-digiSET.TrackerHitCollection = sethitname
-digiSET.TrackerHitAssociationCollection = "OTKBarrelTrackerHitAssociation"
-digiSET.ResolutionU = [0.010]
-digiSET.ResolutionV = [1.000]
-digiSET.UsePlanarTag = True
-digiSET.ParameterizeResolution = False
-digiSET.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452]
-digiSET.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01]
-#digiSET.OutputLevel = DEBUG
-
-digiFTD = PlanarDigiAlg("FTDDigi")
-digiFTD.IsStrip = False
+## OTKBarrel ##
+otkbtool = SmearDigiTool("OTKBarrel")
+otkbtool.ResolutionU = [0.010]
+otkbtool.ResolutionV = [1.000]
+#otkbtool.OutputLevel = DEBUG
+
+digiOTKB = SiTrackerDigiAlg("OTKBarrelDigi")
+digiOTKB.SimTrackHitCollection = "OTKBarrelCollection"
+digiOTKB.TrackerHitCollection = sethitname
+digiOTKB.TrackerHitAssociationCollection = "OTKBarrelTrackerHitAssociation"
+digiOTKB.DigiTool = "SmearDigiTool/OTKBarrel"
+#digiOTKB.OutputLevel = DEBUG
+
+## FTD ##
+ftdtool = SmearDigiTool("FTD")
+ftdtool.ResolutionU = [0.0072]
+ftdtool.ResolutionV = [0.086]
+#ftdtool.OutputLevel = DEBUG
+
+digiFTD = SiTrackerDigiAlg("FTDDigi")
 digiFTD.SimTrackHitCollection = "FTDCollection"
 digiFTD.TrackerHitCollection = ftdhitname
 digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation"
-digiFTD.ResolutionU = [0.0072]
-digiFTD.ResolutionV = [0.086]
-digiFTD.UsePlanarTag = True
-digiFTD.ParameterizeResolution = False
-digiFTD.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452]
-digiFTD.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01]
+digiFTD.DigiTool = "SmearDigiTool/FTD"
 #digiFTD.OutputLevel = DEBUG
 
+## TPC ##
 from Configurables import TPCDigiAlg
 digiTPC = TPCDigiAlg("TPCDigi")
 digiTPC.TPCCollection = "TPCCollection"
@@ -137,6 +140,7 @@ digiTPC.TPCTrackerHitsCol = gashitname
 #digiTPC.N_eff = 30
 #digiTPC.OutputLevel = DEBUG
 
+
 ## Muon Detector ##
 from Configurables import MuonDigiAlg
 digiMuon = MuonDigiAlg("MuonDigiAlg")
@@ -146,8 +150,11 @@ digiMuon.MuonBarrelTrackerHits = "MuonBarrelTrackerHits"
 digiMuon.MuonEndcapTrackerHits = "MuonEndcapTrackerHits"
 digiMuon.WriteNtuple = 0
 digiMuon.OutFileName = "Digi_MUON.root"
+#########################################
 
-# tracking
+################
+# Tracking
+################
 from Configurables import KalTestTool
 # Close multiple scattering and smooth, used by clupatra
 kt010 = KalTestTool("KalTest010")
@@ -266,41 +273,19 @@ tmt.MuonTagEfficiency = 0.95 # muon true tag efficiency, default is 1.0 (100%)
 tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta)
 #tmt.OutputLevel = DEBUG
 
-#from Configurables import ReadDigiAlg
-#readtrk = ReadDigiAlg("ReadDigiAlg")
-#readtrk.SiTracks = "SubsetTracks"
-#readtrk.TPCTracks = "ClupatraTracks"
-#readtrk.FullTracks = "CompleteTracks"
-#readtrk.TPCTracksAssociation = "ClupatraTracksParticleAssociation"
-#readtrk.FullTracksAssociation = "CompleteTracksParticleAssociation"
-#readtrk.OutFileName = "TrackAnaTuple_mu.root"
-
 # output
 from Configurables import PodioOutput
 out = PodioOutput("outputalg")
 out.filename = "Tracking_TDR_o1_v01.root"
-out.outputCommands = ["drop *",
-  "keep ECALBarrel",
-  "keep ECALBarrelParticleAssoCol",
-  "keep HCALBarrel",
-  "keep HCALBarrelParticleAssoCol",
-  "keep ECALEndcaps",
-  "keep HCALEndcaps",
-  "keep ECALEndcapsParticleAssoCol",
-  "keep HCALEndcapsParticleAssoCol",
-  "keep MCParticle",
-  "keep CompleteTracks",
-  "keep CompleteTracksParticleAssociation",
-  "keep RecTofCollection",
-  "keep DndxTracks" ]
+out.outputCommands = ["keep *"]
 
 # ApplicationMgr
 from Configurables import ApplicationMgr
 mgr = ApplicationMgr(
-    TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tof, tmt, out],
+    TopAlg = [podioinput, digiVXD, digiSIT, digiOTKB, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tof, tmt, out],
     EvtSel = 'NONE',
     EvtMax = 10,
     ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc],
-    #HistogramPersistency = 'ROOT',
-    OutputLevel = INFO
+    HistogramPersistency = 'ROOT',
+    OutputLevel = ERROR
 )

Reconstruction/RecPFACyber/src/Algorithm/AxisMergingAlg.cpp

@@ -190,7 +190,6 @@ cout<<endl;
             tmp_goodAxis.push_back( m_newAxisUCol[ic] );
           }
         
-
       else {
         tmp_badAxis.push_back( m_newAxisUCol[ic] );
       }
@@ -573,11 +572,74 @@ StatusCode AxisMergingAlg::BranchMerging( std::vector<Cyber::CaloHalfCluster*>&
       double hough_rho = p_axis->getHoughRho();
       double hough_alpha = p_axis->getHoughAlpha();
 
-      // V plane
-      if(m_axis->getSlayer()==1){
-        double x0 = m_axis->getEnergyCenter().x();
-        double y0 = m_axis->getEnergyCenter().y();
-        double distance = TMath::Abs( x0*TMath::Cos(hough_alpha) + y0*TMath::Sin(hough_alpha) - hough_rho );
+      //Barrel 
+      if( m_axis->getTowerID()[0][0] == Cyber::CaloUnit::System_Barrel ){ 
+        // V plane
+        if(m_axis->getSlayer()==1){
+          double x0 = m_axis->getEnergyCenter().x();
+          double y0 = m_axis->getEnergyCenter().y();
+          double distance = TMath::Abs( x0*TMath::Cos(hough_alpha) + y0*TMath::Sin(hough_alpha) - hough_rho );
+
+          //cout << "    yyy:V rho = " << hough_rho << ", alpha = " << hough_alpha << ", x0 = " << x0 << ", y0 = " << y0 << endl;
+          //cout << "         distanceV = " << distance << endl;
+
+          if (distance<settings.map_floatPars["th_branch_distance"]){
+            m_axisCol[iax]->mergeHalfCluster( m_axisCol[jax] );
+            int axis_type = m_axisCol[iax]->getType() + m_axisCol[jax]->getType();
+            m_axisCol[iax]->setType(axis_type);
+            m_axisCol.erase(m_axisCol.begin()+jax);
+            jax--;
+            if(iax>jax+1) iax--;
+            // cout << "  yyy: axis " << jax << " is merged into axis " << iax << endl;
+          }
+        }
+
+        // U plane
+        else{
+          int m_module = m_axis->getEnergyCenterTower()[1];
+          //cout << "  yyy:branceU: m_module = " << m_module << endl;
+          int p_module = p_axis->getTowerID()[0][1];
+          //cout << "  yyy:branceU: p_module = " << p_module << endl;
+          // Do not merge the two axis in two different modules
+          //if (m_module != p_module) continue;
+
+          TVector3 t_pos = m_axis->getEnergyCenter();
+          // cout << "  yyy:branceU: t_pos = " << t_pos.x() << ", " << t_pos.y() << ", " << t_pos.z() << endl;
+          t_pos.RotateZ( TMath::TwoPi()/Cyber::CaloUnit::Nmodule*(int(Cyber::CaloUnit::Nmodule*3./4.)-m_module) );
+          // cout << "  yyy:branceU: t_pos after rotate to module 6 = " << t_pos.x() << ", " << t_pos.y() << ", " << t_pos.z() << endl;
+          double x0 = t_pos.x();
+          double y0 = t_pos.z();
+          double distance = TMath::Abs( x0*TMath::Cos(hough_alpha) + y0*TMath::Sin(hough_alpha) - hough_rho );
+
+          //cout << "    yyy:U rho = " << hough_rho << ", alpha = " << hough_alpha << ", x0 = " << x0 << ", z0 = " << y0 << endl;
+          //cout << "         distanceU = " << distance << endl;
+
+          if (distance<settings.map_floatPars["th_branch_distance"]){
+            m_axisCol[iax]->mergeHalfCluster( m_axisCol[jax] );
+            int axis_type = m_axisCol[iax]->getType() + m_axisCol[jax]->getType();
+            m_axisCol[iax]->setType(axis_type);
+            m_axisCol.erase(m_axisCol.begin()+jax);
+            //cout << "  yyy: axis " << jax << " is merged into axis " << iax << endl;
+            jax--;
+            if(iax>jax+1) iax--;
+          }
+        }
+      }
+
+      //Endcaps
+      else if(m_axis->getTowerID()[0][0] == Cyber::CaloUnit::System_Endcap) {
+        // V plane
+				double x0, y0, distance;
+        if(m_axis->getSlayer()==1){
+          x0 = m_axis->getEnergyCenter().z();
+          y0 = m_axis->getEnergyCenter().y();
+          distance = TMath::Abs( x0*TMath::Cos(hough_alpha) + y0*TMath::Sin(hough_alpha) - hough_rho );
+				}
+				else{
+          x0 = m_axis->getEnergyCenter().z();
+          y0 = m_axis->getEnergyCenter().x();
+          distance = TMath::Abs( x0*TMath::Cos(hough_alpha) + y0*TMath::Sin(hough_alpha) - hough_rho );
+				}
 
         //cout << "    yyy:V rho = " << hough_rho << ", alpha = " << hough_alpha << ", x0 = " << x0 << ", y0 = " << y0 << endl;
         //cout << "         distanceV = " << distance << endl;
@@ -592,40 +654,10 @@ StatusCode AxisMergingAlg::BranchMerging( std::vector<Cyber::CaloHalfCluster*>&
           // cout << "  yyy: axis " << jax << " is merged into axis " << iax << endl;
         }
       }
-
-      // U plane
       else{
-        int m_module = m_axis->getEnergyCenterTower()[0];
-        //cout << "  yyy:branceU: m_module = " << m_module << endl;
-        int p_module = p_axis->getTowerID()[0][0];
-        //cout << "  yyy:branceU: p_module = " << p_module << endl;
-        // Do not merge the two axis in two different modules
-        //if (m_module != p_module) continue;
-
-        TVector3 t_pos = m_axis->getEnergyCenter();
-        // cout << "  yyy:branceU: t_pos = " << t_pos.x() << ", " << t_pos.y() << ", " << t_pos.z() << endl;
-        t_pos.RotateZ( TMath::TwoPi()/Cyber::CaloUnit::Nmodule*(int(Cyber::CaloUnit::Nmodule*3./4.)-m_module) );
-        // cout << "  yyy:branceU: t_pos after rotate to module 6 = " << t_pos.x() << ", " << t_pos.y() << ", " << t_pos.z() << endl;
-        double x0 = t_pos.x();
-        double y0 = t_pos.z();
-        double distance = TMath::Abs( x0*TMath::Cos(hough_alpha) + y0*TMath::Sin(hough_alpha) - hough_rho );
-
-        //cout << "    yyy:U rho = " << hough_rho << ", alpha = " << hough_alpha << ", x0 = " << x0 << ", z0 = " << y0 << endl;
-        //cout << "         distanceU = " << distance << endl;
-
-        if (distance<settings.map_floatPars["th_branch_distance"]){
-          m_axisCol[iax]->mergeHalfCluster( m_axisCol[jax] );
-          int axis_type = m_axisCol[iax]->getType() + m_axisCol[jax]->getType();
-          m_axisCol[iax]->setType(axis_type);
-          m_axisCol.erase(m_axisCol.begin()+jax);
-          //cout << "  yyy: axis " << jax << " is merged into axis " << iax << endl;
-          jax--;
-          if(iax>jax+1) iax--;
-        }
-
-
+        std::cout<<"ERROR in AxisMergingAlg: Unknown system ID: "<<m_axis->getTowerID()[0][0]<<std::endl;
       }
-
+			
       p_axis=nullptr;
     }
     m_axis=nullptr;
@@ -805,10 +837,7 @@ StatusCode AxisMergingAlg::ConeMerging( std::vector<Cyber::CaloHalfCluster*>& m_
   //cout<<"    Merge! "<<endl;
 
         m_axisCol[iax]->mergeHalfCluster( m_axisCol[jax] );
-        if(m_axisCol[iax]->getType()==1 || m_axisCol[jax]->getType()==1 || m_axisCol[iax]->getType()==5 || m_axisCol[jax]->getType()==5) m_axisCol[iax]->setType(5);
-        else if(m_axisCol[iax]->getType()==3 || m_axisCol[jax]->getType()==3 || m_axisCol[iax]->getType()==6 || m_axisCol[jax]->getType()==6) m_axisCol[iax]->setType(6);
-        else if(m_axisCol[iax]->getType()==4 || m_axisCol[jax]->getType()==4) m_axisCol[iax]->setType(4);
-        else m_axisCol[iax]->setType(7);
+        m_axisCol[iax]->setType( m_axisCol[iax]->getType() + m_axisCol[jax]->getType() );
 
         //delete m_axisCol[jax]; m_axisCol[jax]=nullptr;
         m_axisCol.erase(m_axisCol.begin()+jax);
@@ -852,11 +881,8 @@ StatusCode AxisMergingAlg::ConeMerging( std::vector<Cyber::CaloHalfCluster*>& m_
           relDis <= settings.map_floatPars["relP_Dis"] && relDis>=0 ){
 
         m_axisCol[iax]->mergeHalfCluster( m_axisCol[jax] );
-        if(m_axisCol[iax]->getType()==1 || m_axisCol[jax]->getType()==1 || m_axisCol[iax]->getType()==5 || m_axisCol[jax]->getType()==5) m_axisCol[iax]->setType(5);
-        else if(m_axisCol[iax]->getType()==3 || m_axisCol[jax]->getType()==3 || m_axisCol[iax]->getType()==6 || m_axisCol[jax]->getType()==6) m_axisCol[iax]->setType(6);
-        else if(m_axisCol[iax]->getType()==4 || m_axisCol[jax]->getType()==4) m_axisCol[iax]->setType(4);
-        else m_axisCol[iax]->setType(7);
-
+				m_axisCol[iax]->setType( m_axisCol[iax]->getType()+m_axisCol[jax]->getType() );
+	
         //delete m_axisCol[jax]; m_axisCol[jax]=nullptr;
         m_axisCol.erase(m_axisCol.begin()+jax);
         jax--;

Reconstruction/RecPFACyber/src/Algorithm/ConeClustering2DAlg.cpp

@@ -190,20 +190,38 @@ TVector2 ConeClustering2DAlg::GetProjectedAxis( const Cyber::CaloHalfCluster* m_
 
 TVector2 ConeClustering2DAlg::GetProjectedRelR( const Cyber::Calo1DCluster* m_shower1, const Cyber::Calo1DCluster* m_shower2 ){
   TVector2 paxis1, paxis2;
-  if(m_shower1->getSlayer()==1){ //For V-bars
-    paxis1.Set(m_shower1->getPos().x(), m_shower1->getPos().y());
-    paxis2.Set(m_shower2->getPos().x(), m_shower2->getPos().y());
+  if(m_shower1->getTowerID()[0][0]==Cyber::CaloUnit::System_Barrel){
+
+    if(m_shower1->getSlayer()==1){ //For V-bars
+      paxis1.Set(m_shower1->getPos().x(), m_shower1->getPos().y());
+      paxis2.Set(m_shower2->getPos().x(), m_shower2->getPos().y());
+    }
+    else{
+      float rotAngle = -m_shower1->getTowerID()[0][1]*TMath::TwoPi()/Cyber::CaloUnit::Nmodule;
+      TVector3 raw_pos1 = m_shower1->getPos(); 
+      TVector3 raw_pos2 = m_shower2->getPos(); 
+      raw_pos1.RotateZ(rotAngle); raw_pos2.RotateZ(rotAngle);
+      paxis1.Set(raw_pos1.y(), raw_pos1.z());
+      paxis2.Set(raw_pos2.y(), raw_pos2.z());
+   
+      //paxis1.Set( sqrt(m_shower1->getPos().x()*m_shower1->getPos().x()+m_shower1->getPos().y()*m_shower1->getPos().y()), m_shower1->getPos().z());
+      //paxis2.Set( sqrt(m_shower2->getPos().x()*m_shower2->getPos().x()+m_shower2->getPos().y()*m_shower2->getPos().y()), m_shower2->getPos().z());
+    }
+  }
+
+  else if(m_shower1->getTowerID()[0][0]==Cyber::CaloUnit::System_Endcap){
+    if(m_shower1->getSlayer()==1){
+      paxis1.Set(m_shower1->getPos().z(), m_shower1->getPos().y());
+      paxis2.Set(m_shower2->getPos().z(), m_shower2->getPos().y());
+    }
+    else{
+      paxis1.Set(m_shower1->getPos().z(), m_shower1->getPos().x());
+      paxis2.Set(m_shower2->getPos().z(), m_shower2->getPos().x());
+    }
   }
+
   else{
-    float rotAngle = -m_shower1->getTowerID()[0][0]*TMath::TwoPi()/Cyber::CaloUnit::Nmodule;
-    TVector3 raw_pos1 = m_shower1->getPos(); 
-    TVector3 raw_pos2 = m_shower2->getPos(); 
-    raw_pos1.RotateZ(rotAngle); raw_pos2.RotateZ(rotAngle);
-    paxis1.Set(raw_pos1.y(), raw_pos1.z());
-    paxis2.Set(raw_pos2.y(), raw_pos2.z());
-
-    //paxis1.Set( sqrt(m_shower1->getPos().x()*m_shower1->getPos().x()+m_shower1->getPos().y()*m_shower1->getPos().y()), m_shower1->getPos().z());
-    //paxis2.Set( sqrt(m_shower2->getPos().x()*m_shower2->getPos().x()+m_shower2->getPos().y()*m_shower2->getPos().y()), m_shower2->getPos().z());
+    std::cout<<"Error in ConeClustering2DAlg: Unknown system ID: "<<m_shower1->getTowerID()[0][0]<<std::endl;
   }
 
   return paxis2 - paxis1;

Reconstruction/RecPFACyber/src/Algorithm/EnergyTimeMatchingAlg.cpp

@@ -673,18 +673,24 @@ vector<vector<double>> EnergyTimeMatchingAlg::GetClusterChi2Map( std::vector<std
   double wi_E = settings.map_floatPars["chi2Wi_E"]/(settings.map_floatPars["chi2Wi_E"] + settings.map_floatPars["chi2Wi_T"]);
   double wi_T = settings.map_floatPars["chi2Wi_T"]/(settings.map_floatPars["chi2Wi_E"] + settings.map_floatPars["chi2Wi_T"]);
 
+  //Rotate angle and tower center position for ECAL Barrel
   TVector3 m_vec(0,0,0);
   double rotAngle = -999;
   TVector3 Ctower(0,0,0);
+  int system = -1;
   for(int ish=0; ish<barShowerUCol.size(); ish++){
     if(barShowerUCol[ish].size()==0) continue;
     rotAngle = -(barShowerUCol[ish][0]->getBars())[0]->getModule()*TMath::TwoPi()/Cyber::CaloUnit::Nmodule;
     Ctower.SetX( (barShowerUCol[ish][0]->getBars())[0]->getPosition().x() );
     Ctower.SetY( (barShowerUCol[ish][0]->getBars())[0]->getPosition().y() );
+    system = (barShowerUCol[ish][0]->getBars())[0]->getSystem();
+    break;
   }
   for(int ish=0; ish<barShowerVCol.size(); ish++){
     if(barShowerVCol[ish].size()==0) continue;
     Ctower.SetZ( (barShowerVCol[ish][0]->getBars())[0]->getPosition().z() );
+    system = (barShowerVCol[ish][0]->getBars())[0]->getSystem();
+    break;
   }
   Ctower.RotateZ(rotAngle);
 
@@ -718,13 +724,27 @@ vector<vector<double>> EnergyTimeMatchingAlg::GetClusterChi2Map( std::vector<std
       double Ex = showerX->getEnergy();
       double Ey = showerY->getEnergy();
       double chi2_E = pow(fabs(Ex-Ey)/settings.map_floatPars["sigmaE"], 2);
-      double PosTx = C*(showerY->getT1()-showerY->getT2())/(2*settings.map_floatPars["nMat"]) + showerY->getPos().z();
-      double chi2_tx = pow( fabs(PosTx-showerX->getPos().z())/settings.map_floatPars["sigmaPos"], 2 );
 
-      double PosTy = C*(showerX->getT1()-showerX->getT2())/(2*settings.map_floatPars["nMat"]);
-      m_vec = showerY->getPos();
-      m_vec.RotateZ(rotAngle);
-      double chi2_ty = pow( fabs(PosTy - (m_vec-Ctower).x() )/settings.map_floatPars["sigmaPos"], 2);
+      double PosTx, chi2_tx, PosTy, chi2_ty;
+      if(system == Cyber::CaloUnit::System_Barrel){
+        PosTx = C*(showerY->getT1()-showerY->getT2())/(2*settings.map_floatPars["nMat"]) + showerY->getPos().z();
+        chi2_tx = pow( fabs(PosTx-showerX->getPos().z())/settings.map_floatPars["sigmaPos"], 2 );
+
+        PosTy = C*(showerX->getT1()-showerX->getT2())/(2*settings.map_floatPars["nMat"]);
+        m_vec = showerY->getPos();
+        m_vec.RotateZ(rotAngle);
+        chi2_ty = pow( fabs(PosTy - (m_vec-Ctower).x() )/settings.map_floatPars["sigmaPos"], 2);
+      }
+      else if(system == Cyber::CaloUnit::System_Endcap){
+        PosTx = C*(showerY->getT1()-showerY->getT2())/(2*settings.map_floatPars["nMat"]) + showerY->getPos().x();
+        chi2_tx = pow( fabs(PosTx-showerX->getPos().x())/settings.map_floatPars["sigmaPos"], 2 );
+
+        PosTy = C*(showerX->getT1()-showerX->getT2())/(2*settings.map_floatPars["nMat"]) + showerX->getPos().y();
+        chi2_ty = pow( fabs(PosTy - showerY->getPos().y() )/settings.map_floatPars["sigmaPos"], 2);
+      }
+      else{
+        std::cout<<"ERROR in EnergyTimeMatchingAlg: Unknown system ID: "<<system<<endl;
+      }
 
       if(chi2_E<min_chi2E) min_chi2E=chi2_E;
       if(chi2_tx<min_chi2tx) min_chi2tx=chi2_tx;
@@ -1071,31 +1091,57 @@ StatusCode EnergyTimeMatchingAlg::GetMatchedShowersL0( const Cyber::Calo1DCluste
   //if(barShowerU->getTowerID().size()==0) { barShowerU->setIDInfo(); }
 
   int _layer = barShowerU->getDlayer();
-  int _module = barShowerU->getTowerID()[0][0];
-  float rotAngle = -_module*TMath::TwoPi()/Cyber::CaloUnit::Nmodule;
-
-  for(int ibar=0;ibar<NbarsX;ibar++){
-    double En_x = barShowerU->getBars()[ibar]->getEnergy();
-    TVector3 m_vecx = barShowerU->getBars()[ibar]->getPosition();
-    m_vecx.RotateZ(rotAngle);
-
-    for(int jbar=0;jbar<NbarsY;jbar++){
-      double En_y = barShowerV->getBars()[jbar]->getEnergy();
-      TVector3 m_vecy = barShowerV->getBars()[jbar]->getPosition();
-      m_vecy.RotateZ(rotAngle);
-
-      TVector3 p_hit(m_vecy.x(), (m_vecx.y()+m_vecy.y())/2., m_vecx.z() );
-      p_hit.RotateZ(-rotAngle);
-      double m_Ehit = En_x*En_y/barShowerV->getEnergy() + En_x*En_y/barShowerU->getEnergy();
-      //Create new CaloHit
-      std::shared_ptr<Cyber::CaloHit> hit = std::make_shared<Cyber::CaloHit>();
-      hit->setcellID(_module, _layer);
-      hit->setPosition(p_hit);
-      hit->setEnergy(m_Ehit);
-      m_digiCol.push_back(hit.get());
-      m_bkCol.map_CaloHit["bkHit"].push_back( hit );
+  int _module = barShowerU->getTowerID()[0][1];
+  int _system = barShowerU->getTowerID()[0][0];
+  if(_system == Cyber::CaloUnit::System_Barrel){
+    float rotAngle = -_module*TMath::TwoPi()/Cyber::CaloUnit::Nmodule;
+
+    for(int ibar=0;ibar<NbarsX;ibar++){
+      double En_x = barShowerU->getBars()[ibar]->getEnergy();
+      TVector3 m_vecx = barShowerU->getBars()[ibar]->getPosition();
+      m_vecx.RotateZ(rotAngle);
+
+      for(int jbar=0;jbar<NbarsY;jbar++){
+        double En_y = barShowerV->getBars()[jbar]->getEnergy();
+        TVector3 m_vecy = barShowerV->getBars()[jbar]->getPosition();
+        m_vecy.RotateZ(rotAngle);
+
+        TVector3 p_hit(m_vecy.x(), (m_vecx.y()+m_vecy.y())/2., m_vecx.z() );
+        p_hit.RotateZ(-rotAngle);
+        double m_Ehit = En_x*En_y/barShowerV->getEnergy() + En_x*En_y/barShowerU->getEnergy();
+        //Create new CaloHit
+        std::shared_ptr<Cyber::CaloHit> hit = std::make_shared<Cyber::CaloHit>();
+        hit->setcellID(_system, _module, _layer);
+        hit->setPosition(p_hit);
+        hit->setEnergy(m_Ehit);
+        m_digiCol.push_back(hit.get());
+        m_bkCol.map_CaloHit["bkHit"].push_back( hit );
+      }
     }
   }
+  else{
+
+    for(int ibar=0;ibar<NbarsX;ibar++){
+      double En_x = barShowerU->getBars()[ibar]->getEnergy();
+      TVector3 m_vecx = barShowerU->getBars()[ibar]->getPosition();
+
+      for(int jbar=0;jbar<NbarsY;jbar++){
+        double En_y = barShowerV->getBars()[jbar]->getEnergy();
+        TVector3 m_vecy = barShowerV->getBars()[jbar]->getPosition();
+
+        TVector3 p_hit(m_vecx.x(), m_vecy.y(), (m_vecx.z()+m_vecy.z())/2. );
+        double m_Ehit = En_x*En_y/barShowerV->getEnergy() + En_x*En_y/barShowerU->getEnergy();
+        //Create new CaloHit
+        std::shared_ptr<Cyber::CaloHit> hit = std::make_shared<Cyber::CaloHit>();
+        hit->setcellID(_system, _module, _layer);
+        hit->setPosition(p_hit);
+        hit->setEnergy(m_Ehit);
+        m_digiCol.push_back(hit.get());
+        m_bkCol.map_CaloHit["bkHit"].push_back( hit );
+      }
+    }
+
+  }
 
   outsh->addUnit( barShowerU );
   outsh->addUnit( barShowerV );

Reconstruction/RecPFACyber/src/Algorithm/GlobalClusteringAlg.cpp

@@ -34,17 +34,22 @@ StatusCode GlobalClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
   //cout<<" When begin clustering: "<<ctime(&time_cb)<<endl;
 
   //Threshold and scale factor (Todo) for bars. 
+//double totE = 0.;
+//double totE_rest = 0.;
   for(int ibar=0; ibar<m_bars.size(); ibar++)
   {
     if(m_bars.at(ibar)->getEnergy()>settings.map_floatPars["unit_threshold"])
     {
       m_processbars.push_back(m_bars.at(ibar));
+      //totE += m_bars.at(ibar)->getEnergy();
     }
     else
     {
       m_restbars.push_back(m_bars.at(ibar));
+      //totE_rest += m_bars.at(ibar)->getEnergy();
     }
   }
+//cout<<"Input bar after threshold: "<<m_processbars.size()<<", totE "<<totE<<", rest E "<<totE_rest<<endl;
 
   //Clustering
   Clustering(m_processbars, m_1dclusters);
@@ -54,7 +59,6 @@ StatusCode GlobalClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
   m_datacol.map_1DCluster["bk1DCluster"].insert(m_datacol.map_1DCluster["bk1DCluster"].end(), m_1dclusters.begin(), m_1dclusters.end());
   m_datacol.map_HalfCluster["bkHalfCluster"].insert(m_datacol.map_HalfCluster["bkHalfCluster"].end(), m_halfclusters.begin(), m_halfclusters.end());
 
-
   std::vector<std::shared_ptr<Cyber::CaloHalfCluster>> m_HalfClusterV; m_HalfClusterV.clear();
   std::vector<std::shared_ptr<Cyber::CaloHalfCluster>> m_HalfClusterU; m_HalfClusterU.clear();
   for(int i=0; i<m_halfclusters.size() && m_halfclusters[i]; i++){
@@ -65,11 +69,11 @@ StatusCode GlobalClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
       m_HalfClusterV.push_back(m_halfclusters[i]);
   }
 
-//printf("  GlobalClustering: RestBarCol size %d, 1DCluster size %d, HalfCluster size (%d, %d) \n", m_restbars.size(), m_1dclusters.size(), m_HalfClusterU.size(), m_HalfClusterV.size());
-//for(int ic=0; ic<m_HalfClusterU.size(); ic++) cout<<m_HalfClusterU[ic]->getEnergy()<<'\t';
-//cout<<endl;
-//for(int ic=0; ic<m_HalfClusterV.size(); ic++) cout<<m_HalfClusterV[ic]->getEnergy()<<'\t';
-//cout<<endl;
+  //printf("  GlobalClustering: RestBarCol size %d, 1DCluster size %d, HalfCluster size (%d, %d) \n", m_restbars.size(), m_1dclusters.size(), m_HalfClusterU.size(), m_HalfClusterV.size());
+  //for(int ic=0; ic<m_HalfClusterU.size(); ic++) cout<<m_HalfClusterU[ic]->getEnergy()<<'\t';
+  //cout<<endl;
+  //for(int ic=0; ic<m_HalfClusterV.size(); ic++) cout<<m_HalfClusterV[ic]->getEnergy()<<'\t';
+  //cout<<endl;
 
 
   //Write results into DataCol.

Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp

@@ -7,7 +7,7 @@ using namespace Cyber;
 StatusCode HcalClusteringAlg::ReadSettings(Cyber::Settings& m_settings){
   settings = m_settings;
 
-  if(settings.map_stringPars.find("InputHCALHits")==settings.map_stringPars.end())         settings.map_stringPars["InputHCALHits"] = "HCALBarrel";
+  if(settings.map_stringVecPars.find("InputHCALHits")==settings.map_stringVecPars.end())         settings.map_stringVecPars["InputHCALHits"] = {"HCALBarrel", "HCALEndcaps"};
   if(settings.map_stringPars.find("OutputHCALClusters")==settings.map_stringPars.end())    settings.map_stringPars["OutputHCALClusters"] = "HCALCluster";
 
   //Set initial values
@@ -39,8 +39,10 @@ StatusCode HcalClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
 
   std::vector<Cyber::CaloHit*> m_hcalHits;
   m_hcalHits.clear();
-  for(int ih=0; ih<m_datacol.map_CaloHit[settings.map_stringPars["InputHCALHits"]].size(); ih++)
-    m_hcalHits.push_back( m_datacol.map_CaloHit[settings.map_stringPars["InputHCALHits"]][ih].get() );
+  for(int icl=0; icl<settings.map_stringVecPars["InputHCALHits"].size(); icl++){
+    for(int ih=0; ih<m_datacol.map_CaloHit[settings.map_stringVecPars["InputHCALHits"][icl]].size(); ih++)
+      m_hcalHits.push_back( m_datacol.map_CaloHit[settings.map_stringVecPars["InputHCALHits"][icl]][ih].get() );
+  }
 
   std::vector<std::shared_ptr<Cyber::Calo3DCluster>> m_clusterCol;  
   m_clusterCol.clear();

Reconstruction/RecPFACyber/src/Algorithm/LocalMaxFindingAlg.cpp

@@ -42,6 +42,7 @@ StatusCode LocalMaxFindingAlg::RunAlgorithm( CyberDataCol& m_datacol){
       m_datacol.map_1DCluster["bk1DCluster"].push_back(iter);
     }
     p_HalfClusU->at(iu).get()->setLocalMax(settings.map_stringPars["OutputLocalMaxName"], tmp_localMax);
+//printf("  HalfClusterU #%d: energy %.3f, localMax size %d \n", iu, p_HalfClusU->at(iu).get()->getEnergy(), tmp_localMax.size());
     m_1dClusCol.clear(); 
     ptr_localMax.clear();
   }
@@ -59,6 +60,7 @@ StatusCode LocalMaxFindingAlg::RunAlgorithm( CyberDataCol& m_datacol){
       m_datacol.map_1DCluster["bk1DCluster"].push_back(iter);
     }
     p_HalfClusV->at(iv).get()->setLocalMax(settings.map_stringPars["OutputLocalMaxName"], tmp_localMax);
+//printf("  HalfClusterV #%d: energy %.3f, localMax size %d \n", iv, p_HalfClusV->at(iv).get()->getEnergy(), tmp_localMax.size());
     m_1dClusCol.clear(); 
     ptr_localMax.clear();
   }
@@ -127,22 +129,36 @@ StatusCode LocalMaxFindingAlg::GetLocalMaxBar( std::vector<const Cyber::CaloUnit
 std::vector<const Cyber::CaloUnit*> LocalMaxFindingAlg::getNeighbors( const Cyber::CaloUnit* seed, std::vector<const Cyber::CaloUnit*>& barCol){
   std::vector<const Cyber::CaloUnit*> m_neighbor; m_neighbor.clear();
   for(int i=0;i<barCol.size();i++){
-    bool fl_neighbor = false; 
-    if( seed->getModule()==barCol[i]->getModule() && 
-        seed->getStave()==barCol[i]->getStave() && 
-        seed->getDlayer()==barCol[i]->getDlayer() &&
-        seed->getSlayer()==barCol[i]->getSlayer() &&
-        abs( seed->getBar()-barCol[i]->getBar() )==1 ) fl_neighbor=true;
-    else if( seed->getStave()==barCol[i]->getStave() && 
-             ( ( seed->getModule()-barCol[i]->getModule()==1 && seed->isAtLowerEdgePhi() && barCol[i]->isAtUpperEdgePhi() ) ||
-               ( barCol[i]->getModule()-seed->getModule()==1 && seed->isAtUpperEdgePhi() && barCol[i]->isAtLowerEdgePhi() ) ) ) fl_neighbor=true;
-    else if( seed->getModule()==barCol[i]->getModule() && 
-             ( ( seed->getStave()-barCol[i]->getStave()==1 && seed->isAtLowerEdgeZ() && barCol[i]->isAtUpperEdgeZ() ) || 
-               ( barCol[i]->getStave()-seed->getStave()==1 && seed->isAtUpperEdgeZ() && barCol[i]->isAtLowerEdgeZ() ) ) ) fl_neighbor=true;
-
-    if(fl_neighbor) m_neighbor.push_back(barCol[i]);
+    //bool fl_neighbor = false; 
+    //if( seed->getSystem()==barCol[i]->getSystem() &&
+    //    seed->getModule()==barCol[i]->getModule() && 
+    //    seed->getStave()==barCol[i]->getStave() && 
+    //    seed->getDlayer()==barCol[i]->getDlayer() &&
+    //    seed->getSlayer()==barCol[i]->getSlayer() &&
+    //    abs( seed->getBar()-barCol[i]->getBar() )==1 ) fl_neighbor=true;
+    //else if( seed->getSystem()==barCol[i]->getSystem() && seed->getStave()==barCol[i]->getStave() && 
+    //         ( ( seed->getModule()-barCol[i]->getModule()==1 && seed->isAtLowerEdgePhi() && barCol[i]->isAtUpperEdgePhi() ) ||
+    //           ( barCol[i]->getModule()-seed->getModule()==1 && seed->isAtUpperEdgePhi() && barCol[i]->isAtLowerEdgePhi() ) ) ) fl_neighbor=true;
+    //else if( seed->getSystem()==barCol[i]->getSystem() && seed->getModule()==barCol[i]->getModule() && 
+    //         ( ( seed->getStave()-barCol[i]->getStave()==1 && seed->isAtLowerEdgeZ() && barCol[i]->isAtUpperEdgeZ() ) || 
+    //           ( barCol[i]->getStave()-seed->getStave()==1 && seed->isAtUpperEdgeZ() && barCol[i]->isAtLowerEdgeZ() ) ) ) fl_neighbor=true;
+
+    bool fl_neighbor = seed->isNeighbor( barCol[i] );
+    if(fl_neighbor){ 
+      if(seed->getSystem()==CaloUnit::System_Barrel && seed->getSlayer()==0 && seed->getModule()!=barCol[i]->getModule()) continue;
+      if(seed->getSystem()==CaloUnit::System_Barrel && seed->getSlayer()==1 && seed->getStave()!=barCol[i]->getStave())   continue;
+      if(seed->getSystem()==CaloUnit::System_Endcap && seed->getSlayer()==0 && seed->getPart()!=barCol[i]->getPart() )    continue;
+      if(seed->getSystem()==CaloUnit::System_Endcap && seed->getSlayer()==1 && seed->getStave()!=barCol[i]->getStave() )  continue;
+      m_neighbor.push_back(barCol[i]);
+    }
+  }
+  if(m_neighbor.size()>2){ 
+    std::cout<<"WARNING: "<<m_neighbor.size()<<" hits in neighborCol!!"<<std::endl;
+    printf("Seed cellID: ( %d, %d, %d, %d, %d, %d, %d), position (%.3f, %.3f, %.3f) \n", seed->getSystem(), seed->getModule(), seed->getStave(), seed->getPart(), seed->getDlayer(), seed->getSlayer(), seed->getBar(), seed->getPosition().x(), seed->getPosition().y(), seed->getPosition().z());
+    for(int i=0; i<m_neighbor.size(); i++)
+      printf("Bar #%d cellID: ( %d, %d, %d, %d, %d, %d, %d), position (%.3f, %.3f, %.3f) \n", i, m_neighbor[i]->getSystem(), m_neighbor[i]->getModule(), m_neighbor[i]->getStave(), m_neighbor[i]->getPart(), m_neighbor[i]->getDlayer(), m_neighbor[i]->getSlayer(), m_neighbor[i]->getBar(), m_neighbor[i]->getPosition().x(), m_neighbor[i]->getPosition().y(), m_neighbor[i]->getPosition().z());
+  
   }
-  if(m_neighbor.size()>2) std::cout<<"WARNING: more than 2 hits in neighborCol!!"<<std::endl;
 
   return m_neighbor;
 }

Reconstruction/RecPFACyber/src/Algorithm/PFOReclusteringAlg.cpp

@@ -68,7 +68,6 @@ StatusCode PFOReclusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){
   // }
   // cout<<"-----Charged cluster Ecal total energy: "<<totE_Ecal<<", Hcal total energy: "<<totE_Hcal<<endl;
 
-
   //If P_trk > E_cluster, merge nearby neutral PFO into the charged. 
   ReCluster_MergeToChg(m_chargedPFOs, m_neutralPFOs);
 
@@ -523,7 +522,7 @@ StatusCode PFOReclusteringAlg::ReCluster_SplitFromChg( std::vector< std::shared_
         std::shared_ptr<Cyber::Calo2DCluster> m_2dclus = std::make_shared<Cyber::Calo2DCluster>();
         m_2dclus->addShowerU(m_1dclus_u.get());
         m_2dclus->addShowerV(m_1dclus_v.get());
-        m_2dclus->addTowerID(0., 0., 0.);
+        m_2dclus->addTowerID(0., 0., 0., 0.);
         m_2dclus->setPos(tmp_pos);
    
         // -- Create Half cluster U/V