Digitization/DigiMuon/src/MuonDigiAlg.h

@@ -49,26 +49,26 @@ class MuonDigiAlg : public GaudiAlgorithm
   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 _simhit, dd4hep::DDSegmentation::BitFieldCoder* _m_decoder);
-  double Gethit_sipm_length_Barrel();
-  double Gethit_sipm_length_Endcap();
-  void EdeptoADC();
-  void SaveData_mapcell();
-  bool Mapcell_todata(int _message[6], std::array<unsigned long long, 2> _key);
+  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 Save_trkhit(edm4hep::TrackerHitCollection* _trkhitVec, std::array<unsigned long long, 2> _key, int _pdgid, unsigned long long _cellid, edm4hep::Vector3d _pos);
-  void Renew_strip(std::array<unsigned long long, 2> _key1, std::array<unsigned long long, 2> _key2);
-  void Find_anotherlayer(int _i, edm4hep::Vector3d _pos1, edm4hep::Vector3d _pos2, double _ddposi, int _pdgid);
-  int Get_true_pdgid(int _i, int _pdgid);
-  void Save_onelayer_signal(edm4hep::TrackerHitCollection* _trkhitVec);
-  void Save_pos(int _i);
+  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:
-
 
 
+ protected:
 
   TTree* m_tree;
   TFile* m_wfile;
@@ -81,6 +81,7 @@ class MuonDigiAlg : public GaudiAlgorithm
 
   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;
@@ -130,18 +131,11 @@ class MuonDigiAlg : public GaudiAlgorithm
   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
-  std::array<unsigned long long, 2> key, key1, key2;
-  int pdgid, layer, slayer, strip, Fe, Env, anotherlayer_cell_num, true_pdgid;
-  unsigned long long cellid, abspdgid, cellid1, cellid2; 
-  double Edep, xydist, hit_sipm_length, ADC, hit_strip_length, ADCmean, Hit_max, Hit_min;
-  edm4hep::Vector3d pos, pos1, pos2;
-  edm4hep::Vector3d mcppos;
-  dd4hep::Position ddpos, ddpos1, ddpos2;  
-  int all_message1[6], all_message2[6];  // int layer1, slayer1, strip1, Fe1, Env1, pdgid1;
   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, 

Digitization/DigiSimple/src/SmearDigiTool.cpp

@@ -89,16 +89,16 @@ StatusCode SmearDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHi
   }
 
   auto e = simhit.getEDep();
-  if (e <= m_eThreshold) return StatusCode::SUCCESS;
+  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, "system");
-  int side    = m_decoder->get(cellId, "side"  );
-  int layer   = m_decoder->get(cellId, "layer" );
-  int module  = m_decoder->get(cellId, "module");
-  int sensor  = m_decoder->get(cellId, "sensor");
+  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();
@@ -277,6 +277,7 @@ StatusCode SmearDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHi
 	    << localPointSmeared.u()/dd4hep::mm << " " << localPointSmeared.v()/dd4hep::mm << ")" << endmsg;
 
     auto outhit = hitCol->create();
