Digitization/DigiCalo/src/HcalDigiAlg.cpp

@@ -83,9 +83,9 @@ StatusCode HcalDigiAlg::initialize()
   for(auto& link : name_CaloMCPAsso){
     if(!link.empty())
       _outputCaloMCPAssoCol.push_back( new CaloParticleAssoType(link, Gaudi::DataHandle::Writer, this) );
-  }  
-
-
+  }
+  
+  // --- Ntuple
   if(_writeNtuple){
     std::string s_outfile = _filename;
     m_wfile = new TFile(s_outfile.c_str(), "recreate");
@@ -266,7 +266,7 @@ StatusCode HcalDigiAlg::execute()
       //printf("  Step #%d: En %.2f, rotate angle %.2f, tile pos after rotation (%.2f, %.2f, %.2f), ", iCont, conb.getEnergy(), rotPhi, rot_tilepos.x(), rot_tilepos.y(), rot_tilepos.z());
       //printf("rel pos after rotation (%.2f, %.2f, %.2f) \n", (rot_steppos-rot_tilepos).x(), (rot_steppos-rot_tilepos).y(), (rot_steppos-rot_tilepos).z());
       //printf("  Project to bin (%d, %d), LY %.3f \n", ibinx, ibiny, GSTileResMap->GetBinContent( ibinx, ibiny ));
-
+      
           Npe_att += conb.getEnergy() / _MIPCali * GSTileResMap->GetBinContent( ibinx, ibiny ) * fLY_tempScale;     
           m_step_LY.push_back(GSTileResMap->GetBinContent( ibinx, ibiny ));
         }

Digitization/DigiSimple/CMakeLists.txt

