update DCHDigi to get correct wire position

Digitisers/DCHDigi/CMakeLists.txt

@@ -3,11 +3,10 @@ gaudi_subdir(DCHDigi v0r0)
 find_package(CLHEP REQUIRED;CONFIG)
 find_package(DD4hep COMPONENTS DDG4 REQUIRED)
 find_package(EDM4HEP REQUIRED )
-message("EDM4HEP_INCLUDE_DIRS: ${EDM4HEP_INCLUDE_DIR}")
-message("EDM4HEP_LIB: ${EDM4HEP_LIBRARIES}")
 include_directories(${EDM4HEP_INCLUDE_DIR})
 
 find_package(podio REQUIRED )
+find_package(ROOT COMPONENTS MathCore Physics GenVector Geom REQUIRED)
 
 set(srcs
     src/*.cpp
@@ -15,11 +14,12 @@ set(srcs
 
 gaudi_depends_on_subdirs(
     Detector/DetInterface
+    Detector/DetSegmentation
 )
 ## Modules
 gaudi_add_module(DCHDigi ${srcs}
-    INCLUDE_DIRS FWCore GaudiKernel GaudiAlgLib ${CLHEP_INCLUDE_DIR} DD4hep 
-    LINK_LIBRARIES FWCore GaudiKernel GaudiAlgLib ${CLHEP_LIBRARIES} DD4hep ${DD4hep_COMPONENT_LIBRARIES} DDRec
+    INCLUDE_DIRS FWCore GaudiKernel GaudiAlgLib ${CLHEP_INCLUDE_DIR} DD4hep ROOT 
+    LINK_LIBRARIES FWCore GaudiKernel GaudiAlgLib ${CLHEP_LIBRARIES} DD4hep ${DD4hep_COMPONENT_LIBRARIES} DDRec ROOT DetSegmentation
     -Wl,--no-as-needed 
     EDM4HEP::edm4hep EDM4HEP::edm4hepDict
 )

Digitisers/DCHDigi/src/DCHDigiAlg.cpp

@@ -29,39 +29,83 @@ DCHDigiAlg::DCHDigiAlg(const std::string& name, ISvcLocator* svcLoc)
   // Output collections
   declareProperty("DigiDCHitCollection", w_DigiDCHCol, "Handle of Digi DCHit collection");
   
-  declareProperty("CaloAssociationCollection", w_AssociationCol, "Handle of Association collection");
+  declareProperty("AssociationCollection", w_AssociationCol, "Handle of Association collection");
    
 }
 
 StatusCode DCHDigiAlg::initialize()
 {
-  /*
   m_geosvc = service<IGeomSvc>("GeomSvc");
   if ( !m_geosvc )  throw "DCHDigiAlg :Failed to find GeomSvc ...";
   dd4hep::Detector* m_dd4hep = m_geosvc->lcdd();
   if ( !m_dd4hep )  throw "DCHDigiAlg :Failed to get dd4hep::Detector ...";
-  m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep);
-  */
+  dd4hep::Readout readout = m_dd4hep->readout(m_readout_name);
+  m_segmentation = dynamic_cast<dd4hep::DDSegmentation::GridDriftChamber*>(readout.segmentation().segmentation());
+  m_decoder = m_geosvc->getDecoder(m_readout_name);
+  if (!m_decoder) {
+      error() << "Failed to get the decoder. " << endmsg;
+      return StatusCode::FAILURE;
+  }
+
+  std::string s_output=m_Output; 
+  if(m_WriteAna){
+      m_fout = new TFile(s_output.c_str(),"RECREATE"); 
+      m_tree = new TTree("evt","tree");
+      m_tree->Branch("chamber"  , &m_chamber  );
+      m_tree->Branch("layer"    , &m_layer    );
+      m_tree->Branch("cell"     , &m_cell     );
+      m_tree->Branch("cell_x"   , &m_cell_x   );
+      m_tree->Branch("cell_y"   , &m_cell_y   );
+      m_tree->Branch("simhit_x" , &m_simhit_x );
+      m_tree->Branch("simhit_y" , &m_simhit_y );
+      m_tree->Branch("simhit_z" , &m_simhit_z );
+      m_tree->Branch("hit_x"    , &m_hit_x    );
+      m_tree->Branch("hit_y"    , &m_hit_y    );
+      m_tree->Branch("hit_z"    , &m_hit_z    );
+      m_tree->Branch("dca"      , &m_dca      );
+      m_tree->Branch("hit_dE"   , &m_hit_dE   );
+      m_tree->Branch("hit_dE_dx", &m_hit_dE_dx);
+  }
   std::cout<<"DCHDigiAlg::initialized"<< std::endl;
   return GaudiAlgorithm::initialize();
 }
 
