.build.ci.sh

 #!/bin/bash
-# This is wrapper to run the build.sh on CI
+# This is wrapper to run the build.sh or build-xyz.sh on CI
+
+# The build mode is the suffix in build-xyz.sh
+export BUILD_CI_MODE=${1}
 
 echo "CEPCSW_LCG_RELEASE: ${CEPCSW_LCG_RELEASE}"
 echo "CEPCSW_LCG_PLATFORM: ${CEPCSW_LCG_PLATFORM}"
@@ -36,12 +39,26 @@ function build-with-log() {
 }
 
 function build-with-stdout() {
+    local build_flags=${BUILD_CI_MODE}
+    local source_flag=true
+
+    # Key4hep stack mode
     if [ "$CEPCSW_LCG_RELEASE" = "KEY4HEP_STACK" ]; then
-        ./build-k4.sh
-    else
+        build_flags=k4
+        source_flag=false
+    fi
+
+    # prepend '-' if necessary
+    if [ -n "$build_flags" ]; then
+       build_flags=-${build_flags}
+    fi
+
+    if $source_flag; then
         source setup.sh
-        ./build.sh
     fi
+    
+    ./build${build_flags}.sh
+
 }
 
 if [ -n "${GITHUB_ACTION}" ]; then

.gitignore

@@ -5,3 +5,4 @@ spack*
 ./Generator/options/
 
 InstallArea/
+venv

.gitlab-ci.yml

@@ -26,6 +26,7 @@ workflow:
 
 stages:
   - build
+  - test
 
 ##############################################################################
 # Template for Build and Test
@@ -34,14 +35,22 @@ stages:
 # the test job will be failed if it is executed on a different nodes.
 # Therefore, put the build script and test script together.
 
+.envvar_template:
+  variables:
+    CEPCSW_LCG_RELEASE: LCG
+    CEPCSW_LCG_PLATFORM: x86_64-centos7-gcc11-opt
+    CEPCSW_LCG_VERSION: 103.0.2
+
 .build_template:
+  extends: .envvar_template
   stage: build
-  variables:
-    CEPCSW_LCG_RELEASE:
-    CEPCSW_LCG_PLATFORM:
-    CEPCSW_LCG_VERSION:
   script:
     - bash ./.build.ci.sh
+
+.test_template:
+  extends: .envvar_template
+  stage: test
+  script:
     - bash ./.test.ci.sh
 
 ##############################################################################
@@ -49,16 +58,38 @@ stages:
 ##############################################################################
 build:lcg:el7:
   extends: .build_template
-  variables:
-    CEPCSW_LCG_RELEASE: LCG
-    CEPCSW_LCG_PLATFORM: x86_64-centos7-gcc11-opt
-    CEPCSW_LCG_VERSION: 103.0.2
   tags:
     - centos7
   artifacts:
     paths:
       - InstallArea
+      - build.${CEPCSW_LCG_VERSION}.${CEPCSW_LCG_PLATFORM}
+
+##############################################################################
+# Test CentOS 7 (LCG)
+##############################################################################
+test:lcg:el7:
+  extends: .test_template
+  tags:
+    - centos7
+  artifacts:
+    paths:
       - TDR_o1_v01.tgeo.root
       - TDR_o1_v02.tgeo.root
     reports:
       junit: build.${CEPCSW_LCG_VERSION}.${CEPCSW_LCG_PLATFORM}/cepcsw-ctest-result.xml
+  dependencies:
+    - build:lcg:el7
+
+##############################################################################
+# Build the docs
+##############################################################################
+build:docs:
+  stage: build 
+  tags:
+    - centos7
+  script:
+    - bash ./.build.ci.sh docs
+  artifacts:
+    paths:
+      - docs/build/html/

.readthedocs.yaml

+version: "2"
+
+build:
+  os: "ubuntu-22.04"
+  tools:
+    python: "3.10"
+
+python:
+  install:
+    - requirements: docs/requirements.txt
+
+sphinx:
+  configuration: docs/source/conf.py

.test.ci.sh

@@ -43,4 +43,9 @@ echo "CEPCSW_BLDTOOL: ${CEPCSW_BLDTOOL}"
 
 source setup.sh
 
+# reconfigure if directory change
+pushd $(build-dir)
+cmake ..
+popd
+
 ctest --output-junit $(junit-output) --test-dir $(build-dir)

Analysis/CMakeLists.txt

@@ -2,3 +2,4 @@
 add_subdirectory(TotalInvMass)
 add_subdirectory(TrackInspect)
 add_subdirectory(DumpEvent)
+add_subdirectory(ReadDigi)

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.cpp