@@ -5,6 +5,8 @@ gaudi_add_module(SimpleDigi
                          src/voxel.cpp
                          src/CylinderDigiAlg.cpp
                          src/SiTrackerDigiAlg.cpp
+                         src/TPCPixelDigiAlg.cpp
+                         src/TPCPixelDigiTool.cpp
                          src/SmearDigiTool.cpp
                  LINK GearSvc
                       DigiTool

Digitization/DigiSimple/src/TPCPixelDigiAlg.cpp

+#include "TPCPixelDigiAlg.h"
+
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/Vector3D.h"
+#include "DDRec/ISurface.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+
+DECLARE_COMPONENT(TPCPixelDigiAlg)
+
+TPCPixelDigiAlg::TPCPixelDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc) {
+  // Input collections
+  declareProperty("SimTrackHitCollection", m_inputColHdls, "Handle of the Input SimTrackerHit collection");
+
+
+  // Output collections
+  declareProperty("TPCTrackerHits", m_outputColHdls, "Handle of the output TrackerHit collection");
+  declareProperty("TPCTrackerHitAss", m_assColHdls, "Handle of the Association collection between SimTrackerHit and TrackerHit");
+}
+
+StatusCode TPCPixelDigiAlg::initialize() {
+  m_digiTool = ToolHandle<IDigiTool>(m_digiToolName.value());
+
+  info() << "DigiTool " << m_digiTool.typeAndName() << " found" << endmsg;
+
+  info() << "TPCPixelDigiAlg::initialized" << endmsg;
+  return GaudiAlgorithm::initialize();
+}
+
+
+StatusCode TPCPixelDigiAlg::execute(){
+  auto trkCol = m_outputColHdls.createAndPut();
+  auto assCol = m_assColHdls.createAndPut();
+
+  const edm4hep::SimTrackerHitCollection* simCol = nullptr;
+  try {
+    simCol = m_inputColHdls.get();
+  }
+  catch(GaudiException &e){
+    debug() << "Collection " << m_inputColHdls.fullKey() << " test is unavailable in event " << m_nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  if (simCol->size() == 0) return StatusCode::SUCCESS;
+  debug() << m_inputColHdls.fullKey() << " has cluser "<< simCol->size() << endmsg;
+
+  m_digiTool->Call(simCol, trkCol, assCol);
+
+  debug() << "Created " << trkCol->size() << " hits, "<<simCol->size()<<" clusters"
+          << simCol->size() - trkCol->size() << " hits got dismissed for being out of boundary"
+          << endmsg;
+
+  m_nEvt++;
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiAlg::finalize(){
+  info() << "Processed " << m_nEvt << " events " << endmsg;
+  return GaudiAlgorithm::finalize();
+}

Digitization/DigiSimple/src/TPCPixelDigiAlg.h

+#ifndef TPCPixelDigiAlg_H
+#define TPCPixelDigiAlg_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/NTuple.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/SimPrimaryIonizationClusterCollection.h"
+#include "edm4hep/SimPrimaryIonizationCluster.h"
+
+#include "DigiTool/IDigiTool.h"
+
+class TPCPixelDigiAlg : public GaudiAlgorithm{
+ public:
+  
+  TPCPixelDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+ 
+  virtual StatusCode initialize();
+  virtual StatusCode execute();
+  virtual StatusCode finalize();
+ 
+protected:
+  // Input collections
+  // DataHandle<edm4hep::SimPrimaryIonizationClusterCollection>            m_inputColHdls{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Reader, this};
+  // DataHandle<edm4hep::SimPrimaryIonizationClusterCollection>            m_inputIonClusterCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::SimTrackerHitCollection>            m_inputColHdls{"TPCCollection", Gaudi::DataHandle::Reader, this};
+  // Output collections
+  DataHandle<edm4hep::TrackerHitCollection>               m_outputColHdls{"TPCTrackerHits", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assColHdls{"TPCTrackerHitsAss", Gaudi::DataHandle::Writer, this};
+
+  Gaudi::Property<std::string> m_digiToolName{this, "DigiTool", "TPCPixelDigiTool"};
+
+  ToolHandle<IDigiTool> m_digiTool;
+  int m_nEvt=0;
+};
+#endif

Digitization/DigiSimple/src/TPCPixelDigiTool.cpp

+#include "TPCPixelDigiTool.h"
+#include "k4FWCore/DataHandle.h"
+#include "DataHelper/TrackerHitHelper.h"
+#include "DetIdentifier/CEPCConf.h"
+#include "DetInterface/IGeomSvc.h"
+#include "DataHelper/Circle.h"
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/ISurface.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include <gear/GEAR.h>
+#include <gear/TPCParameters.h>
+
+#include "CLHEP/Vector/ThreeVector.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include <math.h>
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <map>
+#include <cmath>
+#include <algorithm>
+#include <array>
+
+// #include "DataHelper/Circle.h"
+#include "DataHelper/SimpleHelix.h"
+#include "DataHelper/LCCylinder.h"
+#include "EventSeeder/IEventSeeder.h"
+#include <stdexcept>
+#include "GearSvc/IGearSvc.h"
+#include <UTIL/BitField64.h>
+#include <edm4hep/MCParticle.h>
+#include "edm4hep/SimTrackerHit.h"
+
+#include "CLHEP/Vector/TwoVector.h"
+#include "UTIL/ILDConf.h"
+
+#define TRKHITNCOVMATRIX 6
+// using namespace lcio ;
+
+DECLARE_COMPONENT(TPCPixelDigiTool)
+
+StatusCode TPCPixelDigiTool::initialize()
+{
+
+  m_geosvc = service<IGeomSvc>("GeomSvc");
+  if (!m_geosvc)
+  {
+    error() << "Failed to get the GeomSvc" << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+
+  debug() << "Readout name: " << m_readoutName << endmsg;
+  m_decoder = m_geosvc->getDecoder(m_readoutName);
+  if (!m_decoder)
+  {
+    error() << "Failed to get the decoder. " << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  _gear = service<IGearSvc>("GearSvc");
+  if (!_gear)
+  {
+    error() << "Failed to find GearSvc ..." << endmsg;
+    return StatusCode::FAILURE;
+  }
+  _GEAR = _gear->getGearMgr();
+  _f1 = new TF1("fp", "[0]*pow(1+[1],1+[1])*pow(x/[2],[1])*exp(-(1+[1])*x/[2])/TMath::Gamma(1+[1])", 0, 10000);
+  // _f1->SetParameters(3802, 0.487, 2092);
+  _f1->SetParameters(_polyParas[0], _polyParas[1], _polyParas[2]);
+  x_diffusion_p0=0.;
+  y_diffusion_p0=0.;
+  x_diffusion_p1=0.;
+  y_diffusion_p1=0.;
+  t_diffusion_p0 = 0.;
+  t_diffusion_p1 = 0.;
+
+  switch (_magnetic) {
+    case 2:
+        x_diffusion_p0 = 0.001992;
+        x_diffusion_p1 = 0.004387;
+        y_diffusion_p0 = 0.001877;
+        y_diffusion_p1 = 0.004398;
+		t_diffusion_p0 = 0.0002519;
+		t_diffusion_p1 = 0.00332;
+        tpcZRes = tpcTRes * 80; // drift speed 80um/us
+        break;
+
+    case 3:
+        x_diffusion_p0 = 0.001403;
+        x_diffusion_p1 = 0.003011;
+        y_diffusion_p0 = 0.001604;
+        y_diffusion_p1 = 0.003042;
+		t_diffusion_p0 = 5.41279e-05;
+		t_diffusion_p1 = 3.28842e-03;
+        tpcZRes = tpcTRes * 80; // drift speed 80um/us
+        break;
+
+    default:
+        x_diffusion_p0 = 0.0;
+        x_diffusion_p1 = 0.0;
+        y_diffusion_p0 = 0.0;
+        y_diffusion_p1 = 0.0;
+		t_diffusion_p0 = 5.41279e-05;
+		t_diffusion_p1 = 3.28842e-03;
+        tpcZRes = tpcTRes * 80; // drift speed 80um/us
+        break;
+}
+
+ 
+
+  info() << "initialized" << endmsg;
+  return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiTool::Call(const edm4hep::SimTrackerHitCollection *simCol, edm4hep::TrackerHitCollection *hitCol,
+                                  edm4hep::MCRecoTrackerAssociationCollection *assCol)
+{
+
+  edm4hep::SimTrackerHit simhit;
+  StatusCode sc = Call(simhit, hitCol, assCol);
+  if (sc.isFailure())
+    return sc;
+  return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection *hitCol, edm4hep::MCRecoTrackerAssociationCollection *assCol)
+{
+  DataHandle<edm4hep::SimPrimaryIonizationClusterCollection> m_inputIonClusterCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::SimTrackerHitCollection> m_inputSimCol{m_readoutName, Gaudi::DataHandle::Reader, this};
+  UTIL::BitField64 *_cellid_encoder = new UTIL::BitField64(CEPCConf::DetEncoderString::getStringRepresentation());
+  float weight = 1.0;
+
+  // double bReso = _diffRPhi * _diffRPhi;
+  const gear::TPCParameters &gearTPC = _GEAR->getTPCParameters();
+  const gear::PadRowLayout2D &padLayout = gearTPC.getPadLayout();
+  // gear::Vector2D padCoord;
+  const edm4hep::SimTrackerHitCollection *simCol = nullptr;
+  try
+  {
+    simCol = m_inputSimCol.get();
+  }
+  catch (GaudiException &e)
+  {
+    debug() << "Collection " << m_inputSimCol.fullKey() << " is unavailable in event " << endmsg;
+    return StatusCode::SUCCESS;
+  }
+
+  const edm4hep::SimPrimaryIonizationClusterCollection *clusterCol = nullptr;
+  try
+  {
+    clusterCol = m_inputIonClusterCol.get();
+  }
+  catch (GaudiException &e)
+  {
+    debug() << "Collection " << m_inputIonClusterCol.fullKey() << " is unavailable in event " << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  std::map<int, edm4hep::SimTrackerHit> simhitMap;                      // A map to save the sim tracker hit, key is TPC layer, value is 3D position of simhit
+  std::map<int, std::tuple<double, double, double>> digihitMap; // A map to save the XYZ resolutions and  Z position(key is pad index, values are z and rphi resolutions and Z position)
+  std::map<int, int> tpcPixelQ;                                         // A map to save number of electrons in same pixel (key is pixel index, value is number of eles)
+  std::map<int, int> tpcPixelFlag;                                      // A map to mark if a pixel already has a hit (key is pixel index), only used in merging case
+  std::map<int, std::vector<int>> tpcPixelQMap;
+  std::map<int, std::vector<std::tuple<float, float>>> TimeZMap;        // A map to save time and z direction infos
+  
+  for (auto sim : *simCol)
+  { // loop all sim tracker hits, save them in a map, key is TPC layer, value is 3D position of simhit
+    CLHEP::Hep3Vector SimPoint(sim.getPosition()[0], sim.getPosition()[1], sim.getPosition()[2]);
+    int simPadIndex = padLayout.getNearestPad(SimPoint.perp(), SimPoint.phi());
+    int iRowSimHit = padLayout.getRowNumber(simPadIndex);
+    simhitMap[iRowSimHit] = sim;
+  }
+
+  int nclu = 0;
+
+
+  int padIndex;
+  int iRowHit;
+  int poylaRandom;
+  double x_amp_diffusion;
+  double y_amp_diffusion;
+  CLHEP::Hep3Vector point;
+  CLHEP::Hep3Vector thisPoint;
+  gear::Vector2D padCoord;
+  edm4hep::Vector3d posAfterSmear;
+  edm4hep::Vector3d pos;
+
+  int nIonEle;
+  debug() << " Total clusters: " << clusterCol->size() << endmsg;
+  // for (int i=0; i < _nPads;++i){
+  // padCoord=padLayout.getPadCenter(int(i));
+  // std::cout<<" check padxy "<<padCoord[0] * cos(padCoord[1])<<" "<<padCoord[0] * sin(padCoord[1])<<std::endl;
+  // std::cout<<" check padrphi "<<padCoord[0] <<" "<<padCoord[1] <<std::endl;
+  // }
+  
+
+  for (edm4hep::SimPrimaryIonizationCluster cluster : *clusterCol)
+  { // Loop all clusters
+    // break;
+    debug() << " start cluster: " << nclu << " pos:" << cluster.getPosition() << " size: " << cluster.electronPosition_size() << endmsg;
+    nclu += 1;
+    nIonEle = 0;
+
+    for (std::size_t i = 0; i < cluster.electronPosition_size(); ++i)
+    { // loop all electrons in each clusters
+
+      pos = *(cluster.electronPosition_begin() + i);
+      float eleTime = *(cluster.electronTime_begin() + i) / 1000.;//time in us
+
+      if (_MaxIonEles > 0 && nIonEle > _MaxIonEles)
+      {
+        debug() << "Maximum eles achived, break the loop" << endmsg;
+        break;
+      }
+	  driftLength = gearTPC.getMaxDriftLength() - (fabs(pos.z));
+      tpcTRes = (t_diffusion_p0 +t_diffusion_p1 * sqrt(driftLength / 10.)); 
+
+      tpcXRes = x_diffusion_p0 + x_diffusion_p1 * sqrt(driftLength);
+      tpcYRes = y_diffusion_p0 + y_diffusion_p1 * sqrt(driftLength);
+      double randx = gRandom->Gaus(pos.x, tpcXRes);
+      double randy = gRandom->Gaus(pos.y, tpcYRes);
+      // double randz = gRandom->Gaus(pos.z, tpcZRes);
+      double randt = eleTime + gRandom->Gaus(driftLength / 80, tpcTRes); // unit in us
+      double randz = randt*80;
+      //   
+      point.set(randx, randy, randz);
+
+      if (fabs(point.z()) > gearTPC.getMaxDriftLength())
+        point.setZ((fabs(point.z()) / point.z()) * gearTPC.getMaxDriftLength());
+      debug() << "resolution, x: " << tpcXRes << " y: " << tpcYRes << " z: " << tpcZRes << endmsg;
+      debug() << "Max TPC length: " << gearTPC.getMaxDriftLength() << " posz: " << fabs(pos.z) << " driftLength: " << driftLength << endmsg;
+      debug() << "smear:" << pos.x << "," << pos.y << "," << pos.z << " -> " << point << endmsg;
+      poylaRandom = _f1->GetRandom();
+      // poylaRandom = 1;
+      debug() << "polyRandom: " << poylaRandom << endmsg;
+      for (int j = 0; j < poylaRandom; ++j) // Do amplication
+      {
+
+        x_amp_diffusion = gRandom->Gaus(point.x(), _ampDiffusion);
+        y_amp_diffusion = gRandom->Gaus(point.y(), _ampDiffusion);
+        thisPoint.set(x_amp_diffusion, y_amp_diffusion, point.z());
+
+        // following lines are the method to find the pixel center
+        padIndex = padLayout.getNearestPad(thisPoint.perp(), thisPoint.phi());
+        padCoord = padLayout.getPadCenter(int(padIndex));
+        posAfterSmear.x = padCoord[0] * cos(padCoord[1]);
+        posAfterSmear.y = padCoord[0] * sin(padCoord[1]);
+        posAfterSmear.z = point.z();
+
+        digihitMap[padIndex] = std::make_tuple(tpcXRes, tpcYRes, tpcZRes); 
+        if (tpcPixelQ.find(padIndex) != tpcPixelQ.end())
+        {
+          tpcPixelQ[padIndex] += 1;
+        }
+        else
+        {
+          tpcPixelQ[padIndex] = 1;
+        }
+        if (TimeZMap[padIndex].empty())
+        {
+          TimeZMap[padIndex].push_back(std::make_tuple(randt,point.z()));
+        }
+        else
+        {
+          if (!_saveAllTime)
+          {
+
+            if (randt > std::get<0>(TimeZMap[padIndex].back()) + _deadTime)
+            {
+              info() << "now the size of pad " << padIndex << " is " << TimeZMap[padIndex].size() << " Q is  " << tpcPixelQ[padIndex] << endmsg;
+              info() << "now time is: " << randt << " previous time is: " << std::get<0>(TimeZMap[padIndex].back()) << endmsg;
+
+              TimeZMap[padIndex].push_back(std::make_tuple(randt,point.z()));
+              tpcPixelQMap[padIndex].push_back(tpcPixelQ[padIndex]);
+              tpcPixelQ[padIndex]=0;
+            }
+          }
+          else
+          {
+            auto it = TimeZMap[padIndex].rbegin();
+            bool found = false;
+            for (; it != TimeZMap[padIndex].rend(); ++it)
+            {
+              float diff = std::abs(std::get<0>(*it) - randt);
+              if (diff < _minimumTime)
+              {
+                found = true;
+                break;
+              }
+            }
+
+            if (!found)
+            {
+              TimeZMap[padIndex].emplace_back(randt,point.z());
+              tpcPixelQMap[padIndex].push_back(tpcPixelQ[padIndex]);
+              tpcPixelQ[padIndex]=0;
+            }
+          }
+        }
+
+      }
+      nIonEle += 1;
+    }
+  }
+  debug() << "create " << digihitMap.size() << " hits in different pixels before merging" << endmsg;
+  for (auto digiHit = digihitMap.begin(); digiHit != digihitMap.end(); ++digiHit)
+    {
+      padIndex = digiHit->first;
+      std::vector<std::tuple<float, float>> TimeZVec = TimeZMap[padIndex];
+      tpcPixelQMap[padIndex].push_back(tpcPixelQ[padIndex]);
+      std::vector<int> QVec = tpcPixelQMap[padIndex];
+      
+      tpcXRes = std::get<0>(digiHit->second);
+      tpcYRes = std::get<1>(digiHit->second);
+      tpcZRes = std::get<2>(digiHit->second);
+      padCoord = padLayout.getPadCenter(int(digiHit->first));
+
+      double r = padCoord[0];
+      double phi = padCoord[1];
+      padIndex = digiHit->first;
+      iRowHit = padLayout.getRowNumber(padIndex);
+      debug() << "pad center: " << r * cos(phi) << " " << r * sin(phi) << endmsg;
+      debug() << "pad r: " << r << " phi: " << phi << " Q: " << tpcPixelQ[padIndex] << endmsg;
+      posAfterSmear.x = (r * cos(phi));
+      posAfterSmear.y = (r * sin(phi));
+      
+      for (size_t i = 0; i < TimeZVec.size(); ++i) {
+        posAfterSmear.z = std::get<1>(TimeZVec[i]);
+        auto it = simhitMap.find(iRowHit);
+        if (it != simhitMap.end() && !_skipAss)
+      {
+        
+        simhit = it->second;
+        auto outhit = hitCol->create();
+        auto ass = assCol->create();
+        outhit.setPosition(posAfterSmear);
+        (*_cellid_encoder)[lcio::ILDCellID0::subdet] = lcio::ILDDetID::TPC;
+        (*_cellid_encoder)[lcio::ILDCellID0::layer] = iRowHit;
+        (*_cellid_encoder)[lcio::ILDCellID0::module] = 0;
+
+        (*_cellid_encoder)[lcio::ILDCellID0::side] = lcio::ILDDetID::barrel;
+        outhit.setCellID(_cellid_encoder->lowWord());
+        std::array<float, TRKHITNCOVMATRIX> covMat = {float(tpcXRes * tpcXRes),
+                                                      float(tpcYRes * tpcXRes),
+                                                      float(tpcYRes * tpcYRes),
+                                                      float(tpcZRes * tpcXRes),
+                                                      float(tpcZRes * tpcYRes),
+                                                      float(tpcZRes * tpcZRes)};
+
+        outhit.setCovMatrix(covMat);
+        outhit.setEDep(tpcPixelQMap[padIndex][i]);
+        outhit.setTime(std::get<0>(TimeZVec[i]));
+        info() << "out hit, x: " << outhit.getPosition()[0] << " y: " << outhit.getPosition()[1] << " z: " << outhit.getPosition()[2] << " Q: " << outhit.getEDep() << " Time: "<< outhit.getTime() << endmsg;
+        ass.setSim(simhit);
+        ass.setRec(outhit);
+        ass.setWeight(weight);
+      }
+      else if (_skipAss)// do not match the TPCCollection and TPCTrackerHits, all the TPCTrackerHitsAss has been set to 0,0,0
+      {
+        debug() << " skip ass" << endmsg;
+        debug() << " SimTrackerHit " << simhit.getPosition()[0] << endmsg;
+
+        auto outhit = hitCol->create();
+        auto ass = assCol->create();
+        outhit.setPosition(posAfterSmear);
+        (*_cellid_encoder)[lcio::ILDCellID0::subdet] = lcio::ILDDetID::TPC;
+        (*_cellid_encoder)[lcio::ILDCellID0::layer] = iRowHit;
+        (*_cellid_encoder)[lcio::ILDCellID0::module] = 0;
+
+        (*_cellid_encoder)[lcio::ILDCellID0::side] = lcio::ILDDetID::barrel;
+        outhit.setCellID(_cellid_encoder->lowWord());
+        std::array<float, TRKHITNCOVMATRIX> covMat = {float(tpcXRes * tpcXRes),
+                                                      float(tpcYRes * tpcXRes),
+                                                      float(tpcYRes * tpcYRes),
+                                                      float(tpcZRes * tpcXRes),
+                                                      float(tpcZRes * tpcYRes),
+                                                      float(tpcZRes * tpcZRes)};
+
+        outhit.setCovMatrix(covMat);
+        outhit.setEDep(tpcPixelQMap[padIndex][i]);
+        outhit.setTime(std::get<0>(TimeZVec[i]));
+
+        ass.setSim(simhit);
+        ass.setRec(outhit);
+        ass.setWeight(weight);
+        info() << "out hit, x: " << outhit.getPosition()[0] << " y: " << outhit.getPosition()[1] << " z: " << outhit.getPosition()[2] << " Q: " << outhit.getEDep() << " Time: "<< outhit.getTime() << endmsg;
+      }
+      else
+      {
+        debug() << " no matches" << endmsg;
+      }
+      }
+    }
+  return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiTool::finalize()
+{
+  StatusCode sc;
+  return sc;
+}

Digitization/DigiSimple/src/TPCPixelDigiTool.h

+#ifndef TPCPixelDigiTool_h
+#define TPCPixelDigiTool_h
+#include "CLHEP/Vector/TwoVector.h"
+#include "GaudiKernel/AlgTool.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "TF1.h"
+#include "DigiTool/IDigiTool.h"
+#include "DetInterface/IGeomSvc.h"
+#include "GearSvc/IGearSvc.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/SimPrimaryIonizationClusterCollection.h"
+#include "edm4hep/SimPrimaryIonizationCluster.h"
+#include <gear/GEAR.h>
+#include "DDRec/DetectorData.h"
+#include "DDRec/SurfaceManager.h"
+#include "TRandom3.h"
+
+#include <vector>
+
+class TPCPixelDigiTool : public extends<AlgTool, IDigiTool> {
+ public:
+  using extends::extends;
+
+  virtual StatusCode Call(const edm4hep::SimTrackerHitCollection* simCol, edm4hep::TrackerHitCollection* hitCol,
+			  edm4hep::MCRecoTrackerAssociationCollection* assCol) override;
+  virtual StatusCode Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection* hitCol,
+			  edm4hep::MCRecoTrackerAssociationCollection* assCol) override;
+
+  StatusCode initialize() override;
+  StatusCode finalize() override;
+  // double getPadPhi( CLHEP::Hep3Vector* thisPointRPhi, CLHEP::Hep3Vector* firstPointRPhi, CLHEP::Hep3Vector* middlePointRPhi, CLHEP::Hep3Vector* lastPointRPhi);
+  // double getPadTheta( CLHEP::Hep3Vector* firstPointRPhi, CLHEP::Hep3Vector* middlePointRPhi, CLHEP::Hep3Vector* lastPointRPhi );
+
+ private:
+
+  Gaudi::Property<std::string> m_detName{this, "DetName", "TPC"};
+  Gaudi::Property<std::string> m_readoutName{this, "Readout", "TPCCollection"};
+  
+
+  Gaudi::Property<int> _magnetic{this, "magnetic", 3};// Magnet field 3 Tesla
+  Gaudi::Property<double> _minimumTime{this, "minimumTime", 0.001};// if the time between two electrons less than this _minimumTime(us), then the two electrons will be treated as one 
+  Gaudi::Property<double> _deadTime{this, "deadTime", 300};// 
+  Gaudi::Property<int> _saveAllTime{this, "saveAllTime", 1};// 
+  Gaudi::Property<double> _ampDiffusion{this, "ampDiffusion", 0.0003}; // Amplification factor
+  Gaudi::Property<double> _MaxIonEles{this, "MaxIonEles",-1}; // Max allowed Inoized eles, negetive values mean no limits
+  Gaudi::Property<std::vector<double> > _polyParas{this, "polyParas", {3802., 0.487, 2092}};//amplification polynominal paramters
+  Gaudi::Property<double> _skipAss{this, "skipAss",1}; // Skip the mactching between digi hits and sim tracker hits
+
+
+
+
+  TF1 *_f1;
+  SmartIF<IRndmGenSvc> m_randSvc;
+  SmartIF<IGeomSvc> m_geosvc;
+  SmartIF<IGearSvc> _gear;
+  gear::GearMgr* _GEAR;
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+  const dd4hep::rec::SurfaceMap* m_surfaces;
+  double x_diffusion_p0;
+  double x_diffusion_p1;
+  double y_diffusion_p0; 
+  double y_diffusion_p1;
+  double t_diffusion_p0; 
+  double t_diffusion_p1;
+  double tpcXRes;
+  double tpcYRes;
+  double tpcZRes;
+  double tpcTRes;
+  double driftLength;
+  
+
+
+};
+#endif

README.md

@@ -38,3 +38,7 @@ $ ./run.sh Examples/options/helloalg.py
 * Reconstruction: Reconstruction
 
 
+## CyberPFA-5.0.1-dev (developing)
+* Based on CEPCSW tag tdr 24.12.0
+
+

Reconstruction/JetClustering/CMakeLists.txt

+# Set the CMake module path if necessary
+
+set(FASTJET_ROOT "/cvmfs/sft.cern.ch/lcg/views/LCG_103/x86_64-centos7-gcc11-opt/")
+set(INCLUDE_DIRS "${FASTJET_ROOT}/include")
+set(LIBRARY_DIRS "${FASTJET_ROOT}/lib")
+find_library(FASTJET_LIB fastjet PATHS ${LIBRARY_DIRS})
+
+# Modules
+gaudi_add_module(JetClustering
+                 SOURCES src/JetClustering.cpp
+                 LINK GearSvc
+                      Gaudi::GaudiKernel
+                      k4FWCore::k4FWCore
+                      DetInterface
+                      DataHelperLib
+                      ${GEAR_LIBRARIES} 
+                      ${GSL_LIBRARIES} 
+                      ${ROOT_LIBRARIES}
+                      ${FASTJET_LIB}
+                      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+)
+install(TARGETS JetClustering
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)

Reconstruction/JetClustering/src/JetClustering.cpp

+#include "JetClustering.h"
+
+#include "GaudiKernel/DataObject.h"
+#include "GaudiKernel/IHistogramSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/SmartDataPtr.h"
+
+#include "DetInterface/IGeomSvc.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include <math.h>
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/PseudoJet.hh"
+//time measurement
+#include <chrono>
+#include "TLorentzVector.h"
+
+using namespace fastjet;
+using namespace std;
+
+DECLARE_COMPONENT( JetClustering )
+
+//------------------------------------------------------------------------------
+JetClustering::JetClustering( const string& name, ISvcLocator* pSvcLocator )
+    : Algorithm( name, pSvcLocator ) {
+
+  declareProperty("InputPFOs", m_PFOColHdl, "Handle of the Input PFO collection");
+  declareProperty("Algorithm", m_algo, "Jet clustering algorithm");
+  declareProperty("nJets", m_nJets, "Number of jets to be clustered");
+  declareProperty("R", m_R, "Jet clustering radius");
+  declareProperty("OutputFile", m_outputFile, "Output file name");
+
+}
+
+//------------------------------------------------------------------------------
+StatusCode JetClustering::initialize(){
+
+  _nEvt = 0;
+
+  // Create a new TTree
+  _file = TFile::Open(m_outputFile.value().c_str(), "RECREATE");
+  _tree = new TTree("jets", "jets");
+  _tree->Branch("jet1_px", &jet1_px, "jet1_px/D");
+  _tree->Branch("jet1_py", &jet1_py, "jet1_py/D");
+  _tree->Branch("jet1_pz", &jet1_pz, "jet1_pz/D");
+  _tree->Branch("jet1_E", &jet1_E, "jet1_E/D");
+  _tree->Branch("jet1_costheta", &jet1_costheta, "jet1_costheta/D");
+  _tree->Branch("jet1_phi", &jet1_phi, "jet1_phi/D");
+  _tree->Branch("jet1_pt", &jet1_pt, "jet1_pt/D");
+  _tree->Branch("jet1_nconstituents", &jet1_nconstituents, "jet1_nconstituents/I");
+  _tree->Branch("jet2_px", &jet2_px, "jet2_px/D");
+  _tree->Branch("jet2_py", &jet2_py, "jet2_py/D");
+  _tree->Branch("jet2_pz", &jet2_pz, "jet2_pz/D");
+  _tree->Branch("jet2_E", &jet2_E, "jet2_E/D");
+  _tree->Branch("jet2_costheta", &jet2_costheta, "jet2_costheta/D");
+  _tree->Branch("jet2_phi", &jet2_phi, "jet2_phi/D");
+  _tree->Branch("jet2_pt", &jet2_pt, "jet2_pt/D");
+  _tree->Branch("jet2_nconstituents", &jet2_nconstituents, "jet2_nconstituents/I");
+  _tree->Branch("constituents_E1tot", &constituents_E1tot, "constituents_E1tot/D");
+  _tree->Branch("constituents_E2tot", &constituents_E2tot, "constituents_E2tot/D");
+  _tree->Branch("mass", &mass, "mass/D");
+  _tree->Branch("mass_allpfo", &mass_allpfo, "mass_allpfo/D");
+  _tree->Branch("ymerge", ymerge, "ymerge[6]/D");
+  
+  return StatusCode::SUCCESS;
+
+}
+
+//------------------------------------------------------------------------------
+StatusCode JetClustering::execute(){
+  
+  const edm4hep::ReconstructedParticleCollection* PFO = nullptr;
+  try {
+    PFO = m_PFOColHdl.get();
+  }
+  catch ( GaudiException &e ){
+    info() << "Collection " << m_PFOColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  // Check if the collection is empty
+  if ( PFO->size() == 0 ) {
+    info() << "Collection " << m_PFOColHdl.fullKey() << " is empty in event " << _nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  //initial indces
+  CleanVars();
+
+  //resize particles. there are many particles are empty, maybe upstream bugs
+  TLorentzVector p4_allpfo; 
+  vector<PseudoJet> input_particles;
+  for(const auto& pfo : *PFO){
+    if ( isnan(pfo.getMomentum()[0]) || isnan(pfo.getMomentum()[1]) || isnan(pfo.getMomentum()[2]) || pfo.getEnergy() == 0) {
+      info() << "Particle " << _particle_index << " px: " << pfo.getMomentum()[0] << " py: " << pfo.getMomentum()[1] << " pz: " << pfo.getMomentum()[2] << " E: " << pfo.getEnergy()  << endmsg;
+      _particle_index++;
+      continue;
+    }
+    else{
+      input_particles.push_back(PseudoJet(pfo.getMomentum()[0], pfo.getMomentum()[1], pfo.getMomentum()[2], pfo.getEnergy()));
+      TLorentzVector p4_pfo; 
+      p4_pfo.SetPxPyPzE( pfo.getMomentum()[0], pfo.getMomentum()[1], pfo.getMomentum()[2], pfo.getEnergy() );
+      p4_allpfo += p4_pfo;
+      _particle_index++;
+      debug() << "Particle " << _particle_index << " px: "<<pfo.getMomentum()[0]<<" py: "<<pfo.getMomentum()[1]<<" pz: "<<pfo.getMomentum()[2]<<" E: "<<pfo.getEnergy()<<endmsg;
+    }
+  }
+  mass_allpfo = p4_allpfo.M();
+
+  // create a jet definition
+  int nJets = m_nJets;
+  double R = m_R;
+
+  JetDefinition jet_def(ee_kt_algorithm);
+  //JetDefinition jet_def(antikt_algorithm, R);
+
+  // run the clustering, extract the jets
+  //std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
+
+  ClusterSequence clust_seq(input_particles, jet_def);
+  
+  vector<PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(nJets));
+  
+  std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
+  
+  //_time = std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count() / 1000000.0;
+  
+  //info() << "FastJet clustering done in seconds: " << _time << endmsg;
+  
+  //string jet_def_str = jet_def.description();
+  //debug() << "Clustering with " << jet_def_str << " for " << nJets << " jets with R = " << R << endmsg;
+  
+  // print the jets
+  info() << "                     px           py           pz           E           costheta  phi " << endmsg;
+  for (unsigned i = 0; i < jets.size(); i++) {
+    info() << "   jet " << i << "             " << jets[i].px() << " " << jets[i].py() << " " << jets[i].pz() << " " << jets[i].E() << " " << jets[i].cos_theta() << " " << jets[i].phi() << endmsg;
+    vector<PseudoJet> constituents = jets[i].constituents();
+    for (unsigned j = 0; j < constituents.size(); j++) {
+      if ( i == 1 )constituents_E1tot += constituents[j].E();
+      if ( i == 2 )constituents_E2tot += constituents[j].E();
+      info() << "    constituent  " << j << " " << constituents[j].px() << " " << constituents[j].py() << " " << constituents[j].pz() << " " << constituents[j].E() << " " << constituents[j].cos_theta() << " " << constituents[j].phi() << endmsg;
+    }
+  }
+  info() << "Total energy of constituents: " << constituents_E1tot << " " << constituents_E2tot << endmsg;
+
+
+  double _ymin[20];
+  for(int i=1; i<6;i++){
+    _ymin[i-1] = clust_seq.exclusive_ymerge (i);
+    info() << " -log10(y" << i << i+1 << ") = " << -log10(_ymin[i-1]) << endmsg;
+  }
+  
+
+  mass = (jets[0] + jets[1]).m();
+
+  // fill the tree
+  jet1_px = jets[0].px();
+  jet1_py = jets[0].py();
+  jet1_pz = jets[0].pz();
+  jet1_E = jets[0].E();
+  jet1_costheta = jets[0].cos_theta();
+  jet1_phi = jets[0].phi();
+  jet1_pt = jets[0].pt();
+  jet1_nconstituents = jets[0].constituents().size();
+  jet2_px = jets[1].px();
+  jet2_py = jets[1].py();
+  jet2_pz = jets[1].pz();
+  jet2_E = jets[1].E();
+  jet2_costheta = jets[1].cos_theta();
+  jet2_phi = jets[1].phi();
+  jet2_pt = jets[1].pt();
+  jet2_nconstituents = jets[1].constituents().size();
+  for(int i=0; i<6; i++){
+    ymerge[i] = _ymin[i];
+  }
+  _tree->Fill();
+
+  _nEvt++;
+  return StatusCode::SUCCESS;
+
+}
+
+//------------------------------------------------------------------------------
+StatusCode JetClustering::finalize(){
+  debug() << "Finalizing..." << endmsg;
+  _file->Write();
+  return StatusCode::SUCCESS;
+}

Reconstruction/JetClustering/src/JetClustering.h

+#ifndef JetClustering_h
+#define JetClustering_h 1
+
+#include "UTIL/ILDConf.h"
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/Algorithm.h"
+#include <random>
+#include "GaudiKernel/NTuple.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#include "edm4hep/ReconstructedParticleData.h"
+#include "edm4hep/ReconstructedParticleCollection.h"
+#include "edm4hep/ReconstructedParticle.h"
+
+class JetClustering : public Algorithm {
+ public:
+  // Constructor of this form must be provided
+  JetClustering( const std::string& name, ISvcLocator* pSvcLocator );
+
+  // Three mandatory member functions of any algorithm
+  StatusCode initialize() override;
+  StatusCode execute() override;
+  StatusCode finalize() override;
+
+ private:
+  DataHandle<edm4hep::ReconstructedParticleCollection> m_PFOColHdl{"PandoraPFOs", Gaudi::DataHandle::Reader, this};
+  Gaudi::Property<std::string> m_algo{this, "Algorithm", "ee_kt_algorithm"};
+  Gaudi::Property<int> m_nJets{this, "nJets", 2};
+  Gaudi::Property<double> m_R{this, "R", 0.6};
+  Gaudi::Property<std::string> m_outputFile{this, "OutputFile", "JetClustering.root"};
+
+  int _nEvt;
+  int _particle_index;
+  int _jet_index;
+  int _nparticles;
+  double _time;
+
+  TFile* _file;
+  TTree* _tree;
+  double jet1_px, jet1_py, jet1_pz, jet1_E;
+  double jet2_px, jet2_py, jet2_pz, jet2_E;
+  double jet1_costheta, jet1_phi;
+  double jet2_costheta, jet2_phi;
+  double jet1_pt, jet2_pt;
+  int jet1_nconstituents, jet2_nconstituents;
+  double constituents_E1tot;
+  double constituents_E2tot;
+  double ymerge[6];
+  double mass;
+  double mass_allpfo;
+
+  void CleanVars(){
+
+    _particle_index = 0;
+    _jet_index = 0;
+    _nparticles = 0;
+    _time = 0;
+    
+    jet1_px = jet1_py = jet1_pz = jet1_E = 0;
+    jet2_px = jet2_py = jet2_pz = jet2_E = 0;
+    jet1_costheta = jet1_phi = 0;
+    jet2_costheta = jet2_phi = 0;
+    jet1_pt = jet2_pt = 0;
+    jet1_nconstituents = jet2_nconstituents = 0;
+    constituents_E1tot = constituents_E2tot = 0;
+    for(int i=0; i<6; i++){
+      ymerge[i] = 0;
+    }
+    mass = 0;
+    mass_allpfo = 0;
+  }
+
+};
+#endif
+

Reconstruction/RecActsTracking/options/RecActsTracking.py

@@ -57,9 +57,22 @@ cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest"
 
 from Configurables import RecActsTracking
 actstracking = RecActsTracking("RecActsTracking")
-actstracking.TGeoFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr24.12.0/SiTrack-tgeo.root")
-actstracking.TGeoConfigFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr24.12.0/SiTrack-tgeo-config.json")
-actstracking.MaterialMapFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr24.12.0/SiTrack-material-maps.json")
+actstracking.TGeoFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr25.1.0/SiTrack-tgeo.root")
+actstracking.TGeoConfigFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr25.1.0/SiTrack-tgeo-config.json")
+actstracking.MaterialMapFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Reconstruction/RecActsTracking/data/tdr25.1.0/SiTrack-material-maps.json")
+
+# config for acts tracking
+actstracking.onSurfaceTolerance = 2e-1
+# actstracking.ExtendSeedRange = True
+# actstracking.SeedDeltaRMin = 10
+# actstracking.SeedDeltaRMax = 200
+# actstracking.SeedRMax = 600
+# actstracking.SeedRMin = 80
+# actstracking.SeedImpactMax = 4
+# actstracking.SeedRMinMiddle = 340
+# actstracking.SeedRMaxMiddle = 380
+actstracking.numMeasurementsCutOff = 2
+actstracking.CKFchi2Cut = 1
 
 ##############################################################################
 # output

Reconstruction/RecActsTracking/src/RecActsTracking.cpp

@@ -22,7 +22,6 @@ DECLARE_COMPONENT(RecActsTracking)
 RecActsTracking::RecActsTracking(const std::string& name, ISvcLocator* svcLoc)
     : GaudiAlgorithm(name, svcLoc)
 {
-    // _encoder = new UTIL::BitField64(lcio::ILDCellID0::encoder_string);
 }
 
 StatusCode RecActsTracking::initialize()
@@ -102,6 +101,9 @@ StatusCode RecActsTracking::initialize()
     m_sit_surfaces = m_geosvc->getSurfaceMap("SIT");
     debug() << "Surface map size: " << m_sit_surfaces->size() << endmsg;
 
+    m_ftd_surfaces = m_geosvc->getSurfaceMap("FTD");
+    debug() << "Surface map size: " << m_ftd_surfaces->size() << endmsg;
+
     vxd_decoder = m_geosvc->getDecoder("VXDCollection");
     if(!vxd_decoder){
         return StatusCode::FAILURE;
@@ -112,6 +114,11 @@ StatusCode RecActsTracking::initialize()
         return StatusCode::FAILURE;
     }
 
+    ftd_decoder = m_geosvc->getDecoder("FTDCollection");
+    if(!ftd_decoder){
+        return StatusCode::FAILURE;
+    }
+
     info() << "ACTS Tracking Geometry initialized successfully!" << endmsg;
     // initialize tgeo detector
     auto logger = Acts::getDefaultLogger("TGeoDetector", Acts::Logging::INFO);
@@ -180,8 +187,8 @@ StatusCode RecActsTracking::execute()
     Selected_Seeds.clear();
 
     chronoStatSvc->chronoStart("read input hits");
-    int successVTX = InitialiseVTX();
 
+    int successVTX = InitialiseVTX();
     if (successVTX == 0)
     {
         _nEvt++;
@@ -195,6 +202,12 @@ StatusCode RecActsTracking::execute()
         return StatusCode::SUCCESS;
     }
 
+    int successFTD = InitialiseFTD();
+    if(successFTD == 0){
+        _nEvt++;
+        return StatusCode::SUCCESS;
+    }
+
     chronoStatSvc->chronoStop("read input hits");
     // info() << "Generated " << SpacePointPtrs.size() << " spacepoints for event " << _nEvt << "!" << endmsg;
     // info() << "Generated " << measurements.size() << " measurements for event " << _nEvt << "!" << endmsg;
@@ -571,27 +584,68 @@ StatusCode RecActsTracking::execute()
         _nEvt++;
         return StatusCode::SUCCESS;
     }
-
     for (const auto& cur_track : constTracks)
     {
         auto writeout_track = trkCol->create();
-        int nVTXHit = 0, nFTDHit = 0, nSITHit = 0;
+        int nVTXHit = 0, nFTDHit = 0, nSITHit = 0, nlayerHits = 9;
         int getFirstHit = 0;
         
         writeout_track.setChi2(cur_track.chi2());
         writeout_track.setNdf(cur_track.nDoF());
         writeout_track.setDEdx(cur_track.absoluteMomentum() / Acts::UnitConstants::GeV);
-        writeout_track.setDEdxError(cur_track.qOverP());
-        for (auto trackState : cur_track.trackStates()){
-            if (trackState.hasUncalibratedSourceLink()){
+        // writeout_track.setDEdxError(cur_track.qOverP());
+        std::array<int, 6> nlayer_VTX{0, 0, 0, 0, 0, 0};
+        std::array<int, 3> nlayer_SIT{0, 0, 0};
+        std::array<int, 3> nlayer_FTD{0, 0, 0};
+
+        for (auto trackState : cur_track.trackStates())
+        {
+            if (trackState.hasUncalibratedSourceLink())
+            {
                 auto cur_measurement_sl = trackState.getUncalibratedSourceLink();
                 const auto& MeasSourceLink = cur_measurement_sl.get<IndexSourceLink>();
                 auto cur_measurement = measurements[MeasSourceLink.index()];
                 auto cur_measurement_gid = MeasSourceLink.geometryId();
-                if (std::find(VXD_volume_ids.begin(), VXD_volume_ids.end(), cur_measurement_gid.volume()) != VXD_volume_ids.end()){
-                    nVTXHit++;
-                } else if (std::find(SIT_acts_volume_ids.begin(), SIT_acts_volume_ids.end(), cur_measurement_gid.volume()) != SIT_acts_volume_ids.end()){
-                    nSITHit++;
+                
+                for (int i = 0; i < 5; i++)
+                {
+                    if (cur_measurement_gid.volume() == VXD_volume_ids[i])
+                    {
+                        if (i == 4){
+                            if (cur_measurement_gid.sensitive() & 1 == 1)
+                            {
+                                nlayer_VTX[5]++;
+                            } else{
+                                nlayer_VTX[4]++;
+                            }
+                        } else {
+                            nlayer_VTX[i]++;
+                        }
+                        nVTXHit++;
+                        break;
+                    }
+                }
+
+                for (int i = 0; i < 3; i++){
+                    if (cur_measurement_gid.volume() == SIT_volume_ids[i]){
+                        nlayer_SIT[i]++;
+                        nSITHit++;
+                        break;
+                    }
+                }
+
+                for (int i = 0; i < 3; i++){
+                    if (cur_measurement_gid.volume() == FTD_positive_volume_ids[i]){
+                        nlayer_FTD[i]++;
+                        nFTDHit++;
+                        break;
+                    }
+
+                    if (cur_measurement_gid.volume() == FTD_negative_volume_ids[i]){
+                        nlayer_FTD[i]++;
+                        nFTDHit++;
+                        break;
+                    }
                 }
 
                 writeout_track.addToTrackerHits(MeasSourceLink.getTrackerHit());
@@ -609,23 +663,55 @@ StatusCode RecActsTracking::execute()
                 }
             }
         }
+        for (int i = 0; i < 6; i++){
+            m_nRec_VTX[i] += nlayer_VTX[i];
+            if (nlayer_VTX[i] == 0) {
+                m_n0EventHits[i]++;
+                nlayerHits--;
+            } else if (nlayer_VTX[i] == 1) {
+                m_n1EventHits[i]++;
+            } else if (nlayer_VTX[i] == 2) {
+                m_n2EventHits[i]++;
+            } else if (nlayer_VTX[i] >= 3) {
+                m_n3EventHits[i]++;
+            }
+        }
+        for (int i = 0; i < 3; i++){
+            m_nRec_SIT[i] += nlayer_SIT[i];
+            if (nlayer_SIT[i] == 0) {
+                m_n0EventHits[i+6]++;
+                nlayerHits--;
+            } else if (nlayer_SIT[i] == 1) {
+                m_n1EventHits[i+6]++;
+            } else if (nlayer_SIT[i] == 2) {
+                m_n2EventHits[i+6]++;
+            } else if (nlayer_SIT[i] >= 3) {
+                m_n3EventHits[i+6]++;
+            }
+        }
+
+        for (int i = 0; i < 3; i++){
+            m_nRec_FTD[i] += nlayer_FTD[i];
+        }
 
         // SubdetectorHitNumbers: VXD = 0, FTD = 1, SIT = 2
+        writeout_track.setDEdxError(nlayerHits);
         writeout_track.addToSubdetectorHitNumbers(nVTXHit);
         writeout_track.addToSubdetectorHitNumbers(nFTDHit);
         writeout_track.addToSubdetectorHitNumbers(nSITHit);
         
+        // TODO: covmatrix need to be converted
         std::array<float, 21> writeout_covMatrix;
         auto cur_track_covariance = cur_track.covariance();
         for (int i = 0; i < 6; i++) {
             for (int j = 0; j < 6-i; j++) {
-                writeout_covMatrix[i * 6 + j] = cur_track_covariance(writeout_indices[i], writeout_indices[j]);
+                writeout_covMatrix[int((13-i)*i/2 + j)] = cur_track_covariance(writeout_indices[i], writeout_indices[j]);
             }
         }
         // location: At{Other|IP|FirstHit|LastHit|Calorimeter|Vertex}|LastLocation
         // TrackState: location, d0, phi, omega, z0, tanLambda, time, referencePoint, covMatrix
         edm4hep::TrackState writeout_trackState{
-            0, // location: AtOther
+            1, // location: AtOther
             cur_track.loc0() / Acts::UnitConstants::mm, // d0
             cur_track.phi(), // phi
             cur_track.qOverP() * sin(cur_track.theta()) * _FCT * m_field, // omega = qop * sin(theta) * _FCT * bf
@@ -653,7 +739,15 @@ StatusCode RecActsTracking::finalize()
     info() << "Total number of **TotalSeeds** processed: " << m_nTotalSeeds << endmsg;
     info() << "Total number of **FoundTracks** processed: " << m_nFoundTracks << endmsg;
     info() << "Total number of **SelectedTracks** processed: " << m_nSelectedTracks << endmsg;
-    
+    info() << "Total number of **LayerHits** processed: " << m_nLayerHits << endmsg;
+    info() << "Total number of **Rec_VTX** processed: " << m_nRec_VTX << endmsg;
+    info() << "Total number of **Rec_SIT** processed: " << m_nRec_SIT << endmsg;
+    info() << "Total number of **Rec_FTD** processed: " << m_nRec_FTD << endmsg;
+    info() << "Total number of **EventHits0** processed: " << m_n0EventHits << endmsg;
+    info() << "Total number of **EventHits1** processed: " << m_n1EventHits << endmsg;
+    info() << "Total number of **EventHits2** processed: " << m_n2EventHits << endmsg;
+    info() << "Total number of **EventHitsmore** processed: " << m_n3EventHits << endmsg;
+
     return GaudiAlgorithm::finalize();
 }
 
@@ -822,6 +916,8 @@ int RecActsTracking::InitialiseVTX()
                 IndexSourceLink sourceLink{moduleGeoId, measurementIdx, hit};
                 sourceLinks.insert(sourceLinks.end(), sourceLink);
                 Acts::SourceLink sl{sourceLink};
+                boost::container::static_vector<Acts::SourceLink, 2> slinks;
+                slinks.emplace_back(sl);
 
                 // get the local position of the hit
                 IndexSourceLink::SurfaceAccessor surfaceAccessor{*trackingGeometry};
@@ -841,6 +937,11 @@ int RecActsTracking::InitialiseVTX()
                 auto acts_global_postion = acts_surface->localToGlobal(geoContext, par, globalFakeMom);
                 debug() << "debug surface at: x:" << acts_global_postion[0] << ", y:" << acts_global_postion[1] << ", z:" << acts_global_postion[2] << endmsg;
 
+                if (ExtendSeedRange.value()) {
+                    SimSpacePoint *hitExt = new SimSpacePoint(hit, acts_global_postion[0], acts_global_postion[1], acts_global_postion[2], 0.002, slinks);
+                    SpacePointPtrs.push_back(hitExt);
+                }
+                
                 // create and store the measurement
                 // Plane covMatrix[6] = {u_direction[0], u_direction[1], resU, v_direction[0], v_direction[1], resV}
                 Acts::ActsSquareMatrix<2> cov = Acts::ActsSquareMatrix<2>::Identity();
@@ -934,7 +1035,7 @@ int RecActsTracking::InitialiseSIT()
             }
 
             // set acts geometry identifier
-            uint64_t acts_volume = SIT_acts_volume_ids[m_layer];
+            uint64_t acts_volume = SIT_volume_ids[m_layer];
             uint64_t acts_boundary = 0;
             uint64_t acts_layer = 2;
             uint64_t acts_approach = 0;
@@ -952,6 +1053,8 @@ int RecActsTracking::InitialiseSIT()
             IndexSourceLink sourceLink{moduleGeoId, measurementIdx, hit};
             sourceLinks.insert(sourceLinks.end(), sourceLink);
             Acts::SourceLink sl{sourceLink};
+            boost::container::static_vector<Acts::SourceLink, 2> slinks;
+            slinks.emplace_back(sl);
 
             // get the local position of the hit
             IndexSourceLink::SurfaceAccessor surfaceAccessor{*trackingGeometry};
@@ -971,6 +1074,11 @@ int RecActsTracking::InitialiseSIT()
             auto acts_global_postion = acts_surface->localToGlobal(geoContext, par, globalFakeMom);
             debug() << "debug surface at: x:" << acts_global_postion[0] << ", y:" << acts_global_postion[1] << ", z:" << acts_global_postion[2] << endmsg;
 
+            if (ExtendSeedRange.value()) {
+                SimSpacePoint *hitExt = new SimSpacePoint(hit, acts_global_postion[0], acts_global_postion[1], acts_global_postion[2], 0.002, slinks);
+                SpacePointPtrs.push_back(hitExt);
+            }
+
             // create and store the measurement
             Acts::ActsSquareMatrix<2> cov = Acts::ActsSquareMatrix<2>::Identity();
             cov(0, 0) = std::max<double>(double(m_covMatrix[2]), eps.value());
@@ -985,4 +1093,116 @@ int RecActsTracking::InitialiseSIT()
     } else { success = 0; }
 
     return success;
+}
+
+int RecActsTracking::InitialiseFTD()
+{
+    const edm4hep::TrackerHitCollection* hitFTDCol = nullptr;
+    const edm4hep::SimTrackerHitCollection* SimhitFTDCol = nullptr;
+
+    try {
+        hitFTDCol = _inFTDTrackHdl.get();
+    } catch (GaudiException& e) {
+        debug() << "Collection " << _inFTDTrackHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+        return 0;
+    }
+
+    try {
+        SimhitFTDCol = _inFTDColHdl.get();
+    } catch (GaudiException& e) {
+        debug() << "Collection " << _inFTDColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+        return 0;
+    }
+
+    if (hitFTDCol && SimhitFTDCol)
+    {
+        int nelem = hitFTDCol->size();
+        debug() << "Number of FTD hits = " << nelem << endmsg;
+        for (int ielem = 0; ielem < nelem; ++ielem)
+        {
+            auto hit = hitFTDCol->at(ielem);
+            auto simhit = SimhitFTDCol->at(ielem);
+            auto simcellid = simhit.getCellID();
+            uint64_t m_system = ftd_decoder->get(simcellid, "system");
+            uint64_t m_side = ftd_decoder->get(simcellid, "side");
+            uint64_t m_layer = ftd_decoder->get(simcellid, "layer");
+            uint64_t m_module = ftd_decoder->get(simcellid, "module");
+            uint64_t m_sensor = ftd_decoder->get(simcellid, "sensor");
+            double acts_x = simhit.getPosition()[0];
+            double acts_y = simhit.getPosition()[1];
+            double acts_z = simhit.getPosition()[2];
+            double momentum_x = simhit.getMomentum()[0];
+            double momentum_y = simhit.getMomentum()[1];
+            double momentum_z = simhit.getMomentum()[2];
+            const Acts::Vector3 globalmom{momentum_x, momentum_y, momentum_z};
+            std::array<float, 6> m_covMatrix = hit.getCovMatrix();
+            
+            if (m_layer > 2) {
+                continue;
+            }
+
+            dd4hep::rec::ISurface* surface = nullptr;
+            auto it = m_ftd_surfaces->find(simcellid);
+            if (it != m_ftd_surfaces->end()) {
+                surface = it->second;
+                if (!surface) {
+                    fatal() << "found surface for FTD cell id " << simcellid << ", but NULL" << endmsg;
+                    return 0;
+                }
+            }
+
+            // set acts geometry identifier
+            uint64_t acts_volume = (acts_z > 0) ? FTD_positive_volume_ids[m_layer] : FTD_negative_volume_ids[m_layer];
+            uint64_t acts_boundary = 0;
+            uint64_t acts_layer = 2;
+            uint64_t acts_approach = 0;
+            uint64_t acts_sensitive = m_module + 1;
+
+            Acts::GeometryIdentifier moduleGeoId;
+            moduleGeoId.setVolume(acts_volume);
+            moduleGeoId.setBoundary(acts_boundary);
+            moduleGeoId.setLayer(acts_layer);
+            moduleGeoId.setApproach(acts_approach);
+            moduleGeoId.setSensitive(acts_sensitive);   
+
+            // create and store the source link
+            uint32_t measurementIdx = measurements.size();
+            IndexSourceLink sourceLink{moduleGeoId, measurementIdx, hit};
+            sourceLinks.insert(sourceLinks.end(), sourceLink);
+            Acts::SourceLink sl{sourceLink};    
+            boost::container::static_vector<Acts::SourceLink, 2> slinks;
+            slinks.emplace_back(sl);
+
+            // get the local position of the hit
+            IndexSourceLink::SurfaceAccessor surfaceAccessor{*trackingGeometry};
+            const Acts::Surface* acts_surface = surfaceAccessor(sl);
+            const Acts::Vector3 globalPos{acts_x, acts_y, acts_z};
+            auto acts_local_postion = acts_surface->globalToLocal(geoContext, globalPos, globalmom, onSurfaceTolerance);
+            if (!acts_local_postion.ok()){
+                info() << "Error: failed to get local position for FTD layer " << m_layer << " module " << m_module << " sensor " << m_sensor << endmsg;
+                acts_local_postion = acts_surface->globalToLocal(geoContext, globalPos, globalmom, 100*onSurfaceTolerance);
+            }
+            const std::array<Acts::BoundIndices, 2> indices{Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc1};
+            const Acts::Vector2 par{acts_local_postion.value()[0], acts_local_postion.value()[1]};
+
+            // debug
+            debug() << "FTD measurements global position(x,y,z): " << simhit.getPosition()[0] << ", " << simhit.getPosition()[1] << ", " << simhit.getPosition()[2]
+                << "; local position(loc0, loc1): "<< acts_local_postion.value()[0] << ", " << acts_local_postion.value()[1] << endmsg;
+            auto acts_global_postion = acts_surface->localToGlobal(geoContext, par, globalFakeMom);
+            debug() << "debug surface at: x:" << acts_global_postion[0] << ", y:" << acts_global_postion[1] << ", z:" << acts_global_postion[2] << endmsg;
+
+            if (ExtendSeedRange.value()) {
+                SimSpacePoint *hitExt = new SimSpacePoint(hit, acts_global_postion[0], acts_global_postion[1], acts_global_postion[2], 0.002, slinks);
+                SpacePointPtrs.push_back(hitExt);
+            }
+
+            // create and store the measurement
+            Acts::ActsSquareMatrix<2> cov = Acts::ActsSquareMatrix<2>::Identity();
+            cov(0, 0) = std::max<double>(double(m_covMatrix[2]), eps.value());
+            cov(1, 1) = std::max<double>(double(m_covMatrix[5]), eps.value());
+            measurements.emplace_back(Acts::Measurement<Acts::BoundIndices, 2>(sl, indices, par, cov));
+        }
+    }
+
+    return 1;
 }
\ No newline at end of file

Reconstruction/RecActsTracking/src/RecActsTracking.h

@@ -166,9 +166,11 @@ class RecActsTracking : public GaudiAlgorithm
         // Input collections
         DataHandle<edm4hep::TrackerHitCollection> _inVTXTrackHdl{"VXDTrackerHits", Gaudi::DataHandle::Reader, this};
         DataHandle<edm4hep::TrackerHitCollection> _inSITTrackHdl{"SITTrackerHits", Gaudi::DataHandle::Reader, this};
+        DataHandle<edm4hep::TrackerHitCollection> _inFTDTrackHdl{"FTDTrackerHits", Gaudi::DataHandle::Reader, this};
 
         DataHandle<edm4hep::SimTrackerHitCollection> _inVTXColHdl{"VXDCollection", Gaudi::DataHandle::Reader, this};
         DataHandle<edm4hep::SimTrackerHitCollection> _inSITColHdl{"SITCollection", Gaudi::DataHandle::Reader, this};
+        DataHandle<edm4hep::SimTrackerHitCollection> _inFTDColHdl{"FTDCollection", Gaudi::DataHandle::Reader, this};
 
         DataHandle<edm4hep::MCParticleCollection> _inMCColHdl{"MCParticle", Gaudi::DataHandle::Reader, this};
 
@@ -183,6 +185,7 @@ class RecActsTracking : public GaudiAlgorithm
         Gaudi::Property<double> m_field{this, "Field", 3.0}; // tesla
         Gaudi::Property<double> onSurfaceTolerance{this, "onSurfaceTolerance", 1e-2}; // mm
         Gaudi::Property<double> eps{this, "eps", 1e-5}; // mm
+        Gaudi::Property<bool> ExtendSeedRange{this, "ExtendSeedRange", false};
 
         // seed finder config
         Gaudi::Property<double> SeedDeltaRMin{this, "SeedDeltaRMin", 4}; // mm
@@ -194,16 +197,18 @@ class RecActsTracking : public GaudiAlgorithm
         Gaudi::Property<double> SeedRMaxMiddle{this, "SeedRMaxMiddle", 24}; // mm
 
         // CKF config
-        Gaudi::Property<double> CKFchi2Cut{this, "CKFchi2Cut", 20};
-        Gaudi::Property<std::size_t> numMeasurementsCutOff{this, "numMeasurementsCutOff", 1};
+        Gaudi::Property<double> CKFchi2Cut{this, "CKFchi2Cut", std::numeric_limits<double>::max()};
+        Gaudi::Property<std::size_t> numMeasurementsCutOff{this, "numMeasurementsCutOff", 1u};
         Gaudi::Property<bool> CKFtwoWay{this, "CKFtwoWay", true};
         
         SmartIF<IGeomSvc> m_geosvc;
         SmartIF<IChronoStatSvc> chronoStatSvc;
         dd4hep::DDSegmentation::BitFieldCoder *vxd_decoder;
         dd4hep::DDSegmentation::BitFieldCoder *sit_decoder;
+        dd4hep::DDSegmentation::BitFieldCoder *ftd_decoder;
         const dd4hep::rec::SurfaceMap* m_vtx_surfaces;
         const dd4hep::rec::SurfaceMap* m_sit_surfaces;
+        const dd4hep::rec::SurfaceMap* m_ftd_surfaces;
 
         // configs to build acts geometry
         Acts::GeometryContext geoContext;
@@ -231,6 +236,7 @@ class RecActsTracking : public GaudiAlgorithm
 
         int InitialiseVTX();
         int InitialiseSIT();
+        int InitialiseFTD();
         const dd4hep::rec::ISurface* getISurface(edm4hep::TrackerHit* hit);
         const Acts::GeometryIdentifier getVTXGid(uint64_t cellid);
         const Acts::GeometryIdentifier getSITGid(uint64_t cellid);
@@ -272,10 +278,11 @@ class RecActsTracking : public GaudiAlgorithm
         // Acts::ParticleHypothesis particleHypothesis = Acts::ParticleHypothesis::chargedGeantino();
 
         // gid convert configuration
-        std::vector<uint64_t> VXD_volume_ids{16, 17, 18, 19, 20};
-        std::vector<uint64_t> SIT_acts_volume_ids{22, 24, 26};
+        std::vector<uint64_t> VXD_volume_ids{20, 21, 22, 23, 24};
+        std::vector<uint64_t> SIT_volume_ids{28, 31, 34};
+        std::vector<uint64_t> FTD_positive_volume_ids{29, 32, 35};
+        std::vector<uint64_t> FTD_negative_volume_ids{27, 7, 2};
         std::vector<uint64_t> SIT_module_nums{7, 10, 14};
-        uint64_t SIT_sensor_nums = 14;
 
         // CKF configuration
         // Acts::MeasurementSelector::Config measurementSelectorCfg;
@@ -290,7 +297,16 @@ class RecActsTracking : public GaudiAlgorithm
         mutable std::atomic<std::size_t> m_nFoundTracks{0};
         mutable std::atomic<std::size_t> m_nSelectedTracks{0};
         mutable std::atomic<std::size_t> m_nStoppedBranches{0};
-        
+        // layer hits, VXD (0-5) & SIT (6-8)
+        std::array<std::atomic<size_t>, 9> m_nLayerHits{0, 0, 0, 0, 0, 0, 0, 0, 0};
+        std::array<std::atomic<size_t>, 6> m_nRec_VTX{0, 0, 0, 0, 0, 0};
+        std::array<std::atomic<size_t>, 3> m_nRec_SIT{0, 0, 0};
+        std::array<std::atomic<size_t>, 3> m_nRec_FTD{0, 0, 0};
+        std::array<std::atomic<size_t>, 9> m_n0EventHits{0, 0, 0, 0, 0, 0, 0, 0, 0};
+        std::array<std::atomic<size_t>, 9> m_n1EventHits{0, 0, 0, 0, 0, 0, 0, 0, 0};
+        std::array<std::atomic<size_t>, 9> m_n2EventHits{0, 0, 0, 0, 0, 0, 0, 0, 0};
+        std::array<std::atomic<size_t>, 9> m_n3EventHits{0, 0, 0, 0, 0, 0, 0, 0, 0};
+
         mutable tbb::combinable<Acts::VectorMultiTrajectory::Statistics> m_memoryStatistics{[]() {
             auto mtj = std::make_shared<Acts::VectorMultiTrajectory>();
             return mtj->statistics();

Reconstruction/RecPFACyber/CMakeLists.txt

@@ -24,6 +24,8 @@ gaudi_add_library(CrystalCaloRecLib
                           ${ROOT_LIBRARIES}
                           ${TMVA_LIBRARIES}
                           EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+                          ${DD4hep_COMPONENT_LIBRARIES}
+                          DetIdentifier
 
 )
 

Reconstruction/RecPFACyber/include/Algorithm/GlobalClusteringAlg.h

@@ -6,7 +6,6 @@
 #include "time.h"
 #include <TTimeStamp.h> 
 #include <ctime>
-
 #include <cstdlib>
 
 using namespace Cyber;

Reconstruction/RecPFACyber/include/Algorithm/HoughClusteringAlg.h

@@ -36,8 +36,12 @@ public:
   
 
 private:
-	std::vector<Cyber::CaloHalfCluster*> p_HalfClusterV;
-  std::vector<Cyber::CaloHalfCluster*> p_HalfClusterU;
+	std::vector<Cyber::CaloHalfCluster*> barrel_HalfClusterV;  // Barrel ECAL, V: bars parallel to z axis
+  std::vector<Cyber::CaloHalfCluster*> barrel_HalfClusterU;  // Barrel ECAL, U: bars perpendicular to z axis
+  std::vector<Cyber::CaloHalfCluster*> endcap0_HalfClusterV; // Endcap ECAL at z~-2900mm, V: bars parallel to x axis
+  std::vector<Cyber::CaloHalfCluster*> endcap0_HalfClusterU; // Endcap ECAL at z~-2900mm, U: bars parallel to y axis
+  std::vector<Cyber::CaloHalfCluster*> endcap1_HalfClusterV; // Endcap ECAL at z~2900mm, V: bars parallel to x axis
+  std::vector<Cyber::CaloHalfCluster*> endcap1_HalfClusterU; // Endcap ECAL at z~2900mm, U: bars parallel to y axis
 
   std::vector<const Cyber::Calo1DCluster*> m_localMaxVCol;
   std::vector<const Cyber::Calo1DCluster*> m_localMaxUCol;

Reconstruction/RecPFACyber/include/Algorithm/TrackExtrapolatingAlg.h

@@ -27,14 +27,16 @@ public:
 
   // Get points in each plane of layer
   StatusCode GetLayerRadius(std::vector<float> & ECAL_layer_radius, std::vector<float> & HCAL_layer_radius);
+  StatusCode GetLayerZ(std::vector<float> & ECAL_layer_z, std::vector<float> & HCAL_layer_z);
   // If the track reach barrel ECAL
   bool IsReachECAL(Cyber::Track * track);
   // Get track state at calorimeter
   StatusCode GetTrackStateAtCalo(Cyber::Track * track, 
                            Cyber::TrackState & trk_state_at_calo);
   // get extrapolated points
-  StatusCode ExtrapolateByRadius( const std::vector<float> & ECAL_layer_radius, std::vector<float> & HCAL_layer_radius,
-                                  const Cyber::TrackState & CALO_trk_state, Cyber::Track* p_track);
+  StatusCode Extrapolate( const std::vector<float> & ECAL_layer_radius, const std::vector<float> & ECAL_layer_z, 
+                          const std::vector<float> & HCAL_layer_radius, const std::vector<float> & HCAL_layer_z,
+                          const Cyber::TrackState & CALO_trk_state, Cyber::Track* p_track);
   // Get the radius rho 
   float GetRho(const Cyber::TrackState & trk_state);
   // Get coordinates of the center of the circle
@@ -44,11 +46,14 @@ public:
   // If the charged particle return back 
   bool IsReturn(float rho, TVector2 & center);
 
-  std::vector<float> GetDeltaPhi(float rho, TVector2 center, float alpha0, vector<float> layer_radius, const Cyber::TrackState & CALO_trk_state);
+  StatusCode GetDeltaPhi( float rho, TVector2 center, float alpha0, 
+                          const std::vector<float> & layer_radius, const std::vector<float> & layer_z, 
+                          const Cyber::TrackState & CALO_trk_state,
+                          vector<float> & barrel_delta_phi, vector<float> & endcap_delta_phi);
   std::vector<TVector3> GetExtrapoPoints(std::string calo_name, 
                                          float rho, TVector2 center, float alpha0, 
                                          const Cyber::TrackState & CALO_trk_state,
-                                         const std::vector<float>& delta_phi);
+                                         const std::vector<float>& barrel_delta_phi, const std::vector<float>& endcap_delta_phi);
 
     // Get phi0 of extrapolated points. Note that this phi0 is not same as the definition of the phi0 in TrackState, but will be stored in TrackState
   float GetExtrapolatedPhi0(float Kappa, float ECAL_phi0, TVector2 center, TVector3 ext_point);

Reconstruction/RecPFACyber/include/CyberDataCol.h

+//=============================================================
+// CyberPFA: a PFA developed for CEPC referenece detector
+// Ver. CyberPFA-5.0.1(2025.01.09)
+//-------------------------------------------------------------
+// Data Collection with CyberPFA EDM
+//-------------------------------------------------------------
+//  Author: Fangyi Guo, Yang Zhang, Weizheng Song, Shengsen Sun
+//          (IHEP, CAS)
+//  Contact: guofangyi@ihep.ac.cn,
+//           sunss@ihep.ac.cn
+//=============================================================
 #ifndef _PANDORAPLUS_DATA_H
 #define _PANDORAPLUS_DATA_H
 #include <iostream>

Reconstruction/RecPFACyber/include/CyberPFAlg.h

+//=============================================================
+// CyberPFA: a PFA developed for CEPC referenece detector
+// Ver. CyberPFA-5.0.1(2025.01.09)
+//-------------------------------------------------------------
+//  Author: Fangyi Guo, Yang Zhang, Weizheng Song, Shengsen Sun
+//          (IHEP, CAS)
+//  Contact: guofangyi@ihep.ac.cn,
+//           sunss@ihep.ac.cn
+//=============================================================
 #ifndef PANDORAPLUS_ALG_H
 #define PANDORAPLUS_ALG_H
 
@@ -8,6 +17,7 @@
 #include <DDRec/CellIDPositionConverter.h>
 #include <DD4hep/Segmentations.h>
 #include "DetInterface/IGeomSvc.h"
+#include "DetIdentifier/CEPCDetectorData.h"
 #include <CrystalEcalSvc/ICrystalEcalSvc.h>
 
 #include "k4FWCore/PodioDataSvc.h"
@@ -73,7 +83,6 @@ public:
   /** Called after data processing for clean up.
    */
   virtual StatusCode finalize() ;
-  
 
 protected:
 
@@ -86,7 +95,7 @@ protected:
   dd4hep::Detector* m_dd4hep;
   dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
   dd4hep::VolumeManager m_volumeManager;
-
+  std::map<std::tuple<int, int, int, int, int>, int> barNumberMapEndcapMap;
   //DataCollection: moved into execute() to ensure everything can be cleand after one event. 
   //CyberDataCol     m_DataCol; 
 
@@ -191,10 +200,10 @@ protected:
   //Raw bars and hits
   TTree* t_SimBar;
   float m_totE_EcalSim, m_totE_HcalSim;
-  FloatVec m_simBar_x, m_simBar_y, m_simBar_z, m_simBar_T1, m_simBar_T2, m_simBar_Q1, m_simBar_Q2; 
+  FloatVec m_simBar_x, m_simBar_y, m_simBar_z, m_simBar_length, m_simBar_nBarInLayer, m_simBar_T1, m_simBar_T2, m_simBar_Q1, m_simBar_Q2; 
   FloatVec m_simBar_truthMC_tag, m_simBar_truthMC_pid, m_simBar_truthMC_px, m_simBar_truthMC_py, m_simBar_truthMC_pz, m_simBar_truthMC_E, 
-           m_simBar_truthMC_EPx, m_simBar_truthMC_EPy, m_simBar_truthMC_EPz, m_simBar_truthMC_weight;
-  IntVec m_simBar_dlayer, m_simBar_stave, m_simBar_slayer, m_simBar_module, m_simBar_bar;
+           m_simBar_truthMC_EPx, m_simBar_truthMC_EPy, m_simBar_truthMC_EPz, m_simBar_truthMC_weight;   
+  IntVec m_simBar_dlayer, m_simBar_stave, m_simBar_slayer, m_simBar_module, m_simBar_bar, m_simBar_system;
   FloatVec m_HcalHit_x, m_HcalHit_y, m_HcalHit_z, m_HcalHit_E, 
            m_HcalHit_truthMC_tag, m_HcalHit_truthMC_pid, m_HcalHit_truthMC_px, m_HcalHit_truthMC_py, m_HcalHit_truthMC_pz, m_HcalHit_truthMC_E,
            m_HcalHit_truthMC_EPx, m_HcalHit_truthMC_EPy, m_HcalHit_truthMC_EPz, m_HcalHit_truthMC_weight;

Reconstruction/RecPFACyber/include/Objects/Calo1DCluster.h

@@ -28,7 +28,7 @@ namespace Cyber{
 
     bool isNeighbor(const Cyber::CaloUnit* m_bar) const; 
     bool inCluster(const Cyber::CaloUnit* iBar) const;
-    void sortByPos() { std::sort(Bars.begin(), Bars.end()); }
+    void sortByPos() { std::sort(Bars.begin(), Bars.end(), compPos); }
 
     double getEnergy() const; 
     TVector3 getPos() const;
@@ -65,6 +65,7 @@ namespace Cyber{
 
     int getDlayer() const { if(Bars.size()>0) return Bars[0]->getDlayer(); return -99;  }
     int getSlayer() const { if(Bars.size()>0) return Bars[0]->getSlayer(); return -99;  }
+    int getSystem() const { if(Bars.size()>0) return Bars[0]->getSystem(); return -99;  }
     std::vector< std::vector<int> > getTowerID() const { return towerID; }
     
   private: 
@@ -73,12 +74,15 @@ namespace Cyber{
     double Energy;
     TVector3 pos;
 
-    std::vector< std::vector<int> > towerID; //[module, stave]
+    std::vector< std::vector<int> > towerID; //[system, module, stave, part]
 
     std::vector< const Cyber::Calo1DCluster* > CousinClusters;
     std::vector< const Cyber::Calo1DCluster* > ChildClusters;
 
     std::vector< std::pair<edm4hep::MCParticle, float> > MCParticleWeight;
+
+    static bool compPos( const Cyber::CaloUnit* hit1, const Cyber::CaloUnit* hit2 )
+      { return *hit1 < *hit2; }
   };
 }
 #endif