+void DCHDigiAlg::Reset()
+{
+
+  std::vector<int  >().swap( m_chamber   );
+  std::vector<int  >().swap( m_layer     );
+  std::vector<int  >().swap( m_cell      );
+  std::vector<float>().swap( m_cell_x    );
+  std::vector<float>().swap( m_cell_y    );
+  std::vector<float>().swap( m_simhit_x  );
+  std::vector<float>().swap( m_simhit_y  );
+  std::vector<float>().swap( m_simhit_z  );
+  std::vector<float>().swap( m_hit_x     );
+  std::vector<float>().swap( m_hit_y     );
+  std::vector<float>().swap( m_hit_z     );
+  std::vector<float>().swap( m_dca       );
+  std::vector<float>().swap( m_hit_dE    );
+  std::vector<float>().swap( m_hit_dE_dx );
+
+}
+
 StatusCode DCHDigiAlg::execute()
 {
+  info() << "Processing " << _nEvt << " events " << endmsg;
+  Reset();
   std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> > id_hits_map;
   edm4hep::TrackerHitCollection* Vec   = w_DigiDCHCol.createAndPut();
   edm4hep::MCRecoTrackerAssociationCollection* AssoVec   = w_AssociationCol.createAndPut();
   const edm4hep::SimTrackerHitCollection* SimHitCol =  r_SimDCHCol.get();
-  if(SimHitCol == 0) 
-  {
-     std::cout<<"not found SimCalorimeterHitCollection"<< std::endl;
-     return StatusCode::SUCCESS;
-  }
   std::cout<<"input sim hit size="<< SimHitCol->size() <<std::endl;
   for( int i = 0; i < SimHitCol->size(); i++ ) 
   {
       edm4hep::SimTrackerHit SimHit = SimHitCol->at(i);
       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 ( id_hits_map.find(id) != id_hits_map.end()) id_hits_map[id].push_back(SimHit);
       else 
@@ -74,11 +118,11 @@ StatusCode DCHDigiAlg::execute()
 
   for(std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> >::iterator iter = id_hits_map.begin(); iter != id_hits_map.end(); iter++)
   {
+    unsigned long long wcellid = iter->first;
     auto trkHit = Vec->create();
-    trkHit.setCellID(iter->first);
+    trkHit.setCellID(wcellid);
     double tot_edep   = 0 ;
     double tot_length = 0 ;
-    double tot_time = 0 ;
     double pos_x = 0 ;
     double pos_y = 0 ;
     double pos_z = 0 ;
@@ -87,50 +131,81 @@ StatusCode DCHDigiAlg::execute()
     {
         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);
+    Wstart = 10*Wstart;// from DD4HEP cm to mm
+    Wend   = 10*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 ;
     for(unsigned int i=0; i< simhit_size; i++)
     {
-        dd4hep::rec::Vector3D  west(0,0,0);
-        dd4hep::rec::Vector3D  east(0,0,0);
-        dd4hep::rec::Vector3D  pos(iter->second.at(i).getPosition()[0], iter->second.at(i).getPosition()[1], iter->second.at(i).getPosition()[2]);
-        dd4hep::rec::Vector3D  numerator = (east-west).cross(west-pos) ;
-        dd4hep::rec::Vector3D  denominator = (east-west) ;
-        float distance = numerator.r()/denominator.r() ;
-        std::cout<<"distance="<<distance<<std::endl;
-        if(distance < min_distance){
-            min_distance = distance;
+        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], iter->second.at(i).getPosition()[1], iter->second.at(i).getPosition()[2]);
+        TVector3  numerator = denominator.Cross(Wstart-pos) ;
+        float tmp_distance = numerator.Mag()/denominator.Mag() ;
+        //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
-        /*
-        tot_x    += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[0];
-        tot_y    += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[1];
-        tot_z    += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[2];
-        */
  
         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){
+            m_simhit_x.push_back(pos.x());
+            m_simhit_y.push_back(pos.y());
+            m_simhit_z.push_back(pos.z());
+        }
     }
     
     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*1000/(tot_length/10) ); // MeV/cm, need check!