+#ifndef READ_DIGI_C
+#define READ_DIGI_C
+
+#include "ReadDigiAlg.h"
+
+DECLARE_COMPONENT(ReadDigiAlg)
+
+ReadDigiAlg::ReadDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+    : GaudiAlgorithm(name, svcLoc), 
+     _nEvt(0)
+{
+  declareProperty("MCParticle",  m_MCParticleCol, "MCParticle collection (input)");
+
+  //Tracker hit
+  declareProperty("SITRawHits", m_SITRawColHdl, "SIT Raw Hit Collection of SpacePoints");
+  declareProperty("SETRawHits", m_SETRawColHdl, "SET Raw Hit Collection of SpacePoints");
+  declareProperty("FTDRawHits", m_FTDRawColHdl, "FTD Raw Hit Collection of SpacePoints");
+  declareProperty("VTXTrackerHits", m_VTXTrackerHitColHdl, "VTX Hit Collection");
+  declareProperty("SITTrackerHits", m_SITTrackerHitColHdl, "SIT Hit Collection");
+  declareProperty("SETTrackerHits", m_SETTrackerHitColHdl, "SET Hit Collection");
+  declareProperty("TPCTrackerHits", m_TPCTrackerHitColHdl, "TPC Hit Collection");
+  declareProperty("FTDSpacePoints", m_FTDSpacePointColHdl, "FTD FTDSpacePoint Collection");
+  declareProperty("FTDPixelTrackerHits", m_FTDPixelTrackerHitColHdl, "handler of FTD Pixel Hit Collection");
+
+  //Track
+  declareProperty("FullTracks", m_fullTrk, "Full Track Collection");
+  declareProperty("TPCTracks", m_TPCTrk, "TPC Track Collection");
+  declareProperty("SiTracks", m_SiTrk, "Si Track Collection");
+
+  //declareProperty("EcalBarrelCollection",  m_ECalBarrelHitCol, "ECal Barrel");
+  //declareProperty("EcalEndcapsCollection", m_ECalEndcapHitCol, "ECal Endcap");
+  //declareProperty("HcalBarrelCollection",  m_HCalBarrelHitCol, "HCal Barrel");
+  //declareProperty("HcalEndcapsCollection", m_HCalEndcapHitCol, "HCal Endcap");
+
+}
+
+StatusCode ReadDigiAlg::initialize()
+{
+  debug() << "begin initialize ReadDigiAlg" << endmsg;
+  std::string s_outfile = _filename;
+  m_wfile = new TFile(s_outfile.c_str(), "recreate");
+  m_mctree = new TTree("MCParticle", "MCParticle");
+  m_trktree = new TTree("RecoTrk", "RecoTrk");
+
+  m_mctree->Branch("N_MCP", &N_MCP);
+  m_mctree->Branch("MCP_px", &MCP_px);
+  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_endPoint_x", &MCP_endPoint_x);
+  m_mctree->Branch("MCP_endPoint_y", &MCP_endPoint_y);
+  m_mctree->Branch("MCP_endPoint_z", &MCP_endPoint_z);
+  m_mctree->Branch("MCP_pdgid", &MCP_pdgid);
+  m_mctree->Branch("MCP_gStatus", &MCP_gStatus);
+
+  m_trktree->Branch("N_SiTrk", &N_SiTrk);
+  m_trktree->Branch("N_TPCTrk", &N_TPCTrk);
+  m_trktree->Branch("N_fullTrk", &N_fullTrk);
+  m_trktree->Branch("N_VTXhit", &N_VTXhit);
+  m_trktree->Branch("N_SIThit", &N_SIThit);
+  m_trktree->Branch("N_TPChit", &N_TPChit);
+  m_trktree->Branch("N_SEThit", &N_SEThit);
+  m_trktree->Branch("N_FTDhit", &N_FTDhit);
+  m_trktree->Branch("N_SITrawhit", &N_SITrawhit);
+  m_trktree->Branch("N_SETrawhit", &N_SETrawhit);
+  m_trktree->Branch("trk_px", &m_trk_px);
+  m_trktree->Branch("trk_py", &m_trk_py);
+  m_trktree->Branch("trk_pz", &m_trk_pz);
+  m_trktree->Branch("trk_p", &m_trk_p);
+  m_trktree->Branch("trk_startpoint_x", &m_trk_startpoint_x);
+  m_trktree->Branch("trk_startpoint_y", &m_trk_startpoint_y);
+  m_trktree->Branch("trk_startpoint_z", &m_trk_startpoint_z);
+  m_trktree->Branch("trk_endpoint_x", &m_trk_endpoint_x);
+  m_trktree->Branch("trk_endpoint_y", &m_trk_endpoint_y);
+  m_trktree->Branch("trk_endpoint_z", &m_trk_endpoint_z);
+  m_trktree->Branch("trk_type", &m_trk_type);
+  m_trktree->Branch("trk_Nhit", &m_trk_Nhit);
+
+  return GaudiAlgorithm::initialize();
+}
+
+StatusCode ReadDigiAlg::execute()
+{
+  if(_nEvt==0) std::cout<<"ReadDigiAlg::execute Start"<<std::endl;
+  debug() << "Processing event: "<<_nEvt<< endmsg; 
+
+  Clear(); 
+
+  try{
+  const edm4hep::MCParticleCollection* const_MCPCol = m_MCParticleCol.get();
+  if(const_MCPCol){
+    N_MCP = const_MCPCol->size(); 
+    for(int i=0; i<N_MCP; i++){
+      edm4hep::MCParticle m_MCp = const_MCPCol->at(i);
+      MCP_px.push_back(m_MCp.getMomentum().x);
+      MCP_py.push_back(m_MCp.getMomentum().y);
+      MCP_pz.push_back(m_MCp.getMomentum().z);
+      MCP_E.push_back(m_MCp.getEnergy());
+      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); 
+      MCP_pdgid.push_back(m_MCp.getPDG());
+      MCP_gStatus.push_back(m_MCp.getGeneratorStatus());
+    }
+  }
+  }catch(GaudiException &e){
+    debug()<<"MC Particle is not available "<<endmsg;
+  }
+
+  try{
+  if(m_VTXTrackerHitColHdl.get())
+    N_VTXhit = m_VTXTrackerHitColHdl.get()->size();
+  }catch(GaudiException &e){
+    debug()<<"VTX hit is not available "<<endmsg;
+  }
+
+  try{
+  if(m_SITTrackerHitColHdl.get())
+    N_SIThit = m_SITTrackerHitColHdl.get()->size();
+  }catch(GaudiException &e){
+    debug()<<"SIT hit is not available "<<endmsg;
+  }
+
+  try{
+  if(m_TPCTrackerHitColHdl.get())
+    N_TPChit = m_TPCTrackerHitColHdl.get()->size();
+  }catch(GaudiException &e){
+    debug()<<"TPC hit is not available "<<endmsg;
+  }
+
+  try{
+  if(m_SETTrackerHitColHdl.get())
+    N_SEThit = m_SETTrackerHitColHdl.get()->size();
+  }catch(GaudiException &e){
+    debug()<<"SET hit is not available "<<endmsg;
+  }
+
+  try{
+  if(m_FTDSpacePointColHdl.get())
+    N_FTDhit = m_FTDSpacePointColHdl.get()->size();
+  }catch(GaudiException &e){
+    debug()<<"FDT space point is not available "<<endmsg;
+  }
+
+  try{
+  if(m_SITRawColHdl.get())
+    N_SITrawhit = m_SITRawColHdl.get()->size();
+  }catch(GaudiException &e){
+    debug()<<"SIT raw hit is not available "<<endmsg;
+  }
+
+  try{
+  if(m_SETRawColHdl.get())
+    N_SETrawhit = m_SETRawColHdl.get()->size();
+  }catch(GaudiException &e){
+    debug()<<"SET raw hit is not available "<<endmsg;
+  }
+
+  try{
+  auto const_SiTrkCol = m_SiTrk.get();
+  if(const_SiTrkCol){
+    N_SiTrk = const_SiTrkCol->size();
+    for(int i=0; i<N_SiTrk; i++){
+      auto m_trk = const_SiTrkCol->at(i);
+      if(m_trk.trackStates_size()==0) continue;
+      if(m_trk.trackerHits_size()==0) continue;
+      edm4hep::TrackState m_trkstate = m_trk.getTrackStates(0);
+
+      HelixClassD * TrkInit_Helix = new HelixClassD();
+      TrkInit_Helix->Initialize_Canonical(m_trkstate.phi, m_trkstate.D0, m_trkstate.Z0, m_trkstate.omega, m_trkstate.tanLambda, _Bfield);
+      TVector3 TrkMom(TrkInit_Helix->getMomentum()[0],TrkInit_Helix->getMomentum()[1],TrkInit_Helix->getMomentum()[2]);
+
+      int NTrkHit = m_trk.trackerHits_size();
+      TVector3 StartPointPos = TVector3((m_trk.getTrackerHits(0)).getPosition().x,(m_trk.getTrackerHits(0)).getPosition().y,(m_trk.getTrackerHits(0)).getPosition().z);
+      TVector3 EndPointPos = TVector3((m_trk.getTrackerHits(NTrkHit - 1)).getPosition().x,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().y,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().z);
+
+      m_trk_Nhit.push_back(NTrkHit);
+      m_trk_px.push_back(TrkMom.x());
+      m_trk_py.push_back(TrkMom.y());
+      m_trk_pz.push_back(TrkMom.z());
+      m_trk_p.push_back(TrkMom.Mag());
+      m_trk_endpoint_x.push_back(EndPointPos.x());
+      m_trk_endpoint_y.push_back(EndPointPos.y());
+      m_trk_endpoint_z.push_back(EndPointPos.z());
+      m_trk_startpoint_x.push_back(StartPointPos.x());
+      m_trk_startpoint_y.push_back(StartPointPos.y());
+      m_trk_startpoint_z.push_back(StartPointPos.z());
+      m_trk_type.push_back(1);
+
+      delete TrkInit_Helix;
+    }
+  }
+  }catch(GaudiException &e){
+    debug()<<"Si track is not available "<<endmsg;
+  }  
+
+  try{
+  auto const_TPCTrkCol = m_TPCTrk.get();
+  if(const_TPCTrkCol){
+    N_TPCTrk = const_TPCTrkCol->size();
+    for(int i=0; i<N_TPCTrk; i++){
+      auto m_trk = const_TPCTrkCol->at(i);
+      if(m_trk.trackStates_size()==0) continue;
+      if(m_trk.trackerHits_size()==0) continue;
+      edm4hep::TrackState m_trkstate = m_trk.getTrackStates(0);
+
+      HelixClassD * TrkInit_Helix = new HelixClassD();
+      TrkInit_Helix->Initialize_Canonical(m_trkstate.phi, m_trkstate.D0, m_trkstate.Z0, m_trkstate.omega, m_trkstate.tanLambda, _Bfield);
+      TVector3 TrkMom(TrkInit_Helix->getMomentum()[0],TrkInit_Helix->getMomentum()[1],TrkInit_Helix->getMomentum()[2]);
+
+      int NTrkHit = m_trk.trackerHits_size();
+      TVector3 StartPointPos = TVector3((m_trk.getTrackerHits(0)).getPosition().x,(m_trk.getTrackerHits(0)).getPosition().y,(m_trk.getTrackerHits(0)).getPosition().z);
+      TVector3 EndPointPos = TVector3((m_trk.getTrackerHits(NTrkHit - 1)).getPosition().x,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().y,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().z);
+
+      m_trk_Nhit.push_back(NTrkHit);
+      m_trk_px.push_back(TrkMom.x());
+      m_trk_py.push_back(TrkMom.y());
+      m_trk_pz.push_back(TrkMom.z());
+      m_trk_p.push_back(TrkMom.Mag());
+      m_trk_endpoint_x.push_back(EndPointPos.x());
+      m_trk_endpoint_y.push_back(EndPointPos.y());
+      m_trk_endpoint_z.push_back(EndPointPos.z());
+      m_trk_startpoint_x.push_back(StartPointPos.x());
+      m_trk_startpoint_y.push_back(StartPointPos.y());
+      m_trk_startpoint_z.push_back(StartPointPos.z());
+      m_trk_type.push_back(2);
+
+      delete TrkInit_Helix;
+    }
+  }
+  }catch(GaudiException &e){
+    debug()<<"TPC track is not available "<<endmsg;
+  }
+
+
+  try{
+  auto const_fullTrkCol = m_fullTrk.get();
+  if(const_fullTrkCol){
+    N_fullTrk = const_fullTrkCol->size();
+    for(int i=0; i<N_fullTrk; i++){
+      auto m_trk = const_fullTrkCol->at(i);
+      if(m_trk.trackStates_size()==0) continue;
+      if(m_trk.trackerHits_size()==0) continue;
+      edm4hep::TrackState m_trkstate = m_trk.getTrackStates(0);
+
+      HelixClassD * TrkInit_Helix = new HelixClassD();
+      TrkInit_Helix->Initialize_Canonical(m_trkstate.phi, m_trkstate.D0, m_trkstate.Z0, m_trkstate.omega, m_trkstate.tanLambda, _Bfield);
+      TVector3 TrkMom(TrkInit_Helix->getMomentum()[0],TrkInit_Helix->getMomentum()[1],TrkInit_Helix->getMomentum()[2]);   
+
+      int NTrkHit = m_trk.trackerHits_size();
+      TVector3 StartPointPos = TVector3((m_trk.getTrackerHits(0)).getPosition().x,(m_trk.getTrackerHits(0)).getPosition().y,(m_trk.getTrackerHits(0)).getPosition().z);
+      TVector3 EndPointPos = TVector3((m_trk.getTrackerHits(NTrkHit - 1)).getPosition().x,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().y,(m_trk.getTrackerHits(NTrkHit - 1)).getPosition().z);
+
+      m_trk_Nhit.push_back(NTrkHit);
+      m_trk_px.push_back(TrkMom.x());
+      m_trk_py.push_back(TrkMom.y());
+      m_trk_pz.push_back(TrkMom.z());
+      m_trk_p.push_back(TrkMom.Mag());
+      m_trk_endpoint_x.push_back(EndPointPos.x());
+      m_trk_endpoint_y.push_back(EndPointPos.y());
+      m_trk_endpoint_z.push_back(EndPointPos.z());
+      m_trk_startpoint_x.push_back(StartPointPos.x());
+      m_trk_startpoint_y.push_back(StartPointPos.y());
+      m_trk_startpoint_z.push_back(StartPointPos.z());
+      m_trk_type.push_back(3);
+
+      delete TrkInit_Helix;
+    }
+  }
+  }catch(GaudiException &e){
+    debug()<<"Full track is not available "<<endmsg;
+  }
+
+/*  const edm4hep::SimCalorimeterHitCollection* ECalBHitCol = m_ECalBarrelHitCol.get();
+  Nhit_EcalB = ECalBHitCol->size();
+  for(int i=0; i<ECalBHitCol->size(); i++){
+    edm4hep::SimCalorimeterHit CaloHit = ECalBHitCol->at(i);
+
+    CaloHit_type.push_back(0);
+    CaloHit_x.push_back(CaloHit.getPosition().x);
+    CaloHit_y.push_back(CaloHit.getPosition().y);
+    CaloHit_z.push_back(CaloHit.getPosition().z);
+    CaloHit_E.push_back(CaloHit.getEnergy());
+    CaloHit_Eem.push_back(CaloHit.getEnergy());
+    CaloHit_Ehad.push_back(0.);
+
+    Etot += CaloHit.getEnergy();
+    Etot_ecal += CaloHit.getEnergy();
+  }
+*/
+
+/*  const edm4hep::SimCalorimeterHitCollection* ECalEHitCol = m_ECalEndcapHitCol.get();
+  Nhit_EcalE = ECalEHitCol->size();
+  for(int i=0; i<ECalEHitCol->size(); i++){
+    edm4hep::SimCalorimeterHit CaloHit = ECalEHitCol->at(i);
+
+    CaloHit_type.push_back(1);
+    CaloHit_x.push_back(CaloHit.getPosition().x);
+    CaloHit_y.push_back(CaloHit.getPosition().y);
+    CaloHit_z.push_back(CaloHit.getPosition().z);
+    CaloHit_E.push_back(CaloHit.getEnergy());
+    CaloHit_Eem.push_back(CaloHit.getEnergy());
+    CaloHit_Ehad.push_back(0.);
+  }
+*/
+/*  const edm4hep::SimCalorimeterHitCollection* HCalBHitCol = m_HCalBarrelHitCol.get();
+  Nhit_HcalB = HCalBHitCol->size();
+  for(int i=0; i<HCalBHitCol->size(); i++){
+    edm4hep::SimCalorimeterHit CaloHit = HCalBHitCol->at(i);
+    //if(CaloHit.getEnergy()<0.0001) continue;
+
+    CaloHit_x.push_back(CaloHit.getPosition().x);
+    CaloHit_y.push_back(CaloHit.getPosition().y);
+    CaloHit_z.push_back(CaloHit.getPosition().z);
+    CaloHit_E.push_back(CaloHit.getEnergy());
+
+    int Nconb = 0;
+    double Econb = 0.;
+    for(int iCont=0; iCont < CaloHit.contributions_size(); ++iCont){
+      auto conb = CaloHit.getContributions(iCont);
+      float conb_En = conb.getEnergy();
+      if( !conb.isAvailable() ) continue;
+      if(conb_En == 0) continue;
+
+      Nconb++;
+      Econb += conb_En;
+    }
+
+    Etot += CaloHit.getEnergy();
+  }
+*/
+
+  m_mctree->Fill(); 
+  m_trktree->Fill();
+  _nEvt++; 
+  return StatusCode::SUCCESS;
+}
+
+StatusCode ReadDigiAlg::finalize()
+{
+  debug() << "begin finalize ReadDigiAlg" << endmsg;
+  m_wfile->cd();
+  m_mctree->Write();
+  m_trktree->Write();
+  m_wfile->Close(); 
+
+  return GaudiAlgorithm::finalize();
+}
+
+StatusCode ReadDigiAlg::Clear()
+{
+
+  N_MCP = -99;
+  MCP_px.clear();
+  MCP_py.clear();
+  MCP_pz.clear();
+  MCP_E.clear();
+  MCP_endPoint_x.clear();
+  MCP_endPoint_y.clear();
+  MCP_endPoint_z.clear();
+  MCP_pdgid.clear();
+  MCP_gStatus.clear(); 
+
+  N_SiTrk = -99;
+  N_TPCTrk = -99;
+  N_fullTrk = -99;
+  N_VTXhit = -99; 
+  N_SIThit = -99; 
+  N_TPChit = -99; 
+  N_SEThit = -99; 
+  N_FTDhit = -99; 
+  N_SITrawhit = -99;
+  N_SETrawhit = -99;
+  m_trk_px.clear(); 
+  m_trk_py.clear(); 
+  m_trk_pz.clear();
+  m_trk_p.clear(); 
+  m_trk_endpoint_x.clear();
+  m_trk_endpoint_y.clear();
+  m_trk_endpoint_z.clear();
+  m_trk_startpoint_x.clear();
+  m_trk_startpoint_y.clear();
+  m_trk_startpoint_z.clear();
+  m_trk_type.clear(); 
+  m_trk_Nhit.clear();
+
+  return StatusCode::SUCCESS;
+}
+
+#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 "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{"MarlinTrkTracks", 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_trktree;
+
+  //MCParticle
+  int N_MCP;
+  FloatVec MCP_px, MCP_py, MCP_pz, MCP_E, 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_trk_endpoint_x, m_trk_endpoint_y, m_trk_endpoint_z, m_trk_startpoint_x, m_trk_startpoint_y, m_trk_startpoint_z;
+  IntVec m_trk_type, m_trk_Nhit; 
+
+
+  //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

