diff --git a/Examples/options/tut_detsim_SDT.py b/Examples/options/tut_detsim_SDT.py index ca5498b772ef7afdf10146dc2a551cd4c751a277..8d490193c7ab07e3cd9df9bef9875bfbcc3866cc 100644 --- a/Examples/options/tut_detsim_SDT.py +++ b/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 diff --git a/Examples/options/tut_detsim_dedx.py b/Examples/options/tut_detsim_dedx.py deleted file mode 100644 index fa49a00e87193a03bef378137157881db4893541..0000000000000000000000000000000000000000 --- a/Examples/options/tut_detsim_dedx.py +++ /dev/null @@ -1,146 +0,0 @@ -#!/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], -) diff --git a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp index 1983ba72f44eca929a75ac4abfd7a01d65fe3a3e..f8dab312c8eb867116cb537732976803a47a3575 100644 --- a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp +++ b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp @@ -1,6 +1,7 @@ #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 diff --git a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h index 6e80ad7972d1fcac831f709cde20f77f8bf97f5a..dfb671e3d9ddb692bcc4ea5eb63f2d5a9f06266d 100644 --- a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.h +++ b/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