Merge pull request #47 from wenxingfang/master

update BetheBlochEquationDedxSimTool

Examples/options/tut_detsim_SDT.py

@@ -139,6 +139,10 @@ 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
 
 ##############################################################################
 # POD I/O

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],
-)

Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp

 #include "BetheBlochEquationDedxSimTool.h"
 #include "G4Step.hh"
 #include "G4SystemOfUnits.hh"
+#include "CLHEP/Random/RandGauss.h"
 
 // https://folk.uib.no/ruv004/
 DECLARE_COMPONENT(BetheBlochEquationDedxSimTool)
@@ -19,22 +20,20 @@ double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep)
 
     m_I = material_Z*10;  // Approximate
 
-    G4double M = gTrack->GetDefinition()->GetPDGMass();//MeV
-    M = pow(10,6)*M; //to eV
+    G4double M = gTrack->GetDefinition()->GetPDGMass()/(CLHEP::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)*material_Z*(0.5*log(2*m_me*pow(gammabeta,2)*Tmax/pow(m_I,2))-pow(beta,2))/(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*material_density; // MeV / cm 
+    dedx = dedx*CLHEP::RandGauss::shoot(m_scale, m_resolution);
+    double Dedx = dedx * (CLHEP::MeV/CLHEP::cm) / (CLHEP::g/CLHEP::cm3) ;
+    return Dedx*material_density;
 }
 
 StatusCode BetheBlochEquationDedxSimTool::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
 

Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h

@@ -3,7 +3,6 @@
 
 #include "DetSimInterface/IDedxSimTool.h"
 #include <GaudiKernel/AlgTool.h>
-#include <random>
 
 class BetheBlochEquationDedxSimTool: public extends<AlgTool, IDedxSimTool> {
     public:
@@ -23,9 +22,6 @@ class BetheBlochEquationDedxSimTool: public extends<AlgTool, IDedxSimTool> {
         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