Analysis/DumpEvent/src/DumpTrackAlg.cpp

+#include "DumpTrackAlg.h"
+
+#if __has_include("edm4hep/EDM4hepVersion.h")
+#include "edm4hep/EDM4hepVersion.h"
+#else
+// Copy the necessary parts from  the header above to make whatever we need to work here
+#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
+// v00-09 is the last version without the capitalization change of the track vector members
+#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
+#endif
+
+
+#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 <math.h>
+
+DECLARE_COMPONENT( DumpTrackAlg )
+
+//------------------------------------------------------------------------------
+DumpTrackAlg::DumpTrackAlg( const std::string& name, ISvcLocator* pSvcLocator )
+    : Algorithm( name, pSvcLocator ) {
+  declareProperty("MCParticleCollection", _inMCColHdl, "Handle of the Input MCParticle collection");
+  declareProperty("TrackCollection", _inTrackColHdl, "Handle of the Input Track collection from CEPCSW");
+
+  m_thisName = name;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpTrackAlg::initialize(){
+  info() << "Booking Ntuple" << endmsg;
+
+  NTuplePtr nt1(ntupleSvc(), "MyTuples/Track"+m_thisName);
+  if ( !nt1 ) {
+    m_tuple = ntupleSvc()->book("MyTuples/Track"+m_thisName,CLID_ColumnWiseTuple,"Tracking result");
+    if ( 0 != m_tuple ) {
+      m_tuple->addItem ("ntrk",      m_nTracks, 0, 500 ).ignore();
+      m_tuple->addIndexedItem ("x",         m_nTracks, m_x ).ignore();
+      m_tuple->addIndexedItem ("y",         m_nTracks, m_y ).ignore();
+      m_tuple->addIndexedItem ("z",         m_nTracks, m_z ).ignore();
+      m_tuple->addIndexedItem ("px",        m_nTracks, m_px ).ignore();
+      m_tuple->addIndexedItem ("py",        m_nTracks, m_py ).ignore();
+      m_tuple->addIndexedItem ("pz",        m_nTracks, m_pz ).ignore();
+      m_tuple->addIndexedItem ("d0",        m_nTracks, m_d0 ).ignore();
+      m_tuple->addIndexedItem ("phi0",      m_nTracks, m_phi0 ).ignore();
+      m_tuple->addIndexedItem ("omega",     m_nTracks, m_omega ).ignore();
+      m_tuple->addIndexedItem ("z0",        m_nTracks, m_z0 ).ignore();
+      m_tuple->addIndexedItem ("tanLambda", m_nTracks, m_tanLambda ).ignore();
+      m_tuple->addIndexedItem ("sigma_d0",        m_nTracks, m_sigma_d0 ).ignore();
+      m_tuple->addIndexedItem ("sigma_phi0",      m_nTracks, m_sigma_phi0 ).ignore();
+      m_tuple->addIndexedItem ("sigma_omega",     m_nTracks, m_sigma_omega ).ignore();
+      m_tuple->addIndexedItem ("sigma_z0",        m_nTracks, m_sigma_z0 ).ignore();
+      m_tuple->addIndexedItem ("sigma_tanLambda", m_nTracks, m_sigma_tanLambda ).ignore();
+      m_tuple->addIndexedItem ("nvxd", m_nTracks, m_nHitsVXD ).ignore();
+      m_tuple->addIndexedItem ("nftd", m_nTracks, m_nHitsFTD ).ignore();
+      m_tuple->addIndexedItem ("nsit", m_nTracks, m_nHitsSIT ).ignore();
+      m_tuple->addIndexedItem ("ngas", m_nTracks, m_nHitsGAS ).ignore();
+      m_tuple->addIndexedItem ("nset", m_nTracks, m_nHitsSET ).ignore();
+    }
+    else { // did not manage to book the N tuple....
+      fatal() << "Cannot book MyTuples/Track"+m_thisName <<endmsg; 
+      return StatusCode::FAILURE;
+    }
+  }
+  else{
+    m_tuple = nt1;
+  }
+
+  auto geomSvc = service<IGeomSvc>("GeomSvc");
+  if(geomSvc){
+    const dd4hep::Direction& field = geomSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
+    m_field = field.z()/dd4hep::tesla;
+    info() << "Magnetic field will obtain from GeomSvc = " << m_field << " tesla" << endmsg;
+  }
+  else{
+    info() << "Failed to find GeomSvc ..." << endmsg;
+    info() << "Magnetic field will use what input through python option for this algorithm namse as Field, now " << m_field << " tesla" << endmsg;
+  }
+
+  _nEvt = 0;
+  return StatusCode::SUCCESS;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpTrackAlg::execute(){
+  const edm4hep::TrackCollection* trackCols = nullptr;
+  try {
+    trackCols = _inTrackColHdl.get();
+  }
+  catch ( GaudiException &e ) {
+    debug() << "Collection " << _inTrackColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+  }
+
+  if(trackCols){
+    m_nTracks = 0;
+    for(auto track : *trackCols){
+      // since possible more than one location=1 TrackState (not deleted in reconstruction), always use last one 
+      for(std::vector<edm4hep::TrackState>::const_iterator it=track.trackStates_end()-1; it!=track.trackStates_begin()-1; it--){
+	edm4hep::TrackState trackState = *it;
+	if(trackState.location!=1)continue;
+	m_d0[m_nTracks]              = trackState.D0;
+	m_phi0[m_nTracks]            = trackState.phi;
+	m_omega[m_nTracks]           = trackState.omega;
+	m_z0[m_nTracks]              = trackState.Z0;
+	m_tanLambda[m_nTracks]       = trackState.tanLambda;
+	m_sigma_d0[m_nTracks]        = std::sqrt(trackState.covMatrix[0]);
+	m_sigma_phi0[m_nTracks]      = std::sqrt(trackState.covMatrix[2]);
+	m_sigma_omega[m_nTracks]     = std::sqrt(trackState.covMatrix[5]);
+	m_sigma_z0[m_nTracks]        = std::sqrt(trackState.covMatrix[9]);
+        m_sigma_tanLambda[m_nTracks] = std::sqrt(trackState.covMatrix[14]);
+	HelixClass helix_fit;
+        helix_fit.Initialize_Canonical(trackState.phi,trackState.D0,trackState.Z0,trackState.omega,trackState.tanLambda,m_field);
+	m_x[m_nTracks] = helix_fit.getReferencePoint()[0];
+	m_y[m_nTracks] = helix_fit.getReferencePoint()[1];
+	m_z[m_nTracks] = helix_fit.getReferencePoint()[2];
+        m_px[m_nTracks] = helix_fit.getMomentum()[0];
+        m_py[m_nTracks] = helix_fit.getMomentum()[1];
+        m_pz[m_nTracks] = helix_fit.getMomentum()[2];
+	//info() << "size = " << track.subDetectorHitNumbers_size() << endmsg;
+	//for(int ii=0;ii<track.subDetectorHitNumbers_size();ii++){
+	//  std::cout << track.getSubDetectorHitNumbers(ii) << " ";
+	//}
+	//std::cout << std::endl;
+#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0)
+	if(track.subdetectorHitNumbers_size()>=5){
+	  m_nHitsVXD[m_nTracks] = track.getSubdetectorHitNumbers(0);
+	  m_nHitsFTD[m_nTracks] = track.getSubdetectorHitNumbers(1);
+	  m_nHitsSIT[m_nTracks] = track.getSubdetectorHitNumbers(2);
+	  m_nHitsGAS[m_nTracks] = track.getSubdetectorHitNumbers(3);
+	  m_nHitsSET[m_nTracks] = track.getSubdetectorHitNumbers(4);
+	}
+#else
+	if(track.subDetectorHitNumbers_size()>=5){
+	  m_nHitsVXD[m_nTracks] = track.getSubDetectorHitNumbers(0);
+	  m_nHitsFTD[m_nTracks] = track.getSubDetectorHitNumbers(1);
+	  m_nHitsSIT[m_nTracks] = track.getSubDetectorHitNumbers(2);
+	  m_nHitsGAS[m_nTracks] = track.getSubDetectorHitNumbers(3);
+	  m_nHitsSET[m_nTracks] = track.getSubDetectorHitNumbers(4);
+	}
+#endif
+	else{
+	  m_nHitsVXD[m_nTracks] = 0;
+          m_nHitsSIT[m_nTracks] = 0;
+          m_nHitsSET[m_nTracks] = 0;
+          m_nHitsFTD[m_nTracks] = 0;
+          m_nHitsGAS[m_nTracks] = 0;
+	}
+	m_nTracks++;
+	break;
+      }
+    }
+  }
+  m_tuple->write();
+  _nEvt++;
+  return StatusCode::SUCCESS;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpTrackAlg::finalize(){
+  debug() << "Finalizing..." << endmsg;
+
+  return StatusCode::SUCCESS;
+}

Analysis/DumpEvent/src/DumpTrackAlg.h

+#ifndef DumpTrackAlg_h
+#define DumpTrackAlg_h 1
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/Algorithm.h"
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/TrackCollection.h"
+
+#include "GaudiKernel/NTuple.h"
+
+class DumpTrackAlg : public Algorithm {
+ public:
+  // Constructor of this form must be provided
+  DumpTrackAlg( 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::MCParticleCollection> _inMCColHdl{"MCParticle", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::TrackCollection> _inTrackColHdl{"SiTracks", Gaudi::DataHandle::Reader, this};
+
+  Gaudi::Property<double> m_field{this, "Field", 3.0};
+
+  NTuple::Tuple*       m_tuple;
+  NTuple::Item<long>   m_nTracks;
+  NTuple::Array<float> m_x;
+  NTuple::Array<float> m_y;
+  NTuple::Array<float> m_z;
+  NTuple::Array<float> m_px;
+  NTuple::Array<float> m_py;
+  NTuple::Array<float> m_pz;
+  NTuple::Array<float> m_d0;
+  NTuple::Array<float> m_phi0;
+  NTuple::Array<float> m_omega;
+  NTuple::Array<float> m_z0;
+  NTuple::Array<float> m_tanLambda;
+  NTuple::Array<float> m_sigma_d0;
+  NTuple::Array<float> m_sigma_phi0;
+  NTuple::Array<float> m_sigma_omega;
+  NTuple::Array<float> m_sigma_z0;
+  NTuple::Array<float> m_sigma_tanLambda;
+  NTuple::Array<int>   m_nHitsVXD;
+  NTuple::Array<int>   m_nHitsFTD;
+  NTuple::Array<int>   m_nHitsSIT;
+  NTuple::Array<int>   m_nHitsGAS;
+  NTuple::Array<int>   m_nHitsSET;
+  
+  int _nEvt;
+  std::string m_thisName;
+};
+
+#endif

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/JetClustering/CMakeLists.txt

+# Set the CMake module path if necessary
+
+# Modules
+gaudi_add_module(JetClustering
+                 SOURCES src/JetClustering.cpp
+                 LINK GearSvc
+                      Gaudi::GaudiKernel
+                      k4FWCore::k4FWCore
+                      DetInterface
+                      DataHelperLib
+                      ${GEAR_LIBRARIES} 
+                      ${GSL_LIBRARIES} 
+                      ${ROOT_LIBRARIES}
+                      ${FASTJET_LIBRARY}
+                      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)

Analysis/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>
+
+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("ymerge", ymerge, "ymerge[6]/D");
+  //_tree->Branch("px", &px);
+  //_tree->Branch("py", &py);
+  //_tree->Branch("pz", &pz);
+  //_tree->Branch("E", &E);
+  
+  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
+  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 ) {
+      debug() << "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() ) );
+      _particle_index++;
+      debug() << "Particle " << _particle_index << " px: "<<pfo.getMomentum()[0]<<" py: "<<pfo.getMomentum()[1]<<" pz: "<<pfo.getMomentum()[2]<<" E: "<<pfo.getEnergy()<<endmsg;
+    }
+  }
+
+  // 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
+  debug() << "                     px           py           pz           E           costheta  phi " << endmsg;
+  for (unsigned i = 0; i < jets.size(); i++) {
+    
+    debug() << "   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 == 0 )constituents_E1tot += constituents[j].E();
+      
+      if ( i == 1 )constituents_E2tot += constituents[j].E();
+      
+      debug() << "    constituent  " << j << " " << constituents[j].px() << " " << constituents[j].py() << " " << constituents[j].pz() << " " << constituents[j].E() << " " << constituents[j].cos_theta() << " " << constituents[j].phi() << endmsg;
+    
+    }
+  
+  }
+  
+  debug() << "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);
+    debug() << " -log10(y" << i << i+1 << ") = " << -log10(_ymin[i-1]) << endmsg;
+  }
+  mass = (jets[0] + jets[1]).m();
+
+  info() << "Total mass of all particles: " << " mass of jet1+jet2: " << mass << endmsg;
+
+  // 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;
+}

