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/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

Analysis/TotalInvMass/src/TotalInvMass.cc

@@ -11,13 +11,13 @@
 #include <EVENT/LCFloatVec.h>
 #include <EVENT/LCParameters.h>
 #include <stdexcept>
-#include <TFile.h> 
+#include <TFile.h>
 #include <TTree.h>
 #include <TH1F.h>
 #include <TVector3.h>
 #include <TRandom.h>
-#include <Rtypes.h> 
-#include <sstream>		
+#include <Rtypes.h>
+#include <sstream>
 #include <cmath>
 #include <vector>
 #include <TMath.h>
@@ -47,14 +47,14 @@ TotalInvMass::TotalInvMass(const std::string& name, ISvcLocator* svcLoc)
       _colName ,
       _colName);
     */
-	
+
     /*
       _leptonID = 13;
       registerProcessorParameter( "LeptonIDTag" ,
       "Lepton ID that will be used in this analysis." ,
       _leptonID ,
       _leptonID);
-    */	
+    */
 
 }
 
@@ -74,7 +74,7 @@ StatusCode TotalInvMass::initialize() {
     _outputTree->SetAutoSave(32*1024*1024);  // autosave every 32MB
     _outputTree->Branch("EventNr", &_eventNr, "EventNr/I");
     _outputTree->Branch("Num", &_Num, "Num/I");
-	
+
     _outputTree->Branch("OriQuarkID", &_OriQuarkID, "OriQuarkID/I");
 
     _outputTree->Branch("ISREn", &_ISREn, "ISREn/F");
@@ -88,6 +88,9 @@ StatusCode TotalInvMass::initialize() {
     _outputTree->Branch("N3En", &_N3En, "N3En/F");
     _outputTree->Branch("N3Pt", &_N3Pt, "N3Pt/F");
 
+    _outputTree->Branch("ArborPFO_E", &_ArborPFO_E);
+	_outputTree->Branch("ArborPFO_Charge", &_ArborPFO_Charge);
+
 
     _outputTree->Branch("OriJ1CosTheta", &_OriJ1CosTheta, "OriJ1CosTheta/F");
     _outputTree->Branch("OriJ2CosTheta", &_OriJ2CosTheta, "OriJ2CosTheta/F");
@@ -115,7 +118,7 @@ StatusCode TotalInvMass::initialize() {
 
     _outputTree->Branch("FrPh", FrPh, "FrPh[4]/F");
     _outputTree->Branch("FrNe", FrNe, "FrNe[4]/F");
-	
+
     _outputTree->Branch("KPF", KPF, "KPF[4]/F");
 
     _outputTree->Branch("UdP", UdP, "UdP[4]/F");
@@ -172,13 +175,16 @@ StatusCode TotalInvMass::initialize() {
 
 }
 
-StatusCode TotalInvMass::execute() 
-{		
+StatusCode TotalInvMass::execute()
+{
     info() << "TotalInvMass::executing..." << endmsg;
 
     EVENT::LCEvent* evtP = nullptr;
 
-    // if (evtP) {		
+    MCParticleColHandler m_mcParticle{_MCPCollectionName, Gaudi::DataHandle::Reader, this};
+
+
+    // if (evtP) {
 
     TLorentzVector ArborTotalP(0, 0, 0, 0);
     TLorentzVector ArborChP(0, 0, 0, 0);
@@ -194,6 +200,9 @@ StatusCode TotalInvMass::execute()
     TLorentzVector ArborISR(0, 0, 0, 0);
     TLorentzVector PandoraISR(0, 0, 0, 0);
 
+    _ArborPFO_E.clear();
+    _ArborPFO_Charge.clear();
+
     // TODO: using the event header
     // _eventNr = evtP->getEventNumber();
 
@@ -221,23 +230,23 @@ StatusCode TotalInvMass::execute()
         }
     }
 
-    nCHPFO_a = 0; 
+    nCHPFO_a = 0;
     nCHPFO_p = 0;
-    nNEPFO_a = 0; 
+    nNEPFO_a = 0;
     nNEPFO_p = 0;
     Type = -100;
-    Charge = -100; 
+    Charge = -100;
     Energy = 0;
-    CluEn = 0; 
-    TrkSumEn = 0; 
+    CluEn = 0;
+    TrkSumEn = 0;
 
-    _EcalTotalE = 0; _HcalTotalE = 0; _EcalCluE = 0; _HcalCluE = 0; 
+    _EcalTotalE = 0; _HcalTotalE = 0; _EcalCluE = 0; _HcalCluE = 0;
     _EcalCluE_p = 0; _HcalCluE_p = 0;
     _HcalEn1=0;_HcalEn2=0;_HcalEn3=0;_HcalEn4=0;_HcalEn5=0;
     _EcalEn1=0;_EcalEn2=0;_EcalEn3=0;_EcalEn4=0;_EcalEn5=0;
 
-    _HDPID = -1; 
-    _OriQuarkID = 0; 
+    _HDPID = -1;
+    _OriQuarkID = 0;
     _OQDir = -10;
     _HDir = -10;
     _visE=0;
@@ -323,12 +332,12 @@ StatusCode TotalInvMass::execute()
     try{
         auto MCPCol = m_mcParticle.get();
 
-        TVector3 tmpP; 
+        TVector3 tmpP;
         TVector3 ISRP(0, 0, 0);
         _N3En = 0;
         _N3Pt = 0;
- 
-        int NNeutrinoCount = 0; 
+
+        int NNeutrinoCount = 0;
 
         for (int s0 = 0; s0 < MCPCol->size(); ++s0) {
             auto MCP = (*MCPCol)[s0];
@@ -340,10 +349,20 @@ StatusCode TotalInvMass::execute()
             TVector3 VTX(VTX0.x, VTX0.y, VTX0.z);
             TVector3 EndP(EndP0.x, EndP0.y, EndP0.z);
 
+            int GenStatus = MCP.getGeneratorStatus();
+            int SimuStatus = MCP.getSimulatorStatus();
+
+            if(0) debug()<<"MCP "<<s0<<": tmpPID = "<<tmpPID<<", NParent = "<<NParent<<", NDaughter = "<<NDaughter<<", GenStatus = "<<GenStatus<<", SimuStatus = "<<SimuStatus<<endmsg;
+
+            // auto aDau = MCP.getDaughters(0);
+            // if(aDau){
+            //     debug()<<"Dau 0: PDG = "<<aDau->getPDG()<<endmsg;
+
+
             if(tmpPID == 22 && NParent == 0 && s0 < 4) {
                 auto tmpP0 = MCP.getMomentum();
                 tmpP = TVector3(tmpP0.x, tmpP0.y, tmpP0.z);
-                ISRP += tmpP; 
+                ISRP += tmpP;
             }
 
             if( (abs(tmpPID) == 12 || abs(tmpPID) == 14 || abs(tmpPID) == 16) && NParent != 0 && NDaughter == 0 && VTX.Mag() < 100 && EndP.Mag() > 100) {
@@ -355,20 +374,25 @@ StatusCode TotalInvMass::execute()
                     tmpP = TVector3(tmpP0.x, tmpP0.y, tmpP0.z);
                     _N3En += tmpP.Mag();
                     _N3Pt += tmpP.Perp();
-                    cout<<"Found Neutrino: "<<NNeutrinoCount<<" En "<<_N3En<<" Pt "<<_N3Pt<<endl;
+                    debug()<<"Found Neutrino: "<<NNeutrinoCount<<" En "<<_N3En<<" Pt "<<_N3Pt<<endmsg;
                 }
             }
 
             if(tmpPID == 25 && NDaughter > 1 && NParent !=0 ) { //Higgs
-                // cout<<"tmpPID:HDPID"<<tmpPID<<" "<<NDaughter<<" "<<NParent<<endl;
+                // debug()<<"tmpPID:HDPID"<<tmpPID<<" "<<NDaughter<<" "<<NParent<<endmsg;
                 _HDPID = abs(MCP.getDaughters(0).getPDG());
                 _HDir = MCP.getMomentum()[2]/MCP.getEnergy();
 
+                // debug()<<"This is Higgs, has "<<NDaughter<<" daughters"<<endmsg;
+                debug()<<"This is Higgs, has "<<NDaughter<<" daughters"<<endmsg;
+
                 if(NDaughter == 2) {
                     auto D1 = MCP.getDaughters(0);
                     _J1CosTheta = D1.getMomentum()[2]/D1.getEnergy();
                     auto D2 = MCP.getDaughters(1);
                     _J2CosTheta = D2.getMomentum()[2]/D2.getEnergy();
+
+                    debug()<<"---> PDG: "<<D1.getPDG()<<" + "<<D2.getPDG()<<endmsg;
                 }
             }
             if(abs(tmpPID)<7 && NParent == 0) {
@@ -379,9 +403,9 @@ StatusCode TotalInvMass::execute()
             }
 
             if(abs(tmpPID)<7 && NParent == 0) {
-                _OriQuarkID = abs(tmpPID); 
+                _OriQuarkID = abs(tmpPID);
                 _OQDir = MCP.getMomentum()[2]/MCP.getEnergy();
-            }	
+            }
 
 
             if(abs(_J1CosTheta) > abs(_J2CosTheta))
@@ -413,17 +437,21 @@ StatusCode TotalInvMass::execute()
             auto a_RecoP = (*col_RecoNeP)[i0];
             TLorentzVector currP( a_RecoP.getMomentum()[0], a_RecoP.getMomentum()[1], a_RecoP.getMomentum()[2], a_RecoP.getEnergy());
             ArborTotalP += currP;
+
+            _ArborPFO_E.push_back(a_RecoP.getEnergy());
+            _ArborPFO_Charge.push_back(a_RecoP.getCharge());
+
             auto currMom0 = a_RecoP.getMomentum();
             TVector3 currMom(currMom0.x, currMom0.y, currMom0.z);
             //				if(a_RecoP->getType() == 22 && a_RecoP->getEnergy() > 2)	//Compensate...
             //					ArborTotalP += 0.98*currP;
-            //				else 
+            //				else
             //					ArborTotalP += currP;
 
             if(a_RecoP.getCharge() != 0) {
                 ArborChP += currP;
             } else if(a_RecoP.getType() == 310) {
-                ArborKPF += currP; 
+                ArborKPF += currP;
             } else if(a_RecoP.getType() == 22) {
                 if(a_RecoP.getEnergy() < 3.0)
                     ArborFrPh += currP;
@@ -436,21 +464,36 @@ StatusCode TotalInvMass::execute()
                     ArborNeP += currP;
             } else if(a_RecoP.getEnergy() < 3.0) {
                 ArborFrP += currP;
-                cout<<"Undef "<<a_RecoP.getType()<<endl;  
+                if(0) debug()<<"Undef "<<a_RecoP.getType()<<endmsg;
             } else {
                 ArborUdP += currP;
-                cout<<"Undef "<<a_RecoP.getType() << "En "<<a_RecoP.getEnergy()<<endl;
+                if(0) debug()<<"Undef "<<a_RecoP.getType() << "En "<<a_RecoP.getEnergy()<<endmsg;
             }
 
+
+            if(0) debug()<<"[YX debug - BMR] PFO "<<i0<<", E = "<<a_RecoP.getEnergy()<<", Charge = "<<a_RecoP.getCharge()<<endmsg;
+
+
             if(a_RecoP.clusters_size() > 0) {
+
+                int nTrkHit = 0;
+                if(a_RecoP.getCharge()){
+                    auto a_Trk = a_RecoP.getTracks(0);
+                    nTrkHit = a_Trk.trackerHits_size();
+                }
+
                 auto a_clu = a_RecoP.getClusters(0);
                 auto CluPos0 = a_clu.getPosition();
                 TVector3 CluPos(CluPos0.x, CluPos0.y, CluPos0.z);
                 if( CluPos.Perp() < 300 && abs(CluPos.Z()) < 1300 && a_RecoP.getCharge() == 0 )	// 1150-1300
-                    ArborLCAL += currP; 
+                    ArborLCAL += currP;
+
+
+                if(0) debug()<<"[YX debug - BMR] ---> nTrkHit = "<<nTrkHit<<", CluE = "<<a_clu.getEnergy()<<", CluPos = ("<<CluPos.X()<<", "<<CluPos.Y()<<", "<<CluPos.Z()<<")"<<endmsg;
+
 
                 float MinAngleToCH = 1.0E6;
-                float MinAngleToNE = 1.0E6;  
+                float MinAngleToNE = 1.0E6;
 
                 if(a_RecoP.getEnergy() > 5 && a_RecoP.getCharge() == 0) {
                     for(int i1 = 0; i1 < col_RecoNeP->size(); i1++) {
@@ -464,7 +507,7 @@ StatusCode TotalInvMass::execute()
                                     if(tmpAngle < MinAngleToCH) {
                                         MinAngleToCH = tmpAngle;
                                     }
-                                } else {	
+                                } else {
                                     if(tmpAngle < MinAngleToNE) {
                                         MinAngleToNE = tmpAngle;
                                     }
@@ -474,19 +517,19 @@ StatusCode TotalInvMass::execute()
                     }
 
                     if( MinAngleToNE > 0.5 || MinAngleToCH > 0.5 ) {
-                        ArborISR += currP; 
+                        ArborISR += currP;
                     }
                 }
             }
 
-            TrackHit = 0; 
+            TrackHit = 0;
             for(int k = 0; k < 3; k++) {
                 StartPos[k] = 0;
                 EndPos[k] = 0;
-            }		
+            }
 
             Charge = int(a_RecoP.getCharge());
-            CluEn = 0; 
+            CluEn = 0;
             CluEnCom[0] = 0;
             CluEnCom[1] = 0;
 
@@ -510,7 +553,7 @@ StatusCode TotalInvMass::execute()
 
                 if((!Charge) and (currClu.subdetectorEnergies_size() > 2)) {
                     CluEnCom[0] = currClu.getSubdetectorEnergies(0);
-                    CluEnCom[1] = currClu.getSubdetectorEnergies(1);	
+                    CluEnCom[1] = currClu.getSubdetectorEnergies(1);
                     NeCaloE_a[0] += currClu.getSubdetectorEnergies(0);
                     NeCaloE_a[1] += currClu.getSubdetectorEnergies(1);
                 }
@@ -522,7 +565,7 @@ StatusCode TotalInvMass::execute()
 
 
                 if(Charge) {
-                    TrkSumEn += Energy; 
+                    TrkSumEn += Energy;
 
                     if (a_RecoP.tracks_size()>0) {
                         auto a_Trk = a_RecoP.getTracks(0);
@@ -535,19 +578,19 @@ StatusCode TotalInvMass::execute()
 
                             EndPos[0] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[0];
                             EndPos[1] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[1];
-                            EndPos[2] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[2];	
+                            EndPos[2] = (a_Trk.getTrackerHits(TrackHit - 1)).getPosition()[2];
                         }
                     }
 
                     if( Energy > CluEn + sqrt(Energy)) {
-                        ElargeP[0] += Energy; 
-                        ElargeP[1] += CluEn; 
+                        ElargeP[0] += Energy;
+                        ElargeP[1] += CluEn;
                     } else if( fabs(Energy - CluEn) < sqrt(Energy) ) {
                         EequP[0] += Energy;
-                        EequP[1] += CluEn; 
+                        EequP[1] += CluEn;
                     } else {
                         EsmallP[0] += Energy;
-                        EsmallP[1] += CluEn; 
+                        EsmallP[1] += CluEn;
                     }
 
                 }
@@ -559,18 +602,18 @@ StatusCode TotalInvMass::execute()
     }catch (lcio::DataNotAvailableException err) { }
 
     try{
-        //LCCollection* col_RecoPandora = evtP->getCollection( "PandoraPFOs" );			
+        //LCCollection* col_RecoPandora = evtP->getCollection( "PandoraPFOs" );
         //for(int i2 = 0; i2 < col_RecoPandora->getNumberOfElements(); i2++)
         for(int s = 0; s < 1; s++) {
             auto col_PFO_iter = m_arbopfo.get();
-            /* 
+            /*
                if(s==0)
                col_PFO_iter = evtP->getCollection( "ArborChargedCore" );
                else
-               col_PFO_iter = evtP->getCollection( "ArborNeutralCore" );	
+               col_PFO_iter = evtP->getCollection( "ArborNeutralCore" );
             */
 
-            for(int i2 = 0; i2 < col_PFO_iter->size(); i2++) {         
+            for(int i2 = 0; i2 < col_PFO_iter->size(); i2++) {
                 auto a_RecoP = (*col_PFO_iter)[i2];
                 TLorentzVector currP( a_RecoP.getMomentum()[0], a_RecoP.getMomentum()[1], a_RecoP.getMomentum()[2], a_RecoP.getEnergy());
                 PandoraTotalP += currP;
@@ -580,7 +623,7 @@ StatusCode TotalInvMass::execute()
                 } else {
                     nNEPFO_p++;
                 }
-                        
+
                 auto currMom0 = a_RecoP.getMomentum();
                 TVector3 currMom(currMom0.x, currMom0.y, currMom0.z);
 
@@ -592,7 +635,7 @@ StatusCode TotalInvMass::execute()
                     auto currClu = a_RecoP.getClusters(0);
                     CluEn = currClu.getEnergy();
 
-                    _EcalCluE_p += CluEn;	
+                    _EcalCluE_p += CluEn;
                     _HcalCluE_p += currClu.getSubdetectorEnergies(1);
 
                     if(a_RecoP.getEnergy() > 5 && a_RecoP.getCharge() == 0) {
@@ -625,11 +668,11 @@ StatusCode TotalInvMass::execute()
         }
     }catch (lcio::DataNotAvailableException err) { }
 
-    _Mass_a = 0; 
-    _Mass_p = 0; 
+    _Mass_a = 0;
+    _Mass_p = 0;
     _Mass_a_Pisr = 0;
-    _Mass_a_Plcal = 0; 
-    _Mass_p_Pisr = 0; 
+    _Mass_a_Plcal = 0;
+    _Mass_p_Pisr = 0;
 
     _Mass_a = ArborTotalP.M();
     _Mass_p = PandoraTotalP.M();
@@ -656,7 +699,7 @@ StatusCode TotalInvMass::execute()
     PhP[0] = ArborPhP.X();
     PhP[1] = ArborPhP.Y();
     PhP[2] = ArborPhP.Z();
-    PhP[3] = ArborPhP.T();	
+    PhP[3] = ArborPhP.T();
 
     NeP[0] = ArborNeP.X();
     NeP[1] = ArborNeP.Y();
@@ -688,18 +731,18 @@ StatusCode TotalInvMass::execute()
     KPF[2] = ArborKPF.Z();
     KPF[3] = ArborKPF.T();
 
-    cout<<_Mass_a<<" : "<<_Mass_p<<endl;		
+    info()<<_Mass_a<<" : "<<_Mass_p<<endmsg;
 
     _outputTree->Fill();
     _Num++;
-    // }  	  
+    // }
 
 
     info() << "TotalInvMass::execute done" << endmsg;
 
     return StatusCode::SUCCESS;
 
-}	
+}
 
 StatusCode TotalInvMass::finalize()
 {

Analysis/TotalInvMass/src/TotalInvMass.hh

@@ -34,7 +34,11 @@ public:
 
 protected:
     typedef DataHandle<edm4hep::MCParticleCollection> MCParticleColHandler;
-    MCParticleColHandler m_mcParticle{"MCParticle", Gaudi::DataHandle::Reader, this};
+    // MCParticleColHandler m_mcParticle{"MCParticle", Gaudi::DataHandle::Reader, this};
+    // MCParticleColHandler m_mcParticle{"MCParticleGen", Gaudi::DataHandle::Reader, this};
+    Gaudi::Property<std::string> _MCPCollectionName{this,
+            "MCPCollectionName", "MCParticle",
+            "MCPCollectionName"};
 
     typedef DataHandle<edm4hep::CalorimeterHitCollection> CaloHitColHandler;
     CaloHitColHandler m_ecalbarrelhitcol{"ECALBarrel", Gaudi::DataHandle::Reader, this};
@@ -45,8 +49,9 @@ protected:
     CaloHitColHandler m_hcalotherhitcol {"HCALOther", Gaudi::DataHandle::Reader, this};
 
     typedef DataHandle<edm4hep::ReconstructedParticleCollection> RecParticleColHandler;
-    RecParticleColHandler m_reconep{"AncientPFOs", Gaudi::DataHandle::Reader, this};
-    RecParticleColHandler m_arbopfo{"ArborLICHPFOs", Gaudi::DataHandle::Reader, this};
+    RecParticleColHandler m_reconep{"ArborPFO", Gaudi::DataHandle::Reader, this};
+    RecParticleColHandler m_arbopfo{"ArborPFO", Gaudi::DataHandle::Reader, this};
+
 
     Gaudi::Property<std::string> _treeFileName{this,
             "TreeOutputFile", "MCTruth.root",
@@ -61,36 +66,36 @@ protected:
             "OverwriteFile", 0,
             "If zero an already existing file will not be overwritten."};
     int _leptonID;
-    float _cmsE; 
+    float _cmsE;
     TTree *_outputTree, *_outputPFO;
     float _ISRP[3];
-    float _ISREn, _ISRPt; 
-    float _N3En, _N3Pt; 
-    float _NEn, _NPt; 
+    float _ISREn, _ISRPt;
+    float _N3En, _N3Pt;
+    float _NEn, _NPt;
 
-    int _HDPID, _OriQuarkID; 
+    int _HDPID, _OriQuarkID;
     float _OQDir, _HDir;
     int _NMuP, _NMuM, _NChP, _NChM;
     float _P_MuP[4], _P_MuM[4], _P_DL[4];
-    int _EventType; 
-    float _InvMass, _RecoilMass; 
-    float _J1CosTheta, _J2CosTheta, _MaxJetCosTheta; 
-    float _OriJ1CosTheta, _OriJ2CosTheta, _MaxOriJetCosTheta; 
-    float _Mass_a, _Mass_p; 
+    int _EventType;
+    float _InvMass, _RecoilMass;
+    float _J1CosTheta, _J2CosTheta, _MaxJetCosTheta;
+    float _OriJ1CosTheta, _OriJ2CosTheta, _MaxOriJetCosTheta;
+    float _Mass_a, _Mass_p;
 
     int _PID1, _PID2;
     float _PL1[4], _PL2[4], _RPL1[4], _RPL2[4], _SM[4], _P_allCharged[4], _P_allNeutral[4], _P_Higgs[4], _P_allReco[4];
-    float _Hmass; 
+    float _Hmass;
     int _Num;
-    int _NHDaug; 
-    int _HdaughterPID; 
-    int _ZdaughterPID; 
+    int _NHDaug;
+    int _HdaughterPID;
+    int _ZdaughterPID;
     float _Pz[4], _Ph[4], _PzD1[4], _PzD2[4], _PhD1[4], _PhD2[4], _RPzD1[4], _RPzD2[4], _RPhD1[4], _RPhD2[4];
     float _P[4], _SumP[4], _VisP[4], _MissP[4];
-    int _PID, _NFMCP, _MotherFlag, _NNeutrino; 
-    float _ENeutrino, _DiPhMass, _DiPhMassCorr; 
+    int _PID, _NFMCP, _MotherFlag, _NNeutrino;
+    float _ENeutrino, _DiPhMass, _DiPhMassCorr;
     //		float _CosTheta, _Phi, _Charge;
-    float _Mz, _Mrecoil, _MzReco, _MhReco, _MrecoilReco; 
+    float _Mz, _Mrecoil, _MzReco, _MhReco, _MrecoilReco;
     float KthEn[7][9];
     unsigned int _eventNr;
 
@@ -102,12 +107,12 @@ protected:
     float ElargeP[2], EequP[2], EsmallP[2];
 
     float ChP[4], FrP[4], PhP[4], NeP[4], UdP[4], FrPh[4], FrNe[4], KPF[4];
-    float _EcalTotalE, _HcalTotalE, _EcalCluE, _HcalCluE, _EcalCluE_p, _HcalCluE_p; 
+    float _EcalTotalE, _HcalTotalE, _EcalCluE, _HcalCluE, _EcalCluE_p, _HcalCluE_p;
     float _HcalEn1, _HcalEn2,_HcalEn3,_HcalEn4,_HcalEn5;
     float _EcalEn1, _EcalEn2,_EcalEn3,_EcalEn4,_EcalEn5;
 
     int Type, Charge;
-    float Energy, TrkSumEn; 
+    float Energy, TrkSumEn;
     float P[3], CluEnCom[2];
     float CluEn;
 
@@ -115,6 +120,17 @@ protected:
     float StartPos[3], EndPos[3];
     float _visE;
 
+    std::vector<int> _ArborPFO_Charge;
+    std::vector<double> _ArborPFO_E;
+    // std::vector<double> _ArborPFO_CosTheta;
+    // std::vector<int> _ArborPFO_nTrk;
+    // std::vector<int> _ArborPFO_nClu;
+    // std::vector<TVector3> _ArborPFO_TrkP3;
+    // std::vector<double> _ArborPFO_CluE;
+    // std::vector<double> _ArborPFO_CluHitESum;
+    // std::vector<double> _ArborPFO_CluHitESum;
+
+
     std::string _fileName;
     std::ostream *_output;
     std::string _histFileName;

CMakeLists.txt

@@ -6,11 +6,6 @@ project(CEPCSW)
 include(cmake/CEPCSWEnv.cmake)
 #
 
-find_package(ROOT COMPONENTS RIO Tree)
-
-find_package(Gaudi)
-
-
 include(GNUInstallDirs)
 include(CTest)
 
@@ -19,6 +14,17 @@ if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
     "Install path prefix, prepended onto install directories." FORCE )
 endif()
 
