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.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("TPCTracksAssociation", m_TPCTrkAssoHdl, "TPC Track - MCParticle Collection");
+  declareProperty("FullTracksAssociation", m_fullTrkAssoHdl, "Full Track - MCParticle 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_sitrktree = new TTree("RecoSiTrk", "RecoSiTrk");
+  m_tpctrktree = new TTree("RecoTPCTrk", "RecoTPCTrk");
+  m_fulltrktree = new TTree("RecoFullTrk", "RecoFullTrk");
+
+  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_VTX_x", &MCP_VTX_x);
+  m_mctree->Branch("MCP_VTX_y", &MCP_VTX_y);
+  m_mctree->Branch("MCP_VTX_z", &MCP_VTX_z);
+  m_mctree->Branch("MCP_endPoint_x", &MCP_endPoint_x);
+  m_mctree->Branch("MCP_endPoint_y", &MCP_endPoint_y);
+  m_mctree->Branch("MCP_endPoint_z", &MCP_endPoint_z);
+  m_mctree->Branch("MCP_pdgid", &MCP_pdgid);
+  m_mctree->Branch("MCP_gStatus", &MCP_gStatus);
+
+  //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_sitrktree->Branch("N_SiTrk", &N_SiTrk);
+  m_sitrktree->Branch("N_TPCTrk", &N_TPCTrk);
+  m_sitrktree->Branch("N_fullTrk", &N_fullTrk);
+  m_sitrktree->Branch("trk_px", &m_trk_px);
+  m_sitrktree->Branch("trk_py", &m_trk_py);
+  m_sitrktree->Branch("trk_pz", &m_trk_pz);
+  m_sitrktree->Branch("trk_p", &m_trk_p);
+  m_sitrktree->Branch("trk_Nhit", &m_trk_Nhit);
+  m_sitrktree->Branch("trkstate_d0", &m_trkstate_d0 );
+  m_sitrktree->Branch("trkstate_z0", &m_trkstate_z0 );
+  m_sitrktree->Branch("trkstate_phi", &m_trkstate_phi );
+  m_sitrktree->Branch("trkstate_tanL", &m_trkstate_tanL );
+  m_sitrktree->Branch("trkstate_omega", &m_trkstate_omega );
+  m_sitrktree->Branch("trkstate_kappa", &m_trkstate_kappa );
+  m_sitrktree->Branch("trkstate_refx", &m_trkstate_refx );
+  m_sitrktree->Branch("trkstate_refy", &m_trkstate_refy );
+  m_sitrktree->Branch("trkstate_refz", &m_trkstate_refz );
+  m_sitrktree->Branch("trkstate_tag", &m_trkstate_tag );
+  m_sitrktree->Branch("trkstate_location", &m_trkstate_location );
+  m_sitrktree->Branch("truthMC_tag", &m_truthMC_tag);
+  m_sitrktree->Branch("truthMC_pid", &m_truthMC_pid);
+  m_sitrktree->Branch("truthMC_px", &m_truthMC_px);
+  m_sitrktree->Branch("truthMC_py", &m_truthMC_py);
+  m_sitrktree->Branch("truthMC_pz", &m_truthMC_pz);
+  m_sitrktree->Branch("truthMC_E", &m_truthMC_E);
+  m_sitrktree->Branch("truthMC_EPx", &m_truthMC_EPx);
+  m_sitrktree->Branch("truthMC_EPy", &m_truthMC_EPy);
+  m_sitrktree->Branch("truthMC_EPz", &m_truthMC_EPz);
+  m_sitrktree->Branch("truthMC_weight", &m_truthMC_weight);
+
+  m_tpctrktree->Branch("N_SiTrk", &N_SiTrk);
+  m_tpctrktree->Branch("N_TPCTrk", &N_TPCTrk);
+  m_tpctrktree->Branch("N_fullTrk", &N_fullTrk);
+  m_tpctrktree->Branch("trk_px", &m_trk_px);
+  m_tpctrktree->Branch("trk_py", &m_trk_py);
+  m_tpctrktree->Branch("trk_pz", &m_trk_pz);
+  m_tpctrktree->Branch("trk_p", &m_trk_p);
+  m_tpctrktree->Branch("trk_Nhit", &m_trk_Nhit);
+  m_tpctrktree->Branch("trkstate_d0", &m_trkstate_d0 );
+  m_tpctrktree->Branch("trkstate_z0", &m_trkstate_z0 );
+  m_tpctrktree->Branch("trkstate_phi", &m_trkstate_phi );
+  m_tpctrktree->Branch("trkstate_tanL", &m_trkstate_tanL );
+  m_tpctrktree->Branch("trkstate_omega", &m_trkstate_omega );
+  m_tpctrktree->Branch("trkstate_kappa", &m_trkstate_kappa );
+  m_tpctrktree->Branch("trkstate_refx", &m_trkstate_refx );
+  m_tpctrktree->Branch("trkstate_refy", &m_trkstate_refy );
+  m_tpctrktree->Branch("trkstate_refz", &m_trkstate_refz );
+  m_tpctrktree->Branch("trkstate_tag", &m_trkstate_tag );
+  m_tpctrktree->Branch("trkstate_location", &m_trkstate_location );
+  m_tpctrktree->Branch("truthMC_tag", &m_truthMC_tag);
+  m_tpctrktree->Branch("truthMC_pid", &m_truthMC_pid);
+  m_tpctrktree->Branch("truthMC_px", &m_truthMC_px);
+  m_tpctrktree->Branch("truthMC_py", &m_truthMC_py);
+  m_tpctrktree->Branch("truthMC_pz", &m_truthMC_pz);
+  m_tpctrktree->Branch("truthMC_E", &m_truthMC_E);
+  m_tpctrktree->Branch("truthMC_EPx", &m_truthMC_EPx);
+  m_tpctrktree->Branch("truthMC_EPy", &m_truthMC_EPy);
+  m_tpctrktree->Branch("truthMC_EPz", &m_truthMC_EPz);
+  m_tpctrktree->Branch("truthMC_weight", &m_truthMC_weight);
+
+  m_fulltrktree->Branch("N_SiTrk", &N_SiTrk);
+  m_fulltrktree->Branch("N_TPCTrk", &N_TPCTrk);
+  m_fulltrktree->Branch("N_fullTrk", &N_fullTrk);
+  m_fulltrktree->Branch("trk_px", &m_trk_px);
+  m_fulltrktree->Branch("trk_py", &m_trk_py);
+  m_fulltrktree->Branch("trk_pz", &m_trk_pz);
+  m_fulltrktree->Branch("trk_p", &m_trk_p);
+  m_fulltrktree->Branch("trk_Nhit", &m_trk_Nhit);
+  m_fulltrktree->Branch("trk_type", &m_trk_type);
+  m_fulltrktree->Branch("trkstate_d0", &m_trkstate_d0 );
+  m_fulltrktree->Branch("trkstate_z0", &m_trkstate_z0 );
+  m_fulltrktree->Branch("trkstate_phi", &m_trkstate_phi );
+  m_fulltrktree->Branch("trkstate_tanL", &m_trkstate_tanL );
+  m_fulltrktree->Branch("trkstate_omega", &m_trkstate_omega );
+  m_fulltrktree->Branch("trkstate_kappa", &m_trkstate_kappa );
+  m_fulltrktree->Branch("trkstate_refx", &m_trkstate_refx );
+  m_fulltrktree->Branch("trkstate_refy", &m_trkstate_refy );
+  m_fulltrktree->Branch("trkstate_refz", &m_trkstate_refz );
+  m_fulltrktree->Branch("trkstate_tag", &m_trkstate_tag );
+  m_fulltrktree->Branch("trkstate_location", &m_trkstate_location );
+  m_fulltrktree->Branch("truthMC_tag", &m_truthMC_tag);
+  m_fulltrktree->Branch("truthMC_pid", &m_truthMC_pid);
+  m_fulltrktree->Branch("truthMC_px", &m_truthMC_px);
+  m_fulltrktree->Branch("truthMC_py", &m_truthMC_py);
+  m_fulltrktree->Branch("truthMC_pz", &m_truthMC_pz);
+  m_fulltrktree->Branch("truthMC_E", &m_truthMC_E);
+  m_fulltrktree->Branch("truthMC_EPx", &m_truthMC_EPx);
+  m_fulltrktree->Branch("truthMC_EPy", &m_truthMC_EPy);
+  m_fulltrktree->Branch("truthMC_EPz", &m_truthMC_EPz);
+  m_fulltrktree->Branch("truthMC_weight", &m_truthMC_weight);
+
+  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_VTX_x.push_back(m_MCp.getVertex().x);
+      MCP_VTX_y.push_back(m_MCp.getVertex().y);
+      MCP_VTX_z.push_back(m_MCp.getVertex().z);
+      MCP_endPoint_x.push_back(m_MCp.getEndpoint().x); 
+      MCP_endPoint_y.push_back(m_MCp.getEndpoint().y); 
+      MCP_endPoint_z.push_back(m_MCp.getEndpoint().z); 
+      MCP_pdgid.push_back(m_MCp.getPDG());
+      MCP_gStatus.push_back(m_MCp.getGeneratorStatus());
+    }
+  }
+  }catch(GaudiException &e){
+    debug()<<"MC Particle is not available "<<endmsg;
+  }
+  m_mctree->Fill(); 
+
+  //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;
+  //}
+
+  Clear(); 
+  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;
+
+      for(int istat=0; istat<m_trk.trackStates_size(); istat++){
+        edm4hep::TrackState m_trkstate = m_trk.getTrackStates(istat);
+
+        m_trkstate_d0.push_back( m_trkstate.D0 );
+        m_trkstate_z0.push_back( m_trkstate.Z0 );
+        m_trkstate_phi.push_back( m_trkstate.phi );
+        m_trkstate_tanL.push_back( m_trkstate.tanLambda );
+        //m_trkstate_kappa.push_back( m_trkstate.Kappa);
+        m_trkstate_omega.push_back( m_trkstate.omega );
+        m_trkstate_refx.push_back( m_trkstate.referencePoint.x );
+        m_trkstate_refy.push_back( m_trkstate.referencePoint.y );
+        m_trkstate_refz.push_back( m_trkstate.referencePoint.z );
+        m_trkstate_location.push_back( m_trkstate.location );
+        m_trkstate_tag.push_back(i);
+      }
+
+      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();
+      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());
+
+      delete TrkInit_Helix;
+    }
+  }
+  }catch(GaudiException &e){
+    debug()<<"Si track is not available "<<endmsg;
+  }  
+  m_sitrktree->Fill();
+
+  Clear(); 
+  try{
+  auto const_TPCTrkCol = m_TPCTrk.get();
+  auto const_TPCTrkLinkCol = m_TPCTrkAssoHdl.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;
+      int NTrkHit = m_trk.trackerHits_size();
+
+      for(int istat=0; istat<m_trk.trackStates_size(); istat++){
+        edm4hep::TrackState m_trkstate = m_trk.getTrackStates(istat);
+
+        m_trkstate_d0.push_back( m_trkstate.D0 );
+        m_trkstate_z0.push_back( m_trkstate.Z0 );
+        m_trkstate_phi.push_back( m_trkstate.phi );
+        m_trkstate_tanL.push_back( m_trkstate.tanLambda );
+        //m_trkstate_kappa.push_back( m_trkstate.Kappa);
+        m_trkstate_omega.push_back( m_trkstate.omega );
+        m_trkstate_refx.push_back( m_trkstate.referencePoint.x );
+        m_trkstate_refy.push_back( m_trkstate.referencePoint.y );
+        m_trkstate_refz.push_back( m_trkstate.referencePoint.z );
+        m_trkstate_location.push_back( m_trkstate.location );
+        m_trkstate_tag.push_back(i);
+      }
+
+      if(const_TPCTrkLinkCol){
+        for(auto iter: *const_TPCTrkLinkCol){
+          if(iter.getRec()==m_trk){
+            edm4hep::MCParticle m_mcp = iter.getSim();
+            m_truthMC_pid.push_back(m_mcp.getPDG());
+            m_truthMC_px.push_back(m_mcp.getMomentum().x);
+            m_truthMC_py.push_back(m_mcp.getMomentum().y);
+            m_truthMC_pz.push_back(m_mcp.getMomentum().z);
+            m_truthMC_E.push_back(m_mcp.getEnergy());
+            m_truthMC_EPx.push_back(m_mcp.getEndpoint().x);
+            m_truthMC_EPy.push_back(m_mcp.getEndpoint().y);
+            m_truthMC_EPz.push_back(m_mcp.getEndpoint().z);
+            m_truthMC_weight.push_back(iter.getWeight()/NTrkHit);
+            m_truthMC_tag.push_back(i);
+          }
+        }
+      }
+
+      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]);
+
+      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());
+
+      delete TrkInit_Helix;
+    }
+  }
+  }catch(GaudiException &e){
+    debug()<<"TPC track is not available "<<endmsg;
+  }
+  m_tpctrktree->Fill();
+
+
+  Clear(); 
+  try{
+  auto const_fullTrkCol = m_fullTrk.get();
+  auto const_fullTrkLinkCol = m_fullTrkAssoHdl.get();
+  if(const_fullTrkCol){
+    N_fullTrk = const_fullTrkCol->size();
+    for(int i=0; i<N_fullTrk; i++){
+      auto m_trk = const_fullTrkCol->at(i);
+//cout<<"In track #"<<i<<": track state size "<<m_trk.trackStates_size()<<", track hit size "<<m_trk.trackerHits_size()<<endl;
+      if(m_trk.trackStates_size()==0) continue;
+      //if(m_trk.trackerHits_size()==0) continue;
+      int NTrkHit = m_trk.trackerHits_size();
+
+      for(int istat=0; istat<m_trk.trackStates_size(); istat++){
+        edm4hep::TrackState m_trkstate = m_trk.getTrackStates(istat);
+
+        m_trkstate_d0.push_back( m_trkstate.D0 );
+        m_trkstate_z0.push_back( m_trkstate.Z0 );
+        m_trkstate_phi.push_back( m_trkstate.phi );
+        m_trkstate_tanL.push_back( m_trkstate.tanLambda );
+        //m_trkstate_kappa.push_back( m_trkstate.Kappa);
+        m_trkstate_omega.push_back( m_trkstate.omega );
+        m_trkstate_refx.push_back( m_trkstate.referencePoint.x );
+        m_trkstate_refy.push_back( m_trkstate.referencePoint.y );
+        m_trkstate_refz.push_back( m_trkstate.referencePoint.z );
+        m_trkstate_location.push_back( m_trkstate.location );
+        m_trkstate_tag.push_back(i);
+      }
+
+      int NSihit = 0;
+      int NTPChit = 0;
+      int Nelse = 0;
+      for(int ihit=0; ihit<NTrkHit; ihit++){
+        auto m_hit = m_trk.getTrackerHits(ihit);
+        if(m_hit.getType()==8) NSihit++;
+        else if(m_hit.getType()==0) NTPChit++;
+        else Nelse++;
+      }
+      //cout<<"In track #"<<i<<": Nhit "<<NTrkHit<<", Si Hit "<<NSihit<<", TPC hit "<<NTPChit<<", else "<<Nelse<<endl;
+
+      if(const_fullTrkLinkCol){
+        for(auto iter: *const_fullTrkLinkCol){
+          if(iter.getRec()==m_trk){
+            edm4hep::MCParticle m_mcp = iter.getSim();
+            m_truthMC_pid.push_back(m_mcp.getPDG());
+            m_truthMC_px.push_back(m_mcp.getMomentum().x);
+            m_truthMC_py.push_back(m_mcp.getMomentum().y);
+            m_truthMC_pz.push_back(m_mcp.getMomentum().z);
+            m_truthMC_E.push_back(m_mcp.getEnergy());
+            m_truthMC_EPx.push_back(m_mcp.getEndpoint().x);
+            m_truthMC_EPy.push_back(m_mcp.getEndpoint().y);
+            m_truthMC_EPz.push_back(m_mcp.getEndpoint().z);
+            m_truthMC_weight.push_back(iter.getWeight()/NTrkHit);
+            m_truthMC_tag.push_back(i);
+          }
+        }
+      }
+
+      edm4hep::TrackState m_trkstate = m_trk.getTrackStates(0);
+      if(m_trkstate.location!=1){
+        std::cout<<"ERROR: first track state is not IP! Will use this for track momentum "<<std::endl;
+      }
+
+      double m_pt = (0.299792458*_Bfield)/fabs(m_trkstate.omega*1000.);
+      double m_phi = m_trkstate.phi;
+      double m_pz = m_trkstate.tanLambda * m_pt;
+      TVector3 TrkMom(m_pt*cos(m_phi), m_pt*sin(m_phi), m_pz);
+      //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]);
+
+      m_trk_Nhit.push_back(NTrkHit);
+      m_trk_type.push_back(m_trk.getType());
+      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());
+
+      //delete TrkInit_Helix;
+    }
+  }
+  }catch(GaudiException &e){
+    debug()<<"Full track is not available "<<endmsg;
+  }
+  m_fulltrktree->Fill();
+
+/*  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();
+  }
+*/
+
+  _nEvt++; 
+  return StatusCode::SUCCESS;
+}
+
+StatusCode ReadDigiAlg::finalize()
+{
+  debug() << "begin finalize ReadDigiAlg" << endmsg;
+  m_wfile->cd();
+  m_mctree->Write();
+  m_sitrktree->Write();
+  m_tpctrktree->Write();
+  m_fulltrktree->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_VTX_x.clear();
+  MCP_VTX_y.clear();
+  MCP_VTX_z.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_Nhit.clear();
+  m_trk_type.clear();
+  m_trkstate_d0.clear();
+  m_trkstate_z0.clear();
+  m_trkstate_phi.clear();
+  m_trkstate_tanL.clear();
+  m_trkstate_omega.clear();
+  m_trkstate_kappa.clear();
+  m_trkstate_refx.clear();
+  m_trkstate_refy.clear();
+  m_trkstate_refz.clear();
+  m_trkstate_tag.clear();
+  m_trkstate_location.clear();
+  m_truthMC_tag.clear();
+  m_truthMC_pid.clear();
+  m_truthMC_px.clear();
+  m_truthMC_py.clear();
+  m_truthMC_pz.clear();
+  m_truthMC_E.clear();
+  m_truthMC_EPx.clear();
+  m_truthMC_EPy.clear();
+  m_truthMC_EPz.clear();
+  m_truthMC_weight.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 "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/CEPCV4_FullDet_GearOutput.xml

+<gear>
+    <global detectorName="CEPC_v4" />
+    <!--Gear XML file automatically created with GearXML::createXMLFile ....-->
+    <BField type="ConstantBField" x="0.000000000e+00" y="0.000000000e+00" z="3.000000000e+00" />
+    <detectors>
+        <detector geartype="TPCParameters" name="TPC">
+            <maxDriftLength value="2.225000000e+03" />
+            <driftVelocity value="0.000000000e+00" />
+            <coordinateType value="polar" />
+            <modules>
+                <module>
+                    <moduleID value="0" />
+                    <readoutFrequency value="0.000000000e+00" />
+                    <PadRowLayout2D type="FixedPadSizeDiskLayout" rMin="3.840000000e+02" rMax="1.716000000e+03" padHeight="6.000000000e+00" padWidth="1.000000000e+00" maxRow="222" padGap="0.000000000e+00" phiMax="6.283185307e+00" />
+                    <offset x_r="0.000000000e+00" y_phi="0.000000000e+00" />
+                    <angle value="0.000000000e+00" />
+                    <enlargeActiveAreaBy value="0.000000000e+00" />
+                </module>
+            </modules>
+            <parameter name="TPCGasProperties_RadLen" type="double" value="1.155205461e+05" />
+            <parameter name="TPCGasProperties_dEdx" type="double" value="2.668179899e-07" />
+            <parameter name="TPCInnerWallProperties_RadLen" type="double" value="2.740688665e+03" />
+            <parameter name="TPCInnerWallProperties_dEdx" type="double" value="1.241647394e-05" />
+            <parameter name="TPCOuterWallProperties_RadLen" type="double" value="6.495646008e+03" />
+            <parameter name="TPCOuterWallProperties_dEdx" type="double" value="5.300932694e-06" />
+            <parameter name="TPCWallProperties_RadLen" type="double" value="2.740688665e+03" />
+            <parameter name="TPCWallProperties_dEdx" type="double" value="1.241647394e-05" />
+            <parameter name="tpcInnerRadius" type="double" value="3.290000000e+02" />
+            <parameter name="tpcInnerWallThickness" type="double" value="2.500000000e+01" />
+            <parameter name="tpcOuterRadius" type="double" value="1.808000000e+03" />
+            <parameter name="tpcOuterWallThickness" type="double" value="6.000000000e+01" />
+            <parameter name="tpcZAnode" type="double" value="4.600000000e+03" />
+        </detector>
+        <detector name="EcalBarrel" geartype="CalorimeterParameters">
+            <layout type="Barrel" symmetry="8" phi0="0.000000000e+00" />
+            <dimensions inner_r="1.847415655e+03" outer_z="2.350000000e+03" />
+            <layer repeat="19" thickness="5.250000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+            <layer repeat="1" thickness="6.300000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+            <layer repeat="9" thickness="7.350000000e+00" absorberThickness="4.200000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+        </detector>
+        <detector name="EcalEndcap" geartype="CalorimeterParameters">
+            <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" />
+            <dimensions inner_r="4.000000000e+02" outer_r="2.088800000e+03" inner_z="2.450000000e+03" />
+            <layer repeat="19" thickness="5.250000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+            <layer repeat="1" thickness="6.300000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+            <layer repeat="9" thickness="7.350000000e+00" absorberThickness="4.200000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+        </detector>
+        <detector name="EcalPlug" geartype="CalorimeterParameters">
+            <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" />
+            <dimensions inner_r="2.400000000e+02" outer_r="4.000000000e+02" inner_z="2.450000000e+03" />
+            <layer repeat="19" thickness="5.250000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+            <layer repeat="1" thickness="6.300000000e+00" absorberThickness="2.100000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+            <layer repeat="9" thickness="7.350000000e+00" absorberThickness="4.200000000e+00" cellSize0="1.016666667e+01" cellSize1="1.016666667e+01" />
+        </detector>
+        <detector name="YokeBarrel" geartype="CalorimeterParameters">
+            <layout type="Barrel" symmetry="12" phi0="0.000000000e+00" />
+            <dimensions inner_r="4.173929932e+03" outer_z="4.072000000e+03" />
+            <layer repeat="1" thickness="4.000000000e+01" absorberThickness="0.000000000e+00" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" />
+            <layer repeat="9" thickness="1.400000000e+02" absorberThickness="1.000000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" />
+            <layer repeat="3" thickness="6.000000000e+02" absorberThickness="5.600000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" />
+            <layer repeat="1" thickness="4.000000000e+01" absorberThickness="4.027623145e-320" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" />
+        </detector>
+        <detector name="YokeEndcap" geartype="CalorimeterParameters">
+            <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" />
+            <dimensions inner_r="3.200000000e+02" outer_r="7.414929932e+03" inner_z="4.072000000e+03" />
+            <layer repeat="1" thickness="1.000000000e+02" absorberThickness="1.000000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" />
+            <layer repeat="9" thickness="1.400000000e+02" absorberThickness="1.000000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" />
+            <layer repeat="2" thickness="6.000000000e+02" absorberThickness="5.600000000e+02" cellSize0="3.000000000e+01" cellSize1="3.000000000e+01" />
+        </detector>
+        <detector name="YokePlug" geartype="CalorimeterParameters">
+            <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" />
+            <dimensions inner_r="3.200000000e+02" outer_r="2.849254326e+03" inner_z="3.781430000e+03" />
+            <parameter name="YokePlugThickness" type="double" value="2.905700000e+02" />
+        </detector>
+        <detector name="HcalBarrel" geartype="CalorimeterParameters">
+            <layout type="Barrel" symmetry="8" phi0="1.570796327e+00" />
+            <dimensions inner_r="2.058000000e+03" outer_z="2.350000000e+03" />
+            <layer repeat="40" thickness="2.673000000e+01" absorberThickness="2.000000000e+01" cellSize0="1.000000000e+01" cellSize1="1.000000000e+01" />
+            <parameter name="Hcal_barrel_number_modules" type="int" value="5" />
+            <parameter name="N_cells_z" type="int" value="91" />
+            <parameter name="FrameWidth" type="double" value="1.000000000e+00" />
+            <parameter name="Hcal_lateral_structure_thickness" type="double" value="1.000000000e+01" />
+            <parameter name="Hcal_modules_gap" type="double" value="2.000000000e+00" />
+            <parameter name="Hcal_outer_radius" type="double" value="3.144432447e+03" />
+            <parameter name="Hcal_virtual_cell_size" type="double" value="1.000000000e+01" />
+            <parameter name="InnerOctoSize" type="double" value="1.704903012e+03" />
+            <parameter name="RPC_PadSeparation" type="double" value="0.000000000e+00" />
+            <parameter name="TPC_Ecal_Hcal_barrel_halfZ" type="double" value="2.350000000e+03" />
+        </detector>
+        <detector name="HcalEndcap" geartype="CalorimeterParameters">
+            <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" />
+            <dimensions inner_r="3.500000000e+02" outer_r="3.144432447e+03" inner_z="2.650000000e+03" />
+            <layer repeat="40" thickness="2.673000000e+01" absorberThickness="2.000000000e+01" cellSize0="1.000000000e+01" cellSize1="1.000000000e+01" />
+            <parameter name="FrameWidth" type="double" value="0.000000000e+00" />
+            <parameter name="Hcal_virtual_cell_size" type="double" value="1.000000000e+01" />
+        </detector>
+        <detector name="HcalRing" geartype="CalorimeterParameters">
+            <layout type="Endcap" symmetry="2" phi0="0.000000000e+00" />
+            <dimensions inner_r="2.138800000e+03" outer_r="3.144432447e+03" inner_z="2.450000000e+03" />
+            <layer repeat="6" thickness="2.673000000e+01" absorberThickness="2.000000000e+01" cellSize0="1.000000000e+01" cellSize1="1.000000000e+01" />
+            <parameter name="FrameWidth" type="double" value="0.000000000e+00" />
+            <parameter name="Hcal_virtual_cell_size" type="double" value="1.000000000e+01" />
+        </detector>
+        <detector name="Lcal" geartype="CalorimeterParameters">
+            <layout type="Endcap" symmetry="1" phi0="0.000000000e+00" />
+            <dimensions inner_r="3.225828541e+01" outer_r="9.880000000e+01" inner_z="9.519000000e+02" />
+            <layer repeat="30" thickness="4.290000000e+00" absorberThickness="3.500000000e+00" cellSize0="1.039714290e+00" cellSize1="1.308996939e-01" />
+            <parameter name="beam_crossing_angle" type="double" value="0.000000000e+00" />
+        </detector>
+        <detector name="VXD" geartype="ZPlanarParameters">
+            <type technology="HYBRID" />
+            <shell halfLength="1.450000000e+02" gap="0.000000000e+00" innerRadius="6.500000000e+01" outerRadius="6.549392000e+01" radLength="3.527597571e+02" />
+            <layers>
+                <layer nLadders="10" phi0="-1.570796327e+00">
+                    <ladder distance="1.600000000e+01" thickness="1.000000000e+00" width="1.150000000e+01" length="6.250000000e+01" offset="-1.874869853e+00" radLength="1.014262421e+03" />
+                    <sensitive distance="1.595000000e+01" thickness="5.000000000e-02" width="1.100000000e+01" length="6.250000000e+01" offset="-1.624869853e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="10" phi0="-1.570796327e+00">
+                    <ladder distance="1.700000000e+01" thickness="1.000000000e+00" width="1.150000000e+01" length="6.250000000e+01" offset="-1.874869853e+00" radLength="1.014262421e+03" />
+                    <sensitive distance="1.800000000e+01" thickness="5.000000000e-02" width="1.100000000e+01" length="6.250000000e+01" offset="-1.624869853e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="11" phi0="-1.570796327e+00">
+                    <ladder distance="3.700000000e+01" thickness="1.000000000e+00" width="2.250000000e+01" length="1.250000000e+02" offset="-1.837940563e+00" radLength="1.014262421e+03" />
+                    <sensitive distance="3.695000000e+01" thickness="5.000000000e-02" width="2.200000000e+01" length="1.250000000e+02" offset="-1.587940563e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="11" phi0="-1.570796327e+00">
+                    <ladder distance="3.800000000e+01" thickness="1.000000000e+00" width="2.250000000e+01" length="1.250000000e+02" offset="-1.837940563e+00" radLength="1.014262421e+03" />
+                    <sensitive distance="3.900000000e+01" thickness="5.000000000e-02" width="2.200000000e+01" length="1.250000000e+02" offset="-1.587940563e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="17" phi0="-1.570796327e+00">
+                    <ladder distance="5.800000000e+01" thickness="1.000000000e+00" width="2.250000000e+01" length="1.250000000e+02" offset="-2.636744400e+00" radLength="1.014262421e+03" />
+                    <sensitive distance="5.795000000e+01" thickness="5.000000000e-02" width="2.200000000e+01" length="1.250000000e+02" offset="-2.386744400e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="17" phi0="-1.570796327e+00">
+                    <ladder distance="5.900000000e+01" thickness="1.000000000e+00" width="2.250000000e+01" length="1.250000000e+02" offset="-2.636744400e+00" radLength="1.014262421e+03" />
+                    <sensitive distance="6.000000000e+01" thickness="5.000000000e-02" width="2.200000000e+01" length="1.250000000e+02" offset="-2.386744400e+00" radLength="9.366070445e+01" />
+                </layer>
+            </layers>
+        </detector>
+        <detector name="FTD" geartype="FTDParameters">
+            <layers>
+                <layer nPetals="16" nSensors="1" isDoubleSided="0" sensorType="PIXEL" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="2.205200000e+02">
+                    <support thickness="1.000000000e+00" width="1.224000000e+02" lengthMin="1.173582968e+01" lengthMax="6.042957721e+01" rInner="2.950000000e+01" radLength="2.807467352e+02" />
+                    <sensitive thickness="2.000000000e-02" width="1.224000000e+02" lengthMin="1.173582968e+01" lengthMax="6.042957721e+01" rInner="2.950000000e+01" radLength="9.366070445e+01" />
+                </layer>
+                <layer nPetals="16" nSensors="1" isDoubleSided="0" sensorType="PIXEL" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="3.715200000e+02">
+                    <support thickness="1.000000000e+00" width="1.213600000e+02" lengthMin="1.214956740e+01" lengthMax="6.042957721e+01" rInner="3.054000000e+01" radLength="2.807467352e+02" />
+                    <sensitive thickness="2.000000000e-02" width="1.213600000e+02" lengthMin="1.214956740e+01" lengthMax="6.042957721e+01" rInner="3.054000000e+01" radLength="9.366070445e+01" />
+                </layer>
+                <layer nPetals="16" nSensors="2" isDoubleSided="1" sensorType="STRIP" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="6.462000000e+02">
+                    <support thickness="2.000000000e+00" width="2.665000000e+02" lengthMin="1.292930388e+01" lengthMax="1.189495957e+02" rInner="3.250000000e+01" radLength="2.807467352e+02" />
+                    <sensitive thickness="2.000000000e-01" width="2.665000000e+02" lengthMin="1.292930388e+01" lengthMax="1.189495957e+02" rInner="3.250000000e+01" radLength="9.366070445e+01" />
+                </layer>
+                <layer nPetals="16" nSensors="2" isDoubleSided="1" sensorType="STRIP" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="8.472000000e+02">
+                    <support thickness="2.000000000e+00" width="2.750000000e+02" lengthMin="1.352604098e+01" lengthMax="1.229278430e+02" rInner="3.400000000e+01" radLength="2.807467352e+02" />
+                    <sensitive thickness="2.000000000e-01" width="2.750000000e+02" lengthMin="1.352604098e+01" lengthMax="1.229278430e+02" rInner="3.400000000e+01" radLength="9.366070445e+01" />
+                </layer>
+                <layer nPetals="16" nSensors="2" isDoubleSided="1" sensorType="STRIP" petalOpenningAngle="1.963495408e-01" phi0="0.000000000e+00" alpha="0.000000000e+00" zoffset="0.000000000e+00" zsign0="1.000000000e+00" zposition="9.262000000e+02">
+                    <support thickness="2.000000000e+00" width="2.735000000e+02" lengthMin="1.412277808e+01" lengthMax="1.229278430e+02" rInner="3.550000000e+01" radLength="2.807467352e+02" />
+                    <sensitive thickness="2.000000000e-01" width="2.735000000e+02" lengthMin="1.412277808e+01" lengthMax="1.229278430e+02" rInner="3.550000000e+01" radLength="9.366070445e+01" />
+                </layer>
+            </layers>
+            <parameter name="strip_angle_deg" type="double" value="5.000000000e+00" />
+            <parameter name="strip_length_mm" type="double" value="1.600000000e+03" />
+            <parameter name="strip_pitch_mm" type="double" value="1.000000000e-02" />
+            <parameter name="strip_width_mm" type="double" value="1.000000000e-03" />
+        </detector>
+        <detector name="SIT" geartype="ZPlanarParameters">
+            <type technology="CCD" />
+            <shell halfLength="0.000000000e+00" gap="0.000000000e+00" innerRadius="0.000000000e+00" outerRadius="0.000000000e+00" radLength="0.000000000e+00" />
+            <layers>
+                <layer nLadders="10" phi0="0.000000000e+00">
+                    <ladder distance="1.531000000e+02" thickness="1.000000000e+00" width="9.916044311e+01" length="3.680000000e+02" offset="0.000000000e+00" radLength="2.134851878e+02" />
+                    <sensitive distance="1.529000000e+02" thickness="2.000000000e-01" width="9.916044311e+01" length="3.680000000e+02" offset="0.000000000e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="10" phi0="0.000000000e+00">
+                    <ladder distance="1.544000000e+02" thickness="1.000000000e+00" width="1.001352022e+02" length="3.680000000e+02" offset="0.000000000e+00" radLength="2.134851878e+02" />
+                    <sensitive distance="1.554000000e+02" thickness="2.000000000e-01" width="1.001352022e+02" length="3.680000000e+02" offset="0.000000000e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="19" phi0="0.000000000e+00">
+                    <ladder distance="3.001000000e+02" thickness="1.000000000e+00" width="9.988891763e+01" length="6.440000000e+02" offset="0.000000000e+00" radLength="2.134851878e+02" />
+                    <sensitive distance="2.999000000e+02" thickness="2.000000000e-01" width="9.988891763e+01" length="6.440000000e+02" offset="0.000000000e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="19" phi0="0.000000000e+00">
+                    <ladder distance="3.014000000e+02" thickness="1.000000000e+00" width="1.003895291e+02" length="6.440000000e+02" offset="0.000000000e+00" radLength="2.134851878e+02" />
+                    <sensitive distance="3.024000000e+02" thickness="2.000000000e-01" width="1.003895291e+02" length="6.440000000e+02" offset="0.000000000e+00" radLength="9.366070445e+01" />
+                </layer>
+            </layers>
+            <parameter name="sensor_length_mm" type="double" value="9.200000000e+01" />
+            <parameter name="strip_angle_deg" type="double" value="7.000000000e+00" />
+            <parameter name="strip_length_mm" type="double" value="9.200000000e+01" />
+            <parameter name="strip_pitch_mm" type="double" value="5.000000000e-02" />
+            <parameter name="strip_width_mm" type="double" value="1.250000000e-02" />
+            <parameter name="n_sensors_per_ladder" type="IntVec" value="8 8 14 14" />
+        </detector>
+        <detector name="SET" geartype="ZPlanarParameters">
+            <type technology="CCD" />
+            <shell halfLength="0.000000000e+00" gap="0.000000000e+00" innerRadius="0.000000000e+00" outerRadius="0.000000000e+00" radLength="0.000000000e+00" />
+            <layers>
+                <layer nLadders="24" phi0="0.000000000e+00">
+                    <ladder distance="1.811100000e+03" thickness="1.000000000e+00" width="4.766190158e+02" length="2.300000000e+03" offset="0.000000000e+00" radLength="2.134851878e+02" />
+                    <sensitive distance="1.810900000e+03" thickness="2.000000000e-01" width="4.766190158e+02" length="2.300000000e+03" offset="0.000000000e+00" radLength="9.366070445e+01" />
+                </layer>
+                <layer nLadders="24" phi0="0.000000000e+00">
+                    <ladder distance="1.812400000e+03" thickness="1.000000000e+00" width="4.770139733e+02" length="2.300000000e+03" offset="0.000000000e+00" radLength="2.134851878e+02" />
+                    <sensitive distance="1.813400000e+03" thickness="2.000000000e-01" width="4.770139733e+02" length="2.300000000e+03" offset="0.000000000e+00" radLength="9.366070445e+01" />
+                </layer>
+            </layers>
+            <parameter name="sensor_length_mm" type="double" value="9.200000000e+01" />
+            <parameter name="strip_angle_deg" type="double" value="7.000000000e+00" />
+            <parameter name="strip_length_mm" type="double" value="9.200000000e+01" />
+            <parameter name="strip_pitch_mm" type="double" value="5.000000000e-02" />
+            <parameter name="strip_width_mm" type="double" value="1.250000000e-02" />
+            <parameter name="n_sensors_per_ladder" type="IntVec" value="50 50" />
+        </detector>
+        <detector name="BeamPipe" geartype="GearParameters">
+            <parameter name="BeamPipeHalfZ" type="double" value="7.300000000e+02" />
+            <parameter name="BeamPipeProperties_RadLen" type="double" value="3.527597571e+02" />
+            <parameter name="BeamPipeProperties_dEdx" type="double" value="2.941795296e-04" />
+            <parameter name="BeamPipeRadius" type="double" value="1.400000000e+01" />
+            <parameter name="BeamPipeThickness" type="double" value="5.000000000e-01" />
+            <parameter name="RInner" type="DoubleVec" value="1.400000000e+01 1.400000000e+01 2.500000000e+00 1.300000000e+01 1.300000000e+01 1.300000000e+01 1.550000000e+01 1.550000000e+01 1.900000000e+01 1.900000000e+01 2.500000000e+01 2.500000000e+01 1.300000000e+01 1.300000000e+01 2.050000000e+01 2.050000000e+01 2.300000000e+01 2.300000000e+01 2.600000000e+01 2.600000000e+01 3.200000000e+01 3.200000000e+01" />
+            <parameter name="ROuter" type="DoubleVec" value="1.450000000e+01 1.450000000e+01 1.800000000e+01 1.800000000e+01 1.550000000e+01 1.550000000e+01 1.900000000e+01 1.900000000e+01 2.500000000e+01 2.500000000e+01 3.300000000e+01 3.300000000e+01 1.550000000e+01 1.550000000e+01 2.300000000e+01 2.300000000e+01 2.600000000e+01 2.600000000e+01 3.200000000e+01 3.200000000e+01 4.000000000e+01 4.000000000e+01" />
+            <parameter name="Z" type="DoubleVec" value="0.000000000e+00 5.000000000e+02 7.000000000e+02 7.010000000e+02 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 2.200000000e+03 3.950000000e+03 3.950000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03 4.450000000e+03" />
+        </detector>
+        <detector name="CoilParameters" geartype="GearParameters">
+            <parameter name="Coil_cryostat_c_modules_half_z" type="double" value="1.224000000e+03" />
+            <parameter name="Coil_cryostat_c_modules_inner_radius" type="double" value="3.348930000e+03" />
+            <parameter name="Coil_cryostat_c_modules_outer_radius" type="double" value="3.599930000e+03" />
+            <parameter name="Coil_cryostat_half_z" type="double" value="3.872000000e+03" />
+            <parameter name="Coil_cryostat_inner_cyl_half_z" type="double" value="3.872000000e+03" />
+            <parameter name="Coil_cryostat_inner_cyl_inner_radius" type="double" value="3.173930000e+03" />
+            <parameter name="Coil_cryostat_inner_cyl_outer_radius" type="double" value="3.963930000e+03" />
+            <parameter name="Coil_cryostat_inner_radius" type="double" value="3.173930000e+03" />
+            <parameter name="Coil_cryostat_mandrel_half_z" type="double" value="3.675000000e+03" />
+            <parameter name="Coil_cryostat_mandrel_inner_radius" type="double" value="3.599930000e+03" />
+            <parameter name="Coil_cryostat_mandrel_outer_radius" type="double" value="3.827930000e+03" />
+            <parameter name="Coil_cryostat_modules_half_z" type="double" value="7.960000000e+02" />
+            <parameter name="Coil_cryostat_modules_inner_radius" type="double" value="3.348930000e+03" />
+            <parameter name="Coil_cryostat_modules_outer_radius" type="double" value="3.599930000e+03" />
+            <parameter name="Coil_cryostat_outer_cyl_half_z" type="double" value="3.872000000e+03" />
+            <parameter name="Coil_cryostat_outer_cyl_inner_radius" type="double" value="3.893930000e+03" />
+            <parameter name="Coil_cryostat_outer_cyl_outer_radius" type="double" value="3.923930000e+03" />
+            <parameter name="Coil_cryostat_outer_radius" type="double" value="3.923930000e+03" />
+            <parameter name="Coil_cryostat_scint1_inner_radius" type="double" value="3.263930000e+03" />
+            <parameter name="Coil_cryostat_scint1_outer_radius" type="double" value="3.273930000e+03" />
+            <parameter name="Coil_cryostat_scint1_zposend" type="double" value="3.972000000e+03" />
+            <parameter name="Coil_cryostat_scint1_zposin" type="double" value="3.772000000e+03" />
+            <parameter name="Coil_cryostat_scint2_inner_radius" type="double" value="3.278930000e+03" />
+            <parameter name="Coil_cryostat_scint2_outer_radius" type="double" value="3.288930000e+03" />
+            <parameter name="Coil_cryostat_scint2_zposend" type="double" value="3.972000000e+03" />
+            <parameter name="Coil_cryostat_scint2_zposin" type="double" value="3.772000000e+03" />
+            <parameter name="Coil_cryostat_scint3_inner_radius" type="double" value="3.833930000e+03" />
+            <parameter name="Coil_cryostat_scint3_outer_radius" type="double" value="3.843930000e+03" />
+            <parameter name="Coil_cryostat_scint3_zposend" type="double" value="3.972000000e+03" />
+            <parameter name="Coil_cryostat_scint3_zposin" type="double" value="3.772000000e+03" />
+            <parameter name="Coil_cryostat_scint4_inner_radius" type="double" value="3.818930000e+03" />
+            <parameter name="Coil_cryostat_scint4_outer_radius" type="double" value="3.828930000e+03" />
+            <parameter name="Coil_cryostat_scint4_zposend" type="double" value="3.972000000e+03" />
+            <parameter name="Coil_cryostat_scint4_zposin" type="double" value="3.772000000e+03" />
+            <parameter name="Coil_cryostat_side_l_half_z" type="double" value="2.500000000e+01" />
+            <parameter name="Coil_cryostat_side_l_inner_radius" type="double" value="3.213930000e+03" />
+            <parameter name="Coil_cryostat_side_l_outer_radius" type="double" value="3.893930000e+03" />
+            <parameter name="Coil_cryostat_side_r_half_z" type="double" value="2.500000000e+01" />
+            <parameter name="Coil_cryostat_side_r_inner_radius" type="double" value="3.213930000e+03" />
+            <parameter name="Coil_cryostat_side_r_outer_radius" type="double" value="3.893930000e+03" />
+            <parameter name="Coil_material_c_modules" type="string" value="aluminium" />
+            <parameter name="Coil_material_inner_cyl" type="string" value="aluminium" />
+            <parameter name="Coil_material_mandrel" type="string" value="aluminium" />
+            <parameter name="Coil_material_modules" type="string" value="aluminium" />
+            <parameter name="Coil_material_outer_cyl" type="string" value="aluminium" />
+            <parameter name="Coil_material_scint1" type="string" value="polystyrene" />
+            <parameter name="Coil_material_scint2" type="string" value="polystyrene" />
+            <parameter name="Coil_material_scint3" type="string" value="polystyrene" />
+            <parameter name="Coil_material_scint4" type="string" value="polystyrene" />
+            <parameter name="Coil_material_side_l" type="string" value="aluminium" />
+            <parameter name="Coil_material_side_r" type="string" value="aluminium" />
+        </detector>
+        <detector name="MokkaParameters" geartype="GearParameters">
+            <parameter name="Ecal_endcap_outer_radius" type="string" value="2088.8" />
+            <parameter name="Ecal_endcap_plug_rmin" type="string" value="240" />
+            <parameter name="Ecal_endcap_zmax" type="string" value="2635" />
+            <parameter name="Ecal_endcap_zmin" type="string" value="2450" />
+            <parameter name="Ecal_outer_radius" type="string" value="2028" />
+            <parameter name="Hcal_R_max" type="string" value="3144.43" />
+            <parameter name="Hcal_endcap_zmin" type="string" value="2650" />
+            <parameter name="Lcal_z_begin" type="string" value="951.9" />
+            <parameter name="Lcal_z_thickness" type="string" value="128.1" />
+            <parameter name="MokkaModel" type="string" value="CEPC_v4" />
+            <parameter name="MokkaVersion" type="string" value="void" />
+            <parameter name="SIT1_Half_Length_Z" type="string" value="368" />
+            <parameter name="SIT1_Radius" type="string" value="152.9" />
+            <parameter name="SIT2_Half_Length_Z" type="string" value="644" />
+            <parameter name="SIT2_Radius" type="string" value="299.9" />
+            <parameter name="SiTrackerEndcap" type="string" value="FTD_PIXEL,29.5,151.9,220,16;FTD_PIXEL,30.54,151.9,371,16;FTD_STRIP,32.5,299,645,16;FTD_STRIP,34,309,846,16;FTD_STRIP,35.5,309,925,16" />
+            <parameter name="SiTrackerLayerStructure" type="string" value="FTD_PIXEL,Si:-0.02,CarbonFiber:1;FTD_STRIP,Si:-0.2,CarbonFiber:2,Si:-0.2" />
+            <parameter name="TPC_Ecal_Hcal_barrel_halfZ" type="string" value="2350" />
+            <parameter name="Yoke_Z_start_endcaps" type="string" value="4072" />
+            <parameter name="Yoke_barrel_inner_radius" type="string" value="4173.929931640625" />
+            <parameter name="calorimeter_region_rmax" type="string" value="3144.43" />
+            <parameter name="calorimeter_region_zmax" type="string" value="3736.43" />
+            <parameter name="tracker_region_rmax" type="string" value="1842.9" />
+            <parameter name="tracker_region_zmax" type="string" value="2350" />
+            <parameter name="world_box_hx" type="string" value="" />
+            <parameter name="world_box_hy" type="string" value="" />
+            <parameter name="world_box_hz" type="string" value="" />
+        </detector>
+        <detector name="VXDInfra" geartype="GearParameters">
+            <parameter name="ActiveLayerProperties_dEdx" type="double" value="3.870163611e-04" />
+            <parameter name="BeSupportEndplateThickness" type="double" value="2.000000000e+00" />
+            <parameter name="BeSupport_dEdx" type="double" value="2.941795296e-04" />
+            <parameter name="CryostatAlHalfZ" type="double" value="1.766000000e+02" />
+            <parameter name="CryostatAlInnerR" type="double" value="2.420000000e+01" />
+            <parameter name="CryostatAlRadius" type="double" value="1.000000000e+02" />
+            <parameter name="CryostatAlThickness" type="double" value="5.000000000e-01" />
+            <parameter name="CryostatAlZEndCap" type="double" value="1.768500000e+02" />
+            <parameter name="CryostatFoamRadius" type="double" value="9.000000000e+01" />
+            <parameter name="CryostatFoamThickness" type="double" value="1.000000000e+01" />
+            <parameter name="Cryostat_RadLen" type="double" value="8.896317758e+01" />
+            <parameter name="Cryostat_dEdx" type="double" value="4.350185478e-04" />
+            <parameter name="ElectronicEndLength" type="double" value="1.000000000e+01" />
+            <parameter name="ElectronicEndThickness" type="double" value="1.000000000e-01" />
+            <parameter name="StripLineBeamPipeRadius" type="double" value="2.430000000e+01" />
+            <parameter name="VXDEndPlateInnerRadius" type="double" value="3.000000000e+01" />
+            <parameter name="VXDSupport_dEdx" type="double" value="5.431907412e-05" />
+            <parameter name="LadderGaps" type="DoubleVec" value="0.000000000e+00 0.000000000e+00 0.000000000e+00 0.000000000e+00 0.000000000e+00 0.000000000e+00" />
+            <parameter name="StripLineFinalZ" type="DoubleVec" value="1.500000000e+02 1.500000000e+02 1.500000000e+02 1.500000000e+02 1.500000000e+02 1.500000000e+02" />
+        </detector>
+    </detectors>
+    <materials>
+        <material name="VXDFoamShellMaterial" A="1.043890843e+01" Z="5.612886646e+00" density="2.500000000e+01" radLength="1.751650267e+04" intLength="6.594366018e+01" />
+        <material name="VXDSupportMaterial" A="2.075865162e+01" Z="1.039383117e+01" density="2.765900000e+02" radLength="1.014262421e+03" intLength="1.206635688e+02" />
+    </materials>
+</gear>

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/calorimeter/SHcalSc04_Barrel_v04.cpp

@@ -22,10 +22,11 @@
 
 using namespace std;
 
-using dd4hep::BUILD_ENVELOPE;
+using dd4hep::_toString;
 using dd4hep::Box;
-using dd4hep::DetElement;
+using dd4hep::BUILD_ENVELOPE;
 using dd4hep::Detector;
+using dd4hep::DetElement;
 using dd4hep::IntersectionSolid;
 using dd4hep::Material;
 using dd4hep::PlacedVolume;
@@ -40,679 +41,607 @@ using dd4hep::Transform3D;
 using dd4hep::Trapezoid;
 using dd4hep::Tube;
 using dd4hep::Volume;
-using dd4hep::_toString;
 
 using dd4hep::rec::LayeredCalorimeterData;
 
 // After reading in all the necessary parameters.
 // To check the radius range and the space for placing the total layers
-static bool validateEnvelope(double rInner, double rOuter, double radiatorThickness, double layerThickness, int layerNumber, int nsymmetry){
-  
+static bool validateEnvelope(double rInner, double rOuter, double radiatorThickness, double layerThickness, int layerNumber, int nsymmetry)
+{
+
   bool Error = false;
   bool Warning = false;
-  double spaceAllowed = rOuter*cos(M_PI/nsymmetry) - rInner;
-  double spaceNeeded  = (radiatorThickness + layerThickness)* layerNumber;
-  double spaceToleranted  = (radiatorThickness + layerThickness)* (layerNumber+1);
-  double rOuterRecommaned = ( rInner + spaceNeeded )/cos(M_PI/16.);
-  int layerNumberRecommaned = floor( ( spaceAllowed ) / (radiatorThickness + layerThickness) );
-  
-  
-  if( spaceNeeded > spaceAllowed )
-    {
-      printout( dd4hep::ERROR,  "SHcalSc04_Barrel_v04", " Layer number is more than it can be built! "  ) ;
-      Error = true;
-    }
-  else if ( spaceToleranted < spaceAllowed )
-    {
-      printout( dd4hep::WARNING,  "SHcalSc04_Barrel_v04", " Layer number is less than it is able to build!" ) ;
-      Warning = true;
-    }
+  double spaceAllowed = rOuter * cos(M_PI / nsymmetry) - rInner;
+  double spaceNeeded = (radiatorThickness + layerThickness) * layerNumber;
+  double spaceToleranted = (radiatorThickness + layerThickness) * (layerNumber + 1);
+  double rOuterRecommaned = (rInner + spaceNeeded) / cos(M_PI / 16.);
+  int layerNumberRecommaned = floor((spaceAllowed) / (radiatorThickness + layerThickness));
+
+  if (spaceNeeded > spaceAllowed)
+  {
+    printout(dd4hep::ERROR, "SHcalSc04_Barrel_v04", " Layer number is more than it can be built! ");
+    Error = true;
+  }
+  else if (spaceToleranted < spaceAllowed)
+  {
+    printout(dd4hep::WARNING, "SHcalSc04_Barrel_v04", " Layer number is less than it is able to build!");
+    Warning = true;
+  }
   else
-    {
-      printout( dd4hep::DEBUG,  "SHcalSc04_Barrel_v04"," has been validated and start to build it." ) ;
-      Error = false;
-      Warning = false;
-    }
-
-  if( Error )
-    {
-      cout<<"\n ============> First Help Documentation <=============== \n"
-	  <<" When you see this message, that means you are crashing the module. \n"
-	  <<" Please take a cup of cafe, and think about what you want to do! \n"
-	  <<" Here are few FirstAid# for you. Please read them carefully. \n"
-	  <<" \n"
-	  <<" ###  FirstAid 1: ###\n"
-	  <<" If you want to build HCAL within the rInner and rOuter range, \n"
-	  <<" please reduce the layer number to " << layerNumberRecommaned <<" \n"
-	  <<" with the current layer thickness structure. \n"
-	  <<" \n"
-	  <<" You may redisgn the layer structure and thickness, too. \n"
-	  <<" \n"
-	  <<" ###  FirstAid 2: ###\n"
-	  <<" If you want to build HCAL with this layer number and the layer thickness, \n"
-	  <<" you have to update rOuter to "<< rOuterRecommaned*10. <<"*mm \n"
-	  <<" and to inform other subdetector, you need this space for building your design. \n"
-	  <<" \n"
-	  <<" ###  FirstAid 3: ###\n"
-	  <<" Do you think that you are looking for another type of HCAL driver? \n"
-	  <<" \n"
-	  <<endl; 
-      throw lcgeo::GeometryException(  "SHcalSc04_Barrel: Error: Layer number is more than it can be built!"   ) ;
-    } 
-  else if( Warning )
-    {
-      cout<<"\n ============> First Help Documentation <=============== \n"
-	  <<" When you see this warning message, that means you are changing the module. \n"
-	  <<" Please take a cup of cafe, and think about what you want to do! \n"
-	  <<" Here are few FirstAid# for you. Please read them carefully. \n"
-	  <<" \n"
-	  <<" ###  FirstAid 1: ###\n"
-	  <<" If you want to build HCAL within the rInner and rOuter range, \n"
-	  <<" You could build the layer number up to " << layerNumberRecommaned <<" \n"
-	  <<" with the current layer thickness structure. \n"
-	  <<" \n"
-	  <<" You may redisgn the layer structure and thickness, too. \n"
-	  <<" \n"
-	  <<" ###  FirstAid 2: ###\n"
-	  <<" If you want to build HCAL with this layer number and the layer thickness, \n"
-	  <<" you could reduce rOuter to "<< rOuterRecommaned*10. <<"*mm \n"
-	  <<" and to reduce the back plate thickness, which you may not need for placing layer. \n"
-	  <<" \n"
-	  <<" ###  FirstAid 3: ###\n"
-	  <<" Do you think that you are looking for another type of HCAL driver? \n"
-	  <<" \n"
-	  <<endl; 
-      return Warning;
-    }
-  else { return true; }
+  {
+    printout(dd4hep::DEBUG, "SHcalSc04_Barrel_v04", " has been validated and start to build it.");
+    Error = false;
+    Warning = false;
+  }
 
+  if (Error)
+  {
+    cout << "\n ============> First Help Documentation <=============== \n"
+         << " When you see this message, that means you are crashing the module. \n"
+         << " Please take a cup of cafe, and think about what you want to do! \n"
+         << " Here are few FirstAid# for you. Please read them carefully. \n"
+         << " \n"
+         << " ###  FirstAid 1: ###\n"
+         << " If you want to build HCAL within the rInner and rOuter range, \n"
+         << " please reduce the layer number to " << layerNumberRecommaned << " \n"
+         << " with the current layer thickness structure. \n"
+         << " \n"
+         << " You may redisgn the layer structure and thickness, too. \n"
+         << " \n"
+         << " ###  FirstAid 2: ###\n"
+         << " If you want to build HCAL with this layer number and the layer thickness, \n"
+         << " you have to update rOuter to " << rOuterRecommaned * 10. << "*mm \n"
+         << " and to inform other subdetector, you need this space for building your design. \n"
+         << " \n"
+         << " ###  FirstAid 3: ###\n"
+         << " Do you think that you are looking for another type of HCAL driver? \n"
+         << " \n"
+         << endl;
+    throw lcgeo::GeometryException("SHcalSc04_Barrel: Error: Layer number is more than it can be built!");
+  }
+  else if (Warning)
+  {
+    cout << "\n ============> First Help Documentation <=============== \n"
+         << " When you see this warning message, that means you are changing the module. \n"
+         << " Please take a cup of cafe, and think about what you want to do! \n"
+         << " Here are few FirstAid# for you. Please read them carefully. \n"
+         << " \n"
+         << " ###  FirstAid 1: ###\n"
+         << " If you want to build HCAL within the rInner and rOuter range, \n"
+         << " You could build the layer number up to " << layerNumberRecommaned << " \n"
+         << " with the current layer thickness structure. \n"
+         << " \n"
+         << " You may redisgn the layer structure and thickness, too. \n"
+         << " \n"
+         << " ###  FirstAid 2: ###\n"
+         << " If you want to build HCAL with this layer number and the layer thickness, \n"
+         << " you could reduce rOuter to " << rOuterRecommaned * 10. << "*mm \n"
+         << " and to reduce the back plate thickness, which you may not need for placing layer. \n"
+         << " \n"
+         << " ###  FirstAid 3: ###\n"
+         << " Do you think that you are looking for another type of HCAL driver? \n"
+         << " \n"
+         << endl;
+    return Warning;
+  }
+  else
+  {
+    return true;
+  }
 }
 
-static Ref_t create_detector(Detector& theDetector, xml_h element, SensitiveDetector sens)  {
+static Ref_t create_detector(Detector &theDetector, xml_h element, SensitiveDetector sens)
+{
 
   double boundarySafety = 0.0001;
 
-  xml_det_t   x_det       = element;
-  string      det_name    = x_det.nameStr();
-  int           det_id    = x_det.id();
-  DetElement  sdet( det_name, det_id );
+  xml_det_t x_det = element;
+  string det_name = x_det.nameStr();
+  int det_id = x_det.id();
+  DetElement sdet(det_name, det_id);
 
   // --- create an envelope volume and position it into the world ---------------------
-  
-  Volume envelope = dd4hep::xml::createPlacedEnvelope( theDetector,  element , sdet ) ;
-  
-  dd4hep::xml::setDetectorTypeFlag( element, sdet ) ;
-  
-  if( theDetector.buildType() == BUILD_ENVELOPE ) return sdet ;
+
+  Volume envelope = dd4hep::xml::createPlacedEnvelope(theDetector, element, sdet);
+
+  dd4hep::xml::setDetectorTypeFlag(element, sdet);
+
+  if (theDetector.buildType() == BUILD_ENVELOPE)
+    return sdet;
 
   //-----------------------------------------------------------------------------------
 
-  xml_comp_t    x_staves          = x_det.staves();
-  Material      stavesMaterial    = theDetector.material(x_staves.materialStr());
-  Material      air               = theDetector.air();
+  xml_comp_t x_staves = x_det.staves(); // absorber
+  Material stavesMaterial = theDetector.material(x_staves.materialStr());
+  Material air = theDetector.air();
 
   PlacedVolume pv;
 
   sens.setType("calorimeter");
 
-//====================================================================
-//
-// Read all the constant from ILD_o1_v05.xml
-// Use them to build HcalBarrel
-//
-//====================================================================
-  double      Hcal_inner_radius   = theDetector.constant<double>("Hcal_inner_radius");
-  double      Hcal_outer_radius   = theDetector.constant<double>("Hcal_outer_radius");
-  double      Hcal_half_length    = theDetector.constant<double>("Hcal_half_length");
-  int         Hcal_inner_symmetry = theDetector.constant<int>("Hcal_inner_symmetry");
-  int         Hcal_outer_symmetry = 0; // Fixed shape for Tube, and not allow to modify from compact xml.
-
-  double      Hcal_radiator_thickness          = theDetector.constant<double>("Hcal_radiator_thickness");
-  double      Hcal_chamber_thickness           = theDetector.constant<double>("Hcal_chamber_thickness");
-  double      Hcal_back_plate_thickness        = theDetector.constant<double>("Hcal_back_plate_thickness");
-  double      Hcal_lateral_plate_thickness     = theDetector.constant<double>("Hcal_lateral_structure_thickness");
-  double      Hcal_stave_gaps                  = theDetector.constant<double>("Hcal_stave_gaps");
-  double      Hcal_modules_gap                 = theDetector.constant<double>("Hcal_modules_gap"); 
-  double      Hcal_middle_stave_gaps           = theDetector.constant<double>("Hcal_middle_stave_gaps");
-  double      Hcal_layer_air_gap               = theDetector.constant<double>("Hcal_layer_air_gap");
-  //double      Hcal_cells_size                  = theDetector.constant<double>("Hcal_cells_size");
-
-  int         Hcal_nlayers                     = theDetector.constant<int>("Hcal_nlayers");
-
-  printout( dd4hep::INFO, "SHcalSc04_Barrel_v04", "Hcal_inner_radius : %e   - Hcal_outer_radius %e ", Hcal_inner_radius , Hcal_outer_radius);
+  //====================================================================
+  //
+  // Read all the constant from ILD_o1_v05.xml
+  // Use them to build HcalBarrel
+  //
+  //====================================================================
+  double Hcal_inner_radius = theDetector.constant<double>("Hcal_inner_radius");
+  double Hcal_outer_radius = theDetector.constant<double>("Hcal_outer_radius");
+  double Hcal_half_length = theDetector.constant<double>("Hcal_half_length");
+  int Hcal_inner_symmetry = theDetector.constant<int>("Hcal_inner_symmetry");
+  int Hcal_outer_symmetry = 0; // Fixed shape for Tube, and not allow to modify from compact xml
+  double Hcal_cell_size = theDetector.constant<double>("Hcal_cell_size");
+  double Hcal_cell_size_abnormal = theDetector.constant<double>("Hcal_cell_size_abnormal");
+  double Hcal_scintillator_air_gap = theDetector.constant<double>("Hcal_scintillator_air_gap");
+  double Hcal_scintillator_ESR_thickness = theDetector.constant<double>("Hcal_scintillator_ESR_thickness");
+  double Hcal_scintillator_thickness = theDetector.constant<double>("Hcal_scintillator_thickness");
+  double Hcal_radiator_thickness = theDetector.constant<double>("Hcal_radiator_thickness");
+  double Hcal_chamber_thickness = theDetector.constant<double>("Hcal_chamber_thickness");
+  double Hcal_back_plate_thickness = theDetector.constant<double>("Hcal_back_plate_thickness");
+  double Hcal_lateral_plate_thickness = theDetector.constant<double>("Hcal_lateral_structure_thickness");
+  double Hcal_stave_gaps = theDetector.constant<double>("Hcal_stave_gaps");
+  // double Hcal_modules_gap = theDetector.constant<double>("Hcal_modules_gap");
+  double Hcal_middle_stave_gaps = theDetector.constant<double>("Hcal_middle_stave_gaps");
+  double Hcal_layer_air_gap = theDetector.constant<double>("Hcal_layer_air_gap");
+  // double      Hcal_cells_size                  = theDetector.constant<double>("Hcal_cells_size");
+
+  int Hcal_nlayers = theDetector.constant<int>("Hcal_nlayers");
+
+  //=================================================================================
+  //  if Hcal_outer_radius stands for the radius of circumcircle, comment next line 
+  Hcal_outer_radius /= cos(M_PI / Hcal_inner_symmetry);
+  //=================================================================================
+  printout(dd4hep::INFO, "SHcalSc04_Barrel_v04", "Hcal_inner_radius : %e   - Hcal_outer_radius %e ", Hcal_inner_radius, Hcal_outer_radius);
   validateEnvelope(Hcal_inner_radius, Hcal_outer_radius, Hcal_radiator_thickness, Hcal_chamber_thickness, Hcal_nlayers, Hcal_inner_symmetry);
 
-  Readout readout = sens.readout();
-  dd4hep::Segmentation seg = readout.segmentation();
-  
-  dd4hep::DDSegmentation::BitField64 encoder = seg.decoder();
-  encoder.setValue(0) ;
-  
-
-  //fg: this is a bit tricky: we first have to check whether a multi segmentation is used and then, if so, we
-  //    will check the <subsegmentation key="" value=""/> element, for which subsegmentation to use for filling the 
-  //    DDRec:LayeredCalorimeterData information.
-  //    Additionally, we need to figure out if there is a TiledLayerGridXY instance defined -
-  //    and in case, there is more than one defined, we need to pick the correct one specified via the
-  //    <subsegmentation key="" value=""/> element.
-  //    This involves a lot of casting:  review the API in DD4hep !!
+  // fg: this is a bit tricky: we first have to check whether a multi segmentation is used and then, if so, we
+  //     will check the <subsegmentation key="" value=""/> element, for which subsegmentation to use for filling the
+  //     DDRec:LayeredCalorimeterData information.
+  //     Additionally, we need to figure out if there is a TiledLayerGridXY instance defined -
+  //     and in case, there is more than one defined, we need to pick the correct one specified via the
+  //     <subsegmentation key="" value=""/> element.
+  //     This involves a lot of casting:  review the API in DD4hep !!
 
   // check if we have a multi segmentation :
 
-  dd4hep::DDSegmentation::MultiSegmentation* multiSeg = 
-    dynamic_cast< dd4hep::DDSegmentation::MultiSegmentation*>( seg.segmentation() ) ;
-  
-  dd4hep::DDSegmentation::TiledLayerGridXY* tileSeg = 0 ;
-    
-  int sensitive_slice_number = -1 ;
-
-  if( multiSeg ){
-
-    try{ 
-      // check if we have an entry for the subsegmentation to be used 
-      xml_comp_t segxml = x_det.child( _Unicode( subsegmentation ) ) ;
-
-      std::string keyStr = segxml.attr<std::string>( _Unicode(key) ) ;
-      int keyVal = segxml.attr<int>( _Unicode(value) )  ;
-
-      encoder[ keyStr ] =  keyVal ;
-
-      // if we have a multisegmentation that uses the slice as key, we need to know for the
-      // computation of the layer parameters in LayeredCalorimeterData::Layer below
-      if( keyStr == "slice" ){
-	sensitive_slice_number = keyVal ;
-      }
-
-
-    } catch(const std::runtime_error &) {
-      throw lcgeo::GeometryException(  "SHcalSc04_Barrel: Error: MultiSegmentation specified but no "
-                                       " <subsegmentation key="" value=""/> element defined for detector ! " ) ;
-    }
-    
-    // check if we have a TiledLayerGridXY segmentation :
-    const dd4hep::DDSegmentation::TiledLayerGridXY* ts0 =
-      dynamic_cast<const dd4hep::DDSegmentation::TiledLayerGridXY*>(  &multiSeg->subsegmentation( encoder.getValue() ) ) ;
-    
-    tileSeg = const_cast<dd4hep::DDSegmentation::TiledLayerGridXY*>( ts0 ) ;
-    
-    if( ! tileSeg ){ // if the current segmentation is not a tileSeg, we see if there is another one
-      
-      for( auto s : multiSeg->subSegmentations() ){
-	const dd4hep::DDSegmentation::TiledLayerGridXY* ts =
-	  dynamic_cast<const dd4hep::DDSegmentation::TiledLayerGridXY*>( s.segmentation ) ;
-	
-	if( ts ) {
-	  tileSeg = const_cast<dd4hep::DDSegmentation::TiledLayerGridXY*>( ts ) ;
-	  break ;
-	}
-      }
-    }
-    
-  } else {
-    
-    tileSeg = 
-      dynamic_cast< dd4hep::DDSegmentation::TiledLayerGridXY*>( seg.segmentation() ) ;
-  }
-  
-  
-
-  std::vector<double> cellSizeVector = seg.cellDimensions( encoder.getValue() ); //Assume uniform cell sizes, provide dummy cellID
-  double cell_sizeX      = cellSizeVector[0];
-  double cell_sizeY      = cellSizeVector[1];
-
-
   //========== fill data for reconstruction ============================
-  LayeredCalorimeterData* caloData = new LayeredCalorimeterData ;
-  caloData->layoutType = LayeredCalorimeterData::BarrelLayout ;
-  caloData->inner_symmetry = Hcal_inner_symmetry  ;
-  caloData->outer_symmetry = Hcal_outer_symmetry  ;
-  caloData->phi0 = 0 ; // fg: also hardcoded below 
+  LayeredCalorimeterData *caloData = new LayeredCalorimeterData;
+  caloData->layoutType = LayeredCalorimeterData::BarrelLayout;
+  caloData->inner_symmetry = Hcal_inner_symmetry;
+  caloData->outer_symmetry = Hcal_outer_symmetry;
+  caloData->phi0 = 0; // fg: also hardcoded below
 
   /// extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm.
-  caloData->extent[0] = Hcal_inner_radius ;
-  caloData->extent[1] = Hcal_outer_radius ;
-  caloData->extent[2] = 0. ; // Barrel zmin is "0" by default.
-  caloData->extent[3] = Hcal_half_length ;
+  caloData->extent[0] = Hcal_inner_radius;
+  caloData->extent[1] = Hcal_outer_radius;
+  caloData->extent[2] = 0.; // Barrel zmin is "0" by default.
+  caloData->extent[3] = Hcal_half_length;
 
-//====================================================================
-//
-// general calculated parameters
-//
-//====================================================================
+  //====================================================================
+  //
+  // general calculated parameters
+  //
+  //====================================================================
 
-  double Hcal_total_dim_y   = Hcal_nlayers * (Hcal_radiator_thickness + Hcal_chamber_thickness) 
-                            + Hcal_back_plate_thickness;
+  double Hcal_total_dim_y = Hcal_nlayers * (Hcal_radiator_thickness + Hcal_chamber_thickness) + Hcal_back_plate_thickness;
 
-  double Hcal_y_dim1_for_x  = Hcal_outer_radius*cos(M_PI/Hcal_inner_symmetry) - Hcal_inner_radius;
-  double Hcal_bottom_dim_x  = 2.*Hcal_inner_radius*tan(M_PI/Hcal_inner_symmetry)- Hcal_stave_gaps;
-  double Hcal_normal_dim_z  = (2 * Hcal_half_length - Hcal_modules_gap)/2.;
+  // double Hcal_y_dim1_for_x = Hcal_outer_radius * cos(M_PI / Hcal_inner_symmetry) - Hcal_inner_radius;
+  double Hcal_bottom_dim_x = 2. * Hcal_inner_radius * tan(M_PI / Hcal_inner_symmetry) - Hcal_lateral_plate_thickness * 2 - Hcal_stave_gaps;
+  // double Hcal_normal_dim_z = (2 * Hcal_half_length - Hcal_modules_gap) / 2.;
+  double Hcal_normal_dim_z = Hcal_half_length;
 
- //only the middle has the steel plate.
-  double Hcal_regular_chamber_dim_z = Hcal_normal_dim_z - Hcal_lateral_plate_thickness;
+  // only the middle has the steel plate.
+  // double Hcal_regular_chamber_dim_z = Hcal_normal_dim_z - Hcal_lateral_plate_thickness;
+  double Hcal_regular_chamber_dim_z = Hcal_normal_dim_z;
 
-  //double Hcal_cell_dim_x            = Hcal_cells_size;
-  //double Hcal_cell_dim_z            = Hcal_regular_chamber_dim_z / floor (Hcal_regular_chamber_dim_z/Hcal_cell_dim_x);
+  // double Hcal_cell_dim_x            = Hcal_cells_size;
+  // double Hcal_cell_dim_z            = Hcal_regular_chamber_dim_z / floor (Hcal_regular_chamber_dim_z/Hcal_cell_dim_x);
+
+  // ========= Create Hcal Barrel stave   ====================================
+  //  It will be the volume for palcing the Hcal Barrel Chamber(i.e. Layers).
+  //  Itself will be placed into the world volume.
+  // ==========================================================================
 
- 
-// ========= Create Hcal Barrel stave   ====================================
-//  It will be the volume for palcing the Hcal Barrel Chamber(i.e. Layers).
-//  Itself will be placed into the world volume.
-// ==========================================================================
- 
-  double chambers_y_off_correction = 0.;
- 
   // stave modules shaper parameters
-  double BHX  = (Hcal_bottom_dim_x + Hcal_stave_gaps)/2.;
-  double THX  = (Hcal_total_dim_y + Hcal_inner_radius)*tan(M_PI/Hcal_inner_symmetry);
-  double YXH  = Hcal_total_dim_y / 2.;
-  double DHZ  = (Hcal_normal_dim_z - Hcal_lateral_plate_thickness) / 2.;
+  double BHX = (Hcal_bottom_dim_x + Hcal_lateral_plate_thickness * 2 + Hcal_stave_gaps) / 2.;
+  double THX = (Hcal_total_dim_y + Hcal_inner_radius) * tan(M_PI / Hcal_inner_symmetry);
+  double YXH = Hcal_total_dim_y / 2.;
+  // double DHZ = (Hcal_normal_dim_z - Hcal_lateral_plate_thickness) / 2.;
+  double DHZ = Hcal_normal_dim_z;
+
+  Trapezoid stave_shaper(THX, BHX, DHZ, DHZ, YXH);
 
-  Trapezoid stave_shaper(  THX, BHX, DHZ, DHZ, YXH);
+  Tube solidCaloTube(0, Hcal_outer_radius, DHZ + boundarySafety);
 
-  Tube solidCaloTube(0, Hcal_outer_radius, DHZ+boundarySafety);
-  
-  RotationZYX mrot(0,0,M_PI/2.);
+  RotationZYX mrot(0, 0, M_PI / 2.);
 
   Rotation3D mrot3D(mrot);
-  Position mxyzVec(0,0,(Hcal_inner_radius + Hcal_total_dim_y / 2.));
-  Transform3D mtran3D(mrot3D,mxyzVec);
+  Position mxyzVec(0, 0, (Hcal_inner_radius + Hcal_total_dim_y / 2.));
+  Transform3D mtran3D(mrot3D, mxyzVec);
 
   IntersectionSolid barrelModuleSolid(stave_shaper, solidCaloTube, mtran3D);
 
-  Volume  EnvLogHcalModuleBarrel(det_name+"_module",barrelModuleSolid,stavesMaterial);
-
-  EnvLogHcalModuleBarrel.setAttributes(theDetector,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
-
-
-
-
-
-  //stave modules lateral plate shaper parameters
-  double BHX_LP  = BHX;
-  double THX_LP  = THX;
-  double YXH_LP  = YXH;
+  Volume EnvLogHcalModuleBarrel(det_name + "_module", barrelModuleSolid, stavesMaterial);
 
-  //build lateral palte here to simulate lateral plate in the middle of barrel.
-  double DHZ_LP  = Hcal_lateral_plate_thickness/2.0; 
+  EnvLogHcalModuleBarrel.setAttributes(theDetector, x_staves.regionStr(), x_staves.limitsStr(), x_staves.visStr());
 
-  Trapezoid stave_shaper_LP(THX_LP, BHX_LP, DHZ_LP, DHZ_LP, YXH_LP);
+  // stave modules lateral plate shaper parameters
+  // double BHX_LP = BHX;
+  // double THX_LP = THX;
+  // double YXH_LP = YXH;
 
-  Tube solidCaloTube_LP(0, Hcal_outer_radius, DHZ_LP+boundarySafety);
+  // // build lateral palte here to simulate lateral plate in the middle of barrel.
+  // double DHZ_LP = Hcal_lateral_plate_thickness / 2.0;
 
-  IntersectionSolid Module_lateral_plate(stave_shaper_LP, solidCaloTube_LP, mtran3D);
+  // Trapezoid stave_shaper_LP(THX_LP, BHX_LP, DHZ_LP, DHZ_LP, YXH_LP);
 
-  Volume  EnvLogHcalModuleBarrel_LP(det_name+"_Module_lateral_plate",Module_lateral_plate,stavesMaterial);
+  // Tube solidCaloTube_LP(0, Hcal_outer_radius, DHZ_LP + boundarySafety);
 
-  EnvLogHcalModuleBarrel_LP.setAttributes(theDetector,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
+  // IntersectionSolid Module_lateral_plate(stave_shaper_LP, solidCaloTube_LP, mtran3D);
 
+  // Volume EnvLogHcalModuleBarrel_LP(det_name + "_Module_lateral_plate", Module_lateral_plate, stavesMaterial);
 
+  // EnvLogHcalModuleBarrel_LP.setAttributes(theDetector, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
 
 #ifdef SHCALSC04_DEBUG
-  std::cout<< " ==> Hcal_outer_radius: "<<Hcal_outer_radius <<std::endl;
+  std::cout << " ==> Hcal_outer_radius: " << Hcal_outer_radius << std::endl;
 #endif
 
-
-
-
-
-//====================================================================
-//
-// Chambers in the HCAL BARREL
-//
-//====================================================================
+  //====================================================================
+  //
+  // Chambers in the HCAL BARREL
+  //
+  //====================================================================
   // Build Layer Chamber fill with air, which include the tolarance space at the two side
   // place the slice into the Layer Chamber
   // Place the Layer Chamber into the Stave module
   // place the Stave module into the asembly Volume
   // place the module middle lateral plate into asembly Volume
 
-  double x_length; //dimension of an Hcal barrel layer on the x-axis
-  double y_height; //dimension of an Hcal barrel layer on the y-axis
-  double z_width;  //dimension of an Hcal barrel layer on the z-axis
-  x_length = 0.; // Each Layer the x_length has new value.
-  y_height = Hcal_chamber_thickness / 2.;
-  z_width  = Hcal_regular_chamber_dim_z/2.;
+  double x_halflength; // dimension of an Hcal barrel layer on the x-axis
+  double y_halfheight; // dimension of an Hcal barrel layer on the y-axis
+  double z_halfwidth;  // dimension of an Hcal barrel layer on the z-axis
+  x_halflength = 0.;   // Each Layer the x_halflength has new value.
+  y_halfheight = Hcal_chamber_thickness / 2.;
+  z_halfwidth = Hcal_regular_chamber_dim_z;
 
-  double xOffset = 0.;//the x_length of a barrel layer is calculated as a
-  //barrel x-dimension plus (bottom barrel) or minus
+  double xOffset = 0.; // the x_halflength of a barrel layer is calculated as a
+  // barrel x-dimension plus (bottom barrel) or minus
   //(top barrel) an x-offset, which depends on the angle M_PI/Hcal_inner_symmetry
 
-  double xShift = 0.;//Geant4 draws everything in the barrel related to the 
-  //center of the bottom barrel, so we need to shift the layers to
-  //the left (or to the right) with the quantity xShift
+  double xShift = 0.; // Geant4 draws everything in the barrel related to the
+  // center of the bottom barrel, so we need to shift the layers to
+  // the left (or to the right) with the quantity xShift
+
+  // read sensitive cell parameters and create it
+
+  // xml_comp_t x_scintillator = x_det.solid();
+  // Material sensitiveMaterial = theDetector.material(x_scintillator.materialStr());
+  // xml_comp_t x_ESR = x_det.shape();
+  // Material ladderMaterial = theDetector.material(x_ESR.materialStr());
+
+  Volume cell_vol;
+  Volume scintillator_vol;
+  Volume cell_vol_abnormal;
+  Volume scintillator_vol_abnormal;
+  xml_comp_t x_layer_tmp(x_det.child(_U(layer)));
+  for (xml_coll_t k(x_layer_tmp, _U(slice)); k; ++k)
+  {
+    xml_comp_t x_slice = k;
+    if (x_slice.isSensitive())
+    {
+      // child elements: ladder and sensitive
+      xml_comp_t x_sensitive(x_slice.child(_U(sensitive)));
+      xml_comp_t x_ladder(x_slice.child(_U(ladder)));
+      Material sensitiveMaterial = theDetector.material(x_sensitive.materialStr());
+      Material ladderMaterial = theDetector.material(x_ladder.materialStr());
+
+      // Store scintillator thickness
+      double cell_thickness = Hcal_scintillator_thickness + 2 * Hcal_scintillator_ESR_thickness;
+
+      double cell_size = 2 * Hcal_scintillator_ESR_thickness + Hcal_cell_size;
+      // place sensitive cell into ESR cell
+      string cell_name = "cell";
+      string scintillator_name = cell_name + "_scintillator";
+      // create sensitive cell with ESR
+      Box cell_box(cell_size / 2., cell_size / 2., cell_thickness / 2.);
+      // string cell_name=
+      // Volume cell_vol_tmp(cell_name, cell_box, ladderMaterial);
+      cell_vol = Volume(cell_name, cell_box, ladderMaterial);
+      cell_vol.setAttributes(theDetector, x_ladder.regionStr(), x_ladder.limitsStr(), x_ladder.visStr());
+      Box scintillator_box(Hcal_cell_size / 2., Hcal_cell_size / 2., Hcal_scintillator_thickness / 2.);
+      // Volume scintillator_vol_tmp(scintillator_name, scintillator_box, sensitiveMaterial);
+      scintillator_vol = Volume(scintillator_name, scintillator_box, sensitiveMaterial);
+      scintillator_vol.setAttributes(theDetector, x_sensitive.regionStr(), x_sensitive.limitsStr(), x_sensitive.visStr());
+      scintillator_vol.setSensitiveDetector(sens);
+      cell_vol.placeVolume(scintillator_vol, Position(0., 0., 0.));
+      // cout << x_ladder.visStr() << "   " << x_sensitive.visStr() << endl;
+      ///////////////////////////////
+      // abnormal cell             //
+      ///////////////////////////////
+      double cell_size_abnormal = 2 * Hcal_scintillator_ESR_thickness + Hcal_cell_size_abnormal;
+      // place sensitive cell into ESR cell
+      string cell_name_abnormal = "cell_abnormal";
+      string scintillator_name_abnormal = cell_name_abnormal + "_scintillator_abnormal";
+      // create sensitive cell with ESR
+      Box cell_box_abnormal(cell_size_abnormal / 2., cell_size / 2., cell_thickness / 2.);
+      // string cell_name=
+      // Volume cell_vol_tmp_abnormal(cell_name_abnormal, cell_box_abnormal, ladderMaterial);
+      cell_vol_abnormal = Volume(cell_name_abnormal, cell_box_abnormal, ladderMaterial);
+      cell_vol_abnormal.setAttributes(theDetector, x_ladder.regionStr(), x_ladder.limitsStr(), x_ladder.visStr());
+      Box scintillator_box_abnormal(Hcal_cell_size_abnormal / 2., Hcal_cell_size / 2., Hcal_scintillator_thickness / 2.);
+      // Volume scintillator_vol_tmp_abnormal(scintillator_name_abnormal, scintillator_box_abnormal, sensitiveMaterial);
+      scintillator_vol_abnormal = Volume(scintillator_name_abnormal, scintillator_box_abnormal, sensitiveMaterial);
+      scintillator_vol_abnormal.setAttributes(theDetector, x_sensitive.regionStr(), x_sensitive.limitsStr(), x_sensitive.visStr());
+      scintillator_vol_abnormal.setSensitiveDetector(sens);
+      cell_vol_abnormal.placeVolume(scintillator_vol_abnormal, Position(0., 0., 0.));
+    }
+  }
 
+  // sensitive cell creating done
 
   //-------------------- start loop over HCAL layers ----------------------
+  // int n_x_layer_1 = 0;
+  for (int layer_id = 1; layer_id <= (Hcal_nlayers); layer_id++)
+  {
 
-  for (int layer_id = 1; layer_id <= (2*Hcal_nlayers); layer_id++)
-    {
- 
-      double TanPiDiv8 = tan(M_PI/Hcal_inner_symmetry);
-      double x_total   = 0.;
-      double x_halfLength;
-      x_length  = 0.;
-
-      int logical_layer_id = 0;
-
-      if ( (layer_id < Hcal_nlayers)
-	   || (layer_id > Hcal_nlayers && layer_id < (2*Hcal_nlayers)) )
-	logical_layer_id = layer_id % Hcal_nlayers;
-      else if ( (layer_id == Hcal_nlayers) 
-		|| (layer_id == 2*Hcal_nlayers) ) logical_layer_id = Hcal_nlayers;
-
-      //---- bottom barrel------------------------------------------------------------
-      if( logical_layer_id *(Hcal_radiator_thickness + Hcal_chamber_thickness)
-	  < (Hcal_outer_radius * cos(M_PI/Hcal_inner_symmetry) - Hcal_inner_radius ) ) {
-	xOffset = (logical_layer_id * Hcal_radiator_thickness 
-		   + (logical_layer_id -1) * Hcal_chamber_thickness) * TanPiDiv8;
-
-	x_total  = Hcal_bottom_dim_x/2 - Hcal_middle_stave_gaps/2 + xOffset;
-	x_length = x_total - 2*Hcal_layer_air_gap;
-	x_halfLength = x_length/2.;
-
-      } else {//----- top barrel -------------------------------------------------
-	double y_layerID = logical_layer_id * (Hcal_radiator_thickness + Hcal_chamber_thickness) + Hcal_inner_radius;
-	double ro_layer = Hcal_outer_radius - Hcal_radiator_thickness;
-	
-	x_total = sqrt( ro_layer * ro_layer - y_layerID * y_layerID);
-	
-	x_length = x_total - Hcal_middle_stave_gaps;
-	
-	x_halfLength = x_length/2.;
-	
-	xOffset = (logical_layer_id * Hcal_radiator_thickness 
-		   + (logical_layer_id - 1) * Hcal_chamber_thickness - Hcal_y_dim1_for_x) / TanPiDiv8
-	  + Hcal_chamber_thickness / TanPiDiv8;
-	
-      }
+    double TanPiDiv8 = tan(M_PI / Hcal_inner_symmetry);
+    double x_total = 0.;
+    // double x_halfLength;
+    x_halflength = 0.;
 
-      double xAbsShift = (Hcal_middle_stave_gaps/2 + Hcal_layer_air_gap + x_halfLength);
-      
-      if (layer_id <= Hcal_nlayers)     xShift = - xAbsShift;
-      else if (layer_id > Hcal_nlayers) xShift = xAbsShift;
-
-      x_length = x_length/2.;
-
-      
-      //calculate the size of a fractional tile
-      //-> this sets fract_cell_dim_x
-      
-      //double fract_cell_dim_x = 0.;
-      //this->CalculateFractTileSize(2*x_length, Hcal_cell_dim_x, fract_cell_dim_x);
-      
-      //Vector newFractCellDim(fract_cell_dim_x, Hcal_chamber_thickness, Hcal_cell_dim_z);
-      //theBarrilRegSD->SetFractCellDimPerLayer(layer_id, newFractCellDim);
-
-      encoder["layer"] = logical_layer_id ;
-      cellSizeVector = seg.segmentation()->cellDimensions( encoder.getValue() ); 
-      cell_sizeX      = cellSizeVector[0];
-      cell_sizeY      = cellSizeVector[1];
-
-      LayeredCalorimeterData::Layer caloLayer ;
-      caloLayer.cellSize0 = cell_sizeX;
-      caloLayer.cellSize1 = cell_sizeY;
-
-
-      //--------------------------------------------------------------------------------
-      //  build chamber box, with the calculated dimensions 
-      //-------------------------------------------------------------------------------
-      printout( dd4hep::DEBUG,  "SHcalSc04_Barrel_v04", " \n Start to build layer chamber - layer_id: %d", layer_id ) ;
-      printout( dd4hep::DEBUG,  "SHcalSc04_Barrel_v04"," chamber x:y:z:  %e:%e:%e", x_length*2., z_width*2. , y_height*2. );
-
-      //check if x_length (scintillator length) is divisible with x_integerTileSize
-      if( layer_id <= Hcal_nlayers) {
-	double fracPart, intPart;
-	double temp = x_length*2./cell_sizeX;
-	fracPart = modf(temp, &intPart);
-	int noOfIntCells = int(temp);
-
-
-	if( tileSeg !=0 ){
-
-	  tileSeg->setBoundaryLayerX(x_length);
-	  
-	  if (fracPart == 0){ //divisible
-	    if ( noOfIntCells%2 ) {
-	      if( tileSeg !=0 ) tileSeg->setLayerOffsetX(0);
-	    }
-	    else {
-	      if( tileSeg !=0 ) tileSeg->setLayerOffsetX(1);
-	    }
-	    tileSeg->setFractCellSizeXPerLayer(0);
-	  }
-	  else if (fracPart>0){
-	    if ( noOfIntCells%2 ) {
-	      if( tileSeg !=0 ) tileSeg->setLayerOffsetX(1);
-	    }
-	    else {
-	      if( tileSeg !=0 ) tileSeg->setLayerOffsetX(0);
-	    }
-	    tileSeg->setFractCellSizeXPerLayer( (fracPart+1.0)/2.0*cell_sizeX );
-	  }
-	  
-	  if ( (int)( (z_width*2.) / cell_sizeX)%2 ){
-	    if( tileSeg !=0 ) tileSeg->setLayerOffsetY(0);
-	  }
-	  else {
-	    if( tileSeg !=0 ) tileSeg->setLayerOffsetY(1);
-	  }
-
-	}
-      }
-      Box ChamberSolid((x_length + Hcal_layer_air_gap),  //x + air gaps at two side, do not need to build air gaps individualy.
-			 z_width,   //z attention!
-			 y_height); //y attention!
+    // int layer_id = 0;
 
-      string ChamberLogical_name      = det_name+_toString(layer_id,"_layer%d");
+    // if ((layer_id < Hcal_nlayers) || (layer_id > Hcal_nlayers && layer_id < (2 * Hcal_nlayers)))
+    //   layer_id = layer_id % Hcal_nlayers;
+    // else if ((layer_id == Hcal_nlayers) || (layer_id == 2 * Hcal_nlayers))
+    //   layer_id = Hcal_nlayers;
 
-      Volume ChamberLogical(ChamberLogical_name, ChamberSolid, air);   
+    //---- bottom barrel------------------------------------------------------------
+    // if (layer_id * (Hcal_radiator_thickness + Hcal_chamber_thickness) < (Hcal_outer_radius * cos(M_PI / Hcal_inner_symmetry) - Hcal_inner_radius))
+    // {
+    xOffset = (layer_id * Hcal_radiator_thickness + (layer_id - 1) * Hcal_chamber_thickness) * TanPiDiv8;
 
+    x_total = Hcal_bottom_dim_x - Hcal_middle_stave_gaps + xOffset * 2;
+    x_halflength = x_total - 2 * Hcal_layer_air_gap;
+    // x_halfLength = x_halflength / 2.;
+    // }
+    // else
+    // { //----- top barrel -------------------------------------------------
+    //   double y_layerID = layer_id * (Hcal_radiator_thickness + Hcal_chamber_thickness) + Hcal_inner_radius;
+    //   double ro_layer = Hcal_outer_radius - Hcal_radiator_thickness;
 
+    //   x_total = sqrt(ro_layer * ro_layer - y_layerID * y_layerID);
 
-      double layer_thickness = y_height*2.;
+    //   x_halflength = x_total - Hcal_middle_stave_gaps;
 
-      double nRadiationLengths=0.;
-      double nInteractionLengths=0.;
-      double thickness_sum=0;
+    //   x_halfLength = x_halflength / 2.;
 
-      nRadiationLengths   = Hcal_radiator_thickness/(stavesMaterial.radLength());
-      nInteractionLengths = Hcal_radiator_thickness/(stavesMaterial.intLength());
-      thickness_sum       = Hcal_radiator_thickness;
+    //   xOffset = (layer_id * Hcal_radiator_thickness + (layer_id - 1) * Hcal_chamber_thickness - Hcal_y_dim1_for_x) / TanPiDiv8 + Hcal_chamber_thickness / TanPiDiv8;
+    // }
 
-//====================================================================
-// Create Hcal Barrel Chamber without radiator
-// Place into the Hcal Barrel stave, after each radiator 
-//====================================================================
-      xml_coll_t c(x_det,_U(layer));
-      xml_comp_t   x_layer = c;
-      string layer_name      = det_name+_toString(layer_id,"_layer%d");
-      
-      // Create the slices (sublayers) within the Hcal Barrel Chamber.
-      double slice_pos_z = layer_thickness/2.;
-      int slice_number = 0;
-
-      for(xml_coll_t k(x_layer,_U(slice)); k; ++k)  {
-	xml_comp_t x_slice = k;
-	string   slice_name      = layer_name + _toString(slice_number,"_slice%d");
-	double   slice_thickness = x_slice.thickness();
-	Material slice_material  = theDetector.material(x_slice.materialStr());
-	DetElement slice(layer_name,_toString(slice_number,"slice%d"),x_det.id());
-	
-	slice_pos_z -= slice_thickness/2.;
-	
-	// Slice volume & box
-	Volume slice_vol(slice_name,Box(x_length,z_width,slice_thickness/2.),slice_material);
-
-	nRadiationLengths   += slice_thickness/(2.*slice_material.radLength());
-	nInteractionLengths += slice_thickness/(2.*slice_material.intLength());
-	thickness_sum       += slice_thickness/2;
-	
-	if ( x_slice.isSensitive() ) {
-
-	  slice_vol.setSensitiveDetector(sens);
-
-	  // if we have a multisegmentation based on slices, we need to use the correct slice here
-	  if ( sensitive_slice_number<0  || sensitive_slice_number == slice_number ) {
-	
-	    //Store "inner" quantities
-	    caloLayer.inner_nRadiationLengths   = nRadiationLengths;
-	    caloLayer.inner_nInteractionLengths = nInteractionLengths;
-	    caloLayer.inner_thickness = thickness_sum;
-	    //Store scintillator thickness
-	    caloLayer.sensitive_thickness = slice_thickness;
-	    
-	    //Reset counters to measure "outside" quantitites
-	    nRadiationLengths=0.;
-	    nInteractionLengths=0.;
-	    thickness_sum = 0.;
-	  }
-	}
-
-	nRadiationLengths   += slice_thickness/(2.*slice_material.radLength());
-	nInteractionLengths += slice_thickness/(2.*slice_material.intLength());
-	thickness_sum       += slice_thickness/2;
-
-
-	// Set region, limitset, and vis.
-	slice_vol.setAttributes(theDetector,x_slice.regionStr(),x_slice.limitsStr(),x_slice.visStr());
-	// slice PlacedVolume
-	PlacedVolume slice_phv = ChamberLogical.placeVolume(slice_vol,Position(0.,0.,slice_pos_z));
-
-	slice_phv.addPhysVolID("layer",logical_layer_id).addPhysVolID("slice", slice_number );
-
-	
-	if ( x_slice.isSensitive() ) {
-	  int tower_id  = (layer_id > Hcal_nlayers)? 1:-1;
-	  slice_phv.addPhysVolID("tower",tower_id);
-	  printout( dd4hep::DEBUG,  "SHcalSc04_Barrel_v04", "  logical_layer_id:  %d  tower_id:  %d", logical_layer_id, tower_id  ) ;
-	}
-	
-	slice.setPlacement(slice_phv);
-	// Increment x position for next slice.
-	slice_pos_z -= slice_thickness/2.;
-	// Increment slice number.
-	++slice_number;             
-      }
+    x_halflength = x_halflength / 2.;
 
-      //Store "outer" quantities
-      caloLayer.outer_nRadiationLengths = nRadiationLengths;
-      caloLayer.outer_nInteractionLengths = nInteractionLengths;
-      caloLayer.outer_thickness = thickness_sum;
+    LayeredCalorimeterData::Layer caloLayer;
+    caloLayer.cellSize0 = Hcal_cell_size;
+    caloLayer.cellSize1 = Hcal_cell_size;
 
+    //--------------------------------------------------------------------------------
+    //  build chamber box, with the calculated dimensions
+    //-------------------------------------------------------------------------------
+    printout(dd4hep::DEBUG, "SHcalSc04_Barrel_v04", " \n Start to build layer chamber - layer_id: %d", layer_id);
+    printout(dd4hep::DEBUG, "SHcalSc04_Barrel_v04", " chamber x:y:z:  %e:%e:%e", x_halflength * 2., z_halfwidth * 2., y_halfheight * 2.);
 
-//---------------------------  Chamber Placements -----------------------------------------
+    // check if x_halflength (scintillator length) is divisible with x_integerTileSize
 
-      double chamber_x_offset, chamber_y_offset, chamber_z_offset;
-      chamber_x_offset = xShift;
+    Box ChamberSolid((x_halflength + Hcal_layer_air_gap), // x + air gaps at two side, do not need to build air gaps individualy.
+                     z_halfwidth,                         // z attention!
+                     y_halfheight);                       // y attention!
 
-      chamber_z_offset = 0;
+    string ChamberLogical_name = det_name + _toString(layer_id, "_layer%d");
 
+    Volume ChamberLogical(ChamberLogical_name, ChamberSolid, air);
+    ChamberLogical.setAttributes(theDetector, x_layer_tmp.regionStr(), x_layer_tmp.limitsStr(), x_layer_tmp.visStr());
 
-      chamber_y_offset = -(-Hcal_total_dim_y/2. 
-			   + (logical_layer_id-1) *(Hcal_chamber_thickness + Hcal_radiator_thickness)
-			   + Hcal_radiator_thickness + Hcal_chamber_thickness/2.);
+    double layer_thickness = y_halfheight * 2.;
 
+    double nRadiationLengths = 0.;
+    double nInteractionLengths = 0.;
+    double thickness_sum = 0;
 
-      pv =  EnvLogHcalModuleBarrel.placeVolume(ChamberLogical,
-					       Position(chamber_x_offset,
-							chamber_z_offset,
-							chamber_y_offset + chambers_y_off_correction));
-      
+    nRadiationLengths = Hcal_radiator_thickness / (stavesMaterial.radLength());
+    nInteractionLengths = Hcal_radiator_thickness / (stavesMaterial.intLength());
+    thickness_sum = Hcal_radiator_thickness;
 
+    //====================================================================
+    // Create Hcal Barrel Chamber without radiator
+    // Place into the Hcal Barrel stave, after each radiator
+    //====================================================================
+    xml_coll_t c(x_det, _U(layer));
+    xml_comp_t x_layer = c;
+    string layer_name = det_name + _toString(layer_id, "_layer%d");
 
-      //-----------------------------------------------------------------------------------------
-      if( layer_id <= Hcal_nlayers ){ // add the first set of layers to the reconstruction data
-	
-	caloLayer.distance = Hcal_inner_radius + Hcal_total_dim_y/2.0 - (chamber_y_offset + chambers_y_off_correction)
-	  - caloLayer.inner_thickness ; // Will be added later at "DDMarlinPandora/DDGeometryCreator.cc:226" to get center of sensitive element
+    // Create the slices (sublayers) within the Hcal Barrel Chamber.
+    double slice_pos_z = layer_thickness / 2.;
+    int slice_number = 0;
 
-	caloLayer.absorberThickness = Hcal_radiator_thickness ;
-	
-	caloData->layers.push_back( caloLayer ) ;
+    for (xml_coll_t k(x_layer, _U(slice)); k; ++k)
+    {
+      xml_comp_t x_slice = k;
+      string slice_name = layer_name + _toString(slice_number, "_slice%d");
+      double slice_thickness = x_slice.thickness();
+      Material slice_material = theDetector.material(x_slice.materialStr());
+      DetElement slice(layer_name, _toString(slice_number, "slice%d"), x_det.id());
+
+      slice_pos_z -= slice_thickness / 2.;
+
+      // Slice volume & box
+      Volume slice_vol(slice_name, Box(x_halflength, z_halfwidth, slice_thickness / 2.), slice_material);
+      slice_vol.setAttributes(theDetector, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
+
+      nRadiationLengths += slice_thickness / (2. * slice_material.radLength());
+      nInteractionLengths += slice_thickness / (2. * slice_material.intLength());
+      thickness_sum += slice_thickness / 2;
+      if (x_slice.isSensitive())
+      {
+        // Store "inner" quantities
+        caloLayer.inner_nRadiationLengths = nRadiationLengths;
+        caloLayer.inner_nInteractionLengths = nInteractionLengths;
+        caloLayer.inner_thickness = thickness_sum;
+        // Store scintillator thickness
+        caloLayer.sensitive_thickness = Hcal_scintillator_thickness;
+
+        double cell_size = 2 * Hcal_scintillator_ESR_thickness + Hcal_cell_size;
+        double cell_size_abnormal = 2 * Hcal_scintillator_ESR_thickness + Hcal_cell_size_abnormal;
+        double cell_x_offset = 0, cell_z_offset = 0, cell_y_offset = 0;
+        // place cell into slice one by one
+        double total_cell_size = cell_size + Hcal_scintillator_air_gap;                   // include air gap
+        double total_cell_size_abnormal = cell_size_abnormal + Hcal_scintillator_air_gap; // include air gap
+        int n_x = (2 * x_halflength) / total_cell_size;
+        int n_z = (2 * z_halfwidth) / total_cell_size;
+        int nx_abnormal = (2 * x_halflength - n_x * total_cell_size) / (total_cell_size_abnormal - total_cell_size);
+        n_x -= nx_abnormal;
+        // if(layer_id==1){n_x_layer_1=n_x;}
+        // else{
+        //   n_x=(n_x-n_x_layer_1)/2;
+        //   n_x=n_x*2+n_x_layer_1;
+        // }
+        for (int index_cell_z = 0; index_cell_z < n_z; index_cell_z++)
+        {
+          cell_x_offset = (n_x * total_cell_size + nx_abnormal * total_cell_size_abnormal) / 2.;
+          cell_z_offset = (n_z / 2. - index_cell_z - 0.5) * total_cell_size;
+          for (int index_cell_x = 0; index_cell_x < n_x; index_cell_x++)
+          {
+            cell_x_offset -= total_cell_size / 2.;
+            int cell_number = index_cell_x + 100 * index_cell_z;
+            string cell_name = slice_name + _toString(cell_number, "_cell%d");
+            string scintillator_name = cell_name + "_scintillator";
+            // cell_y_offset = slice_thickness / 2.;
+            PlacedVolume cell_phv = slice_vol.placeVolume(cell_vol, Position(cell_x_offset, cell_z_offset, cell_y_offset));
+            cell_phv.addPhysVolID("tile", cell_number);
+
+            DetElement cell(cell_name, _toString(cell_number, "cell%d"), cell_number);
+            cell.setPlacement(cell_phv);
+            cell_x_offset -= total_cell_size / 2.;
+          }
+          for (int index_cell_x = 0; index_cell_x < nx_abnormal; index_cell_x++)
+          {
+            cell_x_offset -= total_cell_size_abnormal / 2.;
+            int cell_number = index_cell_x + n_x + 100 * index_cell_z;
+            string cell_name = slice_name + _toString(cell_number, "_cell%d");
+            string scintillator_name = cell_name + "_scintillator";
+            // cell_y_offset = slice_thickness / 2.;
+            PlacedVolume cell_phv = slice_vol.placeVolume(cell_vol_abnormal, Position(cell_x_offset, cell_z_offset, cell_y_offset));
+            cell_phv.addPhysVolID("tile", cell_number);
+
+            DetElement cell(cell_name, _toString(cell_number, "cell%d"), cell_number);
+            cell.setPlacement(cell_phv);
+            cell_x_offset -= total_cell_size_abnormal / 2.;
+          }
+        }
+        //printf("layerID: %2d, x_length: %3.3lf, x_cell: %3.3lf, x_cell_abnormal: %3.3lf, x_dead: %3.3lf, z_length: %3.3lf,z_cell: %3.3lf, z_dead: %3.3lf,\n", layer_id,  2 * x_halflength,  n_x * total_cell_size,  nx_abnormal * total_cell_size_abnormal,  2 * x_halflength - n_x * total_cell_size - nx_abnormal * total_cell_size_abnormal,  2 * z_halfwidth,  n_z * total_cell_size,  2 * z_halfwidth - n_z * total_cell_size);
+        // Reset counters to measure "outside" quantitites
+        nRadiationLengths = 0.;
+        nInteractionLengths = 0.;
+        thickness_sum = 0.;
+        // }
       }
-      //-----------------------------------------------------------------------------------------
-
-      
-    }//end loop over HCAL nlayers;
-
-  if( tileSeg !=0 ){
-    // check the offsets directly in the TileSeg ...
-    std::vector<double> LOX = tileSeg->layerOffsetX();
-    std::vector<double> LOY = tileSeg->layerOffsetY();
-
-    std::stringstream sts ;
-    sts <<" layerOffsetX(): ";
-    for (std::vector<double>::const_iterator ix = LOX.begin(); ix != LOX.end(); ++ix)
-      sts << *ix << ' ';
-    printout( dd4hep::DEBUG,  "SHcalSc04_Barrel_v04", "%s" , sts.str().c_str() ) ;
-    sts.clear() ; sts.str("") ;
-    sts <<" layerOffsetY(): ";
-    for (std::vector<double>::const_iterator iy = LOY.begin(); iy != LOY.end(); ++iy)
-      sts << *iy << ' ';
-    printout( dd4hep::DEBUG,  "SHcalSc04_Barrel_v04", "%s" , sts.str().c_str() ) ;
+      nRadiationLengths += slice_thickness / (2. * slice_material.radLength());
+      nInteractionLengths += slice_thickness / (2. * slice_material.intLength());
+      thickness_sum += slice_thickness / 2;
+
+      // Set region, limitset, and vis.
+      // slice_vol.setAttributes(theDetector, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
+      // slice PlacedVolume
+      PlacedVolume slice_phv = ChamberLogical.placeVolume(slice_vol, Position(0., 0., slice_pos_z));
+
+      slice_phv.addPhysVolID("layer", layer_id);
+      slice.setPlacement(slice_phv);
+      // Increment z position for next slice.
+      slice_pos_z -= slice_thickness / 2.;
+      // Increment slice number.
+      ++slice_number;
+    }
 
-  }
+    // Store "outer" quantities
+    caloLayer.outer_nRadiationLengths = nRadiationLengths;
+    caloLayer.outer_nInteractionLengths = nInteractionLengths;
+    caloLayer.outer_thickness = thickness_sum;
 
+    //---------------------------  Chamber Placements -----------------------------------------
 
+    double chamber_x_offset, chamber_y_offset, chamber_z_offset;
+    chamber_x_offset = xShift;
 
-//====================================================================
-// Place HCAL Barrel stave module into the envelope
-//====================================================================
-  double stave_phi_offset,  module_z_offset,  lateral_plate_z_offset;
+    chamber_z_offset = 0;
 
-  double Y = Hcal_inner_radius + Hcal_total_dim_y / 2.;
+    chamber_y_offset = -(-Hcal_total_dim_y / 2. + (layer_id - 1) * (Hcal_chamber_thickness + Hcal_radiator_thickness) + Hcal_radiator_thickness + Hcal_chamber_thickness / 2.);
 
-  stave_phi_offset = M_PI/Hcal_inner_symmetry -M_PI/2.;
+    pv = EnvLogHcalModuleBarrel.placeVolume(ChamberLogical,
+                                            Position(chamber_x_offset, chamber_z_offset, chamber_y_offset));
 
+    //-----------------------------------------------------------------------------------------
+    if (layer_id <= Hcal_nlayers)
+    { // add the first set of layers to the reconstruction data
 
-  //-------- start loop over HCAL BARREL staves ----------------------------
+      caloLayer.distance = Hcal_inner_radius + Hcal_total_dim_y / 2.0 - chamber_y_offset - caloLayer.inner_thickness;
+      // Will be added later at "DDMarlinPandora/DDGeometryCreator.cc:226" to get center of sensitive element
 
-  for (int stave_id = 1;
-       stave_id <=Hcal_inner_symmetry;
-       stave_id++)
-    {
-      module_z_offset = - (Hcal_normal_dim_z + Hcal_modules_gap + Hcal_lateral_plate_thickness)/2.;
-      lateral_plate_z_offset = - (Hcal_lateral_plate_thickness + Hcal_modules_gap)/2.;
+      caloLayer.absorberThickness = Hcal_radiator_thickness;
 
-      double phirot = stave_phi_offset;
-      RotationZYX srot(0,phirot,M_PI*0.5);
-      Rotation3D srot3D(srot);
+      caloData->layers.push_back(caloLayer);
+    }
+    //-----------------------------------------------------------------------------------------
 
-      for (int module_id = 1;
-         module_id <=2;
-         module_id++)
-	{
-	  Position sxyzVec(-Y*sin(phirot), Y*cos(phirot), module_z_offset);
-	  Transform3D stran3D(srot3D,sxyzVec);
-	  
-	  // Place Hcal Barrel volume into the envelope volume
-	  pv = envelope.placeVolume(EnvLogHcalModuleBarrel,stran3D);
-	  pv.addPhysVolID("stave",stave_id);
-	  pv.addPhysVolID("module",module_id);
-	  pv.addPhysVolID("system",det_id);
+  } // end loop over HCAL nlayers;
 
-	  const int staveCounter = (stave_id-1)*2+module_id-1;
-	  DetElement stave(sdet, _toString(staveCounter,"stave%d"), staveCounter );
-	  stave.setPlacement(pv);
+  //====================================================================
+  // Place HCAL Barrel stave module into the envelope
+  //====================================================================
+  double stave_phi_offset, module_z_offset = 0;
 
-	  Position xyzVec_LP(-Y*sin(phirot), Y*cos(phirot),lateral_plate_z_offset);
-	  Transform3D tran3D_LP(srot3D,xyzVec_LP);
-	  pv = envelope.placeVolume(EnvLogHcalModuleBarrel_LP,tran3D_LP);
+  double Y = Hcal_inner_radius + Hcal_total_dim_y / 2.;
 
-	  module_z_offset = - module_z_offset;
-	  lateral_plate_z_offset = - lateral_plate_z_offset;
-	}
+  stave_phi_offset = M_PI / Hcal_inner_symmetry - M_PI / 2.;
 
+  //-------- start loop over HCAL BARREL staves ----------------------------
 
-      stave_phi_offset -=  M_PI*2.0/Hcal_inner_symmetry;
-    }  //-------- end loop over HCAL BARREL staves ----------------------------
+  for (int stave_id = 1; stave_id <= Hcal_inner_symmetry; stave_id++)
+  {
+    double phirot = stave_phi_offset;
+    RotationZYX srot(0, phirot, M_PI * 0.5);
+    Rotation3D srot3D(srot);
 
+    // for (int module_id = 1; module_id <= 2; module_id++)
+    // {
+    Position sxyzVec(-Y * sin(phirot), Y * cos(phirot), module_z_offset);
+    Transform3D stran3D(srot3D, sxyzVec);
 
-  sdet.addExtension< LayeredCalorimeterData >( caloData ) ;
+    // Place Hcal Barrel volume into the envelope volume
+    pv = envelope.placeVolume(EnvLogHcalModuleBarrel, stran3D);
+    pv.addPhysVolID("stave", stave_id);
+    // pv.addPhysVolID("module", module_id);
+    pv.addPhysVolID("system", det_id);
 
- 
-  return sdet;
+    // const int staveCounter = (stave_id - 1) * 2 + module_id - 1;
+    const int staveCounter = stave_id - 1;
+    DetElement stave(sdet, _toString(staveCounter, "stave%d"), staveCounter);
+    stave.setPlacement(pv);
 
+    // Position xyzVec_LP(-Y * sin(phirot), Y * cos(phirot), lateral_plate_z_offset);
+    // Transform3D tran3D_LP(srot3D, xyzVec_LP);
+    // pv = envelope.placeVolume(EnvLogHcalModuleBarrel_LP, tran3D_LP);
+
+    // module_z_offset = -module_z_offset;
+    // lateral_plate_z_offset = -lateral_plate_z_offset;
+    // }
+
+    stave_phi_offset -= M_PI * 2.0 / Hcal_inner_symmetry;
+  } //-------- end loop over HCAL BARREL staves ----------------------------
+
+  sdet.addExtension<LayeredCalorimeterData>(caloData);
+
+  return sdet;
 }
 
 DECLARE_DETELEMENT(SHcalSc04_Barrel_v04, create_detector)

Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Endcaps_v02.cpp

+//================================================================================
+// Description — End-Cap AHCAL
+//--------------------------------------------------------------------------------
+// Author: Ji-Yuan CHEN (SJTU; jy_chen@sjtu.edu.cn)
+//--------------------------------------------------------------------------------
+//
+// End-cap AHCAL with 40×40×3 mm³ scintillator tiles.
+// Adding the dead materials (cassette, PCB and electronic components, ESR wrapper and steel absorber) and including the air gaps, the size of each cell becomes 40.3×40.3×27.2 mm³.
+//
+// The inner radius, end-cap thickness, unit size, etc. are directly read from XML file.
+//
+// Structure in a layer: Cassette → ESR → scintillator → ESR → PCB (including electronic components) → cassette → steel absorber.
+//
+// Default layout: The absorber is designed as a whole with some space left for the sensitive units (the cross-section is not a polygon); steel frame on the edges of each module; 16 modules (called 'staves' in this program), 48 layers; 72 rows in r direction, with uniform granularity.  For a schematic diagram, see Page 3 of  <https://indico.ihep.ac.cn/event/22010/contributions/153187/attachments/77666/96430/2024_0324_Calorimeter_Endcaps.pdf>  (the design on the left; the numbers have been modified).
+//================================================================================
+
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/Shapes.h"
+#include "DD4hep/DetType.h"
+#include "XML/Layering.h"
+#include "XML/Utilities.h"
+#include "DDRec/DetectorData.h"
+#include "DDSegmentation/Segmentation.h"
+#include <sstream>
+#include <cassert>
+
+using std::cout;
+using std::endl;
+using std::string;
+using std::vector;
+using std::to_string;
+
+#define MYDEBUG(x) cout << __FILE__ << ":" << __LINE__ << ": " << x << endl;
+#define MYDEBUGVAL(x) cout << __FILE__ << ":" << __LINE__ << ": "<< #x << ": " << x << endl;
+
+using dd4hep::Ref_t;
+using dd4hep::Detector;
+using dd4hep::SensitiveDetector;
+using dd4hep::pi;
+using dd4hep::mm;
+using dd4hep::Material;
+using dd4hep::DetElement;
+using dd4hep::Volume;
+using dd4hep::PolyhedraRegular;
+using dd4hep::Tube;
+using dd4hep::UnionSolid;
+using dd4hep::Trapezoid;
+using dd4hep::Box;
+using dd4hep::SubtractionSolid;
+using dd4hep::PlacedVolume;
+using dd4hep::Position;
+using dd4hep::RotationX;
+using dd4hep::_toString;
+using dd4hep::Transform3D;
+using dd4hep::RotationZ;
+using dd4hep::RotationY;
+
+static Ref_t create_detector(Detector& theDetector, xml_h e, SensitiveDetector sens)
+{
+    xml_det_t x_det = e;
+
+    const string det_name = x_det.nameStr();
+    const string det_type = x_det.typeStr();
+    const int det_id = x_det.id();
+    MYDEBUGVAL(det_name)
+    MYDEBUGVAL(det_type)
+    MYDEBUGVAL(det_id)
+
+    // To prevent overlapping
+    const double boundary_safety = theDetector.constant<double>("boundary_safety");
+
+    // Global geometry
+    const double r_in = theDetector.constant<double>("hcalendcap_inner_radius");
+    const double r_out = theDetector.constant<double>("hcalendcap_outer_radius");
+    const double module_thickness = theDetector.constant<double>("hcalendcap_thickness");
+    const double pos_z = theDetector.constant<double>("hcalendcap_z");
+    const int Nmodules = theDetector.constant<int>("Nmodules");
+    const double angle = 2 * pi / Nmodules;
+
+    // Unit size
+    const double scintillator_xy = theDetector.constant<double>("scintillator_xy");
+    const double scintillator_z = theDetector.constant<double>("scintillator_z");
+    const double wrapped_scintillator_xy = theDetector.constant<double>("wrapped_scintillator_xy");
+    const double wrapped_scintillator_z = theDetector.constant<double>("wrapped_scintillator_z");
+
+    // Dead materials
+    const double cassette_thickness = theDetector.constant<double>("cassette_thickness"); 
+    const double esr_thickness = theDetector.constant<double>("esr_thickness");    // Wrapper
+    [[maybe_unused]] const double sipm_xy = theDetector.constant<double>("sipm_xy");
+    [[maybe_unused]] const double sipm_z = theDetector.constant<double>("sipm_z");
+    const double pcb_thickness = theDetector.constant<double>("pcb_thickness"); 
+    const double absorber_thickness = theDetector.constant<double>("absorber_thickness"); 
+    const double air_gap_xy = 0.5 * (wrapped_scintillator_xy - scintillator_xy) - esr_thickness;
+    [[maybe_unused]] const double air_gap_z = 0.5 * (wrapped_scintillator_z - scintillator_z) - esr_thickness;
+
+    // Odd-shaped cells
+    const int Nodd = theDetector.constant<int>("Nodd");
+    const double short_elongation_1 = theDetector.constant<double>("short_elongation_1");
+    const double short_elongation_2 = theDetector.constant<double>("short_elongation_2");
+    const double short_elongation_3 = theDetector.constant<double>("short_elongation_3");
+    const double short_elongation_4 = theDetector.constant<double>("short_elongation_4");
+    const double short_elongation_5 = theDetector.constant<double>("short_elongation_5");
+
+    const double elongation_angle_esr_out = wrapped_scintillator_xy * tan(0.5 * angle);    // Elongation of the outer edge of ESR: long side - short side
+    const double elongation_angle_esr_in = (wrapped_scintillator_xy - esr_thickness / cos(0.5 * angle)) * tan(0.5 * angle);    // Elongation of the inner edge of ESR: long side - short side
+    const double elongation_angle_sc = scintillator_xy * tan(0.5 * angle);    // Elongation of the scintillator: long side - short side
+    const vector<double> short_elongation_esr_out = { short_elongation_1, short_elongation_2, short_elongation_3, short_elongation_4, short_elongation_5 };
+    assert(Nodd == short_elongation_esr_out.size());
+
+    vector<double> short_elongation_sc(short_elongation_esr_out.size()), long_elongation_sc(short_elongation_esr_out.size()), short_elongation_esr_in(short_elongation_esr_out.size()), long_elongation_esr_in(short_elongation_esr_out.size()), long_elongation_esr_out(short_elongation_esr_out.size());
+    for (int i = 0; i < Nodd; ++i)
+    {
+        short_elongation_sc.at(i) = short_elongation_esr_out.at(i) - (air_gap_xy + esr_thickness) / cos(0.5 * angle);
+        long_elongation_sc.at(i) = short_elongation_sc.at(i) + elongation_angle_sc;
+        short_elongation_esr_in.at(i) = short_elongation_esr_out.at(i) - esr_thickness / cos(0.5 * angle);
+        long_elongation_esr_in.at(i) = short_elongation_esr_in.at(i) + elongation_angle_esr_in;
+        long_elongation_esr_out.at(i) = short_elongation_esr_out.at(i) + elongation_angle_esr_out;
+    }
+
+    // Mechanical structure
+    const double inner_structure_thickness = theDetector.constant<double>("inner_structure_thickness");
+    [[maybe_unused]] const double outer_structure_width = theDetector.constant<double>("outer_structure_width");
+    [[maybe_unused]] const double outer_structure_thickness = theDetector.constant<double>("outer_structure_thickness");
+    const double frame_thickness = theDetector.constant<double>("frame_thickness");
+
+    const double layer_thickness = wrapped_scintillator_z + pcb_thickness + 2 * cassette_thickness + absorber_thickness + 5 * boundary_safety;
+    const double non_absorber_thickness = wrapped_scintillator_z + pcb_thickness + 3 * boundary_safety;
+
+    MYDEBUGVAL(layer_thickness)
+
+    // Module size
+    const double dim_in_frame = (r_in + inner_structure_thickness + frame_thickness) * tan(0.5 * angle);
+    const double dim_out_frame = r_out * sin(0.5 * angle);
+    const double dim_in = dim_in_frame - frame_thickness / cos(0.5 * angle);
+    const double dim_out = dim_out_frame - frame_thickness / cos(0.5 * angle);
+    const double height = r_out * cos(0.5 * angle) - r_in - inner_structure_thickness;
+    const double r0 = r_in + inner_structure_thickness + 0.5 * height;    // Rotation radius for placing the modules
+
+    const int Nrows = (int) (height / (wrapped_scintillator_xy + boundary_safety));
+    const int Nlayers = (int) (module_thickness / layer_thickness);
+    int Ncells_phi;
+
+    MYDEBUGVAL(Nrows)
+    MYDEBUGVAL(Nlayers)
+
+    // Materials
+    Material mat_air(theDetector.material("Air"));
+    Material mat_PCB(theDetector.material("PCB"));
+    [[maybe_unused]] Material mat_SiPM(theDetector.material("G4_Si"));
+    Material mat_sensitive(theDetector.material( x_det.materialStr() ));
+    Material mat_ESR(theDetector.material("G4_ESR"));
+    Material mat_steel(theDetector.material("Steel235"));
+
+    // The detector and mother volumes (world)
+    DetElement AHCAL(det_name, det_id);
+    Volume motherVol = theDetector.pickMotherVolume(AHCAL);
+
+    // Create two tube-like envelopes to represent the end-cap volumes
+//    PolyhedraRegular envelope_side(Nmodules, 0.5 * angle, r_in + inner_structure_thickness, r_out, module_thickness);
+    Tube envelope_side(r_in, r_out, 0.5 * module_thickness);
+    UnionSolid envelope(envelope_side, envelope_side, Position(0, 0, 2 * pos_z));
+    Volume envelopeVol(det_name, envelope, mat_steel);
+    PlacedVolume envelopePlv = motherVol.placeVolume(envelopeVol, Position(0, 0, -pos_z));
+    envelopePlv.addPhysVolID("system", det_id);
+    envelopeVol.setVisAttributes(theDetector, "SeeThrough");
+    AHCAL.setPlacement(envelopePlv);
+    DetElement stave_det(AHCAL, "sector", det_id);
+
+    // Module (stave)
+    Trapezoid stave(dim_in, dim_out, 0.5 * module_thickness, 0.5 * module_thickness, 0.5 * height);
+    Volume stave_vol("stave_vol", stave, mat_steel);
+    stave_vol.setVisAttributes(theDetector, "BlueVis");
+
+    // Layer with only wrapped scintillators, electronic components and PCB
+    Trapezoid layer(dim_in, dim_out, 0.5 * non_absorber_thickness, 0.5 * non_absorber_thickness, 0.5 * height);
+    Volume layer_vol("layer_vol", layer, mat_steel);
+    layer_vol.setVisAttributes(theDetector, "SeeThrough");
+
+    // Scintillator, ESR, SiPM, PCB, cassette, and absorber
+    Box normal_scintillator(0.5 * scintillator_xy, 0.5 * scintillator_z, 0.5 * scintillator_xy);
+    Volume scintillator("scintillator", normal_scintillator, mat_sensitive);    // Order: phi → z → r
+    scintillator.setVisAttributes(theDetector, "SeeThrough");
+    scintillator.setSensitiveDetector(sens);
+
+    [[maybe_unused]] Volume sipm("SiPM", Box(0.5 * sipm_xy, 0.5 * sipm_xy, 0.5 * sipm_z), mat_SiPM);
+    sipm.setVisAttributes(theDetector, "CyanVis");
+
+    Box esr_out(0.5 * wrapped_scintillator_xy, 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_xy);
+    Box esr_in(0.5 * wrapped_scintillator_xy - esr_thickness, 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_xy - esr_thickness);
+    Volume esr("esr", SubtractionSolid(esr_out, esr_in, Position(0, 0, 0)), mat_ESR);
+    esr.setVisAttributes(theDetector, "SeeThrough");
+
+    // Odd-shaped scintillators
+    Trapezoid odd_add_0_scintillator(0, long_elongation_sc.at(0) - short_elongation_sc.at(0), 0.5 * scintillator_z, 0.5 * scintillator_z, 0.5 * scintillator_xy);
+    Trapezoid odd_add_45_scintillator(short_elongation_sc.at(1), long_elongation_sc.at(1), 0.5 * scintillator_z, 0.5 * scintillator_z, 0.5 * scintillator_xy);
+    Trapezoid odd_add_9_scintillator(short_elongation_sc.at(2), long_elongation_sc.at(2), 0.5 * scintillator_z, 0.5 * scintillator_z, 0.5 * scintillator_xy);
+    Trapezoid odd_add_14_scintillator(short_elongation_sc.at(3), long_elongation_sc.at(3), 0.5 * scintillator_z, 0.5 * scintillator_z, 0.5 * scintillator_xy);
+    Trapezoid odd_add_17_scintillator(short_elongation_sc.at(4), long_elongation_sc.at(4), 0.5 * scintillator_z, 0.5 * scintillator_z, 0.5 * scintillator_xy);
+
+    Box odd_scintillator_0_rect(0.5 * (scintillator_xy + short_elongation_sc.at(0)), 0.5 * scintillator_z, 0.5 * scintillator_xy);
+
+    Volume odd_0_scintillator("odd_0_scintillator", UnionSolid(odd_scintillator_0_rect, odd_add_0_scintillator, Position(0.5 * (scintillator_xy + short_elongation_sc.at(0)), 0, 0)), mat_sensitive);
+    odd_0_scintillator.setVisAttributes(theDetector, "SeeThrough");
+    odd_0_scintillator.setSensitiveDetector(sens);
+
+    Volume odd_45_scintillator("odd_45_scintillator", UnionSolid(normal_scintillator, odd_add_45_scintillator, Position(0.5 * scintillator_xy, 0, 0)), mat_sensitive);
+    odd_45_scintillator.setVisAttributes(theDetector, "SeeThrough");
+    odd_45_scintillator.setSensitiveDetector(sens);
+
+    Volume odd_9_scintillator("odd_9_scintillator", UnionSolid(normal_scintillator, odd_add_9_scintillator, Position(0.5 * scintillator_xy, 0, 0)), mat_sensitive);
+    odd_9_scintillator.setVisAttributes(theDetector, "SeeThrough");
+    odd_9_scintillator.setSensitiveDetector(sens);
+
+    Volume odd_14_scintillator("odd_14_scintillator", UnionSolid(normal_scintillator, odd_add_14_scintillator, Position(0.5 * scintillator_xy, 0, 0)), mat_sensitive);
+    odd_14_scintillator.setVisAttributes(theDetector, "SeeThrough");
+    odd_14_scintillator.setSensitiveDetector(sens);
+
+    Volume odd_17_scintillator("odd_17_scintillator", UnionSolid(normal_scintillator, odd_add_17_scintillator, Position(0.5 * scintillator_xy, 0, 0)), mat_sensitive);
+    odd_17_scintillator.setVisAttributes(theDetector, "SeeThrough");
+    odd_17_scintillator.setSensitiveDetector(sens);
+
+    // Odd-shaped ESR
+    // Outer edge
+    Trapezoid odd_add_0_esr_out(short_elongation_esr_out.at(0), long_elongation_esr_out.at(0), 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_xy);
+    Trapezoid odd_add_45_esr_out(short_elongation_esr_out.at(1), long_elongation_esr_out.at(1), 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_xy);
+    Trapezoid odd_add_9_esr_out(short_elongation_esr_out.at(2), long_elongation_esr_out.at(2), 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_xy);
+    Trapezoid odd_add_14_esr_out(short_elongation_esr_out.at(3), long_elongation_esr_out.at(3), 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_xy);
+    Trapezoid odd_add_17_esr_out(short_elongation_esr_out.at(4), long_elongation_esr_out.at(4), 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_z, 0.5 * wrapped_scintillator_xy);
+
+    UnionSolid odd_0_esr_out(esr_out, odd_add_0_esr_out, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+    UnionSolid odd_45_esr_out(esr_out, odd_add_45_esr_out, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+    UnionSolid odd_9_esr_out(esr_out, odd_add_9_esr_out, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+    UnionSolid odd_14_esr_out(esr_out, odd_add_14_esr_out, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+    UnionSolid odd_17_esr_out(esr_out, odd_add_17_esr_out, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+
+    // Inner edge
+    Trapezoid odd_add_0_esr_in(0, long_elongation_esr_in.at(0) - short_elongation_esr_in.at(0), 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_xy - esr_thickness);
+    Trapezoid odd_add_45_esr_in(short_elongation_esr_in.at(1), long_elongation_esr_in.at(1), 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_xy - esr_thickness);
+    Trapezoid odd_add_9_esr_in(short_elongation_esr_in.at(2), long_elongation_esr_in.at(2), 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_xy - esr_thickness);
+    Trapezoid odd_add_14_esr_in(short_elongation_esr_in.at(3), long_elongation_esr_in.at(3), 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_xy - esr_thickness);
+    Trapezoid odd_add_17_esr_in(short_elongation_esr_in.at(4), long_elongation_esr_in.at(4), 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_xy - esr_thickness);
+
+    Box odd_esr_in_0_rect(0.5 * (wrapped_scintillator_xy - 2 * esr_thickness + short_elongation_esr_in.at(0)), 0.5 * wrapped_scintillator_z - esr_thickness, 0.5 * wrapped_scintillator_xy - esr_thickness);
+
+    UnionSolid odd_0_esr_in(odd_esr_in_0_rect, odd_add_0_esr_in, Position(0.5 * (wrapped_scintillator_xy + short_elongation_esr_in.at(0)), 0, 0));
+    UnionSolid odd_45_esr_in(esr_in, odd_add_45_esr_in, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+    UnionSolid odd_9_esr_in(esr_in, odd_add_9_esr_in, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+    UnionSolid odd_14_esr_in(esr_in, odd_add_14_esr_in, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+    UnionSolid odd_17_esr_in(esr_in, odd_add_17_esr_in, Position(0.5 * wrapped_scintillator_xy, 0, 0));
+
+    // Odd-shaped ESR, from the subtraction of outer and inner frames
+    Volume odd_0_esr("odd_0_esr", SubtractionSolid(odd_0_esr_out, odd_0_esr_in, Position(0, 0, 0)), mat_ESR);
+    Volume odd_45_esr("odd_45_esr", SubtractionSolid(odd_45_esr_out, odd_45_esr_in, Position(0, 0, 0)), mat_ESR);
+    Volume odd_9_esr("odd_9_esr", SubtractionSolid(odd_9_esr_out, odd_9_esr_in, Position(0, 0, 0)), mat_ESR);
+    Volume odd_14_esr("odd_14_esr", SubtractionSolid(odd_14_esr_out, odd_14_esr_in, Position(0, 0, 0)), mat_ESR);
+    Volume odd_17_esr("odd_17_esr", SubtractionSolid(odd_17_esr_out, odd_17_esr_in, Position(0, 0, 0)), mat_ESR);
+
+    // Positions
+    const double scintillator_pos_z = -0.5 * non_absorber_thickness + boundary_safety + 0.5 * wrapped_scintillator_z;
+    const double esr_pos_z = scintillator_pos_z;
+    const double pcb_pos_z = esr_pos_z + 0.5 * wrapped_scintillator_z + boundary_safety + 0.5 * pcb_thickness;
+
+    // Loop for placing the units in a layer
+    for (int ir = 1; ir <= Nrows; ++ir)
+    {
+        const double layer_length = dim_in + (ir - 1) * wrapped_scintillator_xy * tan(0.5 * angle);
+        Ncells_phi = (int) (2 * layer_length / wrapped_scintillator_xy);
+        const double layer_phi = Ncells_phi * wrapped_scintillator_xy;
+        const double single_elongation = layer_length - 0.5 * layer_phi;
+
+        Volume slice("slice", Trapezoid(layer_length, layer_length + elongation_angle_esr_out, 0.5 * non_absorber_thickness, 0.5 * non_absorber_thickness, 0.5 * (wrapped_scintillator_xy + boundary_safety)), mat_air);
+        slice.setVisAttributes(theDetector, "SeeThrough");
+        string slicename = "Slice_" + to_string(ir);
+        DetElement sd(stave_det, slicename, det_id);
+
+        Volume slice_pcb("slice_pcb", Box(0.5 * layer_phi, 0.5 * pcb_thickness, 0.5 * wrapped_scintillator_xy), mat_PCB);
+        slice_pcb.setVisAttributes(theDetector, "SeeThrough");
+
+        PlacedVolume pcb_unit = slice.placeVolume(slice_pcb, Position(0, pcb_pos_z, 0));
+        pcb_unit.addPhysVolID("row", ir);
+
+        for (int iphi = 1; iphi <= Ncells_phi; ++iphi)
+        {
+            PlacedVolume scintillator_unit;
+            PlacedVolume esr_unit;
+
+            Position position_scintillator((-0.5 * (Ncells_phi + 1) + iphi) * wrapped_scintillator_xy, scintillator_pos_z, 0);
+            Position position_esr((-0.5 * (Ncells_phi + 1) + iphi) * wrapped_scintillator_xy, esr_pos_z, 0);
+
+            Transform3D transform_scintillator(RotationZ(pi), position_scintillator);
+            Transform3D transform_esr(RotationZ(pi), position_esr);
+
+            if (iphi == 1)
+            {
+                if (single_elongation >= short_elongation_esr_out.at(4))
+                {
+                    scintillator_unit = slice.placeVolume(odd_17_scintillator, transform_scintillator);
+                    esr_unit = slice.placeVolume(odd_17_esr, transform_esr);
+                }
+                else if (single_elongation >= short_elongation_esr_out.at(3))
+                {
+                    scintillator_unit = slice.placeVolume(odd_14_scintillator, transform_scintillator);
+                    esr_unit = slice.placeVolume(odd_14_esr, transform_esr);
+                }
+                else if (single_elongation >= short_elongation_esr_out.at(2))
+                {
+                    scintillator_unit = slice.placeVolume(odd_9_scintillator, transform_scintillator);
+                    esr_unit = slice.placeVolume(odd_9_esr, transform_esr);
+                }
+                else if (single_elongation >= short_elongation_esr_out.at(1))
+                {
+                    scintillator_unit = slice.placeVolume(odd_45_scintillator, transform_scintillator);
+                    esr_unit = slice.placeVolume(odd_45_esr, transform_esr);
+                }
+                else if (single_elongation >= short_elongation_esr_out.at(0))
+                {
+                    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,
+                                         0)));
+                    esr_unit = slice.placeVolume(odd_0_esr, transform_esr);
+                }
+            }
+            else if (iphi == Ncells_phi)
+            {
+                if (single_elongation >= short_elongation_esr_out.at(4))
+                {
+                    scintillator_unit = slice.placeVolume(odd_17_scintillator, position_scintillator);
+                    esr_unit = slice.placeVolume(odd_17_esr, position_esr);
+                }
+                else if (single_elongation >= short_elongation_esr_out.at(3))
+                {
+                    scintillator_unit = slice.placeVolume(odd_14_scintillator, position_scintillator);
+                    esr_unit = slice.placeVolume(odd_14_esr, position_esr);
+                }
+                else if (single_elongation >= short_elongation_esr_out.at(2))
+                {
+                    scintillator_unit = slice.placeVolume(odd_9_scintillator, position_scintillator);
+                    esr_unit = slice.placeVolume(odd_9_esr, position_esr);
+                }
+                else if (single_elongation >= short_elongation_esr_out.at(1))
+                {
+                    scintillator_unit = slice.placeVolume(odd_45_scintillator, position_scintillator);
+                    esr_unit = slice.placeVolume(odd_45_esr, position_esr);
+                }
+                else if (single_elongation >= short_elongation_esr_out.at(0))
+                {
+                    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,
+                                         0));
+                    esr_unit = slice.placeVolume(odd_0_esr, position_esr);
+                }
+            }
+            else
+            {
+                scintillator_unit = slice.placeVolume(scintillator, position_scintillator);
+                esr_unit = slice.placeVolume(esr, position_esr);
+            }
+
+            scintillator_unit.addPhysVolID("row", ir).addPhysVolID("phi", iphi);
+            esr_unit.addPhysVolID("row", ir).addPhysVolID("phi", iphi);
+
+            string scintillator_name = "Scintillator_" + to_string(ir) + "_" + to_string(iphi);
+            DetElement unit(sd, scintillator_name, det_id);
+            unit.setPlacement(scintillator_unit);
+        }
+
+        PlacedVolume plv = layer_vol.placeVolume(slice, Position(0, 0, (-0.5 * (Nrows + 1) + ir) * (wrapped_scintillator_xy + boundary_safety)));
+        plv.addPhysVolID("row", ir);
+        sd.setPlacement(plv);
+    }
+
+    // Loop for placing the layers in a module
+    for (int ilayer = 1; ilayer <= Nlayers; ++ilayer)
+    {
+        PlacedVolume plv = stave_vol.placeVolume(layer_vol, Position(0, (ilayer - 1) * layer_thickness + cassette_thickness + 0.5 * non_absorber_thickness - 0.5 * module_thickness, 0));
+        plv.addPhysVolID("layer", ilayer);
+        DetElement sd(stave_det, _toString(ilayer, "layer_%3d"), det_id);
+        sd.setPlacement(plv);
+    }
+
+    // Loop for placing the modules
+    for (int i = 0; i < Nmodules; ++i)
+    {
+        const double m_rot = i * angle;
+        const double posx = -r0 * sin(m_rot);
+        const double posy = r0 * cos(m_rot);
+
+        Transform3D transform_neg(RotationZ(m_rot) * RotationX(-0.5 * pi), Position(posx, posy, 0));
+        Transform3D transform_pos(RotationZ(m_rot) * RotationY(pi) * RotationX(-0.5 * pi), Position(posx, posy, 2 * pos_z));
+        PlacedVolume plv_neg = envelopeVol.placeVolume(stave_vol, transform_neg);
+        PlacedVolume plv_pos = envelopeVol.placeVolume(stave_vol, transform_pos);
+        plv_neg.addPhysVolID("stave", i);
+        plv_pos.addPhysVolID("stave", i + Nmodules);
+        DetElement sd_neg(AHCAL, _toString(i, "sector%3d"), det_id);
+        DetElement sd_pos(AHCAL, _toString(i + Nmodules, "sector%3d"), det_id);
+        sd_neg.setPlacement(plv_neg);
+        sd_pos.setPlacement(plv_pos);
+    }
+
+    sens.setType("calorimeter");
+
+    MYDEBUG("create_detector FINISHED.");
+    return AHCAL;
+}
+
+DECLARE_DETELEMENT(SHcalSc04_Endcaps_v02, create_detector)

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