Digitization/DigiCalo/src/HcalDigiAlg.h

+#ifndef HCAL_DIGI_ALG_H 
+#define HCAL_DIGI_ALG_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/MutableCaloHitContribution.h"
+#include "edm4hep/MutableSimCalorimeterHit.h"
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/CalorimeterHitCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/MCRecoCaloAssociationCollection.h"
+#include "edm4hep/MCRecoCaloParticleAssociationCollection.h"
+#include "edm4hep/MCParticle.h"
+
+#include <DDRec/DetectorData.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include <DD4hep/Segmentations.h> 
+#include "DetInterface/IGeomSvc.h"
+#include "TVector3.h"
+#include "TRandom3.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TH3.h"
+#include "TH2.h"
+#include "TH1.h"
+#include "TF1.h"
+
+#include <cstdlib>
+#include "time.h"
+#include <TTimeStamp.h> 
+#include <ctime>
+
+#define C 299.79  // unit: mm/ns
+#define PI 3.141592653
+
+class HcalDigiAlg : public GaudiAlgorithm
+{
+ 
+public:
+ 
+  HcalDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+ 
+  /** Called at the begin of the job before anything is read.
+   * Use to initialize the processor, e.g. book histograms.
+   */
+  virtual StatusCode initialize() ;
+ 
+  /** Called for every event - the working horse.
+   */
+  virtual StatusCode execute() ; 
+ 
+  /** Called after data processing for clean up.
+   */
+  virtual StatusCode finalize() ;
+
+  StatusCode MergeHits(const edm4hep::SimCalorimeterHitCollection& m_col, std::vector<edm4hep::SimCalorimeterHit>& m_hits);
+  edm4hep::MutableSimCalorimeterHit find(const std::vector<edm4hep::MutableSimCalorimeterHit>& m_col, unsigned long long& cellid) const;
+
+	void Clear();
+
+protected:
+
+  SmartIF<IGeomSvc> m_geosvc;
+  typedef std::vector<float> FloatVec;
+  typedef std::map<const edm4hep::MCParticle, float> MCParticleToEnergyWeightMap;
+
+	int _nEvt ;
+	TRandom3 rndm;
+	TFile* m_wfile;
+	TTree* t_simHit;
+  TH2D* GSTileResMap; 
+
+  double m_totE, m_totE_truth;
+  double m_MC_EPx, m_MC_EPy, m_MC_EPz;
+	FloatVec m_simHit_x, m_simHit_y, m_simHit_z, m_simHit_E, m_simHit_Etruth, m_simHit_Eatt, 
+           m_simHit_system, m_simHit_stave, m_simHit_layer, m_simHit_tile, m_simHit_idx, m_simHit_idy, m_simHit_row, m_simHit_phi, 
+           m_simHit_HG, m_simHit_LG, m_simHit_rawQ, m_simHit_Npe_scint, m_simHit_Npe_sipm;
+  FloatVec m_step_x, m_step_y, m_step_LY;
+  std::vector<unsigned long long> m_simHit_cellID;
+  std::vector<int> m_simHit_steps;
+
+  dd4hep::Detector* m_dd4hep;
+	dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+	//dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+  std::map<std::string, dd4hep::DDSegmentation::BitFieldCoder*> map_readout_decoder;
+
+  TF1* f_DarkNoise = nullptr;
+
+	Gaudi::Property<float> m_scale{ this, "Scale", 1 };
+
+  // Input collections
+  typedef DataHandle<edm4hep::SimCalorimeterHitCollection>              SimCaloType;
+  typedef DataHandle<edm4hep::CalorimeterHitCollection>                 CaloType;
+  typedef DataHandle<edm4hep::MCRecoCaloAssociationCollection>          CaloSimAssoType;      //Calorimeter - SimCalorimeter
+  typedef DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>  CaloParticleAssoType; //Calorimeter - MCParticle
+
+  Gaudi::Property< std::vector<std::string> > name_SimCaloHit{ this, "SimCaloHitCollection", {"HcalBarrelCollection"} };
+  Gaudi::Property< std::vector<std::string> > name_Readout{ this, "ReadOutName", {"HcalBarrelCollection"} };
+  Gaudi::Property< std::vector<std::string> > name_CaloHit{ this, "CaloHitCollection", {"HCALBarrel"} };
+  Gaudi::Property< std::vector<std::string> > name_CaloAsso{ this, "CaloAssociationCollection", {"HCALBarrelAssoCol"} };
+  Gaudi::Property< std::vector<std::string> > name_CaloMCPAsso{ this, "CaloMCPAssociationCollection", {"HCALBarrelParticleAssoCol"} };
+
+  std::vector<SimCaloType*> _inputSimHitCollection;
+  std::vector<CaloType*> _outputHitCollection;
+  std::vector<CaloSimAssoType*> _outputCaloSimAssoCol;
+  std::vector<CaloParticleAssoType*> _outputCaloMCPAssoCol;
+
+  DataHandle<edm4hep::MCParticleCollection> m_MCParticleCol{"MCParticle", Gaudi::DataHandle::Reader, this};
+  mutable Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"};
+
+  //Input parameters
+  mutable Gaudi::Property<int>   _writeNtuple{this,  "WriteNtuple", 1, "Write ntuple"};
+  mutable Gaudi::Property<int>   _Nskip{this,  "SkipEvt", 0, "Skip event"};
+  mutable Gaudi::Property<float> _seed{this,   "Seed", 2131, "Random Seed"};
+  mutable Gaudi::Property<float> r_cali{this,  "CalibrHCAL", 1, "Global calibration coefficients for HCAL"};
+  mutable Gaudi::Property<int>   _UseRelDigi{this,   "UseRealisticDigi",  1, "If use the realistic model"};
+
+  //add digitization parameters from AHCAL prototype
+  //Scintillation and general
+  mutable Gaudi::Property<float> _MIPCali{this,   "MIPResponse",  0.007126, "MIP response (/GeV)"};
+  mutable Gaudi::Property<float> _MIPLY{this,   "MIPLY",  80, "Detected light yield (p.e./MIP)"};
+  mutable Gaudi::Property<float> _Eth_Mip{this,   "MIPThreshold", 0.1, "Energy Threshold (/MIP)"};
+  mutable Gaudi::Property<float> _TileRes{this,   "TileNonUniformity",  0., "Non-uniformity level of one tile response"};
+
+  mutable Gaudi::Property<int> _UseTileLYMap{this,   "UseTileLYMap",  1., "Use the tile light yield map"};
+  mutable Gaudi::Property<std::string> _TileLYMapFile{this,   "TileLYMapFile",  "", "Use the tile light yield map"};
+  mutable Gaudi::Property<std::string> _EffAttenLength{this,   "EffAttenLength",  "23mm", "Effictive attenuation length (mm)"};
+  mutable Gaudi::Property<float> _TempCoef{this,   "LYTempCoef",  0, "Temperature dependence of light yield (%/K)"};
+
+  //Temperature
+  mutable Gaudi::Property<float> _TempVariation{this,   "TemperatureVariation",  1., "Temperature control variation (K)"};
+
+  //SiPM
+  mutable Gaudi::Property<int>   _Pixel{this,   "SiPMPixel",  57600, "number of SiPM pixels"};
+  mutable Gaudi::Property<float> _SiPMDCR{this,   "SiPMDCR", 1600, "SiPM Dark Count Rate (Hz)"};
+  mutable Gaudi::Property<float> _SiPMXTalk{this,  "SiPMCT", 0.12, "SiPM crosstalk Probability"};
+  mutable Gaudi::Property<float> _TimeInterval{this,  "TimeInterval", 0.000002, "Time interval for one readout (s)"};
+  mutable Gaudi::Property<float> _SiPMGainTempCoef{this,   "SiPMGainTempCoef", -0.03, "Temperature dependence of SiPM gain (%/K)"}; // doi:10.1016/j.nima.2016.09.053
+  mutable Gaudi::Property<float> _SiPMDCRTempCoef{this,   "SiPMDCRTempCoef", 3.34/80, "Temperature dependence of SiPM DCR (10^{k*deltaT})"}; // doi:10.1016/j.nima.2016.09.053
+
+  //ADC
+  mutable Gaudi::Property<int> _ADC{this,   "ADC", 8192, "Total ADC conuts"};
+  mutable Gaudi::Property<int> _ADCSwitch{this,   "ADCSwitch", 8000, "switching point of different gain mode"};
+  mutable Gaudi::Property<float> _GainRatio_12{this,  "GainRatio_12", 50, "Gain-1 over Gain-2"};
+  mutable Gaudi::Property<float> _GainRatio_23{this,  "GainRatio_23", 60, "Gain-2 over Gain-3"};
+  mutable Gaudi::Property<float> _SiPMGainMean{this,    "SiPMGainMean", 2, "SiPM gain: ADC per p.e. for HG (ADC)"};
+  mutable Gaudi::Property<float> _SiPMGainSigma{this,   "SiPMGainSigma", 0.08, "Fluctuation of single photoelctron ADC around the mean value for single device (%)"};
+  mutable Gaudi::Property<float> _SiPMNoiseSigma{this,  "SiPMNoiseSigma", 0, "Sigma of SiPM noise (ADC)"};
+  mutable Gaudi::Property<float> _Pedestal{this,  "Pedestal", 50, "ADC value of pedestal"};
+  mutable Gaudi::Property<float> _PedestalNoiseSigma{this,  "PedestalSigma", 4, "Sigma of electronic noise (ADC)"};
+
+  // Output collections
+  //DataHandle<edm4hep::CalorimeterHitCollection>    w_DigiCaloCol{"DigiCaloCol", Gaudi::DataHandle::Writer, this};
+  //DataHandle<edm4hep::MCRecoCaloAssociationCollection>    w_CaloAssociationCol{"HCALBarrelAssoCol", Gaudi::DataHandle::Writer, this};
+  //DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>    w_MCPCaloAssociationCol{"HCALBarrelParticleAssoCol", Gaudi::DataHandle::Writer, this};
+};
+
+#endif

Digitization/DigiCalo/src/HitStep.h

+#ifndef HIT_STEP_H
+#define HIT_STEP_H
+
+class HitStep{
+
+public: 
+  HitStep (double _Q, double _T): Q(_Q), T(_T) {};
+  //HitStep (double _Q, double _T) { Q=_Q; T=_T; };
+  HitStep() {};
+
+  void setQ(double _Q) { Q =_Q; }
+  void setT(double _T) { T =_T; }
+
+  double getQ() const { return Q; }
+  double getT() const { return T; }
+  inline bool operator < (const HitStep &x) const { 
+    return T <x.T ;
+  }
+
+private: 
+  double Q;
+  double T;
+
+};
+
+#endif

Digitisers/DCHDigi/CMakeLists.txt → Digitization/DigiDCH/CMakeLists.txt

@@ -3,13 +3,13 @@ gaudi_add_module(DCHDigi
                  SOURCES src/DCHDigiAlg.cpp
                  LINK DetInterface
                       DetSegmentation
+                      DataHelperLib
                       k4FWCore::k4FWCore 
                       Gaudi::GaudiKernel
                       Gaudi::GaudiAlgLib
                       ${CLHEP_LIBRARIES}
                       ${DD4hep_COMPONENT_LIBRARIES}
                       ${ROOT_LIBRARIES}
-                      DetSegmentation
                       EDM4HEP::edm4hep EDM4HEP::edm4hepDict
 )
 install(TARGETS DCHDigi

Digitisers/DCHDigi/src/DCHDigiAlg.cpp → Digitization/DigiDCH/src/DCHDigiAlg.cpp

@@ -21,6 +21,9 @@
 #include <algorithm>
 
 #include "time.h"
+#include <unordered_map>
+
+#include "DataHelper/HelixClass.h"
 
 DECLARE_COMPONENT( DCHDigiAlg )
 
@@ -28,15 +31,16 @@ DCHDigiAlg::DCHDigiAlg(const std::string& name, ISvcLocator* svcLoc)
   : GaudiAlgorithm(name, svcLoc),
     _nEvt(0)
 {
-
+  
   // Input collections
   declareProperty("SimDCHitCollection", r_SimDCHCol, "Handle of the Input SimHit collection");
-
+  
   // Output collections
   declareProperty("DigiDCHitCollection", w_DigiDCHCol, "Handle of Digi DCHit collection");
-
+  declareProperty("SignalDigiDCHitCollection", w_SignalDigiDCHCol, "Handle of Signal Digi DCHit collection");
+  
   declareProperty("AssociationCollection", w_AssociationCol, "Handle of Association collection");
-
+   
 }
 
 StatusCode DCHDigiAlg::initialize()
@@ -53,7 +57,10 @@ StatusCode DCHDigiAlg::initialize()
       return StatusCode::FAILURE;
   }
 
