.gitlab-ci.yml

@@ -37,40 +37,47 @@ stages:
 .envvar_template:
   variables:
     CEPCSW_LCG_RELEASE: LCG
-    CEPCSW_LCG_PLATFORM: x86_64-centos7-gcc11-opt
-    CEPCSW_LCG_VERSION: 103.0.2
+    CEPCSW_LCG_PLATFORM: x86_64-el9-gcc11-opt
+    CEPCSW_LCG_VERSION: 105.0.0
 
-.build_template:
+# for k8s
+.build_template_k8s:
   extends: .envvar_template
+  image: cepc/cepcsw-cvmfs:el9
   stage: build
+  tags:
+    - k8s # using k8s as runner
   script:
+    - sed -i 's%^CVMFS_HTTP_PROXY=.*%CVMFS_HTTP_PROXY=http://squid-01.ihep.ac.cn:3128%' /etc/cvmfs/default.local
+    - for repo in sft.cern.ch geant4.cern.ch cepcsw.ihep.ac.cn; do [ -d "/cvmfs/$repo" ] || mkdir /cvmfs/$repo; sudo mount -t cvmfs $repo /cvmfs/$repo; done
     - bash ./.build.ci.sh
     - bash ./.test.ci.sh
 
 ##############################################################################
-# Build & Test in CentOS 7 (LCG)
+# Build & Test in k8s (LCG)
 ##############################################################################
-build:lcg:el7:
-  extends: .build_template
-  tags:
-    - centos7
+build:lcg:el9:k8s:
+  extends: .build_template_k8s
+  rules:
+    - if: $CI_PIPELINE_SOURCE == 'merge_request_event'
+      when: manual
   artifacts:
-    paths:
-      - InstallArea
-      - TDR_o1_v01.tgeo.root
-      - TDR_o1_v02.tgeo.root
     reports:
       junit: build.${CEPCSW_LCG_VERSION}.${CEPCSW_LCG_PLATFORM}/cepcsw-ctest-result.xml
 
 ##############################################################################
 # Build the docs
 ##############################################################################
-build:docs:
-  stage: build 
-  tags:
-    - centos7
+build:docs:k8s:
+  extends: .build_template_k8s
+  image: sphinxdoc/sphinx
   script:
-    - bash ./.build.ci.sh docs
+    - bash build-docs.sh
   artifacts:
     paths:
       - docs/build/html/
+  rules:
+    - if: $CI_PIPELINE_SOURCE == "merge_request_event"
+      changes:
+      - docs/**/*
+

Analysis/AnalysisPID/AnalysisPID.py

+import os, sys
+from Gaudi.Configuration import *
+
+########### k4DataSvc ####################
+from Configurables import k4DataSvc
+podioevent = k4DataSvc("EventDataSvc", input="track.root")
+##########################################
+
+########## CEPCSWData ################# 
+cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest"
+#######################################
+
+
+########## Podio Input ###################
+from Configurables import PodioInput
+inp = PodioInput("InputReader")
+inp.collections = [ "CompleteTracks", 
+                    "CompleteTracksParticleAssociation",
+                    "RecTofCollection",
+                    "DndxTracks" ]
+##########################################
+
+from Configurables import AnalysisPIDAlg
+anaPID = AnalysisPIDAlg("AnalysisPIDAlg")
+anaPID.OutputFile = "./pid.root"
+
+##############################################################################
+# POD I/O
+##############################################################################
+
+
+########################################
+
+from Configurables import ApplicationMgr
+ApplicationMgr( 
+    TopAlg=[inp, anaPID ],
+    EvtSel="NONE",
+    EvtMax=-1,
+    ExtSvc=[podioevent],
+    #OutputLevel=DEBUG
+)

Analysis/AnalysisPID/CMakeLists.txt