+set(CMAKE_SKIP_BUILD_RPATH TRUE)
+set(CMAKE_EXECUTABLE_SUFFIX .exe)
+
+# Put bin/lib/include under CMAKE_BINARY_DIR
+set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
+include_directories(${CMAKE_BINARY_DIR}/include)
+
+# set Install/lib
+set(CMAKE_INSTALL_LIBDIR lib)
+
 # Set up C++ Standard
 # ``-DCMAKE_CXX_STANDARD=<standard>`` when invoking CMake
 set(CMAKE_CXX_STANDARD 17 CACHE STRING "")
@@ -27,13 +33,28 @@ 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)
 include(cmake/CEPCSWDependencies.cmake)
 
 add_subdirectory(Analysis)
 add_subdirectory(Detector)
-add_subdirectory(Digitisers)
+add_subdirectory(Digitization)
 add_subdirectory(Examples)
 add_subdirectory(Generator)
 add_subdirectory(Reconstruction)

Detector/CMakeLists.txt

@@ -5,6 +5,6 @@ add_subdirectory(DetDriftChamber)
 add_subdirectory(DetEcalMatrix)
 add_subdirectory(DetInterface)
 add_subdirectory(DetSegmentation)
-add_subdirectory(GeomSvc)
-add_subdirectory(Identifier)
+add_subdirectory(DetGeomSvc)
+add_subdirectory(DetIdentifier)
 add_subdirectory(MagneticFieldMap)