Analysis/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;
+
+  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;
+
+  }
+
+};
+#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/CMakeLists.txt

+
+gaudi_add_module(ReadDigi
+                 SOURCES src/ReadDigiAlg.cpp
+                         src/HelixClassD.cc
+                 LINK k4FWCore::k4FWCore 
+                      Gaudi::GaudiKernel 
+                      Gaudi::GaudiAlgLib 
+                      ${CLHEP_LIBRARIES} 
+                      ${GSL_LIBRARIES}
+                      ${LCIO_LIBRARIES}
+                      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+                      ${ROOT_LIBRARIES}
+)
+
+target_include_directories(ReadDigi
+  PUBLIC ${LCIO_INCLUDE_DIRS})
+
+install(TARGETS ReadDigi
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)

Analysis/ReadDigi/src/HelixClassD.cc

+#include "HelixClassD.hh"
+#include <math.h>
+#include <stdlib.h>
+#include <iostream>
+//#include "ced_cli.h"
+
+HelixClassD::HelixClassD() {
+  _const_2pi = 2.0*M_PI;
+  _const_pi2 = 0.5*M_PI;
+  _FCT = 2.99792458E-4;
+}
+
+HelixClassD::~HelixClassD() {}
+
+void HelixClassD::Initialize_VP(float * pos, float * mom, float q, float B) {
+  _referencePoint[0] = pos[0];
+  _referencePoint[1] = pos[1];
+  _referencePoint[2] = pos[2];
+  _momentum[0] = mom[0];
+  _momentum[1] = mom[1];
+  _momentum[2] = mom[2];
+  _charge = q;
+  _bField = B;
+  _pxy = sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+  _radius = _pxy / (_FCT*B);
+  _omega = q/_radius;
+  _tanLambda = mom[2]/_pxy;
+  _phiMomRefPoint = atan2(mom[1],mom[0]);
+  _xCentre = pos[0] + _radius*cos(_phiMomRefPoint-_const_pi2*q);
+  _yCentre = pos[1] + _radius*sin(_phiMomRefPoint-_const_pi2*q);
+  _phiRefPoint = atan2(pos[1]-_yCentre,pos[0]-_xCentre);
+  _phiAtPCA = atan2(-_yCentre,-_xCentre);
+  _phi0 = -_const_pi2*q + _phiAtPCA;
+  while (_phi0<0) _phi0+=_const_2pi;
+  while (_phi0>=_const_2pi) _phi0-=_const_2pi;
+  _xAtPCA = _xCentre + _radius*cos(_phiAtPCA);
+  _yAtPCA = _yCentre + _radius*sin(_phiAtPCA);
+  //    _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0);
+  double pxy = double(_pxy);
+  double radius = pxy/double(_FCT*B);
+  double xCentre = double(pos[0]) + radius*double(cos(_phiMomRefPoint-_const_pi2*q));
+  double yCentre = double(pos[1]) + radius*double(sin(_phiMomRefPoint-_const_pi2*q));
+
+  double d0;
+
+  if (q>0) {
+    d0 = double(q)*radius - double(sqrt(xCentre*xCentre+yCentre*yCentre));
+  }
+  else {
+    d0 = double(q)*radius + double(sqrt(xCentre*xCentre+yCentre*yCentre));
+  }
+
+  _d0 = float(d0);
+
+  //     if (fabs(_d0)>0.001 ) {
+  //       std::cout << "New helix : " << std::endl;
+  //       std::cout << " Position : " << pos[0] 
+  // 		<< " " << pos[1]
+  // 		<< " " << pos[2] << std::endl;
+  //       std::cout << " Radius = " << _radius << std::endl;
+  //       std::cout << " RC = " << sqrt(_xCentre*_xCentre+_yCentre*_yCentre) << std::endl;  
+  //       std::cout << " D0 = " << _d0 << std::endl;
+  //     }
+
+  _pxAtPCA = _pxy*cos(_phi0);
+  _pyAtPCA = _pxy*sin(_phi0);
+  float deltaPhi = _phiRefPoint - _phiAtPCA;    
+  float xCircles = -pos[2]*q/(_radius*_tanLambda) - deltaPhi;
+  xCircles = xCircles/_const_2pi;
+  int nCircles;
+  int n1,n2;
+
+  if (xCircles >= 0.) {
+    n1 = int(xCircles);
+    n2 = n1 + 1;
+  }
+  else {
+    n1 = int(xCircles) - 1;
+    n2 = n1 + 1;
+  }
+
+  if (fabs(n1-xCircles) < fabs(n2-xCircles)) {
+    nCircles = n1;
+  }
+  else {
+    nCircles = n2;
+  }
+  _z0 = pos[2] + _radius*_tanLambda*q*(deltaPhi + _const_2pi*nCircles);
+
+}
+
+void HelixClassD::Initialize_Canonical(float phi0, float d0, float z0, 
+    float omega, float tanLambda, float B) {
+  _omega = omega;
+  _d0 = d0;
+  _phi0 = phi0;
+  _z0 = z0;
+  _tanLambda = tanLambda;
+  _charge = omega/fabs(omega);
+  _radius = 1./fabs(omega);
+  _xAtPCA = -_d0*sin(_phi0);
+  _yAtPCA = _d0*cos(_phi0);    
+  _referencePoint[0] = _xAtPCA;
+  _referencePoint[1] = _yAtPCA;
+  _referencePoint[2] = _z0;
+  _pxy = _FCT*B*_radius;
+  _momentum[0] = _pxy*cos(_phi0);
+  _momentum[1] = _pxy*sin(_phi0);
+  _momentum[2] = _tanLambda * _pxy;    
+  _pxAtPCA = _momentum[0];
+  _pyAtPCA = _momentum[1];
+  _phiMomRefPoint = atan2(_momentum[1],_momentum[0]);
+  _xCentre = _referencePoint[0] + 
+    _radius*cos(_phi0-_const_pi2*_charge);
+  _yCentre = _referencePoint[1] + 
+    _radius*sin(_phi0-_const_pi2*_charge);
+  _phiAtPCA = atan2(-_yCentre,-_xCentre);
+  _phiRefPoint =  _phiAtPCA ;
+  _bField = B;
+}
+
+
+void HelixClassD::Initialize_BZ(float xCentre, float yCentre, float radius, 
+    float bZ, float phi0, float B, float signPz,
+    float zBegin) {
+
+  _radius = radius;
+  _pxy = _FCT*B*_radius;
+  _charge = -(bZ*signPz)/fabs(bZ*signPz);
+  _momentum[2] = -_charge*_pxy/(bZ*_radius);
+  _xCentre = xCentre;
+  _yCentre = yCentre;
+  _omega = _charge/radius;
+  _phiAtPCA = atan2(-_yCentre,-_xCentre);
+  _phi0 = -_const_pi2*_charge + _phiAtPCA;
+  while (_phi0<0) _phi0+=_const_2pi;
+  while (_phi0>=_const_2pi) _phi0-=_const_2pi;
+  _xAtPCA = _xCentre + _radius*cos(_phiAtPCA);
+  _yAtPCA = _yCentre + _radius*sin(_phiAtPCA);
+  _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0);
+  _pxAtPCA = _pxy*cos(_phi0);
+  _pyAtPCA = _pxy*sin(_phi0);
+  _referencePoint[2] = zBegin;
+  _referencePoint[0] = xCentre + radius*cos(bZ*zBegin+phi0);
+  _referencePoint[1] = yCentre + radius*sin(bZ*zBegin+phi0);
+  _phiRefPoint = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre);
+  _phiMomRefPoint =  -_const_pi2*_charge + _phiRefPoint;
+  _tanLambda = _momentum[2]/_pxy;
+  _momentum[0] = _pxy*cos(_phiMomRefPoint);
+  _momentum[1] = _pxy*sin(_phiMomRefPoint);
+
+  float deltaPhi = _phiRefPoint - _phiAtPCA;    
+  float xCircles = bZ*_referencePoint[2] - deltaPhi;
+  xCircles = xCircles/_const_2pi;
+  int nCircles;
+  int n1,n2;
+
+  if (xCircles >= 0.) {
+    n1 = int(xCircles);
+    n2 = n1 + 1;
+  }
+  else {
+    n1 = int(xCircles) - 1;
+    n2 = n1 + 1;
+  }
+
+  if (fabs(n1-xCircles) < fabs(n2-xCircles)) {
+    nCircles = n1;
+  }
+  else {
+    nCircles = n2;
+  }  
+  _z0 = _referencePoint[2] - (deltaPhi + _const_2pi*nCircles)/bZ;  
+  _bField = B;
+
+}
+
+const float * HelixClassD::getMomentum() {
+  return _momentum;
+}
+const float * HelixClassD::getReferencePoint() {
+  return _referencePoint;
+}
+float HelixClassD::getPhi0() {
+  if (_phi0<0.0)
+    _phi0 += 2*M_PI;
+  return _phi0;
+}
+float HelixClassD::getD0() {
+  return _d0;
+}
+float HelixClassD::getZ0() {
+  return _z0;
+}
+float HelixClassD::getOmega() {
+  return _omega;
+}
+float HelixClassD::getTanLambda() {
+  return _tanLambda;
+}
+float HelixClassD::getPXY() {
+  return _pxy;
+}
+float HelixClassD::getXC() {
+  return _xCentre;
+}
+
+float HelixClassD::getYC() {
+  return _yCentre;
+}
+
+float HelixClassD::getRadius() {
+  return _radius;
+}
+
+float HelixClassD::getBz() {
+  return _bZ;
+}
+
+float HelixClassD::getPhiZ() {
+  return _phiZ;
+}
+
+float HelixClassD::getCharge() {
+  return _charge;
+}
+
+float HelixClassD::getPointInXY(float x0, float y0, float ax, float ay, 
+    float * ref , float * point) {
+
+  float time;
+
+  float AA = sqrt(ax*ax+ay*ay);
+
+
+  if (AA <= 0) {
+    time = -1.0e+20; 
+    return time;
+  }
+
+
+  float BB = ax*(x0-_xCentre) + ay*(y0-_yCentre);
+  BB = BB / AA;
+
+  float CC = (x0-_xCentre)*(x0-_xCentre) 
+    + (y0-_yCentre)*(y0-_yCentre) - _radius*_radius;
+
+  CC = CC / AA;
+
+  float DET = BB*BB - CC;
+  float tt1 = 0.;
+  float tt2 = 0.;
+  float xx1,xx2,yy1,yy2; 
+
+
+  if (DET < 0 ) {
+    time = 1.0e+10;
+    point[0]=0.0;
+    point[1]=0.0;
+    point[2]=0.0;
+    return time;
+  }
+
+
+  tt1 = - BB + sqrt(DET);
+  tt2 = - BB - sqrt(DET);
+
+  xx1 = x0 + tt1*ax;
+  yy1 = y0 + tt1*ay;
+  xx2 = x0 + tt2*ax;
+  yy2 = y0 + tt2*ay;
+
+  float phi1 = atan2(yy1-_yCentre,xx1-_xCentre);
+  float phi2 = atan2(yy2-_yCentre,xx2-_xCentre);
+  float phi0 = atan2(ref[1]-_yCentre,ref[0]-_xCentre);
+
+  float dphi1 = phi1 - phi0;
+  float dphi2 = phi2 - phi0;
+
+  if (dphi1 < 0 && _charge < 0) {
+    dphi1 = dphi1 + _const_2pi;
+  }
+  else if (dphi1 > 0 && _charge > 0) { 
+    dphi1 = dphi1 - _const_2pi;
+  }
+
+  if (dphi2 < 0 && _charge < 0) {
+    dphi2 = dphi2 + _const_2pi;
+  }
+  else if (dphi2 > 0 && _charge > 0) { 
+    dphi2 = dphi2 - _const_2pi;
+  }
+
+  // Times
+  tt1 = -_charge*dphi1*_radius/_pxy;
+  tt2 = -_charge*dphi2*_radius/_pxy;
+
+  if (tt1 < 0. )
+    std::cout << "WARNING " << tt1 << std::endl;
+  if (tt2 < 0. )
+    std::cout << "WARNING " << tt2 << std::endl;
+
+
+  if (tt1 < tt2) {
+    point[0] = xx1;
+    point[1] = yy1;
+    time = tt1;
+  }
+  else {
+    point[0] = xx2;
+    point[1] = yy2;
+    time = tt2;
+  }
+
+  point[2] = ref[2]+time*_momentum[2];
+
+
+
+  return time;
+
+}
+
+
+float HelixClassD::getPointOnCircle(float Radius, float * ref, float * point) {
+
+  float distCenterToIP = sqrt(_xCentre*_xCentre + _yCentre*_yCentre);
+
+  point[0] = 0.0;
+  point[1] = 0.0;
+  point[2] = 0.0;
+
+  if ((distCenterToIP+_radius)<Radius) {
+    float xx = -1.0e+20;
+    return xx;
+  }
+
+  if ((_radius+Radius)<distCenterToIP) {
+    float xx = -1.0e+20;
+    return xx;
+  }
+
+  float phiCentre = atan2(_yCentre,_xCentre);
+  float phiStar   = Radius*Radius + distCenterToIP*distCenterToIP 
+    - _radius*_radius;
+
+  phiStar = 0.5*phiStar/fmax(1.0e-20,Radius*distCenterToIP);
+
+  if (phiStar > 1.0) 
+    phiStar = 0.9999999;
+
+  if (phiStar < -1.0)
+    phiStar = -0.9999999;
+
+  phiStar = acos(phiStar);
+
+  float tt1,tt2,time;
+
+  float xx1 = Radius*cos(phiCentre+phiStar);
+  float yy1 = Radius*sin(phiCentre+phiStar);
+
+  float xx2 = Radius*cos(phiCentre-phiStar);
+  float yy2 = Radius*sin(phiCentre-phiStar);
+
+
+  float phi1 = atan2(yy1-_yCentre,xx1-_xCentre);
+  float phi2 = atan2(yy2-_yCentre,xx2-_xCentre);
+  float phi0 = atan2(ref[1]-_yCentre,ref[0]-_xCentre);
+
+  float dphi1 = phi1 - phi0;
+  float dphi2 = phi2 - phi0;
+
+  if (dphi1 < 0 && _charge < 0) {
+    dphi1 = dphi1 + _const_2pi;
+  }
+  else if (dphi1 > 0 && _charge > 0) { 
+    dphi1 = dphi1 - _const_2pi;
+  }
+
+  if (dphi2 < 0 && _charge < 0) {
+    dphi2 = dphi2 + _const_2pi;
+  }
+  else if (dphi2 > 0 && _charge > 0) { 
+    dphi2 = dphi2 - _const_2pi;
+  }
+
+  // Times
+  tt1 = -_charge*dphi1*_radius/_pxy;
+  tt2 = -_charge*dphi2*_radius/_pxy;
+
+  if (tt1 < 0. )
+    std::cout << "WARNING " << tt1 << std::endl;
+  if (tt2 < 0. )
+    std::cout << "WARNING " << tt2 << std::endl;
+
+
+  float time2;
+  if (tt1 < tt2) {
+    point[0] = xx1;
+    point[1] = yy1;
+    point[3] = xx2;
+    point[4] = yy2;
+    time = tt1;
+    time2 = tt2;
+  }
+  else {
+    point[0] = xx2;
+    point[1] = yy2;
+    point[3] = xx1;
+    point[4] = yy1;
+    time = tt2;
+    time2 = tt1;
+  }
+
+  point[2] = ref[2]+time*_momentum[2];
+  point[5] = ref[2]+time2*_momentum[2];
+
+
+  return time;
+
+}
+
+
+float HelixClassD::getPointInZ(float zLine, float * ref, float * point) {
+
+  float time = zLine - ref[2];
+
+  if (_momentum[2] == 0.) {
+    time = -1.0e+20;
+    point[0] = 0.;
+    point[1] = 0.;
+    point[2] = 0.;
+    return time;
+  }
+
+  time = time/_momentum[2];
+
+  float phi0 = atan2(ref[1] - _yCentre , ref[0] - _xCentre);
+  float phi = phi0 - _charge*_pxy*time/_radius;
+  float xx = _xCentre + _radius*cos(phi);
+  float yy = _yCentre + _radius*sin(phi);
+
+  point[0] = xx;
+  point[1] = yy;
+  point[2] = zLine;
+
+  return time;
+
+
+}
+
+int HelixClassD::getNCircle(float * xPoint) {
+
+  int nCircles = 0;
+
+  float phi = atan2(xPoint[1]-_yCentre,xPoint[0]-_xCentre);
+  float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre);
+
+  if (fabs(_tanLambda*_radius)>1.0e-20) {
+    float xCircles = phi0 - phi -_charge*(xPoint[2]-_referencePoint[2])/(_tanLambda*_radius);
+    xCircles = xCircles/_const_2pi;
+    int n1,n2;
+
+    if (xCircles >= 0.) {
+      n1 = int(xCircles);
+      n2 = n1 + 1;
+    }
+    else {
+      n1 = int(xCircles) - 1;
+      n2 = n1 + 1;
+    }
+
+    if (fabs(n1-xCircles) < fabs(n2-xCircles)) {
+      nCircles = n1;
+    }
+    else {
+      nCircles = n2;
+    }
+
+  }
+  return nCircles;
+
+}
+
+float HelixClassD::getDistanceToPoint(float * xPoint, float * Distance) {
+
+  float zOnHelix;
+  float phi = atan2(xPoint[1]-_yCentre,xPoint[0]-_xCentre);
+  float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre);
+  //calculate distance to XYprojected centre of Helix, comparing this with distance to radius around centre gives DistXY
+  float DistXY = (_xCentre-xPoint[0])*(_xCentre-xPoint[0]) + (_yCentre-xPoint[1])*(_yCentre-xPoint[1]);
+  DistXY = sqrt(DistXY);
+  DistXY = fabs(DistXY - _radius);
+
+  int nCircles = 0;
+
+  if (fabs(_tanLambda*_radius)>1.0e-20) {
+    float xCircles = phi0 - phi -_charge*(xPoint[2]-_referencePoint[2])/(_tanLambda*_radius);
+    xCircles = xCircles/_const_2pi;
+    int n1,n2;
+
+    if (xCircles >= 0.) {
+      n1 = int(xCircles);
+      n2 = n1 + 1;
+    }
+    else {
+      n1 = int(xCircles) - 1;
+      n2 = n1 + 1;
+    }
+
+    if (fabs(n1-xCircles) < fabs(n2-xCircles)) {
+      nCircles = n1;
+    }
+    else {
+      nCircles = n2;
+    }
+
+  }
+
+  float DPhi = _const_2pi*((float)nCircles) + phi - phi0;
+
+  zOnHelix = _referencePoint[2] - _charge*_radius*_tanLambda*DPhi;
+
+  float DistZ = fabs(zOnHelix - xPoint[2]);
+  float Time;
+
+  if (fabs(_momentum[2]) > 1.0e-20) {
+    Time = (zOnHelix - _referencePoint[2])/_momentum[2];
+  }
+  else {
+    Time = _charge*_radius*DPhi/_pxy;
+  }
+
+  Distance[0] = DistXY;
+  Distance[1] = DistZ;
+  Distance[2] = sqrt(DistXY*DistXY+DistZ*DistZ);
+
+  return Time;
+
+
+}
+
+//When we are not interested in the exact distance, we can check if we are
+//already far enough away in XY, before we start calculating in Z as the
+//distance will only increase
+float HelixClassD::getDistanceToPoint(const std::vector<float>& xPoint, float distCut) {
+  //calculate distance to XYprojected centre of Helix, comparing this with distance to radius around centre gives DistXY
+  float tempx = xPoint[0]-_xCentre;
+  float tempy = xPoint[1]-_yCentre;
+  float tempsq = sqrt(tempx*tempx + tempy*tempy);
+  float tempdf = tempsq - _radius;
+  float DistXY = fabs( tempdf );
+  //If this is bigger than distCut, we dont have to know how much bigger this is
+  if( DistXY > distCut) {
+    return DistXY;
+  }
+
+  int nCircles = 0;
+  float phi = atan2(tempy,tempx);
+  float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre);
+  float phidiff = phi0-phi;
+  float  tempz = xPoint[2] - _referencePoint[2];//Yes referencePoint
+  float tanradius = _tanLambda*_radius;
+  if (fabs(tanradius)>1.0e-20) {
+    float xCircles = (phidiff -_charge*tempz/tanradius)/_const_2pi;
+    int n1,n2;
+    if (xCircles >= 0.) {
+      n1 = int(xCircles);
+      n2 = n1 + 1;
+    }
+    else {
+      n1 = int(xCircles) - 1;
+      n2 = n1 + 1;
+    }
+    if (fabs(n1-xCircles) < fabs(n2-xCircles)) {
+      nCircles = n1;
+    }
+    else {
+      nCircles = n2;
+    }
+  }
+  float DistZ = - tempz - _charge*tanradius*(_const_2pi*((float)nCircles) - phidiff);
+  return sqrt(DistXY*DistXY+DistZ*DistZ);
+}//getDistanceToPoint(vector,float)
+
+float HelixClassD::getDistanceToPoint(const float* xPoint, float distCut) {
+  std::vector<float> xPosition(xPoint, xPoint + 3 );//We are expecting three coordinates, must be +3, last element is excluded!
+  return getDistanceToPoint(xPosition, distCut);
+}//getDistanceToPoint(float*,float)
+
+
+
+void HelixClassD::setHelixEdges(float * xStart, float * xEnd) {
+  for (int i=0; i<3; ++i) {
+    _xStart[i] = xStart[i];
+    _xEnd[i] = xEnd[i];
+  }
+
+}
+
+float HelixClassD::getDistanceToHelix(HelixClassD * helix, float * pos, float * mom) {
+
+  float x01 = getXC();
+  float y01 = getYC();
+
+  float x02 = helix->getXC();
+  float y02 = helix->getYC();
+
+  float rad1 = getRadius();
+  float rad2 = helix->getRadius();
+
+  float distance = sqrt((x01-x02)*(x01-x02)+(y01-y02)*(y01-y02));
+
+  bool singlePoint = true;
+
+  float phi1 = 0;
+  float phi2 = 0;
+
+  if (rad1+rad2<distance) {
+    phi1 = atan2(y02-y01,x02-x01);
+    phi2 = atan2(y01-y02,x01-x02);
+  }
+  else if (distance+rad2<rad1) {
+    phi1 = atan2(y02-y01,x02-x01);
+    phi2 = phi1;
+  }
+  else if (distance+rad1<rad2) {
+    phi1 = atan2(y01-y02,x01-x02);
+    phi2 = phi1;
+  }
+  else {
+    singlePoint = false;
+    float cosAlpha = 0.5*(distance*distance+rad2*rad2-rad1*rad1)/(distance*rad2);
+    float alpha = acos(cosAlpha);
+    float phi0 = atan2(y01-y02,x01-x02);
+    phi1 = phi0 + alpha;
+    phi2 = phi0 - alpha;
+  }
+
+
+  float ref1[3];
+  float ref2[3];
+  for (int i=0;i<3;++i) {
+    ref1[i]=_referencePoint[i];
+    ref2[i]=helix->getReferencePoint()[i];
+  }
+
+  float pos1[3];
+  float pos2[3];
+  float mom1[3];
+  float mom2[3];
+
+
+  if (singlePoint ) {
+
+    float xSect1 = x01 + rad1*cos(phi1);
+    float ySect1 = y01 + rad1*sin(phi1);
+    float xSect2 = x02 + rad2*cos(phi2);
+    float ySect2 = y02 + rad2*sin(phi2);
+    float R1 = sqrt(xSect1*xSect1+ySect1*ySect1);
+    float R2 = sqrt(xSect2*xSect2+ySect2*ySect2);
+
+    getPointOnCircle(R1,ref1,pos1);
+    helix->getPointOnCircle(R2,ref2,pos2);
+
+  }
+  else {    
+
+    float xSect1 = x02 + rad2*cos(phi1);
+    float ySect1 = y02 + rad2*sin(phi1);
+    float xSect2 = x02 + rad2*cos(phi2);
+    float ySect2 = y02 + rad2*sin(phi2);
+
+    //     std::cout << "(xSect1,ySect1)=(" << xSect1 << "," << ySect1 << ")" << std::endl;
+    //     std::cout << "(xSect2,ySect2)=(" << xSect2 << "," << ySect2 << ")" << std::endl;
+
+    float temp21[3];
+    float temp22[3];
+
+    float phiI2  = atan2(ref2[1]-y02,ref2[0]-x02); 
+    float phiF21 = atan2(ySect1-y02,xSect1-x02);
+    float phiF22 = atan2(ySect2-y02,xSect2-x02);
+    float deltaPhi21 = phiF21 - phiI2;
+    float deltaPhi22 = phiF22 - phiI2;
+    float charge2 = helix->getCharge();
+    float pxy2 = helix->getPXY();
+    float pz2   = helix->getMomentum()[2];
+    if (deltaPhi21 < 0 && charge2 < 0) {
+      deltaPhi21 += _const_2pi;
+    }
+    else if (deltaPhi21 > 0 && charge2 > 0) { 
+      deltaPhi21 -= _const_2pi;
+    }
+
+    if (deltaPhi22 < 0 && charge2 < 0) {
+      deltaPhi22 += _const_2pi;
+    }
+    else if (deltaPhi22 > 0 && charge2 > 0) { 
+      deltaPhi22 -= _const_2pi;
+    }
+
+    float time21 = -charge2*deltaPhi21*rad2/pxy2;
+    float time22 = -charge2*deltaPhi22*rad2/pxy2;
+
+    float Z21 = ref2[2] + time21*pz2;
+    float Z22 = ref2[2] + time22*pz2;
+
+    temp21[0] = xSect1; temp21[1] = ySect1; temp21[2] = Z21;
+    temp22[0] = xSect2; temp22[1] = ySect2; temp22[2] = Z22;
+
+
+    //     std::cout << "temp21 = " << temp21[0] << " " << temp21[1] << " " << temp21[2] << std::endl;
+    //     std::cout << "temp22 = " << temp22[0] << " " << temp22[1] << " " << temp22[2] << std::endl;
+
+
+    float temp11[3];
+    float temp12[3];
+
+    float phiI1  = atan2(ref1[1]-y01,ref1[0]-x01); 
+    float phiF11 = atan2(ySect1-y01,xSect1-x01);
+    float phiF12 = atan2(ySect2-y01,xSect2-x01);
+    float deltaPhi11 = phiF11 - phiI1;
+    float deltaPhi12 = phiF12 - phiI1;
+    float charge1 = _charge;
+    float pxy1 = _pxy;
+    float pz1   = _momentum[2];
+    if (deltaPhi11 < 0 && charge1 < 0) {
+      deltaPhi11 += _const_2pi;
+    }
+    else if (deltaPhi11 > 0 && charge1 > 0) { 
+      deltaPhi11 -= _const_2pi;
+    }
+
+    if (deltaPhi12 < 0 && charge1 < 0) {
+      deltaPhi12 += _const_2pi;
+    }
+    else if (deltaPhi12 > 0 && charge1 > 0) { 
+      deltaPhi12 -= _const_2pi;
+    }
+
+    float time11 = -charge1*deltaPhi11*rad1/pxy1;
+    float time12 = -charge1*deltaPhi12*rad1/pxy1;
+
+    float Z11 = ref1[2] + time11*pz1;
+    float Z12 = ref1[2] + time12*pz1;
+
+    temp11[0] = xSect1; temp11[1] = ySect1; temp11[2] = Z11;
+    temp12[0] = xSect2; temp12[1] = ySect2; temp12[2] = Z12;
+
+    //     std::cout << "temp11 = " << temp11[0] << " " << temp11[1] << " " << temp11[2] << std::endl;
+    //     std::cout << "temp12 = " << temp12[0] << " " << temp12[1] << " " << temp12[2] << std::endl;
+
+    float Dist1 = 0;
+    float Dist2 = 0;
+
+    for (int j=0;j<3;++j) {
+      Dist1 += (temp11[j]-temp21[j])*(temp11[j]-temp21[j]);
+      Dist2 += (temp12[j]-temp22[j])*(temp12[j]-temp22[j]);
+    }
+
+    if (Dist1<Dist2) {
+      for (int l=0;l<3;++l) {
+        pos1[l] = temp11[l];
+        pos2[l] = temp21[l];
+      }
+    }
+    else {
+      for (int l=0;l<3;++l) {
+        pos1[l] = temp12[l];
+        pos2[l] = temp22[l];
+      }
+    }
+
+  }
+
+  getExtrapolatedMomentum(pos1,mom1);
+  helix->getExtrapolatedMomentum(pos2,mom2);
+
+  float helixDistance = 0;
+
+  for (int i=0;i<3;++i) {
+    helixDistance += (pos1[i] - pos2[i])*(pos1[i] - pos2[i]);
+    pos[i] = 0.5*(pos1[i]+pos2[i]);
+    mom[i] = mom1[i]+mom2[i];
+  }
+  helixDistance = sqrt(helixDistance);
+
+  return helixDistance;
+
+}
+
+void HelixClassD::getExtrapolatedMomentum(float * pos, float * momentum) {
+
+  float phi = atan2(pos[1]-_yCentre,pos[0]-_xCentre);
+  if (phi <0.) phi += _const_2pi;
+  phi = phi - _phiAtPCA + _phi0;
+  momentum[0] = _pxy*cos(phi);
+  momentum[1] = _pxy*sin(phi);
+  momentum[2] = _momentum[2];
+
+
+}