CMakeLists.txt

@@ -33,6 +33,20 @@ if(NOT CMAKE_CXX_STANDARD MATCHES "17|20")
   message(FATAL_ERROR "Unsupported C++ standard: ${CMAKE_CXX_STANDARD}")
 endif()
 
+# Build type
+#   Use ``-DCMAKE_BUILD_TYPE=Debug`` when invoking CMake.
+#   By default, it is RelWithDebInfo since Sept 6th, 2024.
+if (NOT CMAKE_CONFIGURATION_TYPES)
+  if (NOT CMAKE_BUILD_TYPE)
+    set(CMAKE_BUILD_TYPE Release
+      CACHE STRING "Choose the type of build, options are: None|Release|MinSizeRel|Debug|RelWithDebInfo" FORCE)
+  else()
+    set(CMAKE_BUILD_TYPE ${CMAKE_BUILD_TYPE}
+      CACHE STRING "Choose the type of build, options are: None|Release|MinSizeRel|Debug|RelWithDebInfo" FORCE)
+  endif()
+endif()
+
+
 list(PREPEND CMAKE_MODULE_PATH $ENV{PANDORAPFA}/cmakemodules)
 list(PREPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake") # (Find*.cmake)
 include(cmake/CEPCSWOptions.cmake)

Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Endcaps_v02.cpp

@@ -315,7 +315,6 @@ static Ref_t create_detector(Detector& theDetector, xml_h e, SensitiveDetector s
                 }
                 else if (single_elongation >= short_elongation_esr_out.at(0))
                 {
-//                    scintillator_unit = slice.placeVolume(odd_0_scintillator, transform_scintillator);
                     scintillator_unit = slice.placeVolume(odd_0_scintillator, Transform3D(RotationZ(pi),
                                 Position((-0.5 * (Ncells_phi + 1) + iphi) * wrapped_scintillator_xy - 0.5 * short_elongation_esr_in.at(0),
                                          scintillator_pos_z,
@@ -347,7 +346,6 @@ static Ref_t create_detector(Detector& theDetector, xml_h e, SensitiveDetector s
                 }
                 else if (single_elongation >= short_elongation_esr_out.at(0))
                 {
-//                    scintillator_unit = slice.placeVolume(odd_0_scintillator, position_scintillator);
                     scintillator_unit = slice.placeVolume(odd_0_scintillator,
                                 Position((-0.5 * (Ncells_phi + 1) + iphi) * wrapped_scintillator_xy + 0.5 * short_elongation_esr_in.at(0),
                                          scintillator_pos_z,

Detector/DetCRD/CMakeLists.txt

@@ -15,11 +15,15 @@ gaudi_add_module(DetCRD
                          src/Other/Lumical_v01_geo_beampipe.cpp
                          src/Other/CRDBeamPipe_v01_geo.cpp
                          src/Muon/Muon_Barrel_v01_01.cpp
-                         src/Muon/Muon_Endcap_v01_01.cpp
+                         src/Muon/Muon_Endcap_v01_02.cpp
                          src/Tracker/SiTrackerSkewRing_v01_geo.cpp
                          src/Tracker/SiTrackerStaggeredLadder_v01_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_otkendcap_v01_geo.cpp
 
 		 LINK ${DD4hep_COMPONENT_LIBRARIES}
 )

Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Barrel_Short_v02.xml

@@ -12,20 +12,20 @@
 
         <constant name="crystal_r" value="4.1*cm"/>
         <constant name="crystal_phi" value="1*cm"/>
-        <constant name="crystal_z" value="1*cm"/>
+        <constant name="crystal_z_barrel" value="1*cm"/>
 
-        <constant name="esr_thickness" value="0.1*mm"/>
+        <constant name="esr_thickness_barrel" value="0.1*mm"/>
         <constant name="sipm_r" value="0.8*mm"/>
         <constant name="sipm_phi" value="3*mm"/>
-        <constant name="sipm_z" value="3*mm"/>
-        <constant name="pcb_thickness" value="2.2*mm"/>
-        <constant name="cu_thickness" value="1*mm"/>
-        <constant name="fibre_thickness" value="0.1*mm"/>
+        <constant name="sipm_z_barrel" value="3*mm"/>
+        <constant name="pcb_thickness_barrel" value="2.2*mm"/>
+        <constant name="cu_thickness_barrel" value="1*mm"/>
+        <constant name="fibre_thickness_barrel" value="0.1*mm"/>
 
         <constant name="collection_width" value="300*mm"/>
         <constant name="collection_thickness" value="10*mm"/>
 
-        <constant name="boundary_safety" value="1*nm"/>
+        <constant name="boundary_safety_barrel" value="1*nm"/>
     </define>
 
     <regions>

Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Endcap_Short_v01.xml

@@ -24,16 +24,16 @@
         <!-- CrystalXY = CellXY - 2 * BoundarySafety - 2 * ESRThickness-->
         <!-- CrystalZ = CellZ - 2 * BoundarySafety - 2 * ESRThickness-->
         <constant name="Ncell_xy" value="35"/>
-        <constant name="crystal_z" value="4.1*cm"/>
+        <constant name="crystal_z_endcap" value="4.1*cm"/>
 
-        <constant name="esr_thickness" value="0.1*mm"/>
+        <constant name="esr_thickness_endcap" value="0.1*mm"/>
         <constant name="sipm_x" value="3*mm"/>
         <constant name="sipm_y" value="3*mm"/>
-        <constant name="sipm_z" value="0.8*mm"/>
-        <constant name="pcb_thickness" value="2.2*mm"/>
-        <constant name="cu_thickness" value="1*mm"/>
-        <constant name="fibre_thickness" value="0.1*mm"/>
-        <constant name="boundary_safety" value="1*nm"/>
+        <constant name="sipm_z_endcap" value="0.8*mm"/>
+        <constant name="pcb_thickness_endcap" value="2.2*mm"/>
+        <constant name="cu_thickness_endcap" value="1*mm"/>
+        <constant name="fibre_thickness_endcap" value="0.1*mm"/>
+        <constant name="boundary_safety_endcap" value="1*nm"/>
     </define>
 
     <regions>

Detector/DetCRD/compact/CRD_common_v01/Muon_Barrel_v01_01.xml

@@ -12,39 +12,44 @@
 
   <define>
     <!--Muon Barrel-->
-    <constant name="Muon_barrel_superlayer_num" value="6"/>
-    <constant name="Muon_barrel_strip_num_0" value="40"/>
-    <constant name="Muon_barrel_strip_num_1" value="48"/>
-    <constant name="Muon_barrel_strip_num_2" value="62"/>
-    <constant name="Muon_barrel_strip_num_3" value="74"/>
-    <constant name="Muon_barrel_strip_num_4" value="84"/>
-    <constant name="Muon_barrel_strip_num_5" value="96"/>
-    <constant name="Muon_barrel_strip_num" value="100"/>
+    <constant name="Muon_barrel_superlayer_num" value="8"/>
+    <constant name="Muon_barrel_strip_num_0" value="17"/>
+    <constant name="Muon_barrel_strip_num_1" value="29"/>
+    <constant name="Muon_barrel_strip_num_2" value="41"/>
+    <constant name="Muon_barrel_strip_num_3" value="53"/>
+    <constant name="Muon_barrel_strip_num_4" value="65"/>
+    <constant name="Muon_barrel_strip_num_5" value="77"/>
+    <constant name="Muon_barrel_strip_num_6" value="89"/>
+    <constant name="Muon_barrel_strip_num_7" value="101"/>
+    <constant name="Muon_barrel_strip_num_fixed_0" value="106"/>
+    <constant name="Muon_barrel_strip_num_fixed_1" value="115"/>
 
-    <constant name="Muon_barrel_iron_x1" value="Muon_standard_scale"/>
-    <constant name="Muon_barrel_iron_y" value="Muon_strip_z+2*Muon_strip_surf+Muon_Iron_gap_z"/>
+    <!--constant name="Muon_barrel_iron_x1" value="Muon_standard_scale"/-->
+    <constant name="Muon_barrel_iron_y" value="Muon_total_length"/>
     <constant name="Muon_barrel_iron_z" value="Muon_standard_scale"/>
     <constant name="Muon_barrel_iron_posx" value="-1*Muon_standard_scale"/>
 
     <constant name="Muon_barrel_barrel_y" value="Muon_barrel_iron_y"/>
     <constant name="Muon_barrel_barrel_posy" value="0.5*Muon_barrel_barrel_y"/>
 
-    <constant name="Muon_barrel_superlayer_init" value="-21.5*cm"/>
-    <constant name="Muon_barrel_superlayer_gap" value="12.5*cm"/>
-    <constant name="Muon_barrel_superlayer_air_gap" value="0.5*cm"/>
+    <constant name="Muon_barrel_superlayer_init" value="-48*cm"/>
+    <constant name="Muon_barrel_superlayer_gap" value="14*cm"/>
+    <constant name="Muon_barrel_superlayer_endcap_gap" value="10*cm"/>
+    <constant name="Muon_barrel_superlayer_air_gap" value="1*cm"/>
     <constant name="Muon_barrel_superlayer_aluminum_gap" value="0.5*Muon_barrel_superlayer_air_gap"/>
 
-    <constant name="Muon_barrel_superlayer_y" value="2*Muon_strip_y+2*Muon_barrel_superlayer_air_gap"/>
-    <constant name="Muon_barrel_superlayer_z" value="Muon_strip_z+2*Muon_strip_surf+2*Muon_barrel_superlayer_air_gap"/>
+    <constant name="Muon_barrel_superlayer_y" value="2*Muon_strip_y+Muon_barrel_superlayer_air_gap"/>
+    <!--constant name="Muon_barrel_superlayer_z" value="Muon_strip_z+2*Muon_strip_surf+2*Muon_barrel_superlayer_air_gap"/-->
+
     <!--Checkout-->
-    <constant name="Muon_barrel_inner_radius" value="483*cm"/>
+    <constant name="Muon_barrel_inner_radius" value="4245*mm"/>
     <constant name="Muon_barrel_barrel_num" value="2"/>
     <constant name="Muon_barrel_iron_part_num" value="12"/>
 
   </define>
 
   <detectors>
-    <detector id="201" name="Muon_Barrel_v01_01" type="Muon_Barrel_v01_01" readout="MuonBarrelCollection" vis="WhiteVis">
+    <detector id="201" name="MuonBarrel" type="Muon_Barrel_v01_01" readout="MuonBarrelCollection" vis="WhiteVis">
       <position x="0" y="0" z="0"/>
       <barrel id="Muon_barrel_iron_part_num" name="Muon_barrel_barrel" type="Muon_barrel_barrel" vis="SeeThrough">
         <position x="0" y="Muon_barrel_barrel_posy" z="0"/>
@@ -52,8 +57,8 @@
           <material name="Iron"/>
           <position x="Muon_barrel_iron_posx" y="0"/>
           <dimensions x1="0.5*Muon_barrel_iron_x1" y1="0.5*Muon_barrel_iron_y" y2="0.5*Muon_barrel_iron_y" dz="0.5*Muon_barrel_iron_z"/>
-          <superlayer id="Muon_barrel_superlayer_num" name="Muon_barrel_superlayer" type="Muon_barrel_superlayer" vis="GreenVis" material="Air">
-            <aluminum id="0" name="Muon_barrel_superlayer_aluminum" type="Muon_barrel_superlayer_aluminum" vis="GrayVis" material="Aluminum">
+          <superlayer id="Muon_barrel_superlayer_num" name="Muon_barrel_superlayer" type="Muon_barrel_superlayer" vis="BlueVis" material="Air">
+            <aluminum id="0" name="Muon_barrel_superlayer_aluminum" type="Muon_barrel_superlayer_aluminum" vis="RedVis" material="Aluminum">
               <material name="Aluminum"/>
               <position x="0" y="0" z="0"/>
               <stripe id="0" name="Muon_stripe" type="Muon_stripe" vis="GreenVis" material="Air">

Detector/DetCRD/compact/CRD_common_v01/Muon_Endcap_v01_01.xml

@@ -23,7 +23,7 @@
 
     <constant name="Muon_endcap_iron_gap" value="12.5*cm"/>
     <constant name="Muon_endcap_endcap_z" value="Muon_endcap_iron_gap_num*Muon_endcap_iron_gap+2*Muon_endcap_superlayer_num*Muon_strip_y"/>
-    <constant name="Muon_endcap_endcap_posy" value="Muon_strip_z+2*Muon_strip_surf+Muon_Iron_gap_z+0.5*Muon_endcap_endcap_z"/>
+    <constant name="Muon_endcap_endcap_posy" value="0.5*Muon_total_length+0.5*Muon_endcap_endcap_z"/>
   </define>
 
   <detectors>

Detector/DetCRD/compact/CRD_common_v01/Muon_Endcap_v01_02.xml

+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd>
+ 
+  <info name="Muon_Endcap"
+	title="Test with Two Single Muon Endcap"
+	author="Zibing Bai"
+	url="http://cepcgit.ihep.ac.cn"
+	status="development"
+	version="v01_02">
+    <comment>Test with Two Single Muon Endcap</comment>        
+  </info>
+
+  <define>
+    <!--Muon Endcap-->
+    <constant name="Muon_endcap_part_num" value="4"/>
+    <constant name="Muon_endcap_superlayer_num" value="6"/>
+    <constant name="Muon_endcap_layer_num" value="2"/>
+    <constant name="Muon_endcap_strip_num_1" value="72"/>
+    <constant name="Muon_endcap_strip_num_2" value="128"/>
+    <constant name="Muon_endcap_strip_num_cut_1" value="17"/>
+    <constant name="Muon_endcap_strip_num_cut_2" value="74"/>
+    <constant name="Muon_endcap_endcap_rmin" value="68*cm"/>
+    <constant name="Muon_endcap_length_cut_1" value="288*cm"/>
+    <constant name="Muon_endcap_length_cut_gap" value="296*cm"/>
+    <constant name="Muon_endcap_length_cut_2" value="512*cm"/>
+    <constant name="Muon_endcap_gap" value="5*cm"/>
+
+    <constant name="Muon_endcap_iron_gap" value="14*cm"/>
+    <constant name="Muon_endcap_iron_init" value="-60*cm"/>
+    <constant name="Muon_endcap_endcap_z" value="Muon_standard_scale"/>
+    <constant name="Muon_endcap_endcap_posy" value="0.5*Muon_total_length+0.5*Muon_endcap_endcap_z"/>
+  </define>
+
+  <detectors>
+    <detector id="20211" name="MuonEndcap" type="Muon_Endcap_v01_02" readout="MuonEndcapCollection" vis="GrayVis" material="Iron">
+      <material name="Iron"/>
+      <position x="0" y="Muon_endcap_endcap_posy" z="0"/>
+      <dimensions rmin="Muon_endcap_endcap_rmin" dz="0.5*Muon_endcap_endcap_z"/>
+      <stripe id="0" name="Muon_stripe" type="Muon_stripe" vis="GreenVis" material="Air">
+        <material name="Air"/>
+        <dimensions dx="0.5*Muon_strip_x"  dy="0.5*Muon_strip_y"  dz="0.5*Muon_strip_z+Muon_strip_SiPM_z"/>
+        <component id="0" type="Muon_strip_surface" name="Muon_strip_surface" vis="GreenVis" material="BC420">
+          <position x="0" y="0" z="0"/>
+          <dimensions dx="0.5*Muon_strip_surface_x" dy="0.5*Muon_strip_surface_y" dz="0.5*Muon_strip_surface_z"/>
+          <!--cut name="Muon_strip_cut1" vis="GreenVis" material="Air">
+            <position x="Muon_strip_cut1_posx" y="Muon_strip_cut1_posy" z="Muon_strip_cut1_posz"/>
+            <dimensions dx="0.5*Muon_strip_cut1_x" dy="0.5*Muon_strip_cut1_y" dz="0.5*Muon_strip_cut1_z"/>
+          </cut-->
+        </component>
+        <component id="1" type="Muon_strip_scintillator" name="Muon_strip_scintillator" vis="GreenVis" material="BC420">
+          <position x="0" y="0" z="0"/>
+          <dimensions dx="0.5*Muon_strip_scintillator_x" dy="0.5*Muon_strip_scintillator_y" dz="0.5*Muon_strip_scintillator_z"/>
+          <cut name="Muon_strip_cut3" vis="GreenVis" material="Air">
+            <!--position x="Muon_strip_cut3_posx" y="Muon_strip_cut3_posy" z="Muon_strip_cut3_posz"/-->
+            <position x="0" y="0" z="0"/>
+            <dimensions rmin="0" rmax="Muon_strip_cut3_rmax" dz="0.5*Muon_strip_cut3_z"/>
+            <!--comb name="Muon_strip_cut2" vis="GreenVis" material="Air">
+              <position x="Muon_strip_cut2_posx" y="Muon_strip_cut2_posy" z="Muon_strip_cut2_posz"/>
+              <dimensions dx="0.5*Muon_strip_cut2_x" dy="0.5*Muon_strip_cut2_y" dz="0.5*Muon_strip_cut2_z"/>
+            </comb-->
+          </cut>
+        </component>
+        <fiber id="0" type="Muon_fiber_cladding" name="Muon_fiber_cladding" vis="GreenVis" material="Pethylene1">
+          <position x="0" y="0" z="0"/>
+          <dimensions rmin="Muon_fiber_cladding_rmin" rmax="Muon_fiber_cladding_rmax" dz="0.5*Muon_fiber_cladding_z"/>
+        </fiber>
+        <fiber id="1" type="Muon_fiber_core" name="Muon_fiber_core" vis="GreenVis" material="Pethylene2">
+          <position x="0" y="0" z="0"/>
+          <dimensions rmin="0" rmax="Muon_fiber_core_rmax" dz="0.5*Muon_fiber_core_z"/>
+        </fiber>
+        <component id="2" type="Muon_strip_SiPM" name="Muon_strip_SiPM" vis="GreenVis" material="Air">
+          <position x="Muon_strip_SiPM_posx" y="Muon_strip_SiPM_posy" z="Muon_strip_SiPM_posz"/>
+          <dimensions dx="0.5*Muon_strip_SiPM_x" dy="0.5*Muon_strip_SiPM_y" dz="0.5*Muon_strip_SiPM_z"/>
+        </component>
+      </stripe>
+    </detector>
+  </detectors>
+
+
+  <readouts>
+    <readout name="MuonEndcapCollection">
+      <id>system:5,Env:5,Endcap:2,Superlayer:15:4,Layer:2,Stripe:9,SiPM:2</id>
+      <!--id>Endcap:2,Superlayer:2,Env:2,Layer:2,Stripe:10:3,SiPM:2</id-->
+    </readout>
+  </readouts>
+
+</lccdd>

Detector/DetCRD/compact/CRD_common_v01/OTKBarrel_v01_01.xml

+<lccdd>
+  <info name="OTKBarrel_v01_01"
+        title="CepC OTKBarrel"
+        author="D.Yu, "
+        url="http://cepc.ihep.ac.cn"
+        contact="yudian2002@sjtu.edu.cn"
+        status="developing"
+        version="v01">
+    <comment>CepC vertex detector based on MOST2 project </comment>
+  </info>
+  <define>
+    <constant name="OTKBarrel_total_length" value="2*OTKBarrel_half_length" />
+    <constant name="OTKBarrel_ladder_total_thickness" value="15*mm" />
+    <constant name="OTKBarrel_ladder_total_width" value="160*mm" />
+    <constant name="OTKBarrel_ladder_total_length" value="OTKBarrel_total_length" />
+    <constant name="OTKBarrel_ladder_support_thickness" value="1*mm" />
+    <constant name="OTKBarrel_ladder_support_width" value="160*mm" />
+    <constant name="OTKBarrel_ladder_support_length" value="OTKBarrel_total_length" />
+    <constant name="OTKBarrel_ladder_support_height" value="OTKBarrel_ladder_support_thickness" />
+    <constant name="OTKBarrel_flex_width" value="15*mm" />
+    <constant name="OTKBarrel_sensor_length" value="140*mm" />
+    <constant name="OTKBarrel_sensor_thickness" value="500*um" />
+    <constant name="OTKBarrel_sensor_active_width" value="159.5*mm" />
+    <constant name="OTKBarrel_sensor_dead_width" value="0.5*mm" />
+    <constant name="OTKBarrel_pcb_thickness" value="500*um" />
+    <constant name="OTKBarrel_pcb_width" value="28*mm" />
+    <constant name="OTKBarrel_port_pcb_width" value="28*mm" />
+    <constant name="OTKBarrel_pcb_length" value="10*mm" />
+    <constant name="OTKBarrel_pcb_zgap" value="1*mm" />
+    <constant name="OTKBarrel_asic_thickness" value="500*um" />
+    <constant name="OTKBarrel_asic_width" value="130*mm" />
+    <constant name="OTKBarrel_asic_length" value="5*mm" />
+    <constant name="OTKBarrel_asic_zgap" value="1*mm" />
+    
+  </define>
+
+  <detectors>
+    <detector id="DetID_OTKBarrel" limits="otk_limits" name="OTKBarrel" type="SiTracker_otkbarrel_v01" vis="OTKBarrelVis" readout="OTKBarrelCollection" insideTrackingVolume="true">
+      <envelope>
+      <shape type="BooleanShape" operation="Union" material="Air" >
+        <shape type="Tube" rmin="OTKBarrel_inner_radius" rmax="OTKBarrel_outer_radius" dz="OTKBarrel_half_length" />
+      </shape>
+      </envelope>
+
+      <type_flags type="DetType_TRACKER + DetType_BARREL + DetType_STRIP "/>
+
+      <global sensitive_mat="G4_Si" support_mat="G4_C" sensitive_threshold_KeV="64*keV" ladder_offset="0*mm"/>
+      <display ladder="SeeThrough" support="WhiteVis" flex="VXDFlexVis" sens_env="SeeThrough" sens="GrayVis" deadsensor="GreenVis"
+      pcb="GreenVis" asic="yellowVis"/>
+
+      <layer layer_id="0" ladder_radius="OTKBarrel1_inner_radius" n_ladders="45" ladder_offset="0*mm" ladder_thickness="OTKBarrel_ladder_total_thickness"
+      ladder_width="OTKBarrel_ladder_total_width" ladder_length="OTKBarrel_ladder_total_length">
+        <ladder isDoubleSided="false">
+          <ladderSupport height="OTKBarrel_ladder_support_height" length="OTKBarrel_ladder_support_length" thickness="OTKBarrel_ladder_support_thickness"
+          width="OTKBarrel_ladder_support_width" mat="CarbonFiber"/>
+          <flex n_slices="3">
+            <slice length="OTKBarrel_total_length" thickness="60*um" width="OTKBarrel_flex_width" mat="epoxy"/>
+            <slice length="OTKBarrel_total_length" thickness="74*um" width="OTKBarrel_flex_width" mat="Kapton"/>
+            <slice length="OTKBarrel_total_length" thickness="26.8*um" width="OTKBarrel_flex_width" mat="G4_Al"/>
+          </flex>
+          <sensor n_sensors="42" gap="0*mm" thickness="OTKBarrel_sensor_thickness" length="OTKBarrel_sensor_length" active_width="OTKBarrel_sensor_active_width" sensor_mat="G4_Si"
+          dead_width="OTKBarrel_sensor_dead_width"/>
+          <other pcb_thickness="OTKBarrel_pcb_thickness" pcb_width="OTKBarrel_pcb_width" pcb_length="OTKBarrel_pcb_length" pcb_zgap="OTKBarrel_pcb_zgap"  port_pcb_width="OTKBarrel_port_pcb_width" pcb_mat="epoxy"
+          asic_thickness="OTKBarrel_asic_thickness" asic_width="OTKBarrel_asic_width" asic_length="OTKBarrel_asic_length" asic_zgap="OTKBarrel_asic_zgap" asic_mat="G4_Si"/>
+        </ladder>
+      </layer>
+      <layer layer_id="1" ladder_radius="OTKBarrel2_inner_radius" n_ladders="45" ladder_offset="0*mm" ladder_thickness="OTKBarrel_ladder_total_thickness"
+      ladder_width="OTKBarrel_ladder_total_width" ladder_length="OTKBarrel_ladder_total_length">
+        <ladder isDoubleSided="false">
+          <ladderSupport height="OTKBarrel_ladder_support_height" length="OTKBarrel_ladder_support_length" thickness="OTKBarrel_ladder_support_thickness"
+          width="OTKBarrel_ladder_support_width" mat="CarbonFiber"/>
+          <flex n_slices="3">
+            <slice length="OTKBarrel_total_length" thickness="60*um" width="OTKBarrel_flex_width" mat="epoxy"/>
+            <slice length="OTKBarrel_total_length" thickness="74*um" width="OTKBarrel_flex_width" mat="Kapton"/>
+            <slice length="OTKBarrel_total_length" thickness="26.8*um" width="OTKBarrel_flex_width" mat="G4_Al"/>
+          </flex>
+          <sensor n_sensors="42" gap="0*mm" thickness="OTKBarrel_sensor_thickness" length="OTKBarrel_sensor_length" active_width="OTKBarrel_sensor_active_width" sensor_mat="G4_Si"
+          dead_width="OTKBarrel_sensor_dead_width"/>
+          <other pcb_thickness="OTKBarrel_pcb_thickness" pcb_width="OTKBarrel_pcb_width" pcb_length="OTKBarrel_pcb_length" pcb_zgap="OTKBarrel_pcb_zgap"  port_pcb_width="OTKBarrel_port_pcb_width" pcb_mat="epoxy"
+          asic_thickness="OTKBarrel_asic_thickness" asic_width="OTKBarrel_asic_width" asic_length="OTKBarrel_asic_length" asic_zgap="OTKBarrel_asic_zgap" asic_mat="G4_Si"/>
+        </ladder>
+      </layer>
+    </detector>
+  </detectors>
+  
+  <readouts>
+    <readout name="OTKBarrelCollection">
+      <id>system:5,side:-2,layer:9,module:8,active:8,sensor:8</id>
+    </readout>
+  </readouts>
+</lccdd>
+