+
     outhit.setCellID(cellId);
 
     edm4hep::Vector3d smearedPos(globalPointSmeared.x()/dd4hep::mm,

Digitization/DigiSimple/src/TPCDigiAlg.cpp

@@ -791,10 +791,34 @@ StatusCode TPCDigiAlg::execute()
         // (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
+        //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
+
       }
       else {
         // Calculate Point Resolutions according to Ron's Formula

Digitization/DigiSimple/src/TPCDigiAlg.h

@@ -169,6 +169,22 @@ 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;

README.md

@@ -7,10 +7,7 @@ CEPC offline software prototype based on [Key4hep](https://github.com/key4hep).
 
 ## Quick start
 
-SSH to lxlogin (Alma Linux 9) and start the container CentOS 7:
-```
-$ /cvmfs/container.ihep.ac.cn/bin/hep_container shell CentOS7
-```
+SSH to lxlogin (Alma Linux 9).
 
 Before run following commands, please make sure you setup the CVMFS:
 
@@ -20,6 +17,7 @@ $ cd CEPCSW
 $ git checkout master # branch name
 $ source setup.sh
 $ ./build.sh
+$ source setup.sh
 $ ./run.sh Examples/options/helloalg.py
 ```
 
@@ -37,8 +35,3 @@ $ ./run.sh Examples/options/helloalg.py
 
 * Reconstruction: Reconstruction
 
-
-## CyberPFA-5.0.1-dev (developing)
-* Based on CEPCSW tag tdr 24.12.0
-
-

Reconstruction/ParticleID/src/FinalPIDAlg.cpp

@@ -19,7 +19,7 @@ FinalPIDAlg::FinalPIDAlg( const std::string& name, ISvcLocator* pSvcLocator )
   // output
   declareProperty("OutputPFOName", m_outPFOCol, "Reconstructed particles with PID information");
   // PID method
-  declareProperty("Method", m_method = "TPC+TOF", "PID method: TPC, TPC+TOF, TPC+TOF+CALO"); 
+  declareProperty("PIDMethod", m_method = "TPC+TOF+CALO", "PID method: TPC, TPC+TOF, TPC+TOF+CALO"); 
 }
 
 //------------------------------------------------------------------------------
@@ -86,13 +86,18 @@ StatusCode FinalPIDAlg::execute(){
     auto outpfo = pfo.clone();
     outpfocol->push_back(outpfo);
     if(m_method.value().find("TPC+TOF")!=std::string::npos ){
-      if( _hasTPC && _hasTPC && (tofcol->size()+dqdxcol->size()>0) )
-        FillTPCTOFPID(dqdxcol, tofcol, outpfo);
+      std::array<double, 5> chi2s; chi2s.fill(0);
+      if( _hasTPC && dqdxcol->size()>0){
+        FillTPCPID(dqdxcol, outpfo, chi2s);
+      }
+      if( _hasTOF && tofcol->size()>0){
+        FillTOFPID(tofcol, outpfo, chi2s);
+      }
     }
     if(m_method.value().find("CALO")!=std::string::npos ) 
       FillCaloPID(outpfo);
 
-    debug()<<" cyber pfo pid : " << outpfo.getType() << " cyber pfo charge : " << outpfo.getCharge() << " cyber pfo energy "<< outpfo.getEnergy() << endmsg;
+    debug()<<" cyber pfo pid : " << outpfo.getType() << " cyber pfo charge : " << outpfo.getCharge() << " cyber pfo energy "<< outpfo.getEnergy() << " cyber pfo mass "<< outpfo.getMass() << endmsg;
   }
 
   _nEvt++;
@@ -107,33 +112,10 @@ StatusCode FinalPIDAlg::finalize(){
 }
 
 
-/*
-void FinalPIDAlg::FillTPCPID(const edm4hep::ReconstructedParticleCollection* pfocol, const edm4hep::RecDqdxCollection* dqdxcol, edm4hep::ParticleIDCollection* pidcol){
-  for (auto pfo : *pfocol){
-    for (auto trk : pfo.getTracks()){
-      for (auto dqdx : *dqdxcol){
-        if (dqdx.getTrack() == trk){
-          std::array<double, 5> chi2s = dqdx.getHypotheses();
-          int minchi2idx = std::distance(chi2s.begin(), std::min_element(chi2s.begin(), chi2s.end()));
-          auto pid = pidcol->create();
-          pid.setPDG( pfo.getCharge() * PDGIDs.at(minchi2idx) );
-          pid.setLikelihood( chi2s[minchi2idx] );
-          //pfo.setParticleID(pid);
-          
-        }
-      }
-    }
-  }
-}
-*/
-
-
-StatusCode FinalPIDAlg::FillTPCTOFPID(const edm4hep::RecDqdxCollection* dqdxcol, const edm4hep::RecTofCollection* tofcol, edm4hep::MutableReconstructedParticle& pfo){
+void FinalPIDAlg::FillTPCPID(const edm4hep::RecDqdxCollection* dqdxcol, edm4hep::MutableReconstructedParticle& pfo, std::array<double, 5>& chi2s){
 
     for (auto trk : pfo.getTracks()){
 
-      std::array<double, 5> chi2s; chi2s.fill(0);
-      
       for (auto dqdx : *dqdxcol){
         if (dqdx.getTrack() == trk){
           for (int i = 0; i < 5; i++){
@@ -142,6 +124,25 @@ StatusCode FinalPIDAlg::FillTPCTOFPID(const edm4hep::RecDqdxCollection* dqdxcol,
         }
       }
 
+      int pdgid = 0;
+      if ( std::all_of(chi2s.begin(), chi2s.end(), [](double x){return x == 0;}) ){
+        pdgid = pfo.getCharge() * 211;
+      }else{
+        int minchi2idx = std::distance(chi2s.begin(), std::min_element(chi2s.begin(), chi2s.end()));
+        pdgid = pfo.getCharge() * PDGIDs.at(minchi2idx);
+      }
+      pfo.setType( pdgid );
+      pfo.setMass( ParticleMass.at( abs(pdgid) ) );
+      pfo.setEnergy( sqrt(pfo.getMomentum()[0]*pfo.getMomentum()[0] + pfo.getMomentum()[1]*pfo.getMomentum()[1] + pfo.getMomentum()[2]*pfo.getMomentum()[2] + pfo.getMass()*pfo.getMass()) );
+      debug() << " fill tpc pid:  chi2s : " << chi2s[0]<<" " <<chi2s[1]<<" " <<chi2s[2]<<" "<<chi2s[3]<<" "<<chi2s[4]<<"	id:"<<pfo.getType()<<" 	mass:"<<pfo.getMass()<<endmsg;
+    }
+    
+}
+
+void FinalPIDAlg::FillTOFPID(const edm4hep::RecTofCollection* tofcol, edm4hep::MutableReconstructedParticle& pfo, std::array<double, 5>& chi2s){
+
+    for (auto trk : pfo.getTracks()){
+
       for (auto tof : *tofcol){
         if (tof.getTrack() == trk){
           double toft = tof.getTime();
@@ -155,13 +156,21 @@ StatusCode FinalPIDAlg::FillTPCTOFPID(const edm4hep::RecDqdxCollection* dqdxcol,
         }
       }
 
-      int minchi2idx = std::distance(chi2s.begin(), std::min_element(chi2s.begin(), chi2s.end()));
-      int pdgid = pfo.getCharge() * PDGIDs.at(minchi2idx);
+      int pdgid = 0;
+      if ( std::all_of(chi2s.begin(), chi2s.end(), [](double x){return x == 0;}) ){
+        pdgid = pfo.getCharge() * 211;
+      }else{
+        int minchi2idx = std::distance(chi2s.begin(), std::min_element(chi2s.begin(), chi2s.end()));
+        pdgid = pfo.getCharge() * PDGIDs.at(minchi2idx);
+      }
+
       pfo.setType( pdgid );
       pfo.setMass( ParticleMass.at( abs(pdgid) ) );
       pfo.setEnergy( sqrt(pfo.getMomentum()[0]*pfo.getMomentum()[0] + pfo.getMomentum()[1]*pfo.getMomentum()[1] + pfo.getMomentum()[2]*pfo.getMomentum()[2] + pfo.getMass()*pfo.getMass()) );
+      debug() << " fill tof pid: chi2 : " << chi2s[0]<<" " <<chi2s[1]<<" " <<chi2s[2]<<" "<<chi2s[3]<<" "<<chi2s[4]<<"	id:"<<pfo.getType()<<" 	mass:"<<pfo.getMass()<<endmsg;
+
     }
-    return StatusCode::SUCCESS;
+
 }
 
 
@@ -197,13 +206,15 @@ StatusCode FinalPIDAlg::FillCaloPID(edm4hep::MutableReconstructedParticle& pfo){
             int pdgid = 22;
             pfo.setType( pdgid );
             pfo.setMass( ParticleMass.at( abs(pdgid) ) );
-            pfo.setEnergy( sqrt(pfo.getMomentum()[0]*pfo.getMomentum()[0] + pfo.getMomentum()[1]*pfo.getMomentum()[1] + pfo.getMomentum()[2]*pfo.getMomentum()[2] + pfo.getMass()*pfo.getMass()) );
+            double p_scale = sqrt( pfo.getEnergy()*pfo.getEnergy() - pfo.getMass()*pfo.getMass() ) / sqrt(pfo.getMomentum()[0]*pfo.getMomentum()[0] + pfo.getMomentum()[1]*pfo.getMomentum()[1] + pfo.getMomentum()[2]*pfo.getMomentum()[2] );
+            pfo.setMomentum( Vector3f(pfo.getMomentum()[0]*p_scale, pfo.getMomentum()[1]*p_scale, pfo.getMomentum()[2]*p_scale) );
         }
         else{
             int pdgid = 130;
             pfo.setType( pdgid );
             pfo.setMass( ParticleMass.at( abs(pdgid) ) );
-            pfo.setEnergy( sqrt(pfo.getMomentum()[0]*pfo.getMomentum()[0] + pfo.getMomentum()[1]*pfo.getMomentum()[1] + pfo.getMomentum()[2]*pfo.getMomentum()[2] + pfo.getMass()*pfo.getMass()) );
+            double p_scale = sqrt( pfo.getEnergy()*pfo.getEnergy() - pfo.getMass()*pfo.getMass() ) / sqrt(pfo.getMomentum()[0]*pfo.getMomentum()[0] + pfo.getMomentum()[1]*pfo.getMomentum()[1] + pfo.getMomentum()[2]*pfo.getMomentum()[2] );
+            pfo.setMomentum( Vector3f(pfo.getMomentum()[0]*p_scale, pfo.getMomentum()[1]*p_scale, pfo.getMomentum()[2]*p_scale) );
 
         }
     }