+# gaudi_add_header_only_library(AnalysisPIDLib)
+
+# Modules
+gaudi_add_module(AnalysisPID
+    SOURCES src/AnalysisPIDAlg.cpp
+                 LINK Gaudi::GaudiAlgLib
+                      Gaudi::GaudiKernel
+                      DataHelperLib
+                      DetSegmentation
+                      DetInterface
+                      ${GSL_LIBRARIES}
+                      ${GEAR_LIBRARIES}
+                      ${LCIO_LIBRARIES}
+		      EDM4CEPC::edm4cepc EDM4CEPC::edm4cepcDict
+                      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+                      k4FWCore::k4FWCore
+)
+
+target_include_directories(AnalysisPID PUBLIC
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include
+  $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
+install(TARGETS AnalysisPID
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)

Analysis/AnalysisPID/src/AnalysisPIDAlg.cpp

+#include "AnalysisPIDAlg.h"
+
+#include "GaudiKernel/DataObject.h"
+#include "GaudiKernel/IHistogramSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/SmartDataPtr.h"
+
+#include "DetInterface/IGeomSvc.h"
+#include "DataHelper/HelixClass.h"
+
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include <cmath>
+#include "UTIL/ILDConf.h"
+#include "DetIdentifier/CEPCConf.h"
+#include "UTIL/CellIDEncoder.h"
+
+using namespace edm4hep;
+
+DECLARE_COMPONENT( AnalysisPIDAlg )
+
+//------------------------------------------------------------------------------
+AnalysisPIDAlg::AnalysisPIDAlg( const std::string& name, ISvcLocator* pSvcLocator )
+    : Algorithm( name, pSvcLocator ) {
+  declareProperty("CompleteTracks", _FultrkCol, "handler of the input complete track collection");
+  declareProperty("CompleteTracksParticleAssociation", _FultrkParAssCol, "handler of the input track particle association collection");
+  declareProperty("DndxTracks", _inDndxColHdl, "handler of the collection of dN/dx tracks");
+  declareProperty("RecTofCollection", _inTofColHdl, "handler of the collection of tof tracks");
+  // output
+  declareProperty("OutputFile", m_outputFile = "pid.root", "output file name");
+}
+
+//------------------------------------------------------------------------------
+StatusCode AnalysisPIDAlg::initialize(){
+
+  info() << "Booking Ntuple" << endmsg;
+  m_file = new TFile(m_outputFile.value().c_str(),"RECREATE");
+  m_tree = new TTree("pid","pid");
+
+  m_tree->Branch("Nevt",&_nEvt,"Nevt/I");
+  m_tree->Branch("Ndndxtrk",&Ndndxtrk,"Ndndxtrk/I");
+  m_tree->Branch("Ntoftrk",&Ntoftrk,"Ntoftrk/I");
+  m_tree->Branch("Nfullass",&Nfullass,"Nfullass/I");
+  m_tree->Branch("Nfulltrk",&Nfulltrk,"Nfulltrk/I");
+  m_tree->Branch("tpcidx",&tpcidx);
+  m_tree->Branch("tofidx",&tofidx);
+  m_tree->Branch("matchedtpc",&matchedtpc);
+  m_tree->Branch("matchedtof",&matchedtof);
+  m_tree->Branch("truthidx",&truthidx);
+  
+  m_tree->Branch("tof_chi2s",&tof_chi2s);
+  m_tree->Branch("tof_chis",&tof_chis);
+  m_tree->Branch("tof_expt",&tof_expt);
+  m_tree->Branch("tof_meast",&tof_meast);
+  m_tree->Branch("tof_measterr",&tof_measterr);
+
+  m_tree->Branch("tpc_chi2s",&tpc_chi2s);
+  m_tree->Branch("tpc_chis",&tpc_chis);
+  m_tree->Branch("tpc_expdndxs",&tpc_expdndxs);
+  m_tree->Branch("tpc_measdndx",&tpc_measdndx);
+  m_tree->Branch("tpc_measdndxerr",&tpc_measdndxerr);
+
+  m_tree->Branch("tot_chi2s",&tot_chi2s);
+  m_tree->Branch("recoPDG",&recoPDG);
+  m_tree->Branch("tpcrecoPDG",&tpcrecoPDG);
+  m_tree->Branch("tofrecoPDG",&tofrecoPDG);
+
+  //gen trk parameters
+  m_tree->Branch("genpx",&genpx);
+  m_tree->Branch("genpy",&genpy);
+  m_tree->Branch("genpz",&genpz);
+  m_tree->Branch("genE",&genE);
+  m_tree->Branch("genp",&genp);
+  m_tree->Branch("genM",&genM);
+  m_tree->Branch("genphi",&genphi);
+  m_tree->Branch("gentheta",&gentheta);
+  m_tree->Branch("endx",&endx);
+  m_tree->Branch("endy",&endy);
+  m_tree->Branch("endz",&endz);
+  m_tree->Branch("endr",&endr);
+  m_tree->Branch("PDG",&PDG);
+  m_tree->Branch("genstatus",&genstatus);
+  m_tree->Branch("simstatus",&simstatus);
+  m_tree->Branch("isdecayintrker",&isdecayintrker);
+  m_tree->Branch("iscreatedinsim",&iscreatedinsim);
+  m_tree->Branch("isbackscatter",&isbackscatter);
+  m_tree->Branch("isstopped",&isstopped);
+
+  _nEvt = 0;
+  return StatusCode::SUCCESS;
+}
+
+//------------------------------------------------------------------------------
+StatusCode AnalysisPIDAlg::execute(){
+  
+  const edm4hep::TrackCollection* trkCol = nullptr;
+  const edm4hep::RecDqdxCollection* dndxCols = nullptr;
+  const edm4hep::RecTofCollection* tofCols = nullptr;
+  const edm4hep::MCRecoTrackParticleAssociationCollection* fultrkparassCols = nullptr;
+  
+  ClearVars();
+
+  try {
+    trkCol = _FultrkCol.get();
+  }
+  catch ( GaudiException &e ) {
+    debug() << "Complete track collection " << _FultrkCol.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+    Nfulltrk = -1;
+  }
+  if ( trkCol->size() == 0 ) {
+    debug() << "No full track found in event " << _nEvt << endmsg;
+    Nfulltrk = 0;
+  }
+  else{
+    Nfulltrk = trkCol->size();
+  }
+
+  try {
+    fultrkparassCols = _FultrkParAssCol.get();
+  }
+  catch ( GaudiException &e ) {
+    debug() << "Complete track particle association collection " << _FultrkParAssCol.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+    Nfullass = -1;
+    m_tree->Fill();
+    _nEvt++;
+    return StatusCode::SUCCESS;
+  }
+
+  try {
+    dndxCols = _inDndxColHdl.get();
+  }
+  catch ( GaudiException &e ) {
+    debug() << "DndxTrack collection " << _inDndxColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+    Ndndxtrk = -1;
+  }
+
+  if ( dndxCols->size() == 0 ) {
+    debug() << "No dndx track found in event " << _nEvt << endmsg;
+    Ndndxtrk = 0;
+  }
+  else{
+    Ndndxtrk = dndxCols->size();
+  }
+
+  try {
+    tofCols = _inTofColHdl.get();
+  }
+  catch ( GaudiException &e ) {
+    debug() << "TofTrack collection " << _inTofColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+    Ntoftrk = -1;
+  }
+
+  if ( tofCols->size() == 0 ) {
+    debug() << "No tof track found in event " << _nEvt << endmsg;
+    Ntoftrk = 0;
+  }
+  else{
+    Ntoftrk = tofCols->size();
+  }
+
+  if ( fultrkparassCols->size() == 0 ) {
+    debug() << "No full track particle association found in event " << _nEvt << endmsg;
+    Nfullass = 0;
+    m_tree->Fill();
+    _nEvt++;
+    return StatusCode::SUCCESS;
+  }
+  else{
+    Nfullass = fultrkparassCols->size();
+  }
+  
+  info() << "normal eventID: " << _nEvt << endmsg;
+  for(auto track : *trkCol){
+      //truth association match
+      max_weight = -999.;
+      max_weight_idx = -1;
+      ass_idx = 0;
+      for (auto ass : *fultrkparassCols) {
+          if (ass.getRec() == track) {
+              weight = ass.getWeight();
+              if (weight > max_weight) {
+                  max_weight = weight;
+                  max_weight_idx = ass_idx;
+              }
+          }
+          ass_idx++;
+      }
+      truthidx.push_back(max_weight_idx);
+//      if (max_weight_idx < 0) continue;
+      p1=fultrkparassCols->at(max_weight_idx).getSim().getMomentum()[0];
+      p2=fultrkparassCols->at(max_weight_idx).getSim().getMomentum()[1];
+      p3=fultrkparassCols->at(max_weight_idx).getSim().getMomentum()[2];
+      genpx.push_back(p1);
+      genpy.push_back(p2);
+      genpz.push_back(p3);
+      genp.push_back(std::sqrt(p1*p1 + p2*p2 + p3*p3));
+      genE.push_back(fultrkparassCols->at(max_weight_idx).getSim().getEnergy());
+      genM.push_back(fultrkparassCols->at(max_weight_idx).getSim().getMass());
+      genphi.push_back(std::atan2(p2,p1));
+      gentheta.push_back(std::acos(p3/std::sqrt(p1*p1 + p2*p2 + p3*p3)));
+      x1=fultrkparassCols->at(max_weight_idx).getSim().getEndpoint()[0];
+      y1=fultrkparassCols->at(max_weight_idx).getSim().getEndpoint()[1];
+      endx.push_back(x1);
+      endy.push_back(y1);
+      endz.push_back(fultrkparassCols->at(max_weight_idx).getSim().getEndpoint()[2]);
+      endr.push_back(std::sqrt(x1*x1 + y1*y1));
+      PDG.push_back(fultrkparassCols->at(max_weight_idx).getSim().getPDG());
+      genstatus.push_back(fultrkparassCols->at(max_weight_idx).getSim().getGeneratorStatus());
+      simstatus.push_back(fultrkparassCols->at(max_weight_idx).getSim().getSimulatorStatus());
+      isdecayintrker.push_back(fultrkparassCols->at(max_weight_idx).getSim().isDecayedInTracker());//getSimulatorStatus().isDecayInTracker();
+      iscreatedinsim.push_back(fultrkparassCols->at(max_weight_idx).getSim().isCreatedInSimulation());
+      isbackscatter.push_back(fultrkparassCols->at(max_weight_idx).getSim().isBackscatter());
+      isstopped.push_back(fultrkparassCols->at(max_weight_idx).getSim().isStopped());
+
+      //find corresponding dndx track
+      edm4hep::RecDqdx dndxtrk;
+      dndx_index = -1;
+      matched1 = false;
+      for(int i=0; i<Ndndxtrk; i++){
+          if ( dndxCols->at(i).getTrack() == track ){
+              dndxtrk = dndxCols->at(i);
+              dndx_index = i;
+              matched1 = true;
+              break;
+          }
+      }
+      tpcidx.push_back(dndx_index);
+      tpc_chi2s_1.clear();tpc_expdndxs_1.clear();tpc_chis_1.clear();
+      tpcdndx=-1;tpcdndxerr=-1;
+      if( matched1 ) {
+          tpcdndx = dndxtrk.getDQdx().value;
+          tpcdndxerr = dndxtrk.getDQdx().error;
+          debug()<<"tpc_measdndx = "<<dndxtrk.getDQdx().value<<endmsg;
+          for (int idx=0;idx<5;idx++) {
+              double tpc_chi2 = dndxtrk.getHypotheses(idx).chi2;
+              double tpc_expdndx = dndxtrk.getHypotheses(idx).expected;
+              double tpc_chi = ( tpcdndx - tpc_expdndx ) / tpcdndxerr;
+              tpc_chi2s_1.push_back(tpc_chi2);
+              tpc_chis_1.push_back(tpc_chi);
+              tpc_expdndxs_1.push_back(tpc_expdndx);
+              debug()<<" idx : "<< idx <<" tpc_chi2 : "<< tpc_chi2 <<endmsg;    
+          }
+      }
+      else{
+          tpc_chi2s_1={-999,-999,-999,-999,-999};
+          tpc_chis_1={-999,-999,-999,-999,-999};
+          tpc_expdndxs_1={-1,-1,-1,-1,-1};
+      }
+      tpc_measdndx.push_back(tpcdndx);
+      tpc_measdndxerr.push_back(tpcdndxerr);
+      matchedtpc.push_back(matched1);
+      tpc_chi2s.push_back(tpc_chi2s_1);
+      tpc_chis.push_back(tpc_chis_1);
+      tpc_expdndxs.push_back(tpc_expdndxs_1);
+
+      //find corresponding tof track
+      edm4hep::RecTof toftrk;
+      tof_index = -1;
+      matched2 = false;
+      for(int i=0; i<Ntoftrk; i++){
+          if ( tofCols->at(i).getTrack() == track ){
+              toftrk = tofCols->at(i);
+              tof_index = i;
+              matched2 = true;
+              break;
+          }
+      }
+      tofidx.push_back(tof_index);
+      tof_chi2s_1.clear();tof_expt_1.clear();tof_chis_1.clear();
+      toft=-1;tofterr=-1;
+      if ( matched2 ) {
+          toft = toftrk.getTime();
+          std::array<float, 5> tofexpts = toftrk.getTimeExp();
+          tofterr = toftrk.getSigma();
+          debug()<<"tof_meast = "<<toftrk.getTime()<<endmsg;
+          for (int idx=0;idx<5;idx++){
+              double tof_chi = ( toft - tofexpts[idx] ) / tofterr;
+              double tof_chi2 = tof_chi*tof_chi;
+              tof_chi2s_1.push_back(tof_chi2);
+              tof_chis_1.push_back(tof_chi);
+              tof_expt_1.push_back(tofexpts[idx]);
+              debug()<<" idx : "<< idx <<" tof_chi2 : "<< tof_chi2 <<endmsg;    
+          }//end loop over masses
+      }
+      else{
+          tof_chi2s_1={-999,-999,-999,-999,-999};
+          tof_chis_1={-999,-999,-999,-999,-999};
+          tof_expt_1={-1,-1,-1,-1,-1};
+      }
+      tof_meast.push_back(toft);
+      tof_measterr.push_back(tofterr);
+      matchedtof.push_back(matched2);
+      tof_chi2s.push_back(tof_chi2s_1);
+      tof_chis.push_back(tof_chis_1);
+      tof_expt.push_back(tof_expt_1);
+
+      tot_chi2s_1.clear();
+      recpdg = 9999;tpcrecpdg = 9999;tofrecpdg = 9999;
+      if(matched1){ 
+          minchi2idx = std::distance(tpc_chi2s_1.begin(), std::min_element(tpc_chi2s_1.begin(), tpc_chi2s_1.end()));
+          tpcrecpdg = PDGIDs.at(minchi2idx);
+      }
+      if(matched2){
+          minchi2idx = std::distance(tof_chi2s_1.begin(), std::min_element(tof_chi2s_1.begin(), tof_chi2s_1.end()));
+          tofrecpdg = PDGIDs.at(minchi2idx);
+      }
+      if(matched1 || matched2){
+          if(matched1 && matched2){
+          for(int i=0;i<5;i++){
+              tot_chi2s_1.push_back(tpc_chi2s_1[i]+tof_chi2s_1[i]);
+          }
+          }
+          if(matched1 && !matched2) tot_chi2s_1=tpc_chi2s_1;
+          if(!matched1 && matched2) tot_chi2s_1=tof_chi2s_1;
+          minchi2idx = std::distance(tot_chi2s_1.begin(), std::min_element(tot_chi2s_1.begin(), tot_chi2s_1.end()));
+          recpdg = PDGIDs.at(minchi2idx);
+      }
+      else{ 
+          tot_chi2s_1={-999,-999,-999,-999,-999};
+      }
+      int charge = track.getTrackStates(1).omega/fabs(track.getTrackStates(1).omega);
+      tot_chi2s.push_back(tot_chi2s_1);
+      recoPDG.push_back(charge*recpdg);
+      tpcrecoPDG.push_back(charge*tpcrecpdg);
+      tofrecoPDG.push_back(charge*tofrecpdg);
+  }
+m_tree->Fill();
+_nEvt++;
+return StatusCode::SUCCESS;
+}// end execute
+
+//------------------------------------------------------------------------------
+StatusCode AnalysisPIDAlg::finalize(){
+  debug() << "Finalizing..." << endmsg;
+  m_file->cd();
+  m_tree->Write();
+  return StatusCode::SUCCESS;
+}

Analysis/AnalysisPID/src/AnalysisPIDAlg.h

+#ifndef AnalysisPIDAlg_h
+#define AnalysisPIDAlg_h 1
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/Algorithm.h"
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+
+#include "edm4hep/TrackerHit.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "edm4hep/RecDqdx.h"
+#include "edm4hep/RecDqdxCollection.h"
+#include "edm4cepc/RecTofCollection.h"
+#include "edm4hep/TrackState.h"
+#include "edm4hep/Vector3d.h"
+
+#include "TFile.h"
+#include "TTree.h"
+#include <random>
+#include "GaudiKernel/NTuple.h"
+
+class AnalysisPIDAlg : public Algorithm {
+ public:
+  // Constructor of this form must be provided
+  AnalysisPIDAlg( 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::TrackCollection> _FultrkCol{"CompleteTracks", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> _FultrkParAssCol{"CompleteTracksParticleAssociation", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::RecDqdxCollection> _inDndxColHdl{"DndxTracks", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::RecTofCollection> _inTofColHdl{"RecTofCollection", Gaudi::DataHandle::Reader, this};
+
+  Gaudi::Property<std::string> m_outputFile{this, "OutputFile", "pid.root"};
+
+  std::vector<double> genpx, genpy, genpz, genE, genp, genM, gentheta, genphi, endx, endy, endz, endr;
+  std::vector<int> PDG, genstatus, simstatus, recoPDG, tpcrecoPDG, tofrecoPDG;
+  std::vector<bool> isdecayintrker, iscreatedinsim, isbackscatter, isstopped;
+  std::vector<bool> matchedtpc, matchedtof;
+  std::vector<int> truthidx, tpcidx, tofidx;
+
+  std::vector<std::vector<double>> tof_chi2s, tof_expt;
+  std::vector<std::vector<double>> tpc_chi2s, tpc_expdndxs;
+  std::vector<std::vector<double>> tof_chis, tpc_chis, tot_chi2s;
+  std::vector<double> tof_meast, tof_measterr;
+  std::vector<double> tpc_measdndx, tpc_measdndxerr;
+  std::vector<double> tof_chi2s_1, tof_expt_1;
+  std::vector<double> tpc_chi2s_1, tpc_expdndxs_1;
+  std::vector<double> tof_chis_1, tpc_chis_1, tot_chi2s_1;
+  
+  int _nEvt;
+  const std::map<int, int> PDGIDs = { 
+    {0, -11},
+    {1, -13},
+    {2, 211},
+    {3, 321},
+    {4, 2212},
+  };
+
+  double tpcdndx, tpcdndxerr, toft, tofterr;
+  int Ndndxtrk, Ntoftrk, Nfulltrk, Nfullass;
+  double max_weight, weight;
+  int max_weight_idx, ass_idx, dndx_index, tof_index; 
+  double p1, p2, p3, x1, y1;
+  bool matched1, matched2;
+  int recpdg, tpcrecpdg, tofrecpdg, minchi2idx;
+
+  TFile* m_file;
+  TTree* m_tree;
+
+  void ClearVars(){
+
+    Ndndxtrk = 0;
+    Ntoftrk = 0;
+    Nfullass = 0;
+    Nfulltrk = 0;
+
+    genpx.clear(); genpy.clear(); genpz.clear(); genE.clear();
+    genp.clear(); genM.clear(); gentheta.clear(); genphi.clear();
+    endx.clear(); endy.clear(); endz.clear(); endr.clear();
+    PDG.clear(); genstatus.clear(); simstatus.clear();
+    recoPDG.clear(); tpcrecoPDG.clear(); tofrecoPDG.clear();
+    isdecayintrker.clear(); iscreatedinsim.clear();
+    isbackscatter.clear(); isstopped.clear();
+
+    tof_chi2s.clear();
+    tof_chis.clear();
+    tof_meast.clear();
+    tof_measterr.clear();
+    tof_expt.clear();
+
+    tpc_chi2s.clear();
+    tpc_chis.clear();
+    tpc_expdndxs.clear();
+    tpc_measdndx.clear();
+    tpc_measdndxerr.clear();
+
+    tot_chi2s.clear();
+
+    matchedtpc.clear();
+    matchedtof.clear();
+    truthidx.clear();
+    tpcidx.clear();
+    tofidx.clear();
+  }
+
+};
+#endif

Analysis/CMakeLists.txt

@@ -4,3 +4,5 @@ add_subdirectory(TrackInspect)
 add_subdirectory(DumpEvent)
 add_subdirectory(ReadDigi)
 add_subdirectory(JetClustering)
+add_subdirectory(GenMatch)
+add_subdirectory(AnalysisPID)

Analysis/GenMatch/CMakeLists.txt

+# Set the CMake module path if necessary
+
+
+# Modules
+gaudi_add_module(GenMatch
+                 SOURCES src/GenMatch.cpp
+                 LINK GearSvc
+                      Gaudi::GaudiKernel
+                      k4FWCore::k4FWCore
+                      DetInterface
+                      DataHelperLib
+                      ${GEAR_LIBRARIES} 
+                      ${GSL_LIBRARIES} 
+                      ${ROOT_LIBRARIES}
+                      ${FASTJET_LIBRARY}
+                      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+)
+install(TARGETS GenMatch
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)

Analysis/GenMatch/jcls_genmatch.py

+import os, sys
+from Gaudi.Configuration import *
+
+########### k4DataSvc ####################
+from Configurables import k4DataSvc
+#podioevent = k4DataSvc("EventDataSvc", input="./FullSim_Samples/nnHgg/Rec_TDR_o1_v01_E240_nnh_gg_test_003.root")
+podioevent = k4DataSvc("EventDataSvc", input="/publicfs/cms/user/wangzebing/CEPC/CEPCSW_tdr24.9.1/CEPCSW/FullSim_samples/Rec_TDR_o1_v01_E240_nnh_gg.root")
+##########################################
+
+########## CEPCSWData ################# 
+cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest"
+#######################################
+
+
+########## Podio Input ###################
+from Configurables import PodioInput
+inp = PodioInput("InputReader")
+inp.collections = [ "CyberPFO", "MCParticle" ]
+##########################################
+
+
+
+from Configurables import GenMatch
+genmatch = GenMatch("GenMatch")
+genmatch.nJets = 2
+genmatch.R = 0.6
+genmatch.OutputFile = "./RecJets_TDR_o1_v01_E240_nnh_gg.root"
+#genmatch.OutputFile = "./FullSim_samples/RecJets_TDR_o1_v01_E240_nnh_gg_CalHits.root"
+
+##############################################################################
+# POD I/O
+##############################################################################
+
+
+########################################
+
+from Configurables import ApplicationMgr
+ApplicationMgr( 
+    TopAlg=[inp, genmatch ],
+    EvtSel="NONE",
+    EvtMax=3,
+    ExtSvc=[podioevent],
+    #OutputLevel=DEBUG
+)

Analysis/GenMatch/src/GenMatch.h

+#ifndef GenMatch_h
+#define GenMatch_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"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/MCParticleCollectionData.h"
+
+class GenMatch : public Algorithm {
+ public:
+  // Constructor of this form must be provided
+  GenMatch( 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{"CyberPFO", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::MCParticleCollection> m_MCParticleGenHdl{"MCParticle", 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", "GenMatch.root"};
+
+  int _nEvt;
+  int _particle_index;
+  int _jet_index;
+  int _nparticles;
+  int jet1_nparticles;
+  int jet2_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, nparticles, jet1_GENMatch_id, jet2_GENMatch_id;
+  double constituents_E1tot;
+  double constituents_E2tot;
+  double ymerge[6];
+  double mass;
+  double mass_gen_match, jet1_GENMatch_mindR, jet2_GENMatch_mindR;
+  double mass_allParticles;
+  //double PFO_Energy_muon[50];
+  //double PFO_Energy_Charge[50];
+  //double PFO_Energy_Neutral[50];
+  //double PFO_Energy_Neutral_singleCluster[600];
+  //double PFO_Energy_Neutral_singleCluster_R[600];
+  int _particle_index_muon, _particle_index_Charge, _particle_index_Neutral, _particle_index_Neutral_singleCluster;
+  int particle_index_muon, particle_index_Charge, particle_index_Neutral, particle_index_Neutral_singleCluster;
+  //double GEN_PFO_Energy_muon[50];
+  //double GEN_PFO_Energy_pipm[120];
+  //double GEN_PFO_Energy_pi0[120];
+  //double GEN_PFO_Energy_ppm[120];
+  //double GEN_PFO_Energy_kpm[120];
+  //double GEN_PFO_Energy_kL[120];
+  //double GEN_PFO_Energy_n[50];
+  //double GEN_PFO_Energy_e[50];
+  //double GEN_PFO_Energy_gamma[120];
+  std::vector<double> PFO_Energy_muon, PFO_Energy_muon_GENMatch_dR, PFO_Energy_muon_GENMatch_E, PFO_Energy_Charge, PFO_Energy_Charge_Ecal, PFO_Energy_Charge_Hcal, PFO_Energy_Charge_GENMatch_dR, PFO_Energy_Charge_GENMatch_E, PFO_Energy_Neutral, PFO_Energy_Neutral_singleCluster, PFO_Energy_Neutral_singleCluster_R;
+  std::vector<int> PFO_Energy_Charge_GENMatch_ID, PFO_Energy_muon_GENMatch_ID;
+  std::vector<double> GEN_PFO_Energy_muon, GEN_PFO_Energy_pipm, GEN_PFO_Energy_pi0, GEN_PFO_Energy_ppm, GEN_PFO_Energy_kpm, GEN_PFO_Energy_kL, GEN_PFO_Energy_n, GEN_PFO_Energy_e, GEN_PFO_Energy_gamma;
+  std::vector<std::vector<double>> PFO_Hits_Charge_E, PFO_Hits_Charge_R, PFO_Hits_Charge_theta, PFO_Hits_Charge_phi;
+  std::vector<std::vector<double>> PFO_Hits_Neutral_E, PFO_Hits_Neutral_R, PFO_Hits_Neutral_theta, PFO_Hits_Neutral_phi;
+  //double PFO_Energy_muon[500];
+
+  int _gen_particle_index;
+  double barrelRatio;
+  double ymerge_gen[6];
+  double GEN_mass;
+  double GEN_ISR_E, GEN_ISR_pt;
+  double GEN_inv_E, GEN_inv_pt;
+  double GEN_jet1_px, GEN_jet1_py, GEN_jet1_pz, GEN_jet1_E;
+  double GEN_jet2_px, GEN_jet2_py, GEN_jet2_pz, GEN_jet2_E;
+  double GEN_jet1_costheta, GEN_jet1_phi;
+  double GEN_jet2_costheta, GEN_jet2_phi;
+  double GEN_jet1_pt, GEN_jet2_pt;
+  int GEN_jet1_nconstituents, GEN_jet2_nconstituents, GEN_nparticles;
+  double GEN_constituents_E1tot;
+  double GEN_constituents_E2tot;
+  int _gen_particle_gamma_index, _gen_particle_pipm_index, _gen_particle_pi0_index, _gen_particle_ppm_index, _gen_particle_kpm_index, _gen_particle_kL_index, _gen_particle_n_index, _gen_particle_e_index;
+  int GEN_particle_gamma_index, GEN_particle_pipm_index, GEN_particle_pi0_index, GEN_particle_ppm_index, GEN_particle_kpm_index, GEN_particle_kL_index, GEN_particle_n_index, GEN_particle_e_index;
+
+  void CleanVars();
+
+};
+#endif
+

Analysis/MaterialScan/README.md

+# Howto scan the material
+
+## Introduction
+Use DD4hep tool [materialScan](https://dd4hep.web.cern.ch/dd4hep/reference/materialScan_8cpp_source.html) to get the material passing through for a geatino partical (i.e. a straight line) at a given angle.
+
+## Material Scan
+
+Select the dd4hep geometry files you wish to scan, and run the given script.
+
+```bash
+# Scan the material for a specific detector
+root -l 'src/ScanMaterial.cpp("input_file.xml", "output_file.txt", bins, world_size)'
+
+# example: scan the material for the BeamPipe detector
+root -l 'src/ScanMaterial.cpp("Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyBeamPipe.xml", "BeamPipe.txt", 100, 10000)'
+```
+
+The "FILE* outFile" truncates the "stdout" output stream to the output_file.txt ("MaterialScanLogs.txt" by default), so you can't see any output. Be patient and wait for the file to be created, it may take 1-2 minutes depending on the detectors & the number of bins.
+
+*bins* is the number of bins for theta and phi, *world_size* is the size of the scanning area (mm).
+
+The output_file.txt contains all the material passing through for a straight line at the given angles with the size of bins of theta and phi.
+
+## Draw the material budget for single xml file
+
+We give an example script [options/extract_x0.py](options/extract_x0.py) to draw the material budget from the scan file generated above (e.x. MaterialScanLogs.txt).
+
+```bash
+# Draw the material budget for single xml file
+# --input_file: the scan file generated above
+# --output_file: the output file name (default: material_budget.png)
+# --detector: the detector name to be used as the title of the plot (default: all_world)
+# --bins: the bins of theta and phi (default: 100), the theta-bins is set to "bins"+1, the phi-bins is set to "bins"*2+1
+# --bias: the bias used to cut the theta range to [bias, bins-bias] (default: 0)
+python options/extract_x0.py \
+    --input_file MaterialScanLogs.txt \
+    --output_file material_budget.png \
+    --detector all_world \
+    --bins 100 \
+    --bias 0
+```
+
+Make sure the property "bins" is the same as the one used in the [Material Scan](#material-scan).
+
+The output is a png file shows the material budget at different theta and phi, and the average X0 for theta/phi.
+
+## Draw the material budget for multiple xml files
+
+For multiple xml files, we can use the script [options/extract_x0_multi.py](options/extract_x0_multi.py) to draw the material budget.
+
+```bash
+# Draw the material budget for multiple xml files
+# --input_files: the scan files generated above, files names separated with spaces
+# --labels: the labels name to be used as the legend of the plot
+# --output_file: the output file name (default: accumulated_material_budget.png)
+# --bins: the bins of theta and phi (default: 100), the theta-bins is set to "bins"+1, the phi-bins is set to "bins"*2+1
+# --bias: the bias used to cut the theta range to [bias, bins-bias] (default: 0)
+python options/extract_x0_multi.py \
+    --input_file BeamPipe.txt Lumical.txt VXD.txt FTD.txt SIT.txt \
+    --labels BeamPipe Lumical VXD FTD SIT \
+    --output_file accumulated_material_budget.png \
+    --bins 100 \
+    --bias 0
+```
+
+The output is a png file shows the average X0 for theta/phi with different labels.
\ No newline at end of file

Analysis/MaterialScan/options/extract_multi_x0.py

+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import argparse
+
+def extract_x0_values(filename):
+    """Extract x0 values from a single file"""
+    x0_values = []
+    counter = 0
+    with open(filename, 'r') as file:
+        for line in file:
+            if '|     0 Average Material' in line:
+                counter += 1
+                items = [item.strip() for item in line.split('|')]
+                x0_value = float(items[1].split()[10])
+                x0_values.append(x0_value)
+                print(f"Processing file {filename}: ", counter, " x0: ", x0_value)
+    return x0_values
+
+def process_multiple_files(file_list, bins, bias):
+    """Process multiple files and return their x0 matrices"""
+    theta_bins = bins + 1
+    phi_bins = bins*2 + 1
+    range_theta = [bias, theta_bins-bias]
+    
+    matrices = []
+    for file in file_list:
+        x0_values = extract_x0_values(file)
+        x0_matrix = np.array(x0_values).reshape(theta_bins, phi_bins)[range_theta[0]:range_theta[1], :]
+        matrices.append(x0_matrix)
+    
+    return matrices
+
+def plot_accumulated_projections(matrices, output_file, bins, bias, labels):
+    """Plot accumulated projections"""
+    theta_bins = bins + 1
+    phi_bins = bins*2 + 1
+    range_theta = [bias, theta_bins-bias]
+    
+    # Create angle arrays
+    phi = np.linspace(-180, 180, phi_bins)
+    theta = np.linspace(range_theta[0]*180/theta_bins, range_theta[1]*180/theta_bins, range_theta[1]-range_theta[0])
+    
+    # Create figure
+    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
+    
+    # Accumulation data
+    accumulated_theta = np.zeros_like(theta)
+    accumulated_phi = np.zeros_like(phi)
+    
+    colors = plt.cm.viridis(np.linspace(0, 1, len(matrices)))
+    
+    # Plot phi projection
+    for i, matrix in enumerate(matrices):
+        phi_projection = np.sum(matrix, axis=0) / (range_theta[1]-range_theta[0])
+        accumulated_phi += phi_projection
+        ax1.fill_between(phi, accumulated_phi, accumulated_phi - phi_projection, 
+                        label=labels[i], color=colors[i], alpha=0.6)
+    
+    # Plot theta projection
+    for i, matrix in enumerate(matrices):
+        theta_projection = np.sum(matrix, axis=1) / phi_bins
+        accumulated_theta += theta_projection
+        ax2.fill_between(theta, accumulated_theta, accumulated_theta - theta_projection,
+                        label=labels[i], color=colors[i], alpha=0.6)
+    
+    # Set phi projection plot
+    ax1.set_title('Accumulated Projection along Phi', fontsize=14, pad=10)
+    ax1.set_xlabel('Phi (degree)', fontsize=12)
+    ax1.set_ylabel('Accumulated X/X0', fontsize=12)
+    ax1.grid(True, linestyle='--', alpha=0.7)
+    ax1.legend(fontsize=10)
+    ax1.spines['top'].set_visible(False)
+    ax1.spines['right'].set_visible(False)
+    
+    # Set theta projection plot
+    ax2.set_title('Accumulated Projection along Theta', fontsize=14, pad=10)
+    ax2.set_xlabel('Theta (degree)', fontsize=12)
+    ax2.set_ylabel('Accumulated X/X0', fontsize=12)
+    ax2.grid(True, linestyle='--', alpha=0.7)
+    ax2.legend(fontsize=10)
+    ax2.spines['top'].set_visible(False)
+    ax2.spines['right'].set_visible(False)
+    
+    plt.suptitle('Detector Material Budget Accumulation Analysis', fontsize=16, y=0.95)
+    plt.tight_layout()
+    
+    plt.savefig(output_file, dpi=600, bbox_inches='tight')
+    plt.show()
+
+def main():
+    parser = argparse.ArgumentParser(description='Process multiple material scan data and generate accumulated distribution plots.')
+    parser.add_argument('--input_files', nargs='+', required=True,
+                        help='List of input files in order')
+    parser.add_argument('--labels', nargs='+', required=True,
+                        help='Labels for each input file')
+    parser.add_argument('--output_file', default='accumulated_material_budget.png',
+                        help='Output PNG file path')
+    parser.add_argument('--bins', type=int, default=100,
+                        help='theta bins is "bins"+1, phi bins is "bins"*2+1, (default: 100)')
+    parser.add_argument('--bias', type=int, default=0,
+                        help='Bias value for theta range (default: 0)')
+    
+    args = parser.parse_args()
+    
+    if len(args.input_files) != len(args.labels):
+        raise ValueError("Number of input files must match number of labels!")
+    
+    # Process all input files
+    matrices = process_multiple_files(args.input_files, args.bins, args.bias)
+    
+    # Plot accumulated graphs
+    plot_accumulated_projections(matrices, args.output_file, args.bins, args.bias, args.labels)
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file

Analysis/MaterialScan/options/extract_x0.py

+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import argparse
+
+def extract_x0_values(filename):
+    x0_values = []
+    counter = 0
+    with open(filename, 'r') as file:
+        for line in file:
+            if '|     0 Average Material' in line:
+                counter += 1
+                items = [item.strip() for item in line.split('|')]
+                x0_value = float(items[1].split()[10])
+                x0_values.append(x0_value)
+                print("processing: ", counter, " x0: ", x0_value)
+    return x0_values
+
+def plot_all_in_one(x0_values, output_file, bins, bias, detector):
+    """Plot all charts in one figure"""
+    theta_bins = bins + 1
+    phi_bins = bins*2 + 1
+    range_theta = [bias, theta_bins-bias]
+    # Adjust figure size and ratio, increase height to accommodate title
+    fig = plt.figure(figsize=(28, 8))
+    
+    # Adjust spacing between subplots and top margin
+    gs = fig.add_gridspec(1, 3, width_ratios=[1.2, 1, 1], wspace=0.25, top=0.85)
+    
+    # Add three subplots
+    ax1 = fig.add_subplot(gs[0])  # Heat map
+    ax2 = fig.add_subplot(gs[1])  # Phi projection
+    ax3 = fig.add_subplot(gs[2])  # Theta projection
+    
+    # Draw heat map
+    X0_matrix = np.array(x0_values).reshape(theta_bins, phi_bins)[range_theta[0]:range_theta[1], :]
+    phi = np.linspace(-180, 180, phi_bins)
+    theta = np.linspace(range_theta[0]*180/theta_bins, range_theta[1]*180/theta_bins, range_theta[1]-range_theta[0])
+    THETA, PHI = np.meshgrid(theta, phi)
+    im = ax1.pcolormesh(THETA, PHI, X0_matrix.T, shading='auto', cmap='viridis')
+    cbar = plt.colorbar(im, ax=ax1)
+    cbar.set_label('X/X0', fontsize=12)
+    ax1.set_title('Material Budget Distribution', fontsize=12, pad=10)
+    ax1.set_xlabel('Theta (angle)', fontsize=10)
+    ax1.set_ylabel('Phi (angle)', fontsize=10)
+    
+    # Draw phi projection
+    phi_projection = np.sum(X0_matrix, axis=0) / (range_theta[1]-range_theta[0])
+    ax2.plot(phi, phi_projection, 'b-', linewidth=2)
+    ax2.fill_between(phi, phi_projection, alpha=0.3)
+    ax2.set_title('Projection along Phi', fontsize=12, pad=10)
+    ax2.set_xlabel('Phi (degree)', fontsize=10)
+    ax2.set_ylabel('Average X/X0', fontsize=10)
+    ax2.grid(True, linestyle='--', alpha=0.7)
+    ax2.spines['top'].set_visible(False)
+    ax2.spines['right'].set_visible(False)
+    
+    # Draw theta projection
+    theta_projection = np.sum(X0_matrix, axis=1) / phi_bins
+    ax3.plot(theta, theta_projection, 'r-', linewidth=2)
+    ax3.fill_between(theta, theta_projection, alpha=0.3)
+    ax3.set_title('Projection along Theta', fontsize=12, pad=10)
+    ax3.set_xlabel('Theta (degree)', fontsize=10)
+    ax3.set_ylabel('Average X/X0', fontsize=10)
+    ax3.grid(True, linestyle='--', alpha=0.7)
+    ax3.spines['top'].set_visible(False)
+    ax3.spines['right'].set_visible(False)
+    
+    # Adjust main title position
+    plt.suptitle('Material Budget Analysis for {}'.format(detector), 
+                fontsize=16,  # Increase font size
+                y=0.95)  # Adjust title vertical position
+    
+    plt.tight_layout()
+    plt.savefig(output_file, dpi=600, bbox_inches='tight')
+    plt.show()
+
+def main():
+    parser = argparse.ArgumentParser(description='Process material scan data and generate plots.')
+    parser.add_argument('--input_file', help='Path to the input text file')
+    parser.add_argument('--output_file', default='material_budget.png', help='Path to the output png file')
+    parser.add_argument('--detector', default='all_world', help='Detector name')
+    parser.add_argument('--bins', type=int, default=100, help='theta bins is "bins"+1, phi bins is "bins"*2+1, (default: 100)')
+    parser.add_argument('--bias', type=int, default=0, help='Bias value for theta range (default: 0)')
+    args = parser.parse_args()
+
+    # all_world  Coil  Ecal  FTD  Hcal  Lumical  Muon  OTK  ParaffinEndcap  SIT  TPC  VXD
+    # bias detectors: BeamPipe onlyTracker
+    # Extract detector name from filename
+    x0_values = extract_x0_values(args.input_file)
+    plot_all_in_one(x0_values, args.output_file, args.bins, args.bias, args.detector)
+
+if __name__ == '__main__':
+    main() 
\ No newline at end of file

Analysis/MaterialScan/src/ScanMaterial.cpp

+#include <DD4hep/Detector.h>
+#include <DDRec/MaterialScan.h>
+#include <fstream>
+
+void ScanMaterial(string infile, string outfile="MaterialScanLogs.txt", int bins=100, double world_size=1000){
+    using namespace dd4hep;
+    using namespace dd4hep::rec;
+
+    if(infile.empty()){
+        cout << "Usage: root -l 'ScanMaterial.cpp(\"input_file.xml\", \"output_file.txt\", bins, world_size)' " << endl;
+        return;
+    }
+    
+    Detector& description = Detector::getInstance();
+    description.fromXML(infile.c_str());
+    MaterialScan scan(description);
+
+    std::ofstream clearFile(outfile.c_str(), std::ios::out | std::ios::trunc);
+    clearFile.close();
+    FILE* outFile = freopen(outfile.c_str(), "w", stdout);
+    for(int thetabin=0; thetabin<=bins; thetabin++){
+        double theta = thetabin * M_PI / bins;
+        double z = world_size*cos(theta);
+        double tranverse = world_size*sin(theta);
+        for(int phibin=-bins; phibin<=bins; phibin++){
+            double phi = phibin * M_PI / bins;
+            double x = tranverse*cos(phi);
+            double y = tranverse*sin(phi);
+            scan.print(0,0,0,x,y,z);
+        }
+    }
+    fclose(outFile);
+    return;
+}
\ No newline at end of file

Analysis/ReadDigi/src/ReadDigiAlg.cpp

@@ -52,6 +52,9 @@ StatusCode ReadDigiAlg::initialize()
   m_mctree->Branch("MCP_py", &MCP_py);
   m_mctree->Branch("MCP_pz", &MCP_pz);
   m_mctree->Branch("MCP_E", &MCP_E);
+  m_mctree->Branch("MCP_VTX_x", &MCP_VTX_x);
+  m_mctree->Branch("MCP_VTX_y", &MCP_VTX_y);
+  m_mctree->Branch("MCP_VTX_z", &MCP_VTX_z);
   m_mctree->Branch("MCP_endPoint_x", &MCP_endPoint_x);
   m_mctree->Branch("MCP_endPoint_y", &MCP_endPoint_y);
   m_mctree->Branch("MCP_endPoint_z", &MCP_endPoint_z);
@@ -176,6 +179,9 @@ StatusCode ReadDigiAlg::execute()
       MCP_py.push_back(m_MCp.getMomentum().y);
       MCP_pz.push_back(m_MCp.getMomentum().z);
       MCP_E.push_back(m_MCp.getEnergy());
+      MCP_VTX_x.push_back(m_MCp.getVertex().x);
+      MCP_VTX_y.push_back(m_MCp.getVertex().y);
+      MCP_VTX_z.push_back(m_MCp.getVertex().z);
       MCP_endPoint_x.push_back(m_MCp.getEndpoint().x); 
       MCP_endPoint_y.push_back(m_MCp.getEndpoint().y); 
       MCP_endPoint_z.push_back(m_MCp.getEndpoint().z); 
@@ -359,8 +365,9 @@ StatusCode ReadDigiAlg::execute()
     N_fullTrk = const_fullTrkCol->size();
     for(int i=0; i<N_fullTrk; i++){
       auto m_trk = const_fullTrkCol->at(i);
+//cout<<"In track #"<<i<<": track state size "<<m_trk.trackStates_size()<<", track hit size "<<m_trk.trackerHits_size()<<endl;
       if(m_trk.trackStates_size()==0) continue;
-      if(m_trk.trackerHits_size()==0) continue;
+      //if(m_trk.trackerHits_size()==0) continue;
       int NTrkHit = m_trk.trackerHits_size();
 
       for(int istat=0; istat<m_trk.trackStates_size(); istat++){
@@ -520,6 +527,9 @@ StatusCode ReadDigiAlg::Clear()
   MCP_py.clear();
   MCP_pz.clear();
   MCP_E.clear();
+  MCP_VTX_x.clear();
+  MCP_VTX_y.clear();
+  MCP_VTX_z.clear();
   MCP_endPoint_x.clear();
   MCP_endPoint_y.clear();
   MCP_endPoint_z.clear();

Analysis/ReadDigi/src/ReadDigiAlg.h

@@ -70,7 +70,7 @@ private :
 
   //MCParticle
   int N_MCP;
-  FloatVec MCP_px, MCP_py, MCP_pz, MCP_E, MCP_endPoint_x, MCP_endPoint_y, MCP_endPoint_z;
+  FloatVec MCP_px, MCP_py, MCP_pz, MCP_E, MCP_VTX_x, MCP_VTX_y, MCP_VTX_z, MCP_endPoint_x, MCP_endPoint_y, MCP_endPoint_z;
   IntVec MCP_pdgid, MCP_gStatus; 
 
   //Tracker 

Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Barrel_v04.cpp

@@ -539,7 +539,7 @@ static Ref_t create_detector(Detector &theDetector, xml_h element, SensitiveDete
             cell_x_offset -= total_cell_size_abnormal / 2.;
           }
         }
-        printf("layerID: %2d, x_length: %3.3lf, x_cell: %3.3lf, x_cell_abnormal: %3.3lf, x_dead: %3.3lf, z_length: %3.3lf,z_cell: %3.3lf, z_dead: %3.3lf,\n", layer_id,  2 * x_halflength,  n_x * total_cell_size,  nx_abnormal * total_cell_size_abnormal,  2 * x_halflength - n_x * total_cell_size - nx_abnormal * total_cell_size_abnormal,  2 * z_halfwidth,  n_z * total_cell_size,  2 * z_halfwidth - n_z * total_cell_size);
+        //printf("layerID: %2d, x_length: %3.3lf, x_cell: %3.3lf, x_cell_abnormal: %3.3lf, x_dead: %3.3lf, z_length: %3.3lf,z_cell: %3.3lf, z_dead: %3.3lf,\n", layer_id,  2 * x_halflength,  n_x * total_cell_size,  nx_abnormal * total_cell_size_abnormal,  2 * x_halflength - n_x * total_cell_size - nx_abnormal * total_cell_size_abnormal,  2 * z_halfwidth,  n_z * total_cell_size,  2 * z_halfwidth - n_z * total_cell_size);
         // Reset counters to measure "outside" quantitites
         nRadiationLengths = 0.;
         nInteractionLengths = 0.;

Detector/DetCRD/CMakeLists.txt

@@ -6,7 +6,11 @@
 gaudi_add_module(DetCRD
                  SOURCES src/Calorimeter/CRDEcal_v01.cpp
                          src/Calorimeter/LongCrystalBarBarrelCalorimeter32Polygon_v01.cpp
+                         src/Calorimeter/LongCrystalBarBarrelCalorimeter32Polygon_v02.cpp
                          src/Calorimeter/LongCrystalBarEndcapCalorimeter_v01.cpp
+                         src/Calorimeter/LongCrystalBarEndcapCalorimeter_v02.cpp
+                         src/Calorimeter/LongCrystalBarEndcapCalorimeter_v03.cpp
+                         src/Calorimeter/LongCrystalBarEndcapCalorimeter_v04.cpp
                          src/Calorimeter/CRDEcal_Short_v02.cpp
                          src/Calorimeter/CRDEcal_Endcap_Short_v01.cpp
                          src/Calorimeter/RotatedPolyhedraBarrelCalorimeter_v01_geo.cpp
@@ -14,19 +18,29 @@ gaudi_add_module(DetCRD
                          src/Calorimeter/Lumical_v01_geo.cpp
                          src/Other/Lumical_v01_geo_beampipe.cpp
                          src/Other/CRDBeamPipe_v01_geo.cpp
-                         src/Muon/Muon_Barrel_v01_01.cpp
+                         src/Muon/Muon_Barrel_v01_04.cpp
                          src/Muon/Muon_Endcap_v01_02.cpp
                          src/Tracker/SiTrackerSkewRing_v01_geo.cpp
+                         src/Tracker/ITK_EndCap_v01.cpp
                          src/Tracker/SiTrackerStitching_v01_geo.cpp
                          src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp
                          src/Tracker/SiTrackerStaggeredLadder_v02_geo.cpp
+                         src/Tracker/SiTrackerStaggeredLadder_v03_geo.cpp
                          src/Tracker/SiTrackerComposite_v01_geo.cpp
+                         src/Tracker/SiTrackerComposite_v02_geo.cpp
+                         src/Tracker/SiTrackerComposite_v03_geo.cpp
                          src/Tracker/TPC_Simple_o1_v01.cpp
                          src/Tracker/TPC_ModularEndcap_o1_v01.cpp
 			 src/Tracker/SiTracker_itkbarrel_v01_geo.cpp
 			 src/Tracker/SiTracker_itkbarrel_v02_geo.cpp
                          src/Tracker/SiTracker_otkbarrel_v01_geo.cpp
+                         src/Tracker/SiTracker_otkbarrel_v02_geo.cpp
+                         src/Tracker/SiTracker_otkendcap_v02_geo.cpp
                          src/Tracker/SiTracker_otkendcap_v01_geo.cpp
+                         src/Tracker/SiTracker_otkbarrel_v02_geo.cpp
+                         src/Tracker/SiTracker_otkendcap_v02_geo.cpp
+                         src/Other/ParaffinEndcap_v01.cpp
+                         src/Other/ConcreteWall_v01.cpp
 
 		 LINK ${DD4hep_COMPONENT_LIBRARIES}
                       DetIdentifier
@@ -74,6 +88,12 @@ add_test(
   WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
   )
 
+add_test(
+  NAME Test_TDR_o1_v01_PFA
+  COMMAND gaudirun.py Detector/DetCRD/scripts/TDR_o1_v01/rec.py
+  WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+  )
+
 foreach(detoption TDR_o1_v01)
   add_test(
     NAME Test_${detoption}_Geo

Detector/DetCRD/compact/CRD_common_v01/ConcreteWall_v01_01.xml

+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd>
+ 
+<define>
+    <constant name="concrete_position_z" value="592.5*cm"/>
+    <constant name="concrete_hole_rmin" value="65*cm"/>
+    <constant name="concrete_wall_size" value="520*cm"/>
+    <constant name="concrete_wall_hz" value="50*cm"/>
+</define>
+
+<detectors>
+    <detector name="ConcreteWall" type="ConcreteWall_v01" vis="VisibleBrown" material="Concrete">
+      <!-- Use cm as unit if you want to use Pandora for reconstruction -->
+      <position x="0"  y="0"  z="concrete_position_z"/>
+      <dimensions rmax="concrete_wall_size" rmin="concrete_hole_rmin" dz="concrete_wall_hz"/>
+    </detector>
+</detectors>  
+
+</lccdd>

Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Barrel_v02_01.xml

+<?xml version="1.0" encoding="UTF-8"?> 
+<lccdd> 
+  <define>
+    <constant name="ecalbarrel_inner_radius"        value="Ecal_barrel_inner_radius"/>
+    <constant name="ecalbarrel_outer_radius"        value="Ecal_barrel_outer_radius"/>
+    <constant name="ecalbarrel_zlength"             value="Ecal_barrel_half_length*2"/>
+
+    <constant name="ecalbarrel_layer"               value="18" />
+    <constant name="ecalbarrel_phimodule_number"    value="32" /> 
+    <constant name="ecalbarrel_Zmodule_number"      value="15"/>  
+    <constant name="ecalbarrel_module_ratation"     value="12.*degree" />
+    <constant name="ecalbarrel_width_crystal"    value="15.2*mm"/>
+
+    <constant name="ecalbarrel_crystal_wrapping"    value="0.1*mm"/>
+    <constant name="ecalbarrel_length_photoelectronic"    value="0.7*mm"/> 
+    <constant name="ecalbarrel_width_photoelectronic"    value="3*mm"/> 
+
+    <constant name="ecalbarrel_length_carbon"    value="5.*mm"/> 
+    <constant name="ecalbarrel_length_cable"    value="0.*mm"/> 
+    <constant name="ecalbarrel_length_cooling"    value="1.*mm"/> 
+    <constant name="ecalbarrel_length_pcb"    value="1.2*mm"/> 
+    <constant name="ecalbarrel_length_asic"    value="1.*mm"/>
+    <constant name="ecalbarrel_length_back"    value="10.*mm"/>
+
+  </define> 
+
+  <regions>
+    <region name="EcalBarrelRegion">
+    </region>
+  </regions>
+
+  <detectors>
+    <detector id="DetID_ECAL" name="EcalBarrel" type="LongCrystalBarBarrelCalorimeter32Polygon_v01" readout="EcalBarrelCollection" vis="Invisible" sensitive="true" region="EcalBarrelRegion">
+      <!-- Use cm as unit if you want to use Pandora for reconstruction -->
+    </detector>
+  </detectors>
+  
+  <readouts>
+    <readout name="EcalBarrelCollection">
+      <segmentation type="NoSegmentation"/>
+      <!--segmentation type="CartesianGridXYZ"
+                    grid_size_x="1*cm"
+                    grid_size_y="1*cm"
+                    grid_size_z="1*cm"/-->
+      <id>system:5,module:5,stave:4,dlayer:5,slayer:6,bar:15</id>
+    </readout>
+  </readouts>
+
+</lccdd>