-    //trkHit.setPosition (edm4hep::Vector3d(tot_x/tot_edep, tot_y/tot_edep, tot_z/tot_edep));//center mass
+    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){
+        m_chamber.push_back(chamber);
+        m_layer  .push_back(layer  );
+        m_cell      .push_back(cellID);
+        m_cell_x    .push_back(Wstart.x());
+        m_cell_y    .push_back(Wstart.y());
+        m_hit_x     .push_back(pos_x);
+        m_hit_y     .push_back(pos_y);
+        m_hit_z     .push_back(pos_z);
+        m_dca       .push_back(min_distance);
+        m_hit_dE    .push_back(trkHit.getEDep());
+        m_hit_dE_dx .push_back(trkHit.getEdx() );
+    }
   }
   std::cout<<"output digi DCHhit size="<< Vec->size() <<std::endl;
   _nEvt ++ ;
 
+  if(m_WriteAna) m_tree->Fill();
+
   return StatusCode::SUCCESS;
 }
 
 StatusCode DCHDigiAlg::finalize()
 {
   info() << "Processed " << _nEvt << " events " << endmsg;
+  if(m_WriteAna){
+      m_fout->cd();
+      m_tree->Write();
+      m_fout->Close();
+  }
   return GaudiAlgorithm::finalize();
 }

Digitisers/DCHDigi/src/DCHDigiAlg.h

@@ -10,7 +10,13 @@
 #include <DDRec/DetectorData.h>
 #include <DDRec/CellIDPositionConverter.h>
 #include "DetInterface/IGeomSvc.h"
+#include "DDSegmentation/Segmentation.h"
+#include "DetSegmentation/GridDriftChamber.h"
 
+#include "TVector3.h"
+#include "TROOT.h"
+#include "TTree.h"
+#include "TFile.h"
 
 
 
@@ -33,6 +39,7 @@ public:
   /** Called after data processing for clean up.
    */
   virtual StatusCode finalize() ;
+  void Reset();
  
 protected:
 
@@ -40,13 +47,39 @@ protected:
   typedef std::vector<float> FloatVec;
   int _nEvt ;
 
-  //float m_length;
+  TFile* m_fout;
+  TTree* m_tree;
+  std::vector<int  > m_chamber   ;
+  std::vector<int  > m_layer     ;
+  std::vector<int  > m_cell      ;
+  std::vector<float> m_cell_x    ;
+  std::vector<float> m_cell_y    ;
+  std::vector<float> m_simhit_x  ;
+  std::vector<float> m_simhit_y  ;
+  std::vector<float> m_simhit_z  ;
+  std::vector<float> m_hit_x     ;
+  std::vector<float> m_hit_y     ;
+  std::vector<float> m_hit_z     ;
+  std::vector<float> m_dca       ;
+  std::vector<float> m_hit_dE    ;
+  std::vector<float> m_hit_dE_dx ;
+
+
+
   dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+  dd4hep::DDSegmentation::GridDriftChamber* m_segmentation;
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+  
+  Gaudi::Property<std::string> m_readout_name{ this, "readout", "DriftChamberHitsCollection"};//readout for getting segmentation
  
   Gaudi::Property<float> m_res_x     { this, "res_x", 0.11};//mm
   Gaudi::Property<float> m_res_y     { this, "res_y", 0.11};//mm
   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<bool>  m_WriteAna { this, "WriteAna", false};
+
+  Gaudi::Property<std::string> m_Output { this, "output", "ana_DCH_digi.root"};
 
   // Input collections
   DataHandle<edm4hep::SimTrackerHitCollection> r_SimDCHCol{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};

Examples/options/tut_detsim_digi_SDT.py