+  fRandom.SetSeed(105105);//FIXME: set by users
+
   if(m_WriteAna){
+
       NTuplePtr nt( ntupleSvc(), "MyTuples/DCH_digi_evt" );
       if ( nt ) m_tuple = nt;
       else {
@@ -66,17 +73,34 @@ StatusCode DCHDigiAlg::initialize()
             m_tuple->addItem( "simhit_x", m_n_sim, m_simhit_x).ignore();
             m_tuple->addItem( "simhit_y", m_n_sim, m_simhit_y).ignore();
             m_tuple->addItem( "simhit_z", m_n_sim, m_simhit_z).ignore();
+            m_tuple->addItem( "Simdca"      , m_n_sim, m_Simdca       ).ignore();
+            m_tuple->addItem( "simhitT", m_n_sim, m_simhitT    ).ignore();
+            m_tuple->addItem( "simhitmom", m_n_sim, m_simhitmom    ).ignore();
+            m_tuple->addItem( "simPDG", m_n_sim, m_simPDG    ).ignore();
             m_tuple->addItem( "chamber" , m_n_digi, m_chamber  ).ignore();
             m_tuple->addItem( "layer"   , m_n_digi, m_layer    ).ignore();
             m_tuple->addItem( "cell"    , m_n_digi, m_cell     ).ignore();
             m_tuple->addItem( "cell_x"  , m_n_digi, m_cell_x   ).ignore();
             m_tuple->addItem( "cell_y"  , m_n_digi, m_cell_y   ).ignore();
+            m_tuple->addItem( "cell1_x"  , m_n_digi, m_cell1_x   ).ignore();
+            m_tuple->addItem( "cell1_y"  , m_n_digi, m_cell1_y   ).ignore();
             m_tuple->addItem( "hit_x"    , m_n_digi,m_hit_x     ).ignore();
             m_tuple->addItem( "hit_y"    , m_n_digi,m_hit_y     ).ignore();
             m_tuple->addItem( "hit_z"    , m_n_digi,m_hit_z     ).ignore();
-            m_tuple->addItem( "dca"      , m_n_digi,m_dca       ).ignore();
+            m_tuple->addItem( "mom_x"    , m_n_digi,m_mom_x     ).ignore();
+            m_tuple->addItem( "mom_y"    , m_n_digi,m_mom_y     ).ignore();
+            m_tuple->addItem( "dca"      , m_n_digi, m_dca       ).ignore();
+            m_tuple->addItem( "poca_x"   , m_n_digi, m_poca_x    ).ignore();
+            m_tuple->addItem( "poca_y"   , m_n_digi, m_poca_y    ).ignore();
             m_tuple->addItem( "hit_dE"   , m_n_digi,m_hit_dE    ).ignore();
-            m_tuple->addItem( "hit_dE_dx", m_n_digi,m_hit_dE_dx ).ignore();
+            m_tuple->addItem( "truthlength", m_n_digi,m_truthlength).ignore();
+            m_tuple->addItem( "N_noise", m_n_noise, 0, 1000000 ).ignore();
+            m_tuple->addItem( "layer_noise", m_n_noise,m_noise_layer).ignore();
+            m_tuple->addItem( "cell_noise", m_n_noise,m_noise_cell).ignore();
+            m_tuple->addItem( "NoiseT", m_n_noise,m_noise_time).ignore();
+            m_tuple->addItem( "N_noiseCover", m_n_noiseCover, 0, 1000000 ).ignore();
+            m_tuple->addItem( "layer_noiseCover", m_n_noiseCover,m_noiseCover_layer).ignore();
+            m_tuple->addItem( "cell_noiseCover", m_n_noiseCover,m_noiseCover_cell).ignore();
           } else { // did not manage to book the N tuple....
             info() << "    Cannot book N-tuple:" << long( m_tuple ) << endmsg;
           }
@@ -92,8 +116,9 @@ StatusCode DCHDigiAlg::execute()
   m_start = clock();
 
   info() << "Processing " << _nEvt << " events " << endmsg;
-  m_evt = _nEvt;
+  if(m_WriteAna) m_evt = _nEvt;
   edm4hep::TrackerHitCollection* Vec   = w_DigiDCHCol.createAndPut();
+  edm4hep::TrackerHitCollection* SignalVec   = w_SignalDigiDCHCol.createAndPut();
   edm4hep::MCRecoTrackerAssociationCollection* AssoVec   = w_AssociationCol.createAndPut();
   const edm4hep::SimTrackerHitCollection* SimHitCol =  r_SimDCHCol.get();
   if (SimHitCol->size() == 0) {
@@ -102,8 +127,8 @@ StatusCode DCHDigiAlg::execute()
   debug()<<"input sim hit size="<< SimHitCol->size() <<endmsg;
 
   auto SimHit0 = SimHitCol->at(0);
-  std::map<unsigned long long, std::vector<decltype(SimHit0)> > id_hits_map;
-  //std::map<unsigned long long, std::vector<edm4hep::SimTrackerhit> > id_hits_map;
+  std::map<unsigned long long, std::vector<decltype(SimHit0)>> id_hits_map;
+
 
   for( int i = 0; i < SimHitCol->size(); i++ ) 
   {
@@ -111,7 +136,12 @@ StatusCode DCHDigiAlg::execute()
       unsigned long long id = SimHit.getCellID();
       float sim_hit_mom = sqrt( SimHit.getMomentum()[0]*SimHit.getMomentum()[0] + SimHit.getMomentum()[1]*SimHit.getMomentum()[1] + SimHit.getMomentum()[2]*SimHit.getMomentum()[2] );//GeV
       if(sim_hit_mom < m_mom_threshold) continue; 
-      if(SimHit.getEDep() <= 0) continue;
+      if(sim_hit_mom > m_mom_threshold_high) continue; 
+      if(SimHit.getEDep() <= m_edep_threshold) continue;
+
+      //Wire efficiency
+      float hitProb = fRandom.Uniform(1.);
+      if(hitProb > m_wireEff) continue;
 
       if ( id_hits_map.find(id) != id_hits_map.end()) id_hits_map[id].push_back(SimHit);
       else 
@@ -125,123 +155,219 @@ StatusCode DCHDigiAlg::execute()
       m_n_sim = 0;
       m_n_digi = 0 ;
   }
+
+  std::vector<TVector2> HitPosition2D;
+  bool ismixNoise[55] = {0};
+  std::cout << " id_hits_map  size = " << id_hits_map.size() << std::endl;
+
   for(auto iter = id_hits_map.begin(); iter != id_hits_map.end(); iter++)
   {
-    unsigned long long wcellid = iter->first;
-    auto trkHit = Vec->create();
-    trkHit.setCellID(wcellid);
-    double tot_edep   = 0 ;
-    double tot_length = 0 ;
-    double pos_x = 0 ;
-    double pos_y = 0 ;
-    double pos_z = 0 ;
-    int simhit_size = iter->second.size();
-    for(unsigned int i=0; i< simhit_size; i++)
-    {
-        tot_edep += iter->second.at(i).getEDep();//GeV
-    }
-    int chamber = m_decoder->get(wcellid, "chamber");
-    int layer   = m_decoder->get(wcellid, "layer"  );
-    int cellID  = m_decoder->get(wcellid, "cellID" );
-    TVector3 Wstart(0,0,0);
-    TVector3 Wend  (0,0,0);
-    m_segmentation->cellposition(wcellid, Wstart, Wend);
-    float dd4hep_mm = dd4hep::mm;
-    //std::cout<<"dd4hep_mm="<<dd4hep_mm<<std::endl;
-//    Wstart =(1/dd4hep_mm)* Wstart;// from DD4HEP cm to mm
-//    Wend   =(1/dd4hep_mm)* Wend  ;
-    //std::cout<<"wcellid="<<wcellid<<",chamber="<<chamber<<",layer="<<layer<<",cellID="<<cellID<<",s_x="<<Wstart.x()<<",s_y="<<Wstart.y()<<",s_z="<<Wstart.z()<<",E_x="<<Wend.x()<<",E_y="<<Wend.y()<<",E_z="<<Wend.z()<<std::endl;
-
-    TVector3  denominator = (Wend-Wstart) ;
-    float min_distance = 999 ;
-    float min_line_distance = 999 ;
-    float tmp_distance =0;
-    for(unsigned int i=0; i< simhit_size; i++)
-    {
-        float sim_hit_mom = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] + iter->second.at(i).getMomentum()[2]*iter->second.at(i).getMomentum()[2] );//GeV
-        float sim_hit_pt = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] );//GeV
-        TVector3  pos(iter->second.at(i).getPosition()[0]*dd4hep_mm, iter->second.at(i).getPosition()[1]*dd4hep_mm, iter->second.at(i).getPosition()[2]*dd4hep_mm);
-
-//        TVector3  numerator = denominator.Cross(Wstart-pos) ;
-//        float tmp_distance = numerator.Mag()/denominator.Mag() ;
-
-        TVector3 sim_mon(iter->second.at(i).getMomentum()[0],iter->second.at(i).getMomentum()[1],iter->second.at(i).getMomentum()[2]);
-        float Steplength = iter->second.at(i).getPathLength();
-        TVector3  pos_start = pos - 0.5 * Steplength * sim_mon.Unit();
-        TVector3  pos_end = pos + 0.5 * Steplength * sim_mon.Unit();
-        if(m_Doca) {
-            tmp_distance = m_segmentation->distanceTrackWire(wcellid,pos_start,pos_end);
-            tmp_distance = tmp_distance/dd4hep_mm; //mm
-        } else {
-            tmp_distance = (m_segmentation->distanceClosestApproach(wcellid,pos)).Mag();
-            tmp_distance = tmp_distance/dd4hep_mm; //mm
-        }
-
-
-       // std::cout << " Steplength= " << Steplength << std::endl;
-       // std::cout<<"tmp_distance="<<tmp_distance<<",x="<<pos.x()<<",y="<<pos.y()<<",z="<<pos.z()<<",mom="<<sim_hit_mom<<",pt="<<sim_hit_pt<<std::endl;
-
-        if(tmp_distance < min_distance){
-            min_distance = tmp_distance;
-            pos_x = pos.x();
-            pos_y = pos.y();
-            pos_z = pos.z();
-        }
-        tot_length += iter->second.at(i).getPathLength();//mm
-        auto asso = AssoVec->create();
-        asso.setRec(trkHit);
-        asso.setSim(iter->second.at(i));
-        asso.setWeight(iter->second.at(i).getEDep()/tot_edep);
-
-        if(m_WriteAna && (nullptr!=m_tuple)) { // && min_distance <0.3){
-            m_simhit_x[m_n_sim] = pos.x();
-            m_simhit_y[m_n_sim] = pos.y();
-            m_simhit_z[m_n_sim] = pos.z();
-            m_n_sim ++ ;
-        }
-    }
-
-    trkHit.setTime(min_distance*1e3/m_velocity);//m_velocity is um/ns, drift time in ns
-    trkHit.setEDep(tot_edep);// GeV
-    trkHit.setEdx (tot_edep/tot_length); // GeV/mm
-    trkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit
-    trkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm
-
-    if(m_WriteAna && (nullptr!=m_tuple)) { // && min_distance <0.3){
-        m_chamber  [m_n_digi] = chamber;
-        m_layer    [m_n_digi] = layer  ;
-        m_cell     [m_n_digi] = cellID;
-        m_cell_x   [m_n_digi] = Wstart.x();
-        m_cell_y   [m_n_digi] = Wstart.y();
-        m_hit_x    [m_n_digi] = pos_x;
-        m_hit_y    [m_n_digi] = pos_y;
-        m_hit_z    [m_n_digi] = pos_z;
-        m_dca      [m_n_digi] = min_distance;
-        m_hit_dE   [m_n_digi] = trkHit.getEDep();
-        m_hit_dE_dx[m_n_digi] = trkHit.getEdx() ;
-        m_n_digi ++ ;
-    }
-  }
+      unsigned long long wcellid = iter->first;
+      auto trkHit = Vec->create();
+      auto SignaltrkHit = SignalVec->create();
+      trkHit.setCellID(wcellid);
+      SignaltrkHit.setCellID(wcellid);
+      float tot_edep   = 0 ;
+      float tot_length = 0 ;
+      float pos_x = 0 ;
+      float pos_y = 0 ;
+      float pos_z = 0 ;
+      float momx,momy,momz = 0;
+      const int simhit_size = iter->second.size();
+      for(unsigned int i=0; i< simhit_size; i++)
+      {
+          tot_edep += iter->second.at(i).getEDep();//GeV
+      }
+      int chamber = m_decoder->get(wcellid, "chamber");
+      int layer   = m_decoder->get(wcellid, "layer"  );
+      int cellID  = m_decoder->get(wcellid, "cellID" );
+      TVector3 Wstart(0,0,0);
+      TVector3 Wend  (0,0,0);
+      m_segmentation->cellposition(wcellid, Wstart, Wend);
+      float dd4hep_mm = dd4hep::mm;
+      if(m_debug) std::cout<<"DCHDigi wcellid ="<<wcellid<< ",chamber="<<chamber<<",layer="<<layer<<",cellID="<<cellID<<",s_x="<<Wstart.x()<<",s_y="<<Wstart.y()<<",s_z="<<Wstart.z()<<",E_x="<<Wend.x()<<",E_y="<<Wend.y()<<",E_z="<<Wend.z()<<std::endl;
+
+      TVector3  denominator = (Wend-Wstart) ;
+      float min_distance = 999 ;
+      float sim_distance = 999 ;
+      float tmp_distance =0;
+      TVector3 hitPosition;
+      TVector3 PCA;
+      for(unsigned int i=0; i< simhit_size; i++)
+      {
+          float sim_hit_mom = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] + iter->second.at(i).getMomentum()[2]*iter->second.at(i).getMomentum()[2] );//GeV
+          if(sim_hit_mom<1e-4) continue;
+
+          TVector3  pos(iter->second.at(i).getPosition()[0]*dd4hep_mm, iter->second.at(i).getPosition()[1]*dd4hep_mm, iter->second.at(i).getPosition()[2]*dd4hep_mm);
+
+          float simpos[3] = {iter->second.at(i).getPosition()[0],iter->second.at(i).getPosition()[1],iter->second.at(i).getPosition()[2]};
+          float simmom[3] = {iter->second.at(i).getMomentum()[0],iter->second.at(i).getMomentum()[1],iter->second.at(i).getMomentum()[2]};
+
+          TVector3 sim_mon(iter->second.at(i).getMomentum()[0],iter->second.at(i).getMomentum()[1],iter->second.at(i).getMomentum()[2]);
+          float Steplength = iter->second.at(i).getPathLength();
+          TVector3  pos_start = pos - 0.5 * Steplength * sim_mon.Unit();
+          TVector3  pos_end = pos + 0.5 * Steplength * sim_mon.Unit();
+          tmp_distance = m_segmentation->Distance(wcellid,pos_start,pos_end,hitPosition,PCA);
+          tmp_distance = tmp_distance/dd4hep_mm; //mm
+          sim_distance = tmp_distance;
+
+          if(tmp_distance < min_distance){
+              min_distance = tmp_distance;
+              pos_x = hitPosition.x();     //pos.x();
+              pos_y = hitPosition.y();     //pos.y();
+              pos_z = pos.z();
+              momx = iter->second.at(i).getMomentum()[0];
+              momy = iter->second.at(i).getMomentum()[1];
+          }
+          tot_length += iter->second.at(i).getPathLength();//mm
+          auto asso = AssoVec->create();
+          asso.setRec(trkHit);
+          asso.setSim(iter->second.at(i));
+          asso.setWeight(iter->second.at(i).getEDep()/tot_edep);
+
+          if(m_WriteAna && (nullptr!=m_tuple)) {
+              m_simhit_x[m_n_sim] = pos.x();
+              m_simhit_y[m_n_sim] = pos.y();
+              m_simhit_z[m_n_sim] = pos.z();
+              m_Simdca[m_n_sim] = sim_distance;
+              m_simhitT[m_n_sim] = iter->second.at(i).getTime();
+              m_simhitmom[m_n_sim] = sim_hit_mom;
+              m_simPDG[m_n_sim] = iter->second.at(i).getMCParticle().getPDG();
+              m_n_sim ++ ;
+          }
+      }
+
+      //Smear drift distance
+      if(m_smear){
+          float dd = min_distance + gRandom->Gaus(0,fabs(m_res_x));
+          if(dd>=0) min_distance = dd;
+          else min_distance = -dd;
+      }
 
+      trkHit.setTime(min_distance*1e3/m_velocity);//m_velocity is um/ns, drift time in ns
+      trkHit.setQuality(1);//1:singal   0:noise of overlap singal  -1:noise
+      trkHit.setEDep(tot_edep);// GeV
+      trkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit
+      trkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm
+
+      //add noise every layer
+      if(m_mixNoise && !ismixNoise[layer]) mixNoise(layer,cellID,Vec,&trkHit,ismixNoise);
+
+      SignaltrkHit.setEDep(tot_edep);// GeV
+      SignaltrkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit
+      SignaltrkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm
+
+      if(m_WriteAna && (nullptr!=m_tuple)) {
+          m_chamber  [m_n_digi] = chamber;
+          m_layer    [m_n_digi] = layer  ;
+          m_cell     [m_n_digi] = cellID;
+          m_cell_x   [m_n_digi] = Wstart.x();
+          m_cell_y   [m_n_digi] = Wstart.y();
+          m_cell1_x   [m_n_digi] = Wend.x();
+          m_cell1_y   [m_n_digi] = Wend.y();
+          m_hit_x    [m_n_digi] = pos_x;
+          m_hit_y    [m_n_digi] = pos_y;
+          m_hit_z    [m_n_digi] = pos_z;
+          m_mom_x    [m_n_digi] = momx ;
+          m_mom_y    [m_n_digi] = momy ;
+          m_dca      [m_n_digi] = min_distance;
+          m_poca_x   [m_n_digi] = PCA.x();
+          m_poca_y   [m_n_digi] = PCA.y();
+          m_hit_dE   [m_n_digi] = trkHit.getEDep();
+          m_truthlength[m_n_digi] = tot_length ;
+          m_n_digi ++ ;
+      }
 
+  }
   debug()<<"output digi DCHhit size="<< Vec->size() <<endmsg;
   _nEvt ++ ;
 
   if(m_WriteAna && (nullptr!=m_tuple)){
       StatusCode status = m_tuple->write();
       if ( status.isFailure() ) {
-        error() << "    Cannot fill N-tuple:" << long( m_tuple ) << endmsg;
-        return StatusCode::FAILURE;
+          error() << "    Cannot fill N-tuple:" << long( m_tuple ) << endmsg;
+          return StatusCode::FAILURE;
       }
   }
   m_end = clock();
-  m_time = (m_end - m_start);
+  if(m_WriteAna){
+      m_time = (m_end - m_start);
+  }
 
   return StatusCode::SUCCESS;
 }
 
 StatusCode DCHDigiAlg::finalize()
 {
-  info() << "Processed " << _nEvt << " events " << endmsg;
-  return GaudiAlgorithm::finalize();
+    info() << "Processed " << _nEvt << " events " << endmsg;
+    return GaudiAlgorithm::finalize();
 }
+
+void DCHDigiAlg::mixNoise(int layerID ,int wireID,
+        edm4hep::TrackerHitCollection* Vec,
+        edm4hep::MutableTrackerHit* trackerHitLayer,
+        bool ismixNoise[55]){
+
+    int maxCellID = m_segmentation->maxWireID(0,layerID);
+    int nNoiseHits = (int) maxCellID*m_noiseEff;
+
+    //loop over noise hits in this layer
+    for(int ic=0;ic<nNoiseHits;ic++){
+        int cellid = (int) maxCellID*fRandom.Uniform(1.);
+        float time = m_timeWindow*fRandom.Uniform(1.);
+        TVector3 Wstart,Wend;
+        m_segmentation->cellposition2(0,layerID,cellid,Wstart,Wend);
+        //found noise overlap with signal
+        if(!(cellid != wireID))
+        {
+            if(time<trackerHitLayer->getTime()){
+                //1. overlap with signal and time > signal time
+                if(!m_mixNoiseAndSkipOverlap){
+                    //set time
+                    trackerHitLayer->setTime(time);
+                    trackerHitLayer->setQuality(0);
+
+                    m_noiseCover_layer[m_n_noiseCover] = layerID;
+                    m_noiseCover_cell[m_n_noiseCover] = cellid;
+
+                    m_n_noiseCover++;
+
+                    m_noise_layer[m_n_noise] = layerID;
+                    m_noise_cell[m_n_noise] = cellid;
+                    m_noise_time[m_n_noise] = time;
+
+                    m_n_noise++;
+                } // end fo m_mixNoiseAndSkipOverlap
+            }else{
+                //2. overlap with signal and time > signal time
+                // TODO adc sum with signal
+                continue;
+            }
+        }else{ // not found noise overlap with signal
+            //3. create a noise hit
+            //create another trackerHit
+            if(!m_mixNoiseAndSkipNoOverlap){
+                auto trkNoiseHit = Vec->create();
+
+                dd4hep::rec::CellID ncellid;
+                m_decoder->set(ncellid,"chamber",0);
+                m_decoder->set(ncellid,"layer",layerID);
+                m_decoder->set(ncellid,"cellID",cellid);
+
+                trkNoiseHit.setCellID(ncellid);
+                trkNoiseHit.setTime(time);
+                trkNoiseHit.setQuality(-1);
+
+                m_noise_layer[m_n_noise] = layerID;
+                m_noise_cell[m_n_noise] = cellid;
+                m_noise_time[m_n_noise] = time;
+
+                m_n_noise++;
+
+            }// end of m_mixNoiseAndSkipNoOverlap
+        } // end of if
+
+
+    }// loop over noise hit
+    ismixNoise[layerID] = true;
+}//end of mixNoise

Digitisers/DCHDigi/src/DCHDigiAlg.h → Digitization/DigiDCH/src/DCHDigiAlg.h

@@ -7,6 +7,7 @@
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/TrackerHitCollection.h"
 #include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCParticle.h"
 
 #include <DDRec/DetectorData.h>
 #include <DDRec/CellIDPositionConverter.h>
@@ -15,9 +16,13 @@
 #include "DetSegmentation/GridDriftChamber.h"
 
 #include "TVector3.h"
+#include "TRandom3.h"
 
