add DCHDigi and DedxSvc

Digitisers/DCHDigi/CMakeLists.txt

+gaudi_subdir(DCHDigi v0r0)
+
+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(CLHEP REQUIRED)
+find_package(podio REQUIRED )
+
+set(srcs
+    src/*.cpp
+)
+
+gaudi_depends_on_subdirs(
+    Detector/DetInterface
+)
+## Modules
+gaudi_add_module(DCHDigi ${srcs}
+    INCLUDE_DIRS FWCore GaudiKernel GaudiAlgLib CLHEP DD4hep 
+    LINK_LIBRARIES FWCore GaudiKernel GaudiAlgLib CLHEP DD4hep ${DD4hep_COMPONENT_LIBRARIES} DDRec
+    -Wl,--no-as-needed 
+    EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+)

Digitisers/DCHDigi/src/DCHDigiAlg.cpp

+/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
+#include "DCHDigiAlg.h"
+
+
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+
+
+#include <math.h>
+#include <cmath>
+#include <algorithm>
+
+DECLARE_COMPONENT( DCHDigiAlg )
+
+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("CaloAssociationCollection", w_AssociationCol, "Handle of Association collection");
+   
+}
+
+StatusCode DCHDigiAlg::initialize()
+{
+  /*
+  m_geosvc = service<IGeoSvc>("GeoSvc");
+  if ( !m_geosvc )  throw "DCHDigiAlg :Failed to find GeoSvc ...";
+  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);
+  */
+  std::cout<<"DCHDigiAlg::initialized"<< std::endl;
+  return GaudiAlgorithm::initialize();
+}
+
+StatusCode DCHDigiAlg::execute()
+{
+  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();
+      
+      if ( id_hits_map.find(id) != id_hits_map.end()) id_hits_map[id].push_back(SimHit);
+      else 
+      {
+          std::vector<edm4hep::SimTrackerHit> vhit;
+          vhit.push_back(SimHit);
+          id_hits_map[id] = vhit ;
+      }
+  }
+
+  for(std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> >::iterator iter = id_hits_map.begin(); iter != id_hits_map.end(); iter++)
+  {
+    auto trkHit = Vec->create();
+    trkHit.setCellID(iter->first);
+    double tot_edep   = 0 ;
+    double tot_length = 0 ;
+    double tot_time = 0 ;
+    double tot_x = 0 ;
+    double tot_y = 0 ;
+    double tot_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
+        tot_length += iter->second.at(i).getPathLength();//mm
+        tot_time += iter->second.at(i).getTime();
+        tot_x    += iter->second.at(i).getPosition()[0];
+        tot_y    += iter->second.at(i).getPosition()[1];
+        tot_z    += iter->second.at(i).getPosition()[2];
+ 
+        auto asso = AssoVec->create();
+        asso.setRec(trkHit);
+        asso.setSim(iter->second.at(i));
+        asso.setWeight(1.0/simhit_size);
+    }
+    
+    trkHit.setTime(tot_time/simhit_size);
+    trkHit.setEDep(tot_edep/simhit_size);
+    trkHit.setEdx (tot_edep*1000/(tot_length/10) ); // MeV/cm
+    trkHit.setPosition (edm4hep::Vector3d(tot_x/simhit_size, tot_y/simhit_size, tot_z/simhit_size));
+  }
+  std::cout<<"output digi DCHhit size="<< Vec->size() <<std::endl;
+  _nEvt ++ ;
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode DCHDigiAlg::finalize()
+{
+  info() << "Processed " << _nEvt << " events " << endmsg;
+  return GaudiAlgorithm::finalize();
+}

Digitisers/DCHDigi/src/DCHDigiAlg.h