Detector/DetCEPCv4/CMakeLists.txt

@@ -4,18 +4,6 @@
 # Based on package: lcgeo
 ################################################################################
 
-find_package(DD4hep COMPONENTS DDRec DDG4 DDParsers REQUIRED)
-# find_package(DD4hep)
-find_package(Geant4)
-include(${Geant4_USE_FILE})
-
-set(CMAKE_MODULE_PATH  ${CMAKE_MODULE_PATH}  ${DD4hep_ROOT}/cmake )
-include( DD4hep )
-
-find_package(ROOT COMPONENTS MathCore GenVector Geom REQUIRED)
-
-# install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/compact DESTINATION Detector/DetCEPCv4)
-
 gaudi_add_module(DetCEPCv4
                  SOURCES src/tracker/VXD04_geo.cpp
                          src/tracker/FTD_Simple_Staggered_geo.cpp
@@ -34,6 +22,7 @@ gaudi_add_module(DetCEPCv4
                          src/calorimeter/SHcalRpc01_EndcapRing.cpp
                          src/calorimeter/SHcalSc04_Barrel_v04.cpp
                          src/calorimeter/SHcalSc04_Endcaps_v01.cpp
+                         src/calorimeter/SHcalSc04_Endcaps_v02.cpp
                          src/calorimeter/Yoke05_Barrel.cpp
                          src/calorimeter/Yoke05_Endcaps.cpp
                          src/other/BoxSupport_o1_v01_geo.cpp

Detector/DetCEPCv4/compact/tpc10_01.xml

@@ -6,7 +6,7 @@
 
   <detectors>
 
-    <detector name="TPC" type="TPC10" vis="TPCVis" id="ILDDetID_TPC" limits="Tracker_limits" readout="TPCCollection" insideTrackingVolume="true">
+    <detector name="TPC" type="TPC10" vis="TPCVis" id="ILDDetID_TPC" limits="TPC_limits" readout="TPCCollection" insideTrackingVolume="true">
 
 
       <envelope vis="ILD_TPCVis">
@@ -86,7 +86,7 @@
 
  <readouts>
     <readout name="TPCCollection">
-      <id>system:5,side:-2,layer:9,module:8,sensor:8</id>
+      <id>system:5,side:-2,layer:13,module:8,sensor:8</id>
     </readout>
  </readouts>
 

Detector/DetCEPCv4/src/tracker/FTD_Simple_Staggered_geo.cpp

@@ -327,7 +327,13 @@ static Ref_t create_element(Detector& theDetector, xml_h e, SensitiveDetector se
   
   PlacedVolume pv;
 
-  sens.setType("tracker");
+  if (x_det.hasChild(_U(sensitive))) {
+    xml_dim_t sd_typ = x_det.child(_U(sensitive));
+    sens.setType(sd_typ.typeStr());
+  }
+  else {
+    sens.setType("tracker");
+  }
 
   // --- create assembly and DetElement for support and service volumes 
 
@@ -1016,7 +1022,11 @@ static Ref_t create_element(Detector& theDetector, xml_h e, SensitiveDetector se
     petalSupport(theDetector, ftd, valuesDict, FTDPetalAirLogical ) ; 
 
     VolVec volV = petalSensor( theDetector, ftd, sens, valuesDict, FTDPetalAirLogical );
-
+    if (x_det.hasAttr(_U(limits))) {
+      for (auto vol : volV) {
+        vol.first.setLimitSet(theDetector, x_det.limitsStr());
+      }
+    }
 
     //---- meassurement surface vectors 
 

Detector/DetCEPCv4/src/tracker/FTD_cepc_geo.cpp

@@ -327,7 +327,13 @@ static Ref_t create_element(Detector& theDetector, xml_h e, SensitiveDetector se
   
   PlacedVolume pv;
 
-  sens.setType("tracker");
+  if (x_det.hasChild(_U(sensitive))) {
+    xml_dim_t sd_typ = x_det.child(_U(sensitive));
+    sens.setType(sd_typ.typeStr());
+  }
+  else {
+    sens.setType("tracker");
+  }
 
   // --- create assembly and DetElement for support and service volumes 
 
@@ -1018,7 +1024,11 @@ static Ref_t create_element(Detector& theDetector, xml_h e, SensitiveDetector se
     petalSupport(theDetector, ftd, valuesDict, FTDPetalAirLogical ) ; 
 
     VolVec volV = petalSensor( theDetector, ftd, sens, valuesDict, FTDPetalAirLogical );
-
+    if (x_det.hasAttr(_U(limits))) {
+      for (auto vol : volV) {
+	vol.first.setLimitSet(theDetector, x_det.limitsStr());
+      }
+    }
 
     //---- meassurement surface vectors