-
-
+namespace dd4hep {
+    namespace rec{
+        class CellIDPositionConverter;
+    }
+}
 class DCHDigiAlg : public GaudiAlgorithm
 {
  
@@ -44,6 +49,15 @@ protected:
   typedef std::vector<float> FloatVec;
   int _nEvt ;
 
+  void addNoiseHit(edm4hep::TrackerHitCollection* Vec,double time, int chamber,
+          int layer, int cellID,unsigned long long wcellid);
+  void mixNoise(int layerID,int wireID,
+          edm4hep::TrackerHitCollection* Vec,
+          edm4hep::MutableTrackerHit* trackerHitLayer,
+          bool ismixNoise[55]);
+
+  TRandom3 fRandom;
+
   NTuple::Tuple* m_tuple = nullptr ;
   NTuple::Item<int> m_evt;
   NTuple::Item<long>   m_n_sim;
@@ -52,17 +66,40 @@ protected:
   NTuple::Array<int  > m_chamber   ;
   NTuple::Array<int  > m_layer     ;
   NTuple::Array<int  > m_cell      ;
+  //The position of wire at -z
   NTuple::Array<float> m_cell_x    ;
   NTuple::Array<float> m_cell_y    ;
+  //The position of wire at +z
+  NTuple::Array<float> m_cell1_x   ;
+  NTuple::Array<float> m_cell1_y   ;
+
   NTuple::Array<float> m_simhit_x  ;
   NTuple::Array<float> m_simhit_y  ;
   NTuple::Array<float> m_simhit_z  ;
+  NTuple::Array<float> m_simhitT ;
+  NTuple::Array<float> m_simhitmom  ;
+  NTuple::Array<int> m_simPDG  ;
   NTuple::Array<float> m_hit_x     ;
   NTuple::Array<float> m_hit_y     ;
   NTuple::Array<float> m_hit_z     ;
+  NTuple::Array<float> m_mom_x     ;
+  NTuple::Array<float> m_mom_y     ;
   NTuple::Array<float> m_dca       ;
+  NTuple::Array<float> m_Simdca       ;
+  NTuple::Array<float> m_poca_x    ;
+  NTuple::Array<float> m_poca_y    ;
   NTuple::Array<float> m_hit_dE    ;
   NTuple::Array<float> m_hit_dE_dx ;
+  NTuple::Array<double> m_truthlength ;
+
+  //Noise hit in cellID
+  NTuple::Item<long>   m_n_noise;
+  NTuple::Array<int>  m_noise_layer;
+  NTuple::Array<int>  m_noise_cell;
+  NTuple::Array<float>  m_noise_time;
+  NTuple::Item<long>   m_n_noiseCover;
+  NTuple::Array<int>  m_noiseCover_layer;
+  NTuple::Array<int>  m_noiseCover_cell;
 
   clock_t m_start,m_end;
 
@@ -77,15 +114,30 @@ protected:
   Gaudi::Property<float> m_res_z     { this, "res_z", 1   };//mm
   Gaudi::Property<float> m_velocity  { this, "drift_velocity", 40};// um/ns
   Gaudi::Property<float> m_mom_threshold { this, "mom_threshold", 0};// GeV
+  Gaudi::Property<float> m_mom_threshold_high { this, "mom_threshold_high", 1e9};// GeV
+  Gaudi::Property<float> m_edep_threshold{ this, "edep_threshold", 0};// GeV
   Gaudi::Property<bool>  m_WriteAna { this, "WriteAna", false};
-  Gaudi::Property<bool>  m_Doca { this, "Doca", false};//1:line dca 0:point dca
+  Gaudi::Property<bool>  m_debug{ this, "debug", false};
+  Gaudi::Property<double>  m_wireEff{ this, "wireEff", 1.0};
+
+  //Smear drift distance
+  Gaudi::Property<bool>  m_smear{ this, "smear", false};
 
+  //Mixed Noise
+  Gaudi::Property<bool>  m_mixNoise{ this, "mixNoise", false};
+  Gaudi::Property<bool>  m_mixNoiseAndSkipOverlap{ this, "mixNoiseAndSkipOverlapthis", false};
+  Gaudi::Property<bool>  m_mixNoiseAndSkipNoOverlap{ this, "mixNoiseAndSkipNoOverlapthis", false};
+  Gaudi::Property<double>  m_noiseEff{ this, "noiseEff", 0.0};
+  Gaudi::Property<double>  m_timeWindow{ this, "timeWindow", 2000.0}; //ns
+  Gaudi::Property<int>  m_layerNum{ this, "layerNum", 55}; //ns
 
   // Input collections
   DataHandle<edm4hep::SimTrackerHitCollection> r_SimDCHCol{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};
   // Output collections
   DataHandle<edm4hep::TrackerHitCollection>    w_DigiDCHCol{"DigiDCHitCollection", Gaudi::DataHandle::Writer, this};
   DataHandle<edm4hep::MCRecoTrackerAssociationCollection>    w_AssociationCol{"DCHitAssociationCollection", Gaudi::DataHandle::Writer, this};
-
+  // Output signal digiHits
+  DataHandle<edm4hep::TrackerHitCollection>    w_SignalDigiDCHCol{"SignalDigiDCHitCollection", Gaudi::DataHandle::Writer, this};
 };
+
 #endif

Digitization/DigiMuon/CMakeLists.txt

+# Modules
+gaudi_add_module(MuonDigi
+                 SOURCES src/MuonDigiAlg.cpp
+                 LINK k4FWCore::k4FWCore
+                      GearSvc
+                      DetInterface
+                      Gaudi::GaudiKernel
+                      Gaudi::GaudiAlgLib 
+                      ${CLHEP_LIBRARIES}
+                      ${GEAR_LIBRARIES} 
+                      ${GSL_LIBRARIES} 
+                      ${LCIO_LIBRARIES}
+                      ${ROOT_LIBRARIES}
+                      EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+)
+install(TARGETS MuonDigi
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)
+

Digitization/DigiMuon/src/MuonDigiAlg.cpp