+#!/usr/bin/env python
+
+import os
+print(os.environ["DD4HEP_LIBRARY_PATH"])
+import sys
+# sys.exit(0)
+
+from Gaudi.Configuration import *
+
+##############################################################################
+# Random Number Svc
+##############################################################################
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Geant4
+rndmengine.SetSingleton = True
+rndmengine.Seeds = [42]
+
+# rndmgensvc = RndmGenSvc("RndmGenSvc")
+# rndmgensvc.Engine = rndmengine.name()
+
+
+##############################################################################
+# Event Data Svc
+##############################################################################
+from Configurables import K4DataSvc
+dsvc = K4DataSvc("EventDataSvc")
+
+
+##############################################################################
+# Geometry Svc
+##############################################################################
+
+# geometry_option = "CepC_v4-onlyTracker.xml"
+geometry_option = "det.xml"
+
+if not os.getenv("DETDRIFTCHAMBERROOT"):
+    print("Can't find the geometry. Please setup envvar DETCEPCV4ROOT." )
+    sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETDRIFTCHAMBERROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+    print("Can't find the compact geometry file: %s"%geometry_path)
+    sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+##############################################################################
+# Physics Generator
+##############################################################################
+from Configurables import GenAlgo
+from Configurables import GtGunTool
+from Configurables import StdHepRdr
+from Configurables import SLCIORdr
+from Configurables import HepMCRdr
+from Configurables import GenPrinter
+
+gun = GtGunTool("GtGunTool")
+# gun.Particles = ["pi+"]
+# gun.EnergyMins = [100.] # GeV
+# gun.EnergyMaxs = [100.] # GeV
+
+gun.Particles = ["e-"]
+
+# gun.PositionXs = [100.] # mm
+# gun.PositionYs = [100.] # mm
+# gun.PositionZs = [0.] # mm
+
+
+gun.EnergyMins = [10.] # GeV
+gun.EnergyMaxs = [10.] # GeV
+
+gun.ThetaMins = [90] # rad; 45deg
+gun.ThetaMaxs = [90] # rad; 45deg
+
+gun.PhiMins = [90] # rad; 0deg
+gun.PhiMaxs = [90] # rad; 360deg
+#gun.PhiMins = [0] # rad; 0deg
+#gun.PhiMaxs = [0] # rad; 360deg
+
+# stdheprdr = StdHepRdr("StdHepRdr")
+# stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
+
+# lciordr = SLCIORdr("SLCIORdr")
+# lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio"
+
+# hepmcrdr = HepMCRdr("HepMCRdr")
+# hepmcrdr.Input = "example_UsingIterators.txt"
+
+genprinter = GenPrinter("GenPrinter")
+
+genalg = GenAlgo("GenAlgo")
+genalg.GenTools = ["GtGunTool"]
+# genalg.GenTools = ["StdHepRdr"]
+# genalg.GenTools = ["StdHepRdr", "GenPrinter"]
+# genalg.GenTools = ["SLCIORdr", "GenPrinter"]
+# genalg.GenTools = ["HepMCRdr", "GenPrinter"]
+
+##############################################################################
+# Detector Simulation
+##############################################################################
+from Configurables import DetSimSvc
+
+detsimsvc = DetSimSvc("DetSimSvc")
+
+# from Configurables import ExampleAnaElemTool
+# example_anatool = ExampleAnaElemTool("ExampleAnaElemTool")
+
+from Configurables import DetSimAlg
+
+detsimalg = DetSimAlg("DetSimAlg")
+
+if int(os.environ.get("VIS", 0)):
+    detsimalg.VisMacs = ["vis.mac"]
+
+detsimalg.RunCmds = [
+#    "/tracking/verbose 1",
+]
+detsimalg.AnaElems = [
+    # example_anatool.name()
+    # "ExampleAnaElemTool"
+    "Edm4hepWriterAnaElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+
+from Configurables import AnExampleDetElemTool
+example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
+
+from Configurables import CalorimeterSensDetTool
+from Configurables import DriftChamberSensDetTool
+
+calo_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
+driftchamber_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+
+# dedxoption = "DummyDedxSimTool"
+dedxoption = "BetheBlochEquationDedxSimTool"
+
+driftchamber_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+    dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+    dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+    dedx_simtool.material_Z = 2
+    dedx_simtool.material_A = 4
+    dedx_simtool.scale = 10
+    dedx_simtool.resolution = 0.0001
+
+##############################################################################
+from Configurables import DCHDigiAlg
+dCHDigiAlg = DCHDigiAlg("DCHDigiAlg")
+dCHDigiAlg.readout = "DriftChamberHitsCollection"
+dCHDigiAlg.drift_velocity = 40#um/ns
+dCHDigiAlg.mom_threshold = 0 #GeV
+dCHDigiAlg.SimDCHitCollection = "DriftChamberHitsCollection"
+dCHDigiAlg.DigiDCHitCollection = "DigiDCHitsCollection"
+dCHDigiAlg.AssociationCollection = "DCHAssociationCollectio"
+dCHDigiAlg.WriteAna  = True
+dCHDigiAlg.output  = "ana_DCH_digi_v1.root"
+
+##############################################################################
+# POD I/O
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "detsim_digi_DCH.root"
+out.outputCommands = ["keep *"]
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+
+from Configurables import ApplicationMgr
+#ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, out],
+ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg],
+                EvtSel = 'NONE',
+                EvtMax = 10,
+                ExtSvc = [rndmengine, dsvc, geosvc],
+)