+#ifndef DCH_DIGI_ALG_H
+#define DCH_DIGI_ALG_H
+
+#include "FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+
+#include <DDRec/DetectorData.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include "DetInterface/IGeoSvc.h"
+
+
+
+
+class DCHDigiAlg : public GaudiAlgorithm
+{
+ 
+public:
+ 
+  DCHDigiAlg(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() ;
+ 
+protected:
+
+  SmartIF<IGeoSvc> m_geosvc;
+  typedef std::vector<float> FloatVec;
+  int _nEvt ;
+
+  //float m_length;
+  dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+ 
+  //Gaudi::Property<float> m_scale     { this, "Scale", 1 };
+  //Gaudi::Property<float> m_resolution{ this, "Res", 0.01 };
+
+  // Input collections
+  DataHandle<edm4hep::SimTrackerHitCollection> r_SimDCHCol{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};
+  // Output collections
+  DataHandle<edm4hep::TrackerHitCollection>    w_DigiDCHCol{"DigiDCHitsCollection", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection>    w_AssociationCol{"DCHitAssociationCollection", Gaudi::DataHandle::Writer, this};
+};
+
+#endif

Examples/options/tut_detsim_dedx.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 = "CepC_v4-onlyVXD.xml"
+
+if not os.getenv("DETCEPCV4ROOT"):
+    print("Can't find the geometry. Please setup envvar DETCEPCV4ROOT." )
+    sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCEPCV4ROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+    print("Can't find the compact geometry file: %s"%geometry_path)
+    sys.exit(-1)
+geometry_path = "/junofs/users/wxfang/MyGit/tmp/sdt/CEPCSW/Detector/DetDriftChamber/compact/det.xml" 
+
+from Configurables import GeoSvc
+geosvc = GeoSvc("GeoSvc")
+geosvc.compact = geometry_path
+
+###### Dedx Svc ##########################################################
+from Configurables import DedxSvc
+dedxsvc = DedxSvc("DedxSvc")
+dedxsvc.scale = 1
+dedxsvc.resolution = 0.01
+dedxsvc.material_Z = 26
+dedxsvc.material_A = 53.9396
+dedxsvc.material_density = 7.874 #g/cm^3
+##############################################################################
+# 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 = ["proton"]
+gun.EnergyMins = [0.1] # GeV
+gun.EnergyMaxs = [100] # GeV
+
+gun.ThetaMins = [90] # rad; 45deg
+gun.ThetaMaxs = [90] # rad; 45deg
+
+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")
+
+# 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")
+
+
+##############################################################################
+# POD I/O
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "test_detsim10.root"
+out.outputCommands = ["keep *"]
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [genalg, detsimalg, out],
+                EvtSel = 'NONE',
+                EvtMax = 10,
+                ExtSvc = [rndmengine, dsvc, geosvc],
+)

Service/DedxSvc/CMakeLists.txt

+gaudi_subdir(DedxSvc v0r0)
+
+set(DedxSvc_srcs
+    src/*.cpp
+)
+find_package(Geant4 REQUIRED ui_all vis_all)
+include(${Geant4_USE_FILE})
+gaudi_install_headers(DedxSvc)
+
+gaudi_add_module(DedxSvc ${DedxSvc_srcs}
+    INCLUDE_DIRS GaudiKernel
+    LINK_LIBRARIES GaudiKernel
+)

Service/DedxSvc/DedxSvc/IDedxSvc.h

+#ifndef I_Dedx_SVC_H
+#define I_Dedx_SVC_H
+
+#include "GaudiKernel/IService.h"
+
+#include "G4Step.hh"
+
+class IDedxSvc: virtual public IService {
+public:
+    DeclareInterfaceID(IDedxSvc, 0, 1); // major/minor version
+    
+    virtual ~IDedxSvc() = default;
+    virtual float pred(const G4Step* aStep)=0 ;
+
+};
+
+
+#endif

Service/DedxSvc/src/DedxSvc.cpp

+#include "DedxSvc.h"
+
+//https://folk.uib.no/ruv004/
+DECLARE_COMPONENT(DedxSvc)
+
+DedxSvc::DedxSvc(const std::string& name, ISvcLocator* svc)
+    : base_class(name, svc)
+{
+}
+
+DedxSvc::~DedxSvc()
+{
+}
+
+float DedxSvc::pred(const G4Step* aStep)
+{
+    G4Track* gTrack = aStep->GetTrack() ;
+    G4int z = gTrack->GetDefinition()->GetPDGCharge();
+    if (z == 0) return 0;
+    G4double M = gTrack->GetDefinition()->GetPDGMass();//MeV
+    M = pow(10,6)*M; //to eV
+    G4double gammabeta=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma();
+    if(gammabeta<0.01)return 0;//too low momentum
+    float beta = gammabeta/sqrt(1.0+pow(gammabeta,2));
+    float gamma = gammabeta/beta;
+    float Tmax = 2*m_me*pow(gammabeta,2)/(1+(2*gamma*m_me/M)+pow(m_me/M,2));
+    float dedx = m_K*pow(z,2)*m_material_Z*(0.5*log(2*m_me*pow(gammabeta,2)*Tmax/pow(m_I,2))-pow(beta,2))/(m_material_A*pow(beta,2));    
+    dedx = dedx*m_scale;// the material density can be absorbed in scale 
+    dedx = dedx*(1+((*m_distribution)(m_generator)));
+    return dedx*m_material_density; // MeV / cm 
+}
+
+StatusCode DedxSvc::initialize()
+{
+    m_distribution = new std::normal_distribution<double>(0, m_resolution);
+    m_me = 0.511*pow(10,6);//0.511 MeV to eV
+    m_K = 0.307075;//const
+    m_I = m_material_Z*10;  // Approximate
+    return StatusCode::SUCCESS;
+}
+
+StatusCode DedxSvc::finalize()
+{
+    return StatusCode::SUCCESS;
+}

Service/DedxSvc/src/DedxSvc.h

+#ifndef Dedx_SVC_H
+#define Dedx_SVC_H
+
+#include "DedxSvc/IDedxSvc.h"
+#include <GaudiKernel/Service.h>
+#include "G4Step.hh"
+#include <random>
+
+class DedxSvc : public extends<Service, IDedxSvc>
+{
+    public:
+        DedxSvc(const std::string& name, ISvcLocator* svc);
+        ~DedxSvc();
+
+
+        StatusCode initialize() override;
+        StatusCode finalize() override;
+        float pred(const G4Step* aStep) override;
+
+    private:
+
+        Gaudi::Property<float> m_material_Z{this, "material_Z", 2};//Default is Helium
+        Gaudi::Property<float> m_material_A{this, "material_A", 4};
+        Gaudi::Property<float> m_material_density{this, "material_density", 0.000178};//g/cm^3
+        Gaudi::Property<float> m_scale{this, "scale", 1};
+        Gaudi::Property<float> m_resolution{this, "resolution", 0.01};
+        float m_me;// Here me is the electron rest mass
+        float m_K; // K was set as a constant.
+        float m_I; // Mean excitation energy
+        std::default_random_engine m_generator;
+        std::normal_distribution<double>* m_distribution;
+};
+
+#endif