Analysis/ReadDigi/src/HelixClassD.hh

+#ifndef HELIXARD_HH
+#define HELIXARD_HH 1
+#include <vector>
+/**
+ *    Utility class to manipulate with different parameterisations <br>
+ *    of helix. Helix can be initialized in a three different <br>
+ *    ways using the following public methods : <br>
+ *    1) Initialize_VP(float * pos, float * mom, float q, float B) : <br>
+ *       initialization of helix is done using <br>
+ *       - position of the reference point : pos[3]; <br>
+ *       - momentum vector at the reference point : mom[3];<br>
+ *       - particle charge : q;<br>
+ *       - magnetic field (in Tesla) : B;<br>
+ *    2) Initialize_BZ(float xCentre, float yCentre, float radius, <br>
+ *				   float bZ, float phi0, float B, float signPz,<br>
+ *				   float zBegin):<br>
+ *       initialization of helix is done according to the following<br>
+ *       parameterization<br>
+ *       x = xCentre + radius*cos(bZ*z + phi0)<br>
+ *       y = yCentre + radius*sin(bZ*z + phi0)<br>
+ *       where (x,y,z) - position of point on the helix<br>
+ *       - (xCentre,yCentre) is the centre of circumference in R-Phi plane<br>
+ *       - radius is the radius of circumference<br>
+ *       - bZ is the helix slope parameter<br>
+ *       - phi0 is the initial phase of circumference<br>
+ *       - B is the magnetic field (in Tesla)<br>
+ *       - signPz is the sign of the z component of momentum vector<br>
+ *       - zBegin is the z coordinate of the reference point;<br>
+ *    3) void Initialize_Canonical(float phi0, float d0, float z0, float omega,<br> 
+ *			      float tanlambda, float B) :<br>
+ *    canonical (LEP-wise) parameterisation with the following parameters<br>
+ *       - phi0 - Phi angle of momentum vector at the point of<br>
+ *         closest approach to IP in R-Phi plane;
+ *       - d0 - signed distance of closest approach to IP in R-Phi plane;<br>
+ *       - z0 - z coordinate of the point of closest approach in R-Phi plane;<br>
+ *       - omega - signed curvature;<br>
+ *       - tanlambda - tangent of dip angle;<br>
+ *       - B - magnetic field (in Tesla);<br>
+ *    A set of public methods (getters) provide access to <br>
+ *    various parameters associated with helix. Helix Class contains<br>
+ *    also several utility methods, allowing for calculation of helix<br>
+ *    intersection points with planes parallel and perpendicular to <br>
+ *    z (beam) axis and determination of the distance of closest approach<br>
+ *    from arbitrary 3D point to the helix. <br>
+ *    @author A. Raspereza (DESY)<br>
+ *    @version $Id: HelixClassD.h,v 1.16 2008-06-05 13:47:18 rasp Exp $<br>
+ *
+ */
+
+//#include "LineClass.h"
+class HelixClassD;
+
+class HelixClassD {
+ public:
+
+/**
+ *  Constructor. Initializations of constants which are used
+ *  to calculate various parameters associated with helix.
+ */ 
+    HelixClassD();
+/**
+ *  Destructor. 
+ */
+    ~HelixClassD();
+/**
+ *   Initialization of helix using <br>
+ *     - position of the reference point : pos[3]; <br>
+ *     - momentum vector at the reference point : mom[3];<br>
+ *     - particle charge : q;<br>
+ *     - magnetic field (in Tesla) : B<br>
+ */  
+    void Initialize_VP(float * pos, float * mom, float q, float B);
+
+/**
+ *   Initialization of helix according to the following<br>
+ *   parameterization<br>
+ *   x = xCentre + radius*cos(bZ*z + phi0)<br>
+ *   y = yCentre + radius*sin(bZ*z + phi0)<br>
+ *     where (x,y,z) - position of point on the helix<br>
+ *     - (xCentre,yCentre) is the centre of circumference in R-Phi plane<br>
+ *     - radius is the radius of circumference<br>
+ *     - bZ is the helix slope parameter<br>
+ *     - phi0 is the initial phase of circumference<br>
+ *     - B is the magnetic field (in Tesla)<br>
+ *     - signPz is the sign of the z component of momentum vector<br>
+ *     - zBegin is the z coordinate of the reference point<br>
+ */  
+    void Initialize_BZ(float xCentre, float yCentre, float radius, 
+				   float bZ, float phi0, float B, float signPz,
+				   float zBegin);
+/**
+ *  Canonical (LEP-wise) parameterisation with the following parameters<br>
+ *     - phi0 - Phi angle of momentum vector at the point of<br>
+ *       closest approach to IP in R-Phi plane;
+ *     - d0 - signed distance of closest approach in R-Phi plane;<br>
+ *     - z0 - z coordinate of the point of closest approach to IP 
+ *       in R-Phi plane;<br>
+ *     - omega - signed curvature;<br>
+ *     - tanlambda - tangent of dip angle;<br>
+ *     - B - magnetic field (in Tesla)<br>
+ */  
+    void Initialize_Canonical(float phi0, float d0, float z0, float omega, 
+			      float tanlambda, float B);
+    /**
+     *  Returns momentum of particle at the point of closest approach <br>
+     *  to IP <br>
+     */
+    const float * getMomentum();
+
+    /**
+     *  Returns reference point of track <br>     
+     */
+    const float * getReferencePoint();
+
+    /**
+     *  Returns Phi angle of the momentum vector <br>
+     *  at the point of closest approach to IP <br>     
+     */
+    float getPhi0();
+
+    /**
+     *  Returns signed distance of closest approach <br>
+     *  to IP in the R-Phi plane <br>     
+     */
+    float getD0();
+
+    /**
+     *  Returns z coordinate of the point of closest 
+     *  approach to IP in the R-Phi plane <br>     
+     */
+    float getZ0();
+
+    /**
+     *  Returns signed curvature of the track <br>     
+     */
+    float getOmega();
+
+    /**
+     *  Returns tangent of dip angle of the track <br>     
+     */
+    float getTanLambda();
+
+    /**
+     *  Returns transverse momentum of the track <br>     
+     */
+    float getPXY();
+
+
+    /**
+     *  Returns x coordinate of circumference
+     */
+    float getXC();
+
+    /**
+     *  Returns y coordinate of circumference
+     */
+    float getYC();
+
+
+     /**
+     *  Returns radius of circumference
+     */
+    float getRadius();   
+
+
+    /**
+     *  Returns helix intersection point with the plane <br>
+     *  parallel to z axis. Plane is defined by two coordinates <br>
+     *  of the point belonging to the plane (x0,y0) and normal <br>
+     *  vector (ax,ay).  ref[3] is the reference point of the helix. <br>
+     *  point[3] - returned vector holding the coordinates of <br>
+     *  intersection point <br>     
+     */
+    float getPointInXY(float x0, float y0, float ax, float ay, 
+			      float * ref , float * point);
+
+    /**
+     *  Returns helix intersection point with the plane <br>
+     *  perpendicular to z axis. Plane is defined by z coordinate <br>
+     *  of the plane. ref[3] is the reference point of the helix. <br>
+     *  point[3] - returned vector holding the coordinates of <br>
+     *  intersection point <br>     
+     */
+    float getPointInZ(float zLine, float * ref, float * point);
+
+    /**
+     * Return distance of the closest approach of the helix to <br>
+     * arbitrary 3D point in space. xPoint[3] - coordinates of <br>
+     * space point. Distance[3] - vector of distances of helix to <br> 
+     * a given point in various projections : <br>
+     * Distance[0] - distance in R-Phi plane <br>
+     * Distance[1] - distance along Z axis <br>
+     * Distance[2] - 3D distance <br> 
+     */
+    float getDistanceToPoint(float * xPoint, float * Distance);
+
+    /**
+     * Return distance of the closest approach of the helix to <br>
+     * arbitrary 3D point in space. xPoint[3] - coordinates of <br>
+     * space point. distCut - limit on the distance between helix <br> 
+     * and the point to reduce calculation time <br>
+     * If R-Phi is found to be greater than distCut, rPhi distance is returned <br>
+     * If the R-Phi distance is not too big, than the exact 3D distance is returned <br>
+     * This function can be used, if the exact distance is not always needed <br>
+     */
+    float getDistanceToPoint(const float* xPoint, float distCut);
+    float getDistanceToPoint(const std::vector<float>& xPoint, float distCut);
+
+    /**
+     * This method calculates coordinates of both intersection <br>
+     * of the helix with a cylinder. <br>
+     * Rotational symmetry with respect to z axis is assumed,  <br>
+     * meaning that cylinder axis corresponds to the z axis <br>
+     * of reference frame. <br>
+     * Inputs : <br> 
+     * Radius - radius of cylinder. <br>
+     * ref[3] - point of closest approach to the origin of the helix. <br>
+     * Output : <br>
+     * point[6] - coordinates of intersection point. <br>
+     * Method returns also generic time, defined as the <br>
+     * ratio of helix length from reference point to the intersection <br>
+     * point to the particle momentum <br>
+     */
+    float getPointOnCircle(float Radius, float * ref, float * point);
+
+    /** Returns distance between two helixes <br>
+     * Output : <br> 
+     * pos[3] - position of the point of closest approach <br>
+     * mom[3] - momentum of V0 <br>
+     */
+    float getDistanceToHelix(HelixClassD * helix, float * pos, float * mom);
+   
+    int getNCircle(float * XPoint);
+
+    /**
+     * Set Edges of helix 
+     */
+    void setHelixEdges(float * xStart, float * xEnd);
+
+    /**
+     * Returns starting point of helix
+     */
+    float * getStartingPoint() {return _xStart;}
+
+    /**
+     * Returns endpoint of helix
+     */
+    float * getEndPoint() {return _xEnd;}
+    
+    /**
+     * Returns BZ for the second parameterization
+     */
+    float getBz();
+
+    /**
+     * Returns Phi for the second parameterization
+     */
+    float getPhiZ();
+
+    /**
+     * Returns extrapolated momentum
+     */
+    void getExtrapolatedMomentum(float * pos, float * momentum);
+
+    /**
+     * Returns charge 
+     */
+    float getCharge();
+
+ private:    
+    float _momentum[3]; // momentum @ ref point 
+    float _referencePoint[3]; // coordinates @ ref point
+    float _phi0; // phi0 in canonical parameterization 
+    float _d0;   // d0 in canonical parameterisation
+    float _z0;   // z0 in canonical parameterisation
+    float _omega; // signed curvuture in canonical parameterisation
+    float _tanLambda; // TanLambda 
+    float _pxy; // Transverse momentum
+    float _charge; // Particle Charge
+    float _bField; // Magnetic field (assumed to point to Z>0)
+    float _radius; // radius of circle in XY plane
+    float _xCentre; // X of circle centre
+    float _yCentre; // Y of circle centre
+    float _phiRefPoint; // Phi w.r.t. (X0,Y0) of circle @ ref point
+    float _phiAtPCA; // Phi w.r.t. (X0,Y0) of circle @ PCA 
+    float _xAtPCA; // X @ PCA
+    float _yAtPCA; // Y @ PCA
+    float _pxAtPCA; // PX @ PCA
+    float _pyAtPCA; // PY @ PCA
+    float _phiMomRefPoint; // Phi of Momentum vector @ ref point
+    float _const_pi; // PI
+    float _const_2pi; // 2*PI
+    float _const_pi2; // PI/2    
+    float _FCT; // 2.99792458E-4
+    float _xStart[3]; // Starting point of track segment
+    float _xEnd[3]; // Ending point of track segment
+
+    float _bZ;
+    float _phiZ;
+
+};
+
+
+#endif