+#include "MuonDigiAlg.h"
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/Vector3D.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/Vector3f.h"
+#include "edm4hep/Cluster.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include <DDRec/CellIDPositionConverter.h>
+
+#include "TVector2.h"
+#include <cmath>
+#include <iostream>
+#include <algorithm>
+#include <map>
+#include <random>
+#include <vector>
+
+
+using namespace std ;
+
+DECLARE_COMPONENT( MuonDigiAlg )
+
+MuonDigiAlg::MuonDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc)
+{
+  // Input collections
+  declareProperty("MuonBarrelHitsCollection", m_inputMuonBarrel, "Handle of the Input SimTrackerHit collection");
+  declareProperty("MuonEndcapHitsCollection", m_inputMuonEndcap, "Handle of the Input SimTrackerHit collection");
+  // Output collections
+  declareProperty("MuonBarrelTrackerHits", m_outputMuonBarrel, "Handle of the output TrackerHit collection");
+  declareProperty("MuonEndcapTrackerHits", m_outputMuonEndcap, "Handle of the output TrackerHit collection");
+}
+
+
+StatusCode MuonDigiAlg::initialize()
+{
+  // set rand seed;
+  m_randSvc = service<IRndmGenSvc>("RndmGenSvc");
+
+  // 
+  if(_writeNtuple)
+  {
+    std::string s_outfile = _filename;
+    m_wfile = new TFile(s_outfile.c_str(), "recreate");
+    m_tree = new TTree("tree", "muon digi tree");
+    m_tree->Branch("n_hit", &n_hit, "n_hit/I");
+    m_tree->Branch("hit_cellid", hit_cellid, "hit_cellid[n_hit]/l");
+    m_tree->Branch("hit_posx", hit_posx, "hit_posx[n_hit]/D");
+    m_tree->Branch("hit_posy", hit_posy, "hit_posy[n_hit]/D");
+    m_tree->Branch("hit_posz", hit_posz, "hit_posz[n_hit]/D");
+    m_tree->Branch("hit_edep", hit_edep, "hit_edep[n_hit]/D");
+    m_tree->Branch("hit_layer", hit_layer, "hit_layer[n_hit]/I");
+    m_tree->Branch("hit_slayer", hit_slayer, "hit_slayer[n_hit]/I");
+    m_tree->Branch("hit_strip", hit_strip, "hit_strip[n_hit]/I");
+    m_tree->Branch("hit_fe", hit_fe, "hit_fe[n_hit]/I");
+    m_tree->Branch("hit_env", hit_env, "hit_env[n_hit]/I");
+    m_tree->Branch("n_cell", &n_cell, "n_cell/I");
+    m_tree->Branch("cell_cellid", cell_cellid, "cell_cellid[n_cell]/l");
+    m_tree->Branch("cell_edep", cell_edep, "cell_edep[n_cell]/D");
+    m_tree->Branch("cell_adc", cell_adc, "cell_adc[n_cell]/D");
+    m_tree->Branch("cell_posx", cell_posx, "cell_posx[n_cell]/D");
+    m_tree->Branch("cell_posy", cell_posy, "cell_posy[n_cell]/D");
+    m_tree->Branch("cell_posz", cell_posz, "cell_posz[n_cell]/D");
+    m_tree->Branch("cell_layer", cell_layer, "cell_layer[n_cell]/I");
+    m_tree->Branch("cell_slayer", cell_slayer, "cell_slayer[n_cell]/I");
+    m_tree->Branch("cell_strip", cell_strip, "cell_strip[n_cell]/I");
+    m_tree->Branch("cell_fe", cell_fe, "cell_fe[n_cell]/I");
+    m_tree->Branch("cell_env", cell_env, "cell_env[n_cell]/I");
+    m_tree->Branch("cell_superlayernumber", &cell_superlayernumber, "cell_superlayernumber/I");
+    m_tree->Branch("cell_pt", &cell_pt, "cell_pt/I");
+  }
+
+  m_geosvc = service<IGeomSvc>("GeomSvc");
+  if(!m_geosvc)
+  {
+    error() << "Failed to get the GeomSvc" << endmsg;
+    return StatusCode::FAILURE;
+  }
+  auto m_dd4hep = m_geosvc->lcdd();
+  if(!m_dd4hep)
+  {
+    error() << "Failed to get the Detector from GeomSvc" << endmsg;
+    return StatusCode::FAILURE;
+  }
+  m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep);
+  if(!m_cellIDConverter)
+  {
+    error() << "Failed to get the m_cellIDConverter." << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  std::string readout_name_barrel = m_inputMuonBarrel.objKey(); 
+  debug() << "Readout name: " << readout_name_barrel << endmsg;
+  m_decoder_barrel = m_geosvc->getDecoder(readout_name_barrel);
+  if(!m_decoder_barrel)
+  {
+    error() << "Failed to get the decoder. " << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  std::string readout_name_endcap = m_inputMuonEndcap.objKey(); 
+  debug() << "Readout name: " << readout_name_endcap << endmsg;
+  m_decoder_endcap = m_geosvc->getDecoder(readout_name_endcap);
+  if(!m_decoder_endcap)
+  {
+    error() << "Failed to get the decoder. " << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  debug() << "m_hitEff: " << m_hitEff << endmsg;
+ 
+  info() << "MuonDigiAlg::initialized" << endmsg;
+  return GaudiAlgorithm::initialize();
+}
+
+
+StatusCode MuonDigiAlg::execute()
+{
+  const edm4hep::SimTrackerHitCollection* STHCol;
+  Clear();
+  trkhitVec = m_outputMuonBarrel.createAndPut();
+  STHCol = nullptr;
+  try 
+  {
+    STHCol = m_inputMuonBarrel.get();
+  }
+  catch(GaudiException &e)
+  {
+    debug() << "Collection " << m_inputMuonBarrel.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+  }
+
+  int hit_data[6];
+  double hit_Edep, hit_xydist, hit_sipm_length, hit_mcpz;
+  double hit_time;
+  edm4hep::Vector3d hit_pos;
+  dd4hep::Position hit_ddpos;
+  unsigned long long cellid;
+  std::array<unsigned long long, 2> cell_key;
+
+  if(STHCol->size()>0)
+  {    
+    for(auto simhit : *STHCol)
+    {
+      GetSimHit(simhit, m_decoder_barrel, 2, hit_data, hit_Edep, cellid, hit_xydist, hit_pos, hit_ddpos, cell_key, hit_mcpz, hit_time);
+      debug()<<"Barrel Simulation Hit::::::: "<<hit_data[5]<<" "<<hit_data[2]<<" "<<hit_data[1]<<" "<<hit_data[0]<<" "<<hit_Edep<<endmsg;
+      SaveData_mapcell(cell_key, hit_Edep, hit_data, hit_pos, hit_time);
+    }
+    for (const auto& item : map_cell_edep) 
+    {
+      cell_key = item.first;
+      hit_data[0] = map_cell_pdgid[cell_key];
+      hit_data[1] = map_cell_layer[cell_key];
+      hit_data[2] = map_cell_slayer[cell_key];
+      hit_data[3] = map_cell_strip[cell_key];
+      hit_data[4] = map_cell_env[cell_key];
+      hit_data[5] = map_cell_fe[cell_key];
+      hit_sipm_length = Gethit_sipm_length_Barrel(m_cellIDConverter->position(cell_key[0]), map_cell_pos[cell_key], hit_data);
+      map_cell_adc[cell_key] = EdeptoADC(item.second, hit_sipm_length);
+      debug() <<"ADC::: "<<map_cell_adc[cell_key]<<endmsg;
+    }
+    // 2 for barrel
+    MuonHandle(2);
+  }
+
+  // Endcaps:
+  Clear();
+  trkhitVec = m_outputMuonEndcap.createAndPut();
+  STHCol = nullptr;
+  try 
+  {
+    STHCol = m_inputMuonEndcap.get();
+  }
+  catch(GaudiException &e)
+  {
+    debug() << "Collection " << m_inputMuonEndcap.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+  }
+  if(STHCol->size()==0) return StatusCode::SUCCESS;
+  for(auto simhit : *STHCol)
+  {
+    GetSimHit(simhit, m_decoder_endcap, 0,hit_data, hit_Edep, cellid, hit_xydist, hit_pos, hit_ddpos, cell_key, hit_mcpz, hit_time);
+    debug() <<"Endcap Simulation Hit::::::: "<<hit_data[4]<<" "<<hit_data[2]<<" "<<hit_data[1]<<" "<<hit_data[0]<<endmsg;
+    SaveData_mapcell(cell_key, hit_Edep, hit_data, hit_pos, hit_time);
+  }
+  for (const auto& item : map_cell_edep) 
+  {
+    cell_key = item.first;
+    hit_data[0] = map_cell_pdgid[cell_key];
+    hit_data[1] = map_cell_layer[cell_key];
+    hit_data[2] = map_cell_slayer[cell_key];
+    hit_data[3] = map_cell_strip[cell_key];
+    hit_data[4] = map_cell_env[cell_key];
+    hit_data[5] = map_cell_fe[cell_key];
+    hit_sipm_length = Gethit_sipm_length_Endcap(m_cellIDConverter->position(cell_key[0]), map_cell_pos[cell_key], hit_data);
+    map_cell_adc[cell_key] = EdeptoADC(item.second, hit_sipm_length);
+    debug()<<"ADC::: "<<map_cell_adc[cell_key]<<endmsg;
+  }
+  // 0 for endcap
+  MuonHandle(0);
+  Clear();
+  m_nEvt++;
+  return StatusCode::SUCCESS;
+}
+
+StatusCode MuonDigiAlg::finalize()
+{
+  if(_writeNtuple)
+  {
+    m_wfile->cd();
+    m_tree->Write();
+    m_wfile->Close();
+    delete m_tree;
+    delete m_wfile;
+  } 
+  info() << "Processed " << m_nEvt << " events " << endmsg;
+  return GaudiAlgorithm::finalize();
+}
+
+
+void MuonDigiAlg::MuonHandle(int barrel_or_endcap)
+{
+  //barrel_or_endcap: 2 is barrel, 0 is endcap
+
+  if(map_cell_edep.size()==0) return;
+  Cut3();
+  std::array<unsigned long long, 2> key1, key2;
+  dd4hep::Position ddpos1, ddpos2;
+  int anotherlayer_cell_num, true_pdgid;
+  for (const auto& item1 : map_cell_edep) 
+  {
+    key1 = item1.first;
+    if(map_cell_edep[key1] <= m_EdepMin) continue;
+    ddpos1 = m_cellIDConverter->position(key1[0]);
+    anotherlayer_cell_num = 0;   
+    for (const auto& item2 : map_cell_edep) 
+    {
+      key2 = item2.first;
+      if(map_cell_edep[key2] <= m_EdepMin) continue;
+      if(barrel_or_endcap == 2) // barrel
+      {
+        if(map_cell_slayer[key1]!=map_cell_slayer[key2] || 
+           map_cell_env[key1]!=map_cell_env[key2] || 
+           map_cell_fe[key1]!=map_cell_fe[key2] || 
+           map_cell_layer[key1]==map_cell_layer[key2] || 
+           map_cell_layer[key1]==2) continue;
+      }
+      if(barrel_or_endcap == 0) // endcap
+      {
+        if(map_cell_slayer[key1]!=map_cell_slayer[key2] || 
+           map_cell_fe[key1]!=map_cell_fe[key2] || 
+           map_cell_layer[key1]==map_cell_layer[key2] || 
+           map_cell_layer[key1]==2) continue;
+        if(std::abs(map_cell_env[key1] - map_cell_env[key2]) == 2) continue;        
+      }
+      ddpos2 = m_cellIDConverter->position(key2[0]);
+      double ddposi = (barrel_or_endcap == 2) ? ddpos2.z() : ddpos2.x();
+      Find_anotherlayer(barrel_or_endcap, key1, key2, ddposi, anotherlayer_cell_num);
+    }
+    if(anotherlayer_cell_num == 0) continue;
+    for(int i=1; i<=anotherlayer_cell_num; i++)
+    {
+      true_pdgid = (std::abs(map_cell_pdgid[key1]) == 13 && std::abs(map_pdgid[i]) == 13) ? 1 : (map_pdgid[i] == -1) ? 2 : 0;
+      edm4hep::Vector3d pos = (barrel_or_endcap == 2) ? edm4hep::Vector3d(ddpos1.x() * 10, ddpos1.y() * 10, map_muonhit[i] * 10) : edm4hep::Vector3d(map_muonhit[i] * 10, ddpos1.y() * 10, ddpos1.z() * 10);
+      Save_trkhit(trkhitVec, key1, true_pdgid, pos);
+    }
+  }
+  Save_onelayer_signal(trkhitVec);
+}
+
+void MuonDigiAlg::Clear()
+{
+  n_hit = 0;
+  n_cell = 0;
+  map_cell_edep.clear();
+  map_cell_adc.clear();
+  map_cell_layer.clear();
+  map_cell_slayer.clear();
+  map_cell_strip.clear();
+  map_cell_fe.clear();
+  map_cell_env.clear();
+  map_cell_pos.clear();
+  map_cell_pdgid.clear();
+  map_cell_time.clear();
+  map_muonhit.clear();
+  map_pdgid.clear();
+}
+
+void MuonDigiAlg::GetSimHit(edm4hep::SimTrackerHit hit_simhit, 
+                            dd4hep::DDSegmentation::BitFieldCoder* hit_m_decoder, 
+                            int hit_i, 
+                            int hit_data[6], 
+                            double & hit_Edep, 
+                            unsigned long long & hit_cellid, 
+                            double & hit_xydist, 
+                            edm4hep::Vector3d & hit_pos, 
+                            dd4hep::Position & hit_ddpos, 
+                            std::array<unsigned long long, 2> & cell_key, 
+                            double & hit_mcpz,
+                            double & hit_time)
+{
+  hit_cellid = hit_simhit.getCellID();
+  auto mcp = hit_simhit.getMCParticle();
+  auto mcppos = mcp.getEndpoint();
+  auto pdgid = mcp.getPDG();    
+  hit_data[0] = std::abs(mcp.getPDG()); //abspdgid    
+  hit_data[1] = hit_m_decoder->get(hit_cellid, "Layer");
+  hit_data[2] = hit_m_decoder->get(hit_cellid, "Superlayer");
+  hit_data[3] = hit_m_decoder->get(hit_cellid, "Stripe");
+  hit_data[4] = hit_m_decoder->get(hit_cellid, "Env");
+  // hit_i = 2 means barrel, hit_i !=2 means endcap
+  hit_data[5] = (hit_i == 2) ? hit_m_decoder->get(hit_cellid, "Fe") : hit_m_decoder->get(hit_cellid, "Endcap");
+  hit_Edep   = (double)hit_simhit.getEDep();//GeV
+  hit_xydist = std::sqrt(mcppos[0]*mcppos[0]+mcppos[1]*mcppos[1]);
+  hit_pos = hit_simhit.getPosition();
+  hit_ddpos = m_cellIDConverter->position(hit_cellid);  
+  cell_key = std::array<unsigned long long, 2>{hit_cellid, (unsigned long long)hit_data[0]};
+  hit_mcpz = std::abs(mcppos[2]);
+  hit_time = (double)hit_simhit.getTime(); // ns
+}
+
+double MuonDigiAlg::Gethit_sipm_length_Barrel(dd4hep::Position hit_ddpos, edm4hep::Vector3d hit_pos, int hit_data[6])
+{
+  // layer 1 is for strips parallel to the beam z-axis
+  if(hit_data[1] == 1)
+  {
+    // We count the number of strips
+    // parpendicular to the beam axis, 
+    // this number times the strip width 4 cm to get 
+    // the length of the strips parallel to the beam axis
+    // hit_data[4] is Envolop, left or right
+    // hit_data[2] is superlayer
+    int tempnum = hit_data[4] + hit_data[2];
+    double strips_count = (tempnum % 2 == 0) ? 115.5 : 106.5; // 115 for long-half-barrel, 106 for short-half-barrel
+    // strip length
+    double length = (strips_count * 4.0) ; // cm
+    // z axis positive half
+    if (hit_ddpos.z()>=0) 
+    {
+
+      // sipm is at the z ~ 0 side
+      // (z<0) =================sipm|z=0|sipm=========c======x=====  (z>0)
+      //  ----->z direction                   center  hit
+
+      return length/2.0 + (hit_pos[2] - hit_ddpos.z());
+    }
+    // z axis negative half
+    else 
+    {
+      // sipm is at the z ~ 0 side
+      // (z<0) =========c====x====sipm|z=0|sipm====================  (z>0)
+      //             center hit      ----->z direction    
+
+      return length/2.0 - (hit_pos[2] - hit_ddpos.z());
+    }
+  } 
+  // layer 2 is for strips perpendicular to the beam axis
+  else if(hit_data[1] == 2)
+  {
+    // We count the number of strips 
+    // parallel to beam direction in each slayer,
+    // this number times the strip width 4cm 
+    // to give the length of strips in each slayer 
+    // parpendicular to the beam axis
+    // strip length
+    double length = strip_length[hit_data[2]-1] * 4.0;
+    
+    // define two vectors
+    TVector2 hit_pos_2vec(hit_pos[0]*0.1, hit_pos[1]*0.1);
+    TVector2 strip_pos_2vec(hit_ddpos.x(), hit_ddpos.y());
+    // define distance is hit - strip center
+    TVector2 distance_2vec = hit_pos_2vec - strip_pos_2vec; 
+    // assume sipm is on the clockwise side
+    //
+    //    ______sipm
+    //   /      \
+    //  /        \sipm
+    // |          |
+    // |     .    |       z>0 point to the out side of the screen
+    // |          |sipm
+    //  \        /
+    //   \______/sipm
+    
+    // define delta_phi as hit - strip_center
+    double dphi = TVector2::Phi_mpi_pi(hit_pos_2vec.Phi()-strip_pos_2vec.Phi());
+    //check delta_phi if it is positive, 
+    // if yes, the hit is on the half of the strip near the sipm
+    // if no, the hit is on the half of the strip farther away from the sipm
+    if (dphi>=0.0) 
+    {
+      return length/2.0 - distance_2vec.Mod();
+    }
+    else 
+    {
+      return length/2.0 + distance_2vec.Mod();
+    }
+  }  
+  // error
+  else 
+  {
+    error() << "Error:: hit_data[1] can only be 1 for parallel to z axis or 2 for perpendicular to z-axis" << endmsg;
+    throw GaudiException("Error:: hit_data[1] can only be 1 for parallel to z axis or 2 for perpendicular to z-axis", "MuonDigiAlg", StatusCode::FAILURE);
+  }
+}
+
+double MuonDigiAlg::Gethit_sipm_length_Endcap(dd4hep::Position hit_ddpos, edm4hep::Vector3d hit_pos, int hit_data[6])
+{
+  // strip length
+  double length = endcap_strip_length[hit_data[3]-1] * 100; // cm
+   
+  // define two vectors
+  TVector2 hit_pos_2vec(hit_pos[0]*0.1, hit_pos[1]*0.1);
+  TVector2 strip_pos_2vec(hit_ddpos.x(), hit_ddpos.y());
+  // define distance is hit - strip center
+  TVector2 distance_2vec = hit_pos_2vec - strip_pos_2vec;
+  // assume sipm is on the clockwise side
+  // define delta_phi as hit - strip_center
+  double dphi = TVector2::Phi_mpi_pi(hit_pos_2vec.Phi()-strip_pos_2vec.Phi());
+  //check delta_phi if it is positive, 
+  // if yes, the hit is on the half of the strip near the sipm
+  // if no, the hit is on the half of the strip farther away from the sipm
+  if (dphi>=0.0)
+  {
+    return length/2.0 - distance_2vec.Mod();
+  }
+  else
+  {
+    return length/2.0 + distance_2vec.Mod();
+  }
+}
+
+double MuonDigiAlg::EdeptoADC(double hit_Edep, double hit_sipm_length)
+{
+  // digitize to ADC 
+  double ADCmean = 34*hit_Edep*47.09/(23*0.00141);//mV
+  if(hit_sipm_length>10)
+  {
+    ADCmean = (16.0813*std::exp(-1*hit_sipm_length/50.8147)+19.5474)*hit_Edep*47.09/(23*0.00141);//mV
+  }
+  return std::max(m_randSvc->generator(Rndm::Landau(ADCmean,7.922))->shoot(), 0.0);
+}
+
+void MuonDigiAlg::SaveData_mapcell(std::array<unsigned long long, 2> cell_key, double hit_Edep, int hit_data[6], edm4hep::Vector3d hit_pos, double hit_time)
+{
+  // check if the key exists, if not initialize be zero
+  if ( map_cell_edep.find(cell_key)==map_cell_edep.end() )
+  {
+    map_cell_edep[cell_key] = 0.0;
+  }
+  // fill the cell loop
+  map_cell_edep[cell_key] += hit_Edep;
+  map_cell_layer[cell_key] = hit_data[1];
+  map_cell_slayer[cell_key] = hit_data[2];
+  map_cell_strip[cell_key] = hit_data[3];
+  map_cell_fe[cell_key] = hit_data[5];
+  map_cell_env[cell_key] = hit_data[4];
+  map_cell_pos[cell_key] = hit_pos;
+  map_cell_pdgid[cell_key] = hit_data[0];
+  // smear the time:
+  map_cell_time[cell_key] = std::max(m_randSvc->generator(Rndm::Gauss(hit_time, m_time_resolution))->shoot(), 0.0);
+}
+
+void MuonDigiAlg::Cut3()
+{
+  //some criteria
+  for(const auto& item : map_cell_edep)
+  {
+    if (map_cell_edep[item.first] <= m_EdepMin || 
+        m_randSvc->generator(Rndm::Flat(0, 1))->shoot() > m_hitEff || 
+        map_cell_adc[item.first] <= m_hit_Edep_min || 
+        (map_cell_adc[item.first] > m_hit_Edep_max && m_hit_Edep_max>0) )
+    {
+      map_cell_edep[item.first] = 0;
+    }
+  }
+}
+
+void MuonDigiAlg::Find_anotherlayer(int barrel_or_endcap, std::array<unsigned long long, 2> key1, std::array<unsigned long long, 2> key2, double ddposi, int & anotherlayer_cell_num)
+{
+  map_cell_strip[key1] = -1;
+  map_cell_strip[key2] = -1;
+  edm4hep::Vector3d pos1 = map_cell_pos[key1];
+  edm4hep::Vector3d pos2 = map_cell_pos[key2];
+  // barrel_or_endcap: 2 is barrel, 0 is endcap
+  // barrel is [2] (z), endcap is [0] (x)
+  double Hit_min = pos1[barrel_or_endcap]*0.1-12;//10cm+2cm
+  double Hit_max = pos1[barrel_or_endcap]*0.1+12;
+  if(ddposi<Hit_min || ddposi>Hit_max) return;
+  map_muonhit[++anotherlayer_cell_num] = ddposi; 
+  map_pdgid[anotherlayer_cell_num] = -1;      
+  if(std::abs(pos1[0]-pos2[0])>40 && std::abs(pos1[1]-pos2[1])>40 && std::abs(pos1[2]-pos2[2])>40) return;
+  map_pdgid[anotherlayer_cell_num] = map_cell_pdgid[key2];
+}
+
+void MuonDigiAlg::Save_trkhit(edm4hep::TrackerHitCollection* trkhitVec, std::array<unsigned long long, 2> key, int pdgid, edm4hep::Vector3d pos)
+{
+  edm4hep::MutableTrackerHit trkHit;
+  trkHit.setEDep((float)map_cell_adc[key]);
+  trkHit.setQuality(pdgid);
+  trkHit.setCellID(key[0]);
+  trkHit.setPosition(pos);
+  trkHit.setTime((float)map_cell_time[key]);
+  trkhitVec->push_back( trkHit );
+}
+
+void MuonDigiAlg::Save_onelayer_signal(edm4hep::TrackerHitCollection* trkhitVec)
+{
+  for (const auto& item : map_cell_strip) 
+  {
+    if(item.second == -1) continue;
+    std::array<unsigned long long, 2> key = item.first;
+    if ( map_cell_edep[key] <= m_EdepMin ) continue;  
+    dd4hep::Position ddpos = m_cellIDConverter->position(key[0]);
+    edm4hep::Vector3d pos = edm4hep::Vector3d(ddpos.x()*10,ddpos.y()*10,ddpos.z()*10);//mm
+    Save_trkhit(trkhitVec, key, 3, pos);  //just one layer signal in a superlayer
+  }
+}

Digitization/DigiMuon/src/MuonDigiAlg.h

+#ifndef MuonDigiAlg_H
+#define MuonDigiAlg_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/NTuple.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+
+#include <DDRec/DetectorData.h>
+#include "DetInterface/IGeomSvc.h"
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/MutableCaloHitContribution.h"
+#include "edm4hep/MutableSimCalorimeterHit.h"
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/CalorimeterHitCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/MCRecoCaloAssociationCollection.h"
+#include "edm4hep/MCRecoCaloParticleAssociationCollection.h"
+
+#include <DDRec/DetectorData.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include <DD4hep/Segmentations.h> 
+#include "DetInterface/IGeomSvc.h"
+#include "TVector3.h"
+#include "TRandom3.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TH3.h"
+#include "TH1.h"
+
+#include <cstdlib>
+#include "time.h"
+#include <TTimeStamp.h> 
+#include <ctime>
+#include <iostream>
+#include <TMath.h>
+#define PI 3.141592653
+class MuonDigiAlg : public GaudiAlgorithm
+{
+ public:
+  
+  MuonDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+ 
+  virtual StatusCode initialize() ;
+  virtual StatusCode execute() ; 
+  virtual StatusCode finalize() ;
+  void MuonHandle(int barrel_or_endcap); // 2->Barrel  0->Endcap
+  void Clear();
+  void GetSimHit(edm4hep::SimTrackerHit hit_simhit, dd4hep::DDSegmentation::BitFieldCoder* hit_m_decoder, int hit_i, 
+                 int hit_data[6], double & hit_Edep, unsigned long long & hit_cellid, double & hit_xydist, 
+                 edm4hep::Vector3d & hit_pos, dd4hep::Position & hit_ddpos, 
+                 std::array<unsigned long long, 2> & cell_key, double & cell_mcpz, double & hit_time);
+  double Gethit_sipm_length_Barrel(dd4hep::Position hit_ddpos, edm4hep::Vector3d hit_pos, int hit_data[6]);
+  double Gethit_sipm_length_Endcap(dd4hep::Position hit_ddpos, edm4hep::Vector3d hit_pos, int hit_data[6]);
+  double EdeptoADC(double hit_Edep, double hit_sipm_length);
+  void SaveData_mapcell(std::array<unsigned long long, 2> cell_key, double hit_Edep, int hit_data[6], edm4hep::Vector3d hit_pos, double hit_time);
+  void Cut3();
+  void Find_anotherlayer(int barrel_or_endcap, std::array<unsigned long long, 2> key1, 
+                         std::array<unsigned long long, 2> key2, double ddposi, int & anotherlayer_cell_num);
+  void Save_trkhit(edm4hep::TrackerHitCollection* trkhitVec, std::array<unsigned long long, 2> key, int pdgid, edm4hep::Vector3d pos);
+  void Save_onelayer_signal(edm4hep::TrackerHitCollection* trkhitVec);
+
+
+
+
+ protected:
+
+  TTree* m_tree;
+  TFile* m_wfile;
+    
+  SmartIF<IGeomSvc> m_geosvc;
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder_barrel;
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder_endcap;
+  dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+  dd4hep::Detector* m_dd4hep;
+
+  std::map<std::array<unsigned long long, 2>, double> map_cell_edep;
+  std::map<std::array<unsigned long long, 2>, double> map_cell_adc;
+  std::map<std::array<unsigned long long, 2>, double> map_cell_time;
+  std::map<std::array<unsigned long long, 2>, int> map_cell_layer;
+  std::map<std::array<unsigned long long, 2>, int> map_cell_slayer;
+  std::map<std::array<unsigned long long, 2>, int> map_cell_strip;
+  std::map<std::array<unsigned long long, 2>, int> map_cell_fe;
+  std::map<std::array<unsigned long long, 2>, int> map_cell_env;
+  std::map<std::array<unsigned long long, 2>, int> map_cell_pdgid;
+  std::map<std::array<unsigned long long, 2>, edm4hep::Vector3d> map_cell_pos;
+  std::map<int, double> map_muonhit;
+  std::map<int, int> map_pdgid;
+
+ 
+#define MAX_SIZE 10000 
+
+  int n_hit, n_cell;
+
+  unsigned long long hit_cellid[MAX_SIZE];
+  double hit_posx[MAX_SIZE];
+  double hit_posy[MAX_SIZE];
+  double hit_posz[MAX_SIZE];
+  double hit_edep[MAX_SIZE];
+  int hit_layer[MAX_SIZE];
+  int hit_slayer[MAX_SIZE];
+  int hit_strip[MAX_SIZE];
+  int hit_fe[MAX_SIZE];
+  int hit_env[MAX_SIZE];
+
+  unsigned long long cell_cellid[MAX_SIZE];
+  double cell_edep[MAX_SIZE];
+  double cell_adc[MAX_SIZE];
+  double cell_posx[MAX_SIZE];
+  double cell_posy[MAX_SIZE];
+  double cell_posz[MAX_SIZE];
+  int cell_layer[MAX_SIZE];
+  int cell_slayer[MAX_SIZE];
+  int cell_strip[MAX_SIZE];
+  int cell_fe[MAX_SIZE];
+  int cell_env[MAX_SIZE];
+  int cell_superlayernumber;
+  double cell_pt;
+
+
+  Gaudi::Property<int>   _writeNtuple{this,  "WriteNtuple", 1, "Write ntuple"};
+  Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"};
+
+
+  Gaudi::Property<double>  m_hitEff{ this, "SiPMEff", 1, "Efficiency of a single hit on a Strip" };
+  Gaudi::Property<double>  m_EdepMin{ this, "EdepMin", 0.0001, "Minimum Edep of a mip" };
+  Gaudi::Property<double>  m_hit_Edep_min{ this, "HitEdepMin", 0.000001, "Minimum Edep of a mip"};  
+  Gaudi::Property<double>  m_hit_Edep_max{ this, "HitEdepMax", 0.1, "Maximum Edep of a mip"};  
+  Gaudi::Property<double>  m_time_resolution{ this, "TimeResolution", 2.0, "Time resolution of a hit, in nano-second ns." };
+
+  // number of strips parallel to beam direction 
+  // in each slayer, each strip width=4cm, so 
+  int strip_length[6] = {26, 38, 50, 62, 74, 86}; //for barrel
+  double endcap_strip_length[193] = {2.12, 2.12, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.19, 2.22,
+                                     2.24, 2.28, 2.32, 2.37, 2.45, 2.65, 2.64, 2.63, 2.62, 2.61,
+                                     2.60, 2.59, 2.57, 2.56, 2.54, 2.53, 2.51, 2.50, 2.48, 2.46, 
+                                     2.44, 2.42, 2.40, 2.38, 2.36, 2.33, 2.31, 2.28, 2.26, 2.23,
+                                     2.20, 2.17, 2.14, 2.11, 2.07, 2.04, 2.00, 1.97, 1.93, 1.89, 
+                                     1.84, 1.80, 1.75, 1.70, 1.65, 1.60, 1.54, 1.48, 1.42, 1.35, 
+                                     1.28, 1.20, 1.12, 1.02, 0.92, 0.80, 0.65, 0.46, 2.20, 2.20,
+                                     2.20, 2.20, 2.20, 2.20, 2.20, 2.21, 2.21, 2.21, 2.21, 2.22,
+                                     2.22, 2.22, 2.23, 2.23, 2.23, 2.24, 2.24, 2.25, 2.25, 2.26,
+                                     2.26, 2.27, 2.28, 2.28, 2.29, 2.30, 2.31, 2.32, 2.32, 2.33,
+                                     2.34, 2.35, 2.36, 2.38, 2.39, 2.40, 2.41, 2.43, 2.44, 2.45,
+                                     2.47, 2.49, 2.50, 2.52, 2.54, 2.56, 2.58, 2.60, 2.62, 2.65, 
+                                     2.67, 2.70, 2.73, 2.76, 2.79, 2.82, 2.86, 2.90, 2.94, 2.99,
+                                     3.04, 3.10, 3.16, 3.23, 3.31, 3.41, 3.53, 3.70, 4.14, 4.12, 
+                                     4.09, 4.06, 4.03, 4.00, 3.97, 3.94, 3.91, 3.87, 3.84, 3.81, 
+                                     3.77, 3.74, 3.70, 3.67, 3.63, 3.59, 3.55, 3.51, 3.47, 3.43, 
+                                     3.38, 3.34, 3.30, 3.25, 3.20, 3.15, 3.10, 3.05, 3.00, 2.95, 
+                                     2.89, 2.83, 2.77, 2.71, 2.65, 2.58, 2.52, 2.45, 2.37, 2.30, 
+                                     2.22, 2.14, 2.05, 1.96, 1.86, 1.76, 1.65, 1.53, 1.40, 1.25, 
+                                     1.09, 0.89, 0.63};  
+
+
+  // Input collections
+  DataHandle<edm4hep::SimTrackerHitCollection>            m_inputMuonBarrel{"MuonBarrelHitsCollection", Gaudi::DataHandle::Reader, this};
+  DataHandle<edm4hep::SimTrackerHitCollection>            m_inputMuonEndcap{"MuonEndcapHitsCollection", Gaudi::DataHandle::Reader, this};
+  // Output collections
+  DataHandle<edm4hep::TrackerHitCollection>               m_outputMuonBarrel{"MuonBarrelTrackerHits", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::TrackerHitCollection>               m_outputMuonEndcap{"MuonEndcapTrackerHits", Gaudi::DataHandle::Writer, this};
+
+  //DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assMuonBarrel{"MuonBarrelTrackerHitAssociationCollection", Gaudi::DataHandle::Writer, this};
+  SmartIF<IRndmGenSvc> m_randSvc;
+  edm4hep::TrackerHitCollection* trkhitVec;
+  int m_nEvt=0;
+};
+#endif
+
+
+

Digitisers/SimpleDigi/CMakeLists.txt → Digitization/DigiSimple/CMakeLists.txt

@@ -3,7 +3,13 @@ gaudi_add_module(SimpleDigi
                  SOURCES src/PlanarDigiAlg.cpp
                          src/TPCDigiAlg.cpp
                          src/voxel.cpp
+                         src/CylinderDigiAlg.cpp
+                         src/SiTrackerDigiAlg.cpp
+                         src/TPCPixelDigiAlg.cpp
+                         src/TPCPixelDigiTool.cpp
+                         src/SmearDigiTool.cpp
                  LINK GearSvc
+                      DigiTool
                       EventSeeder
                       TrackSystemSvcLib
                       DataHelperLib

Digitization/DigiSimple/src/CylinderDigiAlg.cpp

+#include "CylinderDigiAlg.h"
+
+#include "DetIdentifier/CEPCConf.h"
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/Vector3D.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+
+DECLARE_COMPONENT( CylinderDigiAlg )
+
+CylinderDigiAlg::CylinderDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc){
+  // Input collections
+  declareProperty("SimTrackHitCollection", m_inputColHdls, "Handle of the Input SimTrackerHit collection");
+
+  // Output collections
+  declareProperty("TrackerHitCollection", m_outputColHdls, "Handle of the output TrackerHit collection");
+  declareProperty("TrackerHitAssociationCollection", m_assColHdls, "Handle of the Association collection between SimTrackerHit and TrackerHit");
+}
+
+StatusCode CylinderDigiAlg::initialize(){
+  m_geosvc = service<IGeomSvc>("GeomSvc");
+  if(!m_geosvc){
+    error() << "Failed to get the GeomSvc" << endmsg;
+    return StatusCode::FAILURE;
+  }
+  auto detector = m_geosvc->lcdd();
+  if(!detector){
+    error() << "Failed to get the Detector from GeomSvc" << endmsg;
+    return StatusCode::FAILURE;
+  }
+  std::string name = m_inputColHdls.objKey(); 
+  debug() << "Readout name: " << name << endmsg;
+  m_decoder = m_geosvc->getDecoder(name);
+  if(!m_decoder){
+    error() << "Failed to get the decoder. " << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  info() << "CylinderDigiAlg::initialized" << endmsg;
+  return GaudiAlgorithm::initialize();
+}
+
+
+StatusCode CylinderDigiAlg::execute(){
+  auto trkhitVec = m_outputColHdls.createAndPut();
+  auto assVec = m_assColHdls.createAndPut();
+
+  const edm4hep::SimTrackerHitCollection* STHCol = nullptr;
+  try {
+    STHCol = m_inputColHdls.get();
+  }
+  catch(GaudiException &e){
+    debug() << "Collection " << m_inputColHdls.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  if(STHCol->size()==0) return StatusCode::SUCCESS;
+  debug() << m_inputColHdls.fullKey() << " has SimTrackerHit "<< STHCol->size() << endmsg;
+  
+  for(auto simhit : *STHCol){
+    auto particle = simhit.getMCParticle();
+    if(!particle.isAvailable()) continue;
+    
+    auto& mom0 = particle.getMomentum();
+    double pt = sqrt(mom0.x*mom0.x+mom0.y*mom0.y);
+    if(particle.isCreatedInSimulation()&&pt<0.01&&particle.isStopped()) continue;
+
+    auto cellId = simhit.getCellID();
+    int system  = m_decoder->get(cellId, "system");
+    int layer   = m_decoder->get(cellId, "layer"  );
+    int module  = m_decoder->get(cellId, "module");
+    int sensor  = m_decoder->get(cellId, "sensor"  );
+    auto& pos   = simhit.getPosition();
+    auto& mom   = simhit.getMomentum();
+    
+    double phi = atan2(pos[1], pos[0]);
+    double r   = sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
+    double dphi = m_resRPhi/r;
+    phi += randSvc()->generator(Rndm::Gauss(0, dphi))->shoot();
+    double smearedX = r*cos(phi);
+    double smearedY = r*sin(phi);
+    double smearedZ = pos[2] + randSvc()->generator(Rndm::Gauss(0, m_resZ))->shoot();
+
+    auto trkHit = trkhitVec->create();
+    trkHit.setCellID(cellId);
+    trkHit.setTime(simhit.getTime());
+    trkHit.setEDep(simhit.getEDep());
+    trkHit.setPosition (edm4hep::Vector3d(smearedX, smearedY, smearedZ));
+    trkHit.setCovMatrix(std::array<float, 6>{m_resRPhi*m_resRPhi/2, 0, m_resRPhi*m_resRPhi/2, 0, 0, m_resZ*m_resZ});
+    trkHit.setType(1<<CEPCConf::TrkHitTypeBit::CYLINDER);
+    trkHit.addToRawHits(simhit.getObjectID());
+    debug() << "Hit " << simhit.id() << ": " << pos << " -> " << trkHit.getPosition() << "s:" << system << " l:" << layer << " m:" << module << " s:" << sensor
+	    << " pt = " << pt << " " << mom.x << " " << mom.y << " " << mom.z << endmsg;
+
+    auto ass = assVec->create();
+
+    float weight = 1.0;
+
+    debug() <<" Set relation between " << " sim hit " << simhit.id() << " to tracker hit " << trkHit.id() << " with a weight of " << weight << endmsg;
+    ass.setSim(simhit);
+    ass.setRec(trkHit);
+    ass.setWeight(weight);
+  }
+
+  m_nEvt++;
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode CylinderDigiAlg::finalize(){
+  info() << "Processed " << m_nEvt << " events " << endmsg;
+  return GaudiAlgorithm::finalize();
+}

Digitization/DigiSimple/src/CylinderDigiAlg.h

+#ifndef CylinderDigiAlg_H
+#define CylinderDigiAlg_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/NTuple.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+
+#include <DDRec/DetectorData.h>
+#include "DetInterface/IGeomSvc.h"
+
+class CylinderDigiAlg : public GaudiAlgorithm{
+ public:
+  
+  CylinderDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+ 
+  virtual StatusCode initialize() ;
+  virtual StatusCode execute() ; 
+  virtual StatusCode finalize() ;
+ 
+protected:
+  
+  SmartIF<IGeomSvc> m_geosvc;
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+
+  // Input collections
+  DataHandle<edm4hep::SimTrackerHitCollection>            m_inputColHdls{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};
+  // Output collections
+  DataHandle<edm4hep::TrackerHitCollection>               m_outputColHdls{"DCTrackerHits", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assColHdls{"DCTrackerHitAssociationCollection", Gaudi::DataHandle::Writer, this};
+
+  Gaudi::Property<float> m_resRPhi{this, "ResRPhi", 0.1};
+  Gaudi::Property<float> m_resZ   {this, "ResZ", 2.828};
+
+  int m_nEvt=0;
+};
+#endif

Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp → Digitization/DigiSimple/src/PlanarDigiAlg.cpp

@@ -17,7 +17,7 @@
 #include "UTIL/CellIDEncoder.h"
 #include <UTIL/Operators.h>
 */
-#include "Identifier/CEPCConf.h"
+#include "DetIdentifier/CEPCConf.h"
 #include "UTIL/ILDConf.h"
 
 // STUFF needed for GEAR
@@ -54,13 +54,20 @@ PlanarDigiAlg::PlanarDigiAlg(const std::string& name, ISvcLocator* svcLoc)
 
 StatusCode PlanarDigiAlg::initialize()
 {
+  if (_parameterize) {
+    if (_parU.size()!=10 || _parV.size()!=10) {
+      fatal() << "parameters number must be 10! now " << _parU.size() << " for U and " << _parV.size() << " for V" << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+  else {
+    if (_resU.size() !=  _resV.size()) {
+      fatal() << "Inconsistent number of resolutions given for U and V coordinate: " 
+              << "ResolutionU  :" <<   _resU.size() << " != ResolutionV : " <<  _resV.size()
+              << endmsg;
 
-  if( _resU.size() !=  _resV.size() ) {
-    fatal() << "Inconsistent number of resolutions given for U and V coordinate: " 
-      << "ResolutionU  :" <<   _resU.size() << " != ResolutionV : " <<  _resV.size()
-      << endmsg;
-
-    return StatusCode::FAILURE;
+      return StatusCode::FAILURE;
+    }
   }
 
   // get the GEAR manager
@@ -85,16 +92,16 @@ StatusCode PlanarDigiAlg::initialize()
   //TODO:trksystem->init() ;
   //FIXME:SJA gear surface store has now been filled so we can dispose of the MarlinTrkSystem
   //TODO:delete trksystem;
-  
+  // fucd: fix TODO - surface store is needed to fill once always if does not handle tracking algorithm in job
   auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc");
   if ( !_trackSystemSvc ) {
     error() << "Failed to find TrackSystemSvc ..." << endmsg;
     return StatusCode::FAILURE;
   }
-  
+
   MarlinTrk::IMarlinTrkSystem* _trksystem =  _trackSystemSvc->getTrackSystem(this);
   _trksystem->init();
-  
+
   _trackSystemSvc->removeTrackSystem(this);
   
   return GaudiAlgorithm::initialize();
@@ -148,10 +155,11 @@ StatusCode PlanarDigiAlg::execute()
     encoder.setValue(SimTHit.getCellID()) ;
     det_id  = encoder[lcio::ILDCellID0::subdet] ;
   
-    if     ( det_id == lcio::ILDDetID::VXD ){}
-    else if( det_id == lcio::ILDDetID::SIT ){}
-    else if( det_id == lcio::ILDDetID::SET ){}
-    else if( det_id == lcio::ILDDetID::FTD ){}
+    if     ( det_id == CEPCConf::DetID::VXD ){}
+    else if( det_id == CEPCConf::DetID::SIT ){}
+    else if( det_id == CEPCConf::DetID::SET ){}
+    else if( det_id == CEPCConf::DetID::FTD ){}
+    else if( det_id == CEPCConf::DetID::ETD ){}
     else {
       fatal() << "unsupported detector ID = " << det_id << " CellID = " << SimTHit.getCellID()
               << ": file " << __FILE__ << " line " << __LINE__
@@ -166,6 +174,9 @@ StatusCode PlanarDigiAlg::execute()
 
   int i = 0;
   for( auto SimTHit : *STHcol ) {
+    if (SimTHit.getEDep()<=_eThreshold) continue;
+    if (gsl_ran_flat(_rng, 0, 1)>_efficiency) continue;
+    debug() << "MCParticle id " << SimTHit.getMCParticle().id() << endmsg;
 
     const int celId = SimTHit.getCellID() ;
 
@@ -205,28 +216,64 @@ StatusCode PlanarDigiAlg::execute()
     //                << endmsg;
     //        continue;
     //      }
-    if( (_resU.size() > 1 && layer > _resU.size()-1) || (_resV.size() > 1 && layer > _resV.size()-1) ) {
-      fatal() << "layer exceeds resolution vector, please check input parameters ResolutionU and ResolutionV" << endmsg;
-      return StatusCode::FAILURE;
-    }
+    gear::MeasurementSurface const* ms = _GEAR->getMeasurementSurfaceStore().GetMeasurementSurface( encoder.lowWord() );
+    gear::CartesianCoordinateSystem* cartesian = dynamic_cast< gear::CartesianCoordinateSystem* >( ms->getCoordinateSystem() );
+    CLHEP::Hep3Vector uVec = cartesian->getLocalXAxis();
+    CLHEP::Hep3Vector vVec = cartesian->getLocalYAxis();
+
+    float u_direction[2] ;
+    u_direction[0] = uVec.theta();
+    u_direction[1] = uVec.phi();
+    
+    float v_direction[2] ;
+    v_direction[0] = vVec.theta();
+    v_direction[1] = vVec.phi();
+    
+    debug() << " U[0] = "<< u_direction[0] << " U[1] = "<< u_direction[1]
+            << " V[0] = "<< v_direction[0] << " V[1] = "<< v_direction[1]
+            << endmsg ;
+    
+    float resU(0), resV(0), resT(0);
 
-    float resU = ( _resU.size() > 1 ?   _resU.value().at(  layer )     : _resU.value().at(0)   )  ;
-    float resV = ( _resV.size() > 1 ?   _resV.value().at(  layer )     : _resV.value().at(0)   )  ; 
+    if (!_parameterize) {
+      if( (_resU.size() > 1 && layer > _resU.size()-1) || (_resV.size() > 1 && layer > _resV.size()-1) ) {
+        fatal() << "layer exceeds resolution vector, please check input parameters ResolutionU and ResolutionV" << endmsg;
+        return StatusCode::FAILURE;
+      }
+
+      resU = ( _resU.size() > 1 ?   _resU.value().at(  layer )     : _resU.value().at(0)   )  ;
+      resV = ( _resV.size() > 1 ?   _resV.value().at(  layer )     : _resV.value().at(0)   )  ; 
+    }
+    else { // Riccardo's parameterized model
+      auto mom = SimTHit.getMomentum();
+      CLHEP::Hep3Vector momVec(mom[0], mom[1], mom[2]);
+      const double alpha = uVec.azimAngle(momVec, vVec);
+      const double cotanAlpha = 1./tan(alpha);
+      // TODO: title angle (PI/2), magnetic field (3) 
+      const double tanLorentzAngle = (side==0) ? 0. : 0.053 * 3 * cos(M_PI/2.);
+      const double x = fabs(-cotanAlpha - tanLorentzAngle);
+      resU = _parU[0] + _parU[1] * x + _parU[2] * exp(-_parU[9] * x) * cos(_parU[3] * x + _parU[4])
+        + _parU[5] * exp(-0.5 * pow(((x - _parU[6]) / _parU[7]), 2)) + _parU[8] * pow(x, 0.5);
+      
+      const double beta = vVec.azimAngle(momVec, uVec);
+      const double cotanBeta = 1./tan(beta);
+      const double y = fabs(-cotanBeta);
+      resV = _parV[0] + _parV[1] * y + _parV[2] * exp(-_parV[9] * y) * cos(_parV[3] * y + _parV[4])
+        + _parV[5] * exp(-0.5 * pow(((y - _parV[6]) / _parV[7]), 2)) + _parV[8] * pow(y, 0.5);
+    }
+    // parameterize only for position now, todo
+    resT = ( _resT.size() > 1 ? _resT.value().at(layer) : _resT.value().at(0) );
+    // from ps (input unit) to ns (record unit, Geant4)
+    resT *= CLHEP::ps/CLHEP::ns;
 
     debug() << " --- will smear hit with resU = " << resU << " and resV = " << resV << endmsg;
 
     auto& pos =  SimTHit.getPosition() ;  
 
-    //edm4hep::Vector3d smearedPos;
-
-    //GearSurfaces::MeasurementSurface* ms = _GEAR->getMeasurementSurfaceStore().GetMeasurementSurface( SimTHit->getCellID0() );
-    
-    gear::MeasurementSurface const* ms = _GEAR->getMeasurementSurfaceStore().GetMeasurementSurface( encoder.lowWord() );;
     CLHEP::Hep3Vector globalPoint(pos[0],pos[1],pos[2]);
     CLHEP::Hep3Vector localPoint = ms->getCoordinateSystem()->getLocalPoint(globalPoint);
     CLHEP::Hep3Vector localPointSmeared = localPoint;
 
-
     debug() << std::setprecision(8) << "Position of hit before smearing global: ( "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<< " ) "
             << "local: ( " << localPoint.x() << " " << localPoint.y() << " " << localPoint.z() << " )" << endmsg;
 
@@ -252,11 +299,13 @@ StatusCode PlanarDigiAlg::execute()
                 << " sensor" << sensor << " : retries " << tries << endmsg;
       }
 
-      localPointSmeared.setX( localPoint.x() + gsl_ran_gaussian(_rng, resU) );
-      localPointSmeared.setY( localPoint.y() + gsl_ran_gaussian(_rng, resV) );
+      double dx = gsl_ran_gaussian(_rng, resU);
+      double dy = gsl_ran_gaussian(_rng, resV);
+      localPointSmeared.setX( localPoint.x() + dx );
+      localPointSmeared.setY( localPoint.y() + dy );
 
       //check if hit is in boundaries
-      if ( ms->isLocalInBoundary( localPointSmeared ) ) {
+      if ( ms->isLocalInBoundary( localPointSmeared ) && fabs(dx)<=_maxPull*resU && fabs(dy)<=_maxPull*resV ) {
         accept_hit = true;
         break;
       }
@@ -281,10 +330,6 @@ StatusCode PlanarDigiAlg::execute()
             << localPointSmeared.x() << " " << localPointSmeared.y() << " " << localPointSmeared.z() << " )"
             << endmsg;
     
-    //smearedPos[0] = globalPointSmeared.x();
-    //smearedPos[1] = globalPointSmeared.y();
-    //smearedPos[2] = globalPointSmeared.z();
-
     //make the TrackerHitPlaneImpl
     auto trkHit = trkhitVec->create();
 
@@ -292,7 +337,7 @@ StatusCode PlanarDigiAlg::execute()
 
     edm4hep::Vector3d smearedPos(globalPointSmeared.x(), globalPointSmeared.y(), globalPointSmeared.z());
     trkHit.setPosition( smearedPos ) ;
-
+    /*
     gear::CartesianCoordinateSystem* cartesian = dynamic_cast< gear::CartesianCoordinateSystem* >( ms->getCoordinateSystem() ); 
     CLHEP::Hep3Vector uVec = cartesian->getLocalXAxis();
     CLHEP::Hep3Vector vVec = cartesian->getLocalYAxis();
@@ -308,6 +353,7 @@ StatusCode PlanarDigiAlg::execute()
     debug() << " U[0] = "<< u_direction[0] << " U[1] = "<< u_direction[1]
             << " V[0] = "<< v_direction[0] << " V[1] = "<< v_direction[1]
             << endmsg ;
+    */
     // fucd: next TODO: cov[0] = resU*reU, cov[2] = resV*resV, cov[5] = 0
     if(_usePlanarTag){
       std::array<float, 6> cov;
@@ -336,9 +382,10 @@ StatusCode PlanarDigiAlg::execute()
     }
 
     if( _isStrip || (resU!=0&&resV==0) ){
-        trkHit.setType( UTIL::set_bit( trkHit.getType() , UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ) ) ;
+        trkHit.setType( UTIL::set_bit( trkHit.getType() , CEPCConf::TrkHitTypeBit::ONE_DIMENSIONAL ) ) ;
     }
     trkHit.setEDep( SimTHit.getEDep() );
+    trkHit.setTime( SimTHit.getTime() + gsl_ran_gaussian(_rng, resT) );
 
     // make the relation
     auto rel = relCol->create();
@@ -374,5 +421,8 @@ StatusCode PlanarDigiAlg::execute()
 StatusCode PlanarDigiAlg::finalize()
 {
   info() << "Processed " << _nEvt << " events " << endmsg;
+
+  if(_rng) gsl_rng_free(_rng);
+
   return GaudiAlgorithm::finalize();
 }

Digitisers/SimpleDigi/src/PlanarDigiAlg.h → Digitization/DigiSimple/src/PlanarDigiAlg.h

@@ -77,18 +77,26 @@ protected:
   Gaudi::Property<FloatVec> _resU{ this, "ResolutionU", {0.0040} };
   // resolution in direction of v - either one per layer or one for all layers
   Gaudi::Property<FloatVec> _resV{ this, "ResolutionV", {0.0040} };
+  // resolution of t - either one per layer or one for all layers, unit - ps
+  Gaudi::Property<FloatVec> _resT{ this, "ResolutionT", {0.} };
   // whether hits are 1D strip hits
   Gaudi::Property<bool> _isStrip{ this, "IsStrip", false };
   // whether use Planar tag for type and cov, if true, CEPCConf::TrkHitTypeBit::PLANAR bit is set as true
   // cov[0]=thetaU, cov[1]=phiU, cov[2]=resU, cov[0]=thetaV, cov[1]=phiV, cov[2]=resV
   Gaudi::Property<bool> _usePlanarTag{ this, "UsePlanarTag", true };
+  Gaudi::Property<float> _eThreshold{ this, "EnergyThreshold", 0 };
+  Gaudi::Property<float> _efficiency{ this, "Efficiency", 1 };
+  Gaudi::Property<float> _maxPull{ this, "PullCutToRetry", 1000. };
+  Gaudi::Property<bool> _parameterize{ this, "ParameterizeResolution", false};
+  Gaudi::Property<FloatVec> _parU{ this, "ParametersU", {0} };
+  Gaudi::Property<FloatVec> _parV{ this, "ParametersV", {0} };
 
   // Input collections
   DataHandle<edm4hep::EventHeaderCollection> _headerCol{"EventHeaderCol", Gaudi::DataHandle::Reader, this};
   DataHandle<edm4hep::SimTrackerHitCollection> _inColHdl{"VXDCollection", Gaudi::DataHandle::Reader, this};
   // Output collections
   DataHandle<edm4hep::TrackerHitCollection> _outColHdl{"VXDTrackerHits", Gaudi::DataHandle::Writer, this};
-  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _outRelColHdl{"VTXTrackerHitRelations", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _outRelColHdl{"VXDTrackerHitRelations", Gaudi::DataHandle::Writer, this};
 };
 
 #endif

Digitization/DigiSimple/src/SiTrackerDigiAlg.cpp

+#include "SiTrackerDigiAlg.h"
+
+#include "DetIdentifier/CEPCConf.h"
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/Vector3D.h"
+#include "DDRec/ISurface.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+
+DECLARE_COMPONENT(SiTrackerDigiAlg)
+
+SiTrackerDigiAlg::SiTrackerDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc) {
+  // Input collections
+  declareProperty("SimTrackHitCollection", m_inputColHdls, "Handle of the Input SimTrackerHit collection");
+
+  // Output collections
+  declareProperty("TrackerHitCollection", m_outputColHdls, "Handle of the output TrackerHit collection");
+  declareProperty("TrackerHitAssociationCollection", m_assColHdls, "Handle of the Association collection between SimTrackerHit and TrackerHit");
+}
+
+StatusCode SiTrackerDigiAlg::initialize() {
+  m_digiTool = ToolHandle<IDigiTool>(m_digiToolName.value());
+
+  info() << "DigiTool " << m_digiTool.typeAndName() << " found" << endmsg;
+
+  info() << "SiTrackerDigiAlg::initialized" << endmsg;
+  return GaudiAlgorithm::initialize();
+}
+
+
+StatusCode SiTrackerDigiAlg::execute(){
+  auto trkCol = m_outputColHdls.createAndPut();
+  auto assCol = m_assColHdls.createAndPut();
+
+  const edm4hep::SimTrackerHitCollection* simCol = nullptr;
+  try {
+    simCol = m_inputColHdls.get();
+  }
+  catch(GaudiException &e){
+    debug() << "Collection " << m_inputColHdls.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  if (simCol->size() == 0) return StatusCode::SUCCESS;
+  debug() << m_inputColHdls.fullKey() << " has SimTrackerHit "<< simCol->size() << endmsg;
+
+  m_digiTool->Call(simCol, trkCol, assCol);
+
+  debug() << "Created " << trkCol->size() << " hits, "
+          << simCol->size() - trkCol->size() << " hits got dismissed for being out of boundary"
+          << endmsg;
+
+  m_nEvt++;
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode SiTrackerDigiAlg::finalize(){
+  info() << "Processed " << m_nEvt << " events " << endmsg;
+  return GaudiAlgorithm::finalize();
+}

Digitization/DigiSimple/src/SiTrackerDigiAlg.h

+#ifndef SiTrackerDigiAlg_H
+#define SiTrackerDigiAlg_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/NTuple.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+
+#include "DigiTool/IDigiTool.h"
+
+class SiTrackerDigiAlg : public GaudiAlgorithm{
+ public:
+  
+  SiTrackerDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+ 
+  virtual StatusCode initialize();
+  virtual StatusCode execute();
+  virtual StatusCode finalize();
+ 
+protected:
+  // Input collections
+  DataHandle<edm4hep::SimTrackerHitCollection>            m_inputColHdls{"VXDCollection", Gaudi::DataHandle::Reader, this};
+  // Output collections
+  DataHandle<edm4hep::TrackerHitCollection>               m_outputColHdls{"VXDTrackerHits", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assColHdls{"VXDTrackerHitAssociationCollection", Gaudi::DataHandle::Writer, this};
+
+  Gaudi::Property<std::string> m_digiToolName{this, "DigiTool", "SmearDigiTool/VXD"};
+
+  ToolHandle<IDigiTool> m_digiTool;
+  int m_nEvt=0;
+};
+#endif

Digitization/DigiSimple/src/SmearDigiTool.cpp

+#include "SmearDigiTool.h"
+
+#include "DataHelper/TrackerHitHelper.h"
+#include "DetIdentifier/CEPCConf.h"
+#include "DetInterface/IGeomSvc.h"
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/ISurface.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/RndmGenerators.h"
+
+#include "CLHEP/Vector/ThreeVector.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+
+DECLARE_COMPONENT(SmearDigiTool)
+
+StatusCode SmearDigiTool::initialize() {
+  if (m_parameterize) {
+    if (m_parU.size()!=10 || m_parV.size()!=10) {
+      fatal() << "parameters number must be 10! now " << m_parU.size() << " for U and " << m_parV.size() << " for V" << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+  else {
+    if (m_resU.size() != m_resV.size()) {
+      fatal() << "Inconsistent number of resolutions given for U and V coordinate: "
+              << "ResolutionU  :" << m_resU.size() << " != ResolutionV : " << m_resV.size()
+              << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+
+  m_geosvc = service<IGeomSvc>("GeomSvc");
+  if(!m_geosvc){
+    error() << "Failed to get the GeomSvc" << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  if (m_detName=="") {
+    std::string toolName = name();
+    m_detName = toolName.substr(toolName.find(".")+1);
+    debug() << toolName << " --> " << m_detName.value() << endmsg;
+  }
+  m_surfaces = m_geosvc->getSurfaceMap(m_detName.value());
+  debug() << "Surface map size of " << m_detName.value() << ": " << m_surfaces->size() << endmsg;
+  if (msgLevel(MSG::VERBOSE)) {
+    unsigned long old = 0; 
+    for (const auto& pair : *m_surfaces) {
+      verbose() << pair.first << " : " << pair.second << endmsg;
+      if (old==pair.first) error() << old << " repeat!" << endmsg;
+      old = pair.first;
+    }
+  }
+
+  if (m_readoutName=="") m_readoutName = m_detName.value() + "Collection";
+  debug() << "Readout name: " << m_readoutName.value() << endmsg;
+  m_decoder = m_geosvc->getDecoder(m_readoutName.value());
+  if(!m_decoder){
+    error() << "Failed to get the decoder. " << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  m_randSvc = service<IRndmGenSvc>("RndmGenSvc");
+
+  info() << "initialized" << endmsg;
+  return StatusCode::SUCCESS;
+}
+
+StatusCode SmearDigiTool::Call(const edm4hep::SimTrackerHitCollection* simCol, edm4hep::TrackerHitCollection* hitCol,
+			       edm4hep::MCRecoTrackerAssociationCollection* assCol) {
+  for (auto simhit : *simCol) {
+    StatusCode sc = Call(simhit, hitCol, assCol);
+    if (sc.isFailure()) return sc;
+  }
+  return StatusCode::SUCCESS;
+}
+
+StatusCode SmearDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection* hitCol, edm4hep::MCRecoTrackerAssociationCollection* assCol) {
+  if (!simhit.isAvailable()) {
+    error() << "input SimTrackerHit not available!" << endmsg;
+    return StatusCode::SUCCESS;
+  }
+
+  auto e = simhit.getEDep();
+  if (e < m_eThreshold) return StatusCode::SUCCESS;
+  if (m_randSvc->generator(Rndm::Flat(0, 1))->shoot() > m_efficiency) return StatusCode::SUCCESS;
+  auto t = simhit.getTime();
+
+  auto cellId = simhit.getCellID();
+  int system  = m_decoder->get(cellId, CEPCConf::DetCellID::system);
+  int side    = m_decoder->get(cellId, CEPCConf::DetCellID::side);
+  int layer   = m_decoder->get(cellId, CEPCConf::DetCellID::layer);
+  int module  = m_decoder->get(cellId, CEPCConf::DetCellID::module);
+  int sensor  = m_decoder->get(cellId, CEPCConf::DetCellID::sensor);
+
+  auto& pos   = simhit.getPosition();
+  auto& mom   = simhit.getMomentum();
+
+  debug() << "Hit " << simhit.id() << " cell: " << cellId << " d:" << system << " l:" << layer << " m:" << module << " s:" << sensor
+	  << " p = " << mom.x << ", " << mom.y << ", " << mom.z << " e:" << e << " t:" << t << endmsg;
+
+  dd4hep::rec::ISurface* surface = nullptr;
+  auto it = m_surfaces->find(cellId);
+  if (it != m_surfaces->end()) {
+    surface = it->second;
+    if (!surface) {
+      fatal() << "found surface for cell id " << cellId << ", but NULL" << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+  else {
+    fatal() << "not found surface for cell id " << cellId << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  // CLHEP::mm is divided while Edm4hepWriterAnaElemTool write, so pos without unit
+  dd4hep::rec::Vector3D oldPos(pos.x*dd4hep::mm, pos.y*dd4hep::mm, pos.z*dd4hep::mm);
+  dd4hep::rec::Vector3D uVec = surface->u(oldPos);
+  dd4hep::rec::Vector3D vVec = surface->v(oldPos);
+  float u_direction[2];
+  u_direction[0] = uVec.theta();
+  u_direction[1] = uVec.phi();
+
+  float v_direction[2];
+  v_direction[0] = vVec.theta();
+  v_direction[1] = vVec.phi();
+
+  debug() << " U: " << uVec << endmsg;
+  debug() << " V: " << vVec << endmsg;
+  debug() << " N: " << surface->normal() << endmsg;
+  debug() << " O: " << surface->origin() << endmsg;
+
+  float resU(0), resV(0), resT(0);
+  if (!m_parameterize) {
+    if ((m_resU.size() > 1 && layer >= (int)m_resU.size()) || (m_resV.size() > 1 && layer >= (int)m_resV.size())) {
+      fatal() << "layer exceeds resolution vector, please check input parameters ResolutionU and ResolutionV" << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+    resU = (m_resU.size() > 1 ? m_resU.value().at(layer) : m_resU.value().at(0));
+    resV = (m_resV.size() > 1 ? m_resV.value().at(layer) : m_resV.value().at(0));
+  }
+  else { // Riccardo's parameterized model
+    CLHEP::Hep3Vector momVec(mom[0], mom[1], mom[2]);
+    CLHEP::Hep3Vector uVecCLHEP(uVec[0], uVec[1], uVec[2]);
+    CLHEP::Hep3Vector vVecCLHEP(vVec[0], vVec[1], vVec[2]);
+    const double alpha = uVecCLHEP.azimAngle(momVec, vVecCLHEP);
+    const double cotanAlpha = 1./tan(alpha);
+    // TODO: title angle (PI/2), magnetic field (3)
+    const double tanLorentzAngle = (side==0) ? 0. : 0.053 * 3 * cos(M_PI/2.);
+    const double x = fabs(-cotanAlpha - tanLorentzAngle);
+    resU = m_parU[0] + m_parU[1] * x + m_parU[2] * exp(-m_parU[9] * x) * cos(m_parU[3] * x + m_parU[4])
+      + m_parU[5] * exp(-0.5 * pow(((x - m_parU[6]) / m_parU[7]), 2)) + m_parU[8] * pow(x, 0.5);
+
+    const double beta = vVecCLHEP.azimAngle(momVec, uVecCLHEP);
+    const double cotanBeta = 1./tan(beta);
+    const double y = fabs(-cotanBeta);
+    resV = m_parV[0] + m_parV[1] * y + m_parV[2] * exp(-m_parV[9] * y) * cos(m_parV[3] * y + m_parV[4])
+      + m_parV[5] * exp(-0.5 * pow(((y - m_parV[6]) / m_parV[7]), 2)) + m_parV[8] * pow(y, 0.5);
+  }
+  resU *= dd4hep::mm;
+  resV *= dd4hep::mm;
+  // parameterize only for position now, todo
+  resT = (m_resT.size() > 1 ? m_resT.value().at(layer) : m_resT.value().at(0));
+  // from ps (input unit) to ns (record unit, Geant4)
+  resT *= CLHEP::ps/CLHEP::ns;
+  debug() << " --- will smear hit with resU = " << resU/dd4hep::mm << " resV = " << resV/dd4hep::mm << " resT = " << resT << endmsg;
+
+  auto& typeSurface = surface->type();
+  debug() << " Surface id: " << surface->id() << " type:" << endmsg;
+  debug() << " " << typeSurface << endmsg;
+  verbose() << "      count: " << m_surfaces->count(cellId) << " inner: " << surface->innerMaterial().name()
+	    << " outer: " << surface->outerMaterial().name() << endmsg;
+
+  if (typeSurface.isPlane() || typeSurface.isCylinder()) {
+    // scale to cov at edge
+    double scale = 1.0;
+
+    dd4hep::rec::Vector2D localPoint = surface->globalToLocal(oldPos);
+
+    // for planar, same calculation by local is ok, but reduce repeat
+    if (typeSurface.isCylinder()) {
+      dd4hep::rec::Vector3D mom_ddrec(mom.x*dd4hep::GeV, mom.y*dd4hep::GeV, mom.z*dd4hep::GeV);
+      double                length    = simhit.getPathLength()*dd4hep::mm;
+      dd4hep::rec::Vector3D pre       = oldPos - (0.5*length)*mom_ddrec.unit();
+      dd4hep::rec::Vector3D post      = oldPos + (0.5*length)*mom_ddrec.unit();
+      dd4hep::rec::Vector2D localPre  = surface->globalToLocal(pre);
+      dd4hep::rec::Vector2D localPost = surface->globalToLocal(post);
+      localPoint = dd4hep::rec::Vector2D(0.5*(localPre.u()+localPost.u()), 0.5*(localPre.v()+localPost.v()));
+      debug() << "pre: (" << pre.x() << " " << pre.y() << " " << pre.z() << " ) local (" << localPre.u() << ", " << localPre.v() << " ) "
+	      << "post: (" << post.x() << " " << post.y() << " " << post.z() << " ) local (" << localPost.u() << ", " << localPost.v() << " ) " << endmsg;
+    }
+    //dd4hep::rec::Vector3D local3D(localPoint.u(), localPoint.v(), 0);
+    // A small check, if the hit is in the boundaries:
+    if (!surface->insideBounds(oldPos)) {
+      double dSToHit = surface->distance(oldPos);
+      debug() << " global: (" << oldPos.x()/dd4hep::mm << " " << oldPos.y()/dd4hep::mm << " " << oldPos.z()/dd4hep::mm
+	      << ") local: (" << localPoint.u()/dd4hep::mm << ", " << localPoint.v()/dd4hep::mm << " )"
+	      << " distance: " << dSToHit/dd4hep::mm
+	      << " is not within boundaries." << endmsg;
+
+      // FIXME: change position to path at center plane? or enlarge cov? other tracker
+      if (system == CEPCConf::DetID::OTKBarrel || system == CEPCConf::DetID::OTKEndcap) {
+#if 0
+	dd4hep::rec::Vector3D global3D = oldPos;
+	dd4hep::rec::Vector3D mom_ddrec(mom.x*dd4hep::GeV, mom.y*dd4hep::GeV, mom.z*dd4hep::GeV);
+	double                length = simhit.getPathLength()*dd4hep::mm;
+	dd4hep::rec::Vector3D pre    = oldPos - (0.5*length)*mom_ddrec.unit();
+	dd4hep::rec::Vector3D post   = oldPos + (0.5*length)*mom_ddrec.unit();
+	double dSToPre   = surface->distance(pre);
+	double dPreToHit = dSToPre - dSToHit;
+	global3D         = pre + (length/dPreToHit*dSToPre)*mom_ddrec.unit();
+	localPoint       = surface->globalToLocal(global3D);
+	scale = 1;
+
+	debug() << " global: (" << global3D.x()/dd4hep::mm << " " << global3D.y()/dd4hep::mm << " " << global3D.z()/dd4hep::mm
+		<< ") local: (" << localPoint.u()/dd4hep::mm << ", " << localPoint.v()/dd4hep::mm << ", 0 )"
+		<< " distance: " << surface->distance(global3D)/dd4hep::mm
+		<< " is not within boundaries." << endmsg;
+#endif
+      }
+      //else return StatusCode::SUCCESS;
+    }
+    dd4hep::rec::Vector3D globalPointSmeared;//CLHEP::Hep3Vector globalPoint(pos[0],pos[1],pos[2]);
+    dd4hep::rec::Vector2D localPointSmeared;
+
+    debug() << std::setprecision(8) << " Before smearing global: (" << pos[0] << " " << pos[1] << " " << pos[2] << ") "
+	    << "local: (" << localPoint.u()/dd4hep::mm << " " << localPoint.v()/dd4hep::mm << ")" << endmsg;
+
+    unsigned tries = 0;
+
+    bool accept_hit = false;
+
+    // Now try to smear the hit and make sure it is still on the surface
+    while (tries < 100) {
+      if (tries > 0) {
+	debug() << "retry smearing for side" << side << " layer"<< layer<< " module" << module
+		<< " sensor" << sensor << " : retries " << tries << endmsg;
+      }
+
+      double du = m_randSvc->generator(Rndm::Gauss(0, resU))->shoot();
+      double dv = m_randSvc->generator(Rndm::Gauss(0, resV))->shoot();
+      localPointSmeared.u() = localPoint.u() + du;
+      localPointSmeared.v() = localPoint.v() + dv;
+
+      dd4hep::rec::Vector3D local3DSmeared(localPointSmeared.u(), localPointSmeared.v(), 0);
+      globalPointSmeared = surface->localToGlobal(localPointSmeared);
+
+      //check if hit is in boundaries
+      if (surface->insideBounds(globalPointSmeared) && fabs(du) <= m_maxPull*resU && fabs(dv) <= m_maxPull*resV) {
+	accept_hit = true;
+	break;
+      }
+      tries++;
+    }
+
+    if (accept_hit == false) {
+      debug() << "hit could not be smeared within ladder after 100 tries: hit dropped"  << endmsg;
+      return StatusCode::SUCCESS;
+    }
+
+    // for 1D strip measurements: set v to 0! Only the measurement in u counts!
+    if(m_isStrip || (resU != 0 && resV == 0)) localPointSmeared.v() = 0. ;
+    // convert back to global position for TrackerHit
+    globalPointSmeared = surface->localToGlobal(localPointSmeared);
+
+    debug() << "  After smearing global: ("
+	    << globalPointSmeared.x()/dd4hep::mm <<" "<< globalPointSmeared.y()/dd4hep::mm <<" "<< globalPointSmeared.z()/dd4hep::mm << ") "
+	    << "local: ("
+	    << localPointSmeared.u()/dd4hep::mm << " " << localPointSmeared.v()/dd4hep::mm << ")" << endmsg;
+
+    auto outhit = hitCol->create();
+
+    outhit.setCellID(cellId);
+
+    edm4hep::Vector3d smearedPos(globalPointSmeared.x()/dd4hep::mm,
+				 globalPointSmeared.y()/dd4hep::mm,
+				 globalPointSmeared.z()/dd4hep::mm);
+    outhit.setPosition(smearedPos);
+    // recover CLHEP/Geant4 unit
+    resU /= dd4hep::mm;
+    resV /= dd4hep::mm;
+
+    std::bitset<32> type;
+    if (typeSurface.isPlane() && m_usePlanarTag) {
+      std::array<float, 6> cov;
+      cov[0] = u_direction[0];
+      cov[1] = u_direction[1];
+      cov[2] = resU*scale;
+      cov[3] = v_direction[0];
+      cov[4] = v_direction[1];
+      cov[5] = resV*scale;
+      outhit.setCovMatrix(cov);
+
+      type.set(CEPCConf::TrkHitTypeBit::PLANAR);
+    }
+    else if (typeSurface.isPlane()) {
+      outhit.setCovMatrix(CEPC::ConvertToCovXYZ(resU*scale, u_direction[0], u_direction[1], resV*scale, v_direction[0], v_direction[1]));
+    }
+    else if (typeSurface.isCylinder()) {
+      outhit.setCovMatrix(std::array<float, 6>{resU*resU*scale*scale/2., 0, resU*resU*scale*scale/2, 0, 0, resV*resV*scale*scale});
+      type.set(CEPCConf::TrkHitTypeBit::CYLINDER);
+    }
+
+    if(m_isStrip || (resU != 0 && resV == 0)){
+      type.set(CEPCConf::TrkHitTypeBit::ONE_DIMENSIONAL);
+    }
+    outhit.setType((int)type.to_ulong());
+    outhit.setEDep(e);
+    float dt = m_randSvc->generator(Rndm::Gauss(0, resT))->shoot();
+    outhit.setTime(simhit.getTime() + dt);
+    // make the relation
+    auto ass = assCol->create();
+
+    float weight = 1.0;
+
+    debug() <<" Set relation between "
+	    << " sim hit " << simhit.id()
+	    << " to tracker hit " << outhit.id()
+	    << " with a weight of " << weight
+	    << endmsg;
+
+    outhit.addToRawHits(simhit.getObjectID());
+    ass.setSim(simhit);
+    ass.setRec(outhit);
+    ass.setWeight(weight);
+
+    debug() << "-------------------------------------------------------" << endmsg;
+  }
+  else {
+    fatal() << "Not plane and cylinder: " << typeSurface << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode SmearDigiTool::finalize(){
+  StatusCode sc;
+  return sc;
+}

Digitization/DigiSimple/src/SmearDigiTool.h

+#ifndef SmearDigiTool_h
+#define SmearDigiTool_h
+
+#include "GaudiKernel/AlgTool.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+
+#include "DigiTool/IDigiTool.h"
+#include "DetInterface/IGeomSvc.h"
+
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+
+#include "DDRec/DetectorData.h"
+#include "DDRec/SurfaceManager.h"
+
+#include <vector>
+
+class SmearDigiTool : public extends<AlgTool, IDigiTool> {
+ public:
+  using extends::extends;
+  //SmearDigiTool(void* p) { m_pAlgUsing=p; };                                                                                                                         
+
+  virtual StatusCode Call(const edm4hep::SimTrackerHitCollection* simCol, edm4hep::TrackerHitCollection* hitCol,
+			  edm4hep::MCRecoTrackerAssociationCollection* assCol) override;
+  virtual StatusCode Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection* hitCol,
+			  edm4hep::MCRecoTrackerAssociationCollection* assCol) override;
+
+  StatusCode initialize() override;
+  StatusCode finalize() override;
+
+ private:
+  Gaudi::Property<std::string> m_detName{this, "DetName", ""};
+  Gaudi::Property<std::string> m_readoutName{this, "Readout", ""};
+
+  Gaudi::Property<float> m_eThreshold{this, "EnergyThreshold", 0};
+  Gaudi::Property<float> m_efficiency{this, "Efficiency", 1};
+
+  // resolution in direction of u - either one per layer or one for all layers
+  Gaudi::Property<std::vector<float> > m_resU{this, "ResolutionU", {0.004}};
+  // resolution in direction of v - either one per layer or one for all layers
+  Gaudi::Property<std::vector<float> > m_resV{this, "ResolutionV", {0.004}};
+  // resolution of time - either one per layer or one for all layers, ps as unit
+  Gaudi::Property<std::vector<float> > m_resT{this, "ResolutionT", {0.0}};
+  // whether hits are 1D strip hits
+  Gaudi::Property<bool>  m_isStrip{this, "IsStrip", false};
+  // whether use Planar tag for type and cov, if true, CEPCConf::TrkHitTypeBit::PLANAR bit is set as true
+  // cov[0]=thetaU, cov[1]=phiU, cov[2]=resU, cov[0]=thetaV, cov[1]=phiV, cov[2]=resV
+  Gaudi::Property<bool>  m_usePlanarTag{this, "UsePlanarTag", true};
+  Gaudi::Property<float> m_maxPull{this, "PullCutToRetry", 1000.};
+  Gaudi::Property<bool>  m_parameterize{this, "ParameterizeResolution", false};
+  Gaudi::Property<std::vector<float> > m_parU{this, "ParametersU", {0}};
+  Gaudi::Property<std::vector<float> > m_parV{this, "ParametersV", {0}};
+
+  SmartIF<IRndmGenSvc> m_randSvc;
+  SmartIF<IGeomSvc>    m_geosvc;
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder = nullptr;
+  const dd4hep::rec::SurfaceMap* m_surfaces;
+  //void* m_pAlgUsing = nullptr;
+};
+#endif

Digitisers/SimpleDigi/src/TPCDigiAlg.cpp → Digitization/DigiSimple/src/TPCDigiAlg.cpp

@@ -33,6 +33,7 @@
 #include <gear/BField.h>
 //
 #include "UTIL/ILDConf.h"
+#include "DetIdentifier/CEPCConf.h"
 
 #define TRKHITNCOVMATRIX 6
 
@@ -109,8 +110,8 @@ TPCDigiAlg::TPCDigiAlg(const std::string& name, ISvcLocator* svcLoc)
   //                           "Store the pointer to the SimTrackerHits in RawHits (deprecated) ",
   //                           _use_raw_hits_to_store_simhit_pointer,
   //                           bool(false));
-
-  declareProperty("PointResolutionPadPhi",_pointResoPadPhi=0.900,
+  // not used in pixel clustering mode
+  declareProperty("PointResolutionPadPhi",_pointResoPadPhi=0.9,
                   "Pad Phi Resolution constant in TPC");
   //registerProcessorParameter( "PointResolutionPadPhi" ,
   //                           "Pad Phi Resolution constant in TPC"  ,
@@ -124,21 +125,21 @@ TPCDigiAlg::TPCDigiAlg(const std::string& name, ISvcLocator* svcLoc)
   //                           _rejectCellID0 ,
   //                           (int)1) ;
 
-  declareProperty("PointResolutionRPhi",_pointResoRPhi0=0.050,
+  declareProperty("PointResolutionRPhi",_pointResoRPhi0=0.144,
                   "R-Phi Resolution constant in TPC");
   //registerProcessorParameter( "PointResolutionRPhi" ,
   //                           "R-Phi Resolution constant in TPC"  ,
   //                           _pointResoRPhi0 ,
   //                           (float)0.050) ;
 
-  declareProperty("DiffusionCoeffRPhi",_diffRPhi=0.025,
+  declareProperty("DiffusionCoeffRPhi",_diffRPhi=0.0323,
                   "R-Phi Diffusion Coefficent in TPC");
   //registerProcessorParameter( "DiffusionCoeffRPhi" ,
   //                           "R-Phi Diffusion Coefficent in TPC"  ,
   //                           _diffRPhi ,
   //                           (float)0.025) ;
-
-  declareProperty("N_eff",_nEff=22,
+  // CDR: 22, relative with gas
+  declareProperty("N_eff",_nEff=30,
                   "Number of Effective electrons per pad in TPC");
   //registerProcessorParameter( "N_eff" ,
   //                           "Number of Effective electrons per pad in TPC"  ,
@@ -152,7 +153,7 @@ TPCDigiAlg::TPCDigiAlg(const std::string& name, ISvcLocator* svcLoc)
   //                           _pointResoZ0 ,
   //                           (float)0.4) ;
 
-  declareProperty("DiffusionCoeffZ",_diffZ=0.08,
+  declareProperty("DiffusionCoeffZ",_diffZ=0.23,
                   "Z Diffusion Coefficent in TPC");
   //registerProcessorParameter( "DiffusionCoeffZ" ,
   //                           "Z Diffusion Coefficent in TPC"  ,
@@ -194,6 +195,8 @@ TPCDigiAlg::TPCDigiAlg(const std::string& name, ISvcLocator* svcLoc)
   //                           "Defines the maximum number of adjacent hits which can be merged"  ,
   //                           _maxMerge ,
   //                           (int)3) ;
+
+  declareProperty("PixelClustering", _pixelClustering=bool(true), "pixel clustering mode or pad readout mode");
 }
 
 
@@ -512,7 +515,9 @@ StatusCode TPCDigiAlg::execute()
     return StatusCode::SUCCESS;
   }
 
-  _cellid_encoder = new UTIL::BitField64( lcio::ILDCellID0::encoder_string ) ;
+  // _cellid_encoder = new UTIL::BitField64( lcio::ILDCellID0::encoder_string ) ;
+  // _cellid_encoder = new UTIL::BitField64( "system:5,side:-2,layer:13,module:6,sensor:6" ) ;
+  _cellid_encoder = new UTIL::BitField64( CEPCConf::DetEncoderString::getStringRepresentation() ) ;
   //int det_id = 0 ;
   //if ( (STHcol != nullptr) && (STHcol->size()>0) ) {
   //  auto SimTHit = STHcol->at( 0 ) ;
@@ -590,7 +595,7 @@ StatusCode TPCDigiAlg::execute()
       if(mcp.isAvailable()){
 
         // get the pt of the MCParticle, this will used later to uses nominal smearing for low momentum rubish
-        const edm4hep::Vector3f momentumMC= mcp.getMomentum() ;
+        const auto& momentumMC= mcp.getMomentum() ;
         ptSqrdMC = momentumMC[0]*momentumMC[0]+momentumMC[1]*momentumMC[1] ;
         
         debug() << " mcp id = " << mcp.id() 
@@ -761,20 +766,7 @@ StatusCode TPCDigiAlg::execute()
 
       edep = SimTHit.getEDep();
 
-      // Calculate Point Resolutions according to Ron's Formula
-
-      // sigma_{RPhi}^2 = sigma_0^2 + Cd^2/N_{eff} * L_{drift}
-
-      // sigma_0^2 = (50micron)^2 + (900micron*sin(phi))^2
-      // Cd^2/N_{eff}} = 25^2/(22/sin(theta)*h/6mm)
-      // Cd = 25 ( microns / cm^(1/2) )
-      // (this is for B=4T, h is the pad height = pad-row pitch in mm,
-      // theta is the polar angle)
-
-      // sigma_{z}^2 = (400microns)^2 + L_{drift}cm * (80micron/sqrt(cm))^2
-
-      double aReso =_pointResoRPhi0*_pointResoRPhi0 + (_pointResoPadPhi*_pointResoPadPhi * sin(padPhi)*sin(padPhi)) ;
-      double driftLength = gearTPC.getMaxDriftLength() - (fabs(thisPoint.z()));
+      double driftLength = gearTPC.getMaxDriftLength() - (fabs(thisPoint.z())); // mm unit
 
       if (driftLength <0) {
         debug() << " TPCDigiAlg : Warning! driftLength < 0 " << driftLength << " --> Check out your GEAR file!!!!" << endmsg;
@@ -782,16 +774,69 @@ StatusCode TPCDigiAlg::execute()
         debug() << "gearTPC.getMaxDriftLength() = " << gearTPC.getMaxDriftLength() << endmsg;
         driftLength = 0.10;
       }
+      driftLength /= 10; // to cm unit
 
       padheight = padLayout.getPadHeight(padLayout.getNearestPad(thisPoint.perp(),thisPoint.phi()));
+      double ne = _nEff * (padheight/6.0);
+
+      double tpcRPhiRes, tpcZRes;
+      if (_pixelClustering) {
+        // Calculate Point Resolutions according to Chang Yue and ZHAO Guang's work
+
+        // sigma_{RPhi}^2 = (sigma_0^2 + Cd^2 * L_{drift})/N_{eff}
+
+        // sigma_0 = 0.5mm/sqrt(12)
+        // N_{eff} = 30*h/6mm
+        // Cd = 32.3 ( microns / cm^(1/2) )
+        // (this is for B=3T, h is the pad height = pad-row pitch in mm)
+
+        // sigma_{z}^2 = (400microns)^2 + L_{drift}cm * (80micron/sqrt(cm))^2
+
+        // OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO
+        // Pixel readout TPC spatial resolution (preliminary) update, Xin She, 20250228 
+        // + Take "hodoscope" effect into consideration
+        // + The parameters are based on Neff=30 and the pad pitch=0.5mm 
+        //   @T2K gas, B=3T, D_{T} is close to 32.3 um/sqrt(cm) @230 V/cm Drift field
+        // + sigma_x = 0.006001*sqrt(z/10.)+0.1175*exp(-0.09018*z/10.), z>=5.1mm
+        // + sigma_x = 0.1443-0.0047*z, z<5.1mm 
+        //   where sigma_x is the r-phi resolution, unit [mm], 
+        //   z is the driftlength, unit [mm]
+        // OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO
+        double driftLength_in_mm =  driftLength*10.; 
+        if(driftLength_in_mm>=5.1) // using fitted curves
+        {
+            tpcRPhiRes = _fittedRPhiResoParas[0]*std::sqrt(0.1*driftLength_in_mm) + 
+                         _fittedRPhiResoParas[1]*std::exp(-1*_fittedRPhiResoParas[2]*driftLength_in_mm);
+        }
+        else // using linear interpolation
+        {
+            tpcRPhiRes = _fittedRPhiResoParas[3] + _fittedRPhiResoParas[4]*driftLength_in_mm;
+        }
+        // without "hodoscope" effect
+        double aReso = _pointResoRPhi0*_pointResoRPhi0 ;
+        double bReso = _diffRPhi * _diffRPhi ;
+        //tpcRPhiRes   = sqrt(aReso + bReso * driftLength) / sqrt(ne); // driftLength in cm
+        //tpcZRes      = sqrt( _pointResoZ0 * _pointResoZ0 + _diffZ * _diffZ * driftLength ) / sqrt(ne); // driftLength in cm
+        tpcZRes      = sqrt( _pointResoZ0 * _pointResoZ0 + _diffZ * _diffZ * driftLength ) / sqrt(_nEff); // driftLength in cm
 
-      double bReso = ( (_diffRPhi * _diffRPhi) / _nEff ) * sin(padTheta) * ( 6.0 / (padheight) )  * ( 4.0 / bField  ) ;
-
-      double tpcRPhiRes = sqrt( aReso + bReso * (driftLength / 10.0) ); // driftLength in cm
-
-      double tpcZRes  = sqrt(( _pointResoZ0 * _pointResoZ0 )
-                             +
-                             ( _diffZ * _diffZ ) * (driftLength / 10.0) ); // driftLength in cm
+      }
+      else {
+        // Calculate Point Resolutions according to Ron's Formula
+
+        // sigma_{RPhi}^2 = sigma_0^2 + Cd^2/N_{eff} * L_{drift}
+
+        // sigma_0^2 = (50micron)^2 + (900micron*sin(phi))^2
+        // Cd^2/N_{eff}} = 25^2/(22/sin(theta)*h/6mm)
+        // Cd = 25 ( microns / cm^(1/2) )
+        // (this is for B=4T, h is the pad height = pad-row pitch in mm,
+        // theta is the polar angle)
+
+        // sigma_{z}^2 = (400microns)^2 + L_{drift}cm * (80micron/sqrt(cm))^2
+        double aReso = _pointResoRPhi0*_pointResoRPhi0 + (_pointResoPadPhi*_pointResoPadPhi * sin(padPhi)*sin(padPhi)) ;
+        double bReso = ( (_diffRPhi * _diffRPhi) / ne ) * sin(padTheta) * ( 4.0 / bField  ) ;
+        tpcRPhiRes   = sqrt( aReso + bReso * driftLength ); // driftLength in cm
+        tpcZRes      = sqrt( _pointResoZ0 * _pointResoZ0 + _diffZ * _diffZ * driftLength ); // driftLength in cm
+      }
 
       int padIndex = padLayout.getNearestPad(thisPoint.perp(),thisPoint.phi());
 
@@ -829,7 +874,8 @@ StatusCode TPCDigiAlg::execute()
 
       // create a tpc voxel hit and store it for this row
       Voxel_tpc * atpcVoxel = new Voxel_tpc(iRowHit,iPhiHit,iZHit, thisPoint, edep, tpcRPhiRes, tpcZRes);
-
+      debug() << "to seed hit position: " << atpcVoxel->getX() << "," << atpcVoxel->getY() << "," << atpcVoxel->getZ()
+              << " iRow=" << iRowHit << " iPhi=" << iPhiHit << " padIndex=" << padIndex << endmsg;
       _tpcRowHits.at(iRowHit).push_back(atpcVoxel);
       ++numberOfVoxelsCreated;
 
@@ -986,7 +1032,8 @@ StatusCode TPCDigiAlg::execute()
             if( momentum_set ){
 
               const edm4hep::Vector3f Momentum1 = Hit1.getMomentum() ;
-              const edm4hep::Vector3f Momentum2 = Hit1.getMomentum() ;
+              //fucd: Hit1.getMomentum() -> Hit2.getMomentum()
+              const edm4hep::Vector3f Momentum2 = Hit2.getMomentum() ;
 
               CLHEP::Hep3Vector mom1(Momentum1[0],Momentum1[1],Momentum1[2]);
               CLHEP::Hep3Vector mom2(Momentum2[0],Momentum2[1],Momentum2[2]);
@@ -1057,7 +1104,7 @@ StatusCode TPCDigiAlg::execute()
         auto mcp = (_tpcHitMap[ seed_hit ]).getMCParticle() ;
         if(mcp.isAvailable()) {
           ++_NLostPhysicsTPCHits;
-          const edm4hep::Vector3f mom= mcp.getMomentum() ;
+          const auto& mom= mcp.getMomentum() ;
           double ptSQRD = mom[0]*mom[0]+mom[1]*mom[1] ;
           if( ptSQRD > (0.2*0.2) ) ++_NLostPhysicsAbove02GeVPtTPCHits ;
           if( ptSQRD > 1.0 )  ++_NLostPhysicsAbove1GeVPtTPCHits ;
@@ -1164,7 +1211,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
   Voxel_tpc* seed_hit  = aVoxel;
 
   //  if( seed_hit->getRowIndex() > 5 ) return ;
-  debug() << "==============" << endmsg;
+
   //store hit variables
   edm4hep::MutableTrackerHit trkHit;// = _trkhitVec->create();
   //now the hit pos has to be smeared
@@ -1184,7 +1231,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
 
   // make sure the hit is not smeared beyond the TPC Max DriftLength
   if( fabs(point.z()) > gearTPC.getMaxDriftLength() ) point.setZ( (fabs(point.z()) / point.z() ) * gearTPC.getMaxDriftLength() );
-  debug() << "==============" << endmsg;
+  debug() << seed_hit->getX() << "," << seed_hit->getY() << "," << seed_hit->getZ() << " -> " << point << endmsg;
   edm4hep::Vector3d pos(point.x(),point.y(),point.z());
   trkHit.setPosition(pos);
   trkHit.setEDep(seed_hit->getEDep());
@@ -1235,7 +1282,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
     << "\n" ;
     throw errorMsg.str();
   }
-  debug() << "==============" << endmsg;
+
   // For no error in R
   std::array<float,TRKHITNCOVMATRIX> covMat={sin(unsmearedPhi)*sin(unsmearedPhi)*tpcRPhiRes*tpcRPhiRes,
     -cos(unsmearedPhi)*sin(unsmearedPhi)*tpcRPhiRes*tpcRPhiRes,
@@ -1253,7 +1300,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
     << "\n" ;
     throw errorMsg.str();
   }
-  debug() << "==============" << endmsg;
+
   if(pos[0]*pos[0]+pos[1]*pos[1]>0.0){
     //    push back the SimTHit for this TrackerHit
 
@@ -1270,7 +1317,7 @@ void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
     _NRecTPCHits++;
   }
 
-  debug() << "==============" << endmsg;
+
 #ifdef DIGIPLOTS
 //  edm4hep::SimTrackerHit* theSimHit = _tpcHitMap[seed_hit];
 //  double rSimSqrd = theSimHit->getPosition()[0]*theSimHit->getPosition()[0] + theSimHit->getPosition()[1]*theSimHit->getPosition()[1];
@@ -1312,7 +1359,7 @@ void TPCDigiAlg::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge){
 
   unsigned number_of_hits_to_merge = hitsToMerge->size();
 
-
+  debug() << "number_of_hits_to_merge = " << number_of_hits_to_merge << endmsg;
   for(unsigned int ihitCluster = 0; ihitCluster < number_of_hits_to_merge; ++ihitCluster){
 
     sumZ += hitsToMerge->at(ihitCluster)->getZ();
@@ -1324,6 +1371,7 @@ void TPCDigiAlg::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge){
     if (_use_raw_hits_to_store_simhit_pointer) {
       trkHit.addToRawHits(_tpcHitMap[hitsToMerge->at(ihitCluster)].getObjectID());
     }
+    debug() << "raw hit: " << _tpcHitMap[hitsToMerge->at(ihitCluster)].getPosition() << endmsg;
 
     auto rel = _relCol->create();
     rel.setRec (trkHit);
@@ -1363,7 +1411,7 @@ void TPCDigiAlg::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge){
   if( fabs(point.z()) > gearTPC.getMaxDriftLength() ) point.setZ( (fabs(point.z()) / point.z() ) * gearTPC.getMaxDriftLength() );
 
   double pos[3] = {point.x(),point.y(),point.z()};
-
+  debug() << "to hit: " << pos[0] << "," << pos[1] << "," << pos[2] << endmsg;
   //---------------------------------------------------------------------------------
   trkHit.setPosition(pos);
   trkHit.setEDep(sumEDep);

Digitisers/SimpleDigi/src/TPCDigiAlg.h → Digitization/DigiSimple/src/TPCDigiAlg.h

@@ -169,6 +169,24 @@ protected:
   edm4hep::TrackerHitCollection* _trkhitVec;
   edm4hep::MCRecoTrackerAssociationCollection* _relCol;
 
+  // OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO
+  // Pixel readout TPC spatial resolution (preliminary) update, Xin She, 20250228 
+  // + Take "hodoscope" effect into consideration
+  // + The parameters are based on Neff=30 and the pad pitch=0.5mm 
+  //   @T2K gas, B=3T, D_{T} is close to 32.3 um/sqrt(cm) @230 V/cm Drift field
+  // + sigma_x = 0.006001*sqrt(z/10.)+0.1175*exp(-0.09018*z/10.), z>=5.1mm
+  // + sigma_x = 0.1443-0.0047*z, z<5.1mm 
+  //   where sigma_x is the r-phi resolution, unit [mm], 
+  //   z is the driftlength, unit [mm]
+  // Input fitted parameters will be set as optional parameters though Gaudi algo 
+  // interface. If some conditions (e.g. pixel size) had changed, it can be modified
+  // OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO OOOoooOOOoooOOO
+  Gaudi::Property<std::vector<double>> _fittedRPhiResoParas{this, "fittedRPhiResoParas",{0.006001,0.1175,0.009018,0.1443,-0.0047}};
+
+ 
+
+  bool _pixelClustering;
+
   bool _use_raw_hits_to_store_simhit_pointer;
 
   int _rejectCellID0;

Digitization/DigiSimple/src/TPCPixelDigiAlg.cpp

+#include "TPCPixelDigiAlg.h"
+
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/Vector3D.h"
+#include "DDRec/ISurface.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+
+DECLARE_COMPONENT(TPCPixelDigiAlg)
+
+TPCPixelDigiAlg::TPCPixelDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc) {
+  // Input collections
+  declareProperty("SimTrackHitCollection", m_inputColHdls, "Handle of the Input SimTrackerHit collection");
+
+
+  // Output collections
+  declareProperty("TPCTrackerHits", m_outputColHdls, "Handle of the output TrackerHit collection");
+  declareProperty("TPCTrackerHitAss", m_assColHdls, "Handle of the Association collection between SimTrackerHit and TrackerHit");
+}
+
+StatusCode TPCPixelDigiAlg::initialize() {
+  m_digiTool = ToolHandle<IDigiTool>(m_digiToolName.value());
+
+  info() << "DigiTool " << m_digiTool.typeAndName() << " found" << endmsg;
+
+  info() << "TPCPixelDigiAlg::initialized" << endmsg;
+  return GaudiAlgorithm::initialize();
+}
+
+
+StatusCode TPCPixelDigiAlg::execute(){
+  auto trkCol = m_outputColHdls.createAndPut();
+  auto assCol = m_assColHdls.createAndPut();
+
+  const edm4hep::SimTrackerHitCollection* simCol = nullptr;
+  try {
+    simCol = m_inputColHdls.get();
+  }
+  catch(GaudiException &e){
+    debug() << "Collection " << m_inputColHdls.fullKey() << " test is unavailable in event " << m_nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  if (simCol->size() == 0) return StatusCode::SUCCESS;
+  debug() << m_inputColHdls.fullKey() << " has cluser "<< simCol->size() << endmsg;
+
+  m_digiTool->Call(simCol, trkCol, assCol);
+
+  debug() << "Created " << trkCol->size() << " hits, "<<simCol->size()<<" clusters"
+          << simCol->size() - trkCol->size() << " hits got dismissed for being out of boundary"
+          << endmsg;
+
+  m_nEvt++;
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiAlg::finalize(){
+  info() << "Processed " << m_nEvt << " events " << endmsg;
+  return GaudiAlgorithm::finalize();
+}