Analysis/ReadDigi/src/ReadDigiAlg.h

+#ifndef READ_DIGI_H
+#define READ_DIGI_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/TrackerHit.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "HelixClassD.hh"
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TVector3.h"
+using namespace std;
+
+class ReadDigiAlg : public GaudiAlgorithm
+{
+
+public :
+
+  ReadDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+
+  virtual StatusCode initialize();
+  virtual StatusCode execute();
+  virtual StatusCode finalize();
+
+  StatusCode Clear(); 
+
+private :
+  DataHandle<edm4hep::MCParticleCollection> m_MCParticleCol{"MCParticle", Gaudi::DataHandle::Reader, this};
+
+  //DataHandle<edm4hep::TrackerHitCollection> m_SITRawColHdl{"SITTrackerHits", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::TrackerHitCollection> m_SETRawColHdl{"SETTrackerHits", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::TrackerHitCollection> m_FTDRawColHdl{"FTDStripTrackerHits", Gaudi::DataHandle::Reader, this};  
+
+  //DataHandle<edm4hep::TrackerHitCollection> m_VTXTrackerHitColHdl{"VTXTrackerHits", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::TrackerHitCollection> m_SITTrackerHitColHdl{"SITSpacePoints", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::TrackerHitCollection> m_TPCTrackerHitColHdl{"TPCTrackerHits", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::TrackerHitCollection> m_SETTrackerHitColHdl{"SETSpacePoints", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::TrackerHitCollection> m_FTDSpacePointColHdl{"FTDSpacePoints", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::TrackerHitCollection> m_FTDPixelTrackerHitColHdl{"FTDPixelTrackerHits", Gaudi::DataHandle::Reader, this};
+
+  DataHandle<edm4hep::TrackCollection>      m_SiTrk{"SiTracks", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::TrackCollection>      m_TPCTrk{"ClupatraTracks", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::TrackCollection>      m_fullTrk{"CompleteTracks", Gaudi::DataHandle::Reader, this};
+
+  DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> m_TPCTrkAssoHdl{"ClupatraTracksParticleAssociation",  Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> m_fullTrkAssoHdl{"CompleteTracksParticleAssociation",  Gaudi::DataHandle::Reader, this};
+
+  //DataHandle<edm4hep::SimCalorimeterHitCollection> m_ECalBarrelHitCol{"EcalBarrelCollection", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::SimCalorimeterHitCollection> m_ECalEndcapHitCol{"EcalEndcapCollection", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::SimCalorimeterHitCollection> m_HCalBarrelHitCol{"HcalBarrelCollection", Gaudi::DataHandle::Reader, this};
+  //DataHandle<edm4hep::SimCalorimeterHitCollection> m_HCalEndcapHitCol{"HcalEndcapsCollection", Gaudi::DataHandle::Reader, this};
+
+  mutable Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"};
+  mutable Gaudi::Property<double> _Bfield{this, "BField", 3., "Magnet field"};
+
+
+  typedef std::vector<float> FloatVec;
+  typedef std::vector<int>   IntVec;
+
+
+  int _nEvt; 
+  TFile *m_wfile;
+  TTree *m_mctree, *m_sitrktree, *m_tpctrktree, *m_fulltrktree;
+
+  //MCParticle
+  int N_MCP;
+  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 
+  int N_SiTrk, N_TPCTrk, N_fullTrk;
+  //int N_VTXhit, N_SIThit, N_TPChit, N_SEThit, N_FTDhit, N_SITrawhit, N_SETrawhit;
+  FloatVec m_trk_px, m_trk_py, m_trk_pz, m_trk_p;
+  FloatVec m_trkstate_d0, m_trkstate_z0, m_trkstate_phi, m_trkstate_tanL, m_trkstate_omega, m_trkstate_kappa;
+  FloatVec m_trkstate_refx, m_trkstate_refy, m_trkstate_refz;
+  IntVec m_trkstate_tag, m_trkstate_location, m_trk_Nhit, m_trk_type;
+  FloatVec m_truthMC_tag, m_truthMC_pid, m_truthMC_px, m_truthMC_py, m_truthMC_pz, m_truthMC_E;
+  FloatVec m_truthMC_EPx, m_truthMC_EPy, m_truthMC_EPz, m_truthMC_weight;
+
+
+  //ECal
+  //int Nhit_EcalB;
+  //int Nhit_EcalE;
+  //int Nhit_HcalB; 
+  //int Nhit_HcalE;
+  //IntVec CaloHit_type; //0: EM hit.   1: Had hit.
+  //FloatVec CaloHit_x, CaloHit_y, CaloHit_z, CaloHit_E, CaloHit_Eem, CaloHit_Ehad, CaloHit_aveT, CaloHit_Tmin, CaloHit_TmaxE;
+
+
+};
+
+#endif