From e185ef0a930ca2482a2282652cb4215984d976ec Mon Sep 17 00:00:00 2001
From: lintao <lintao51@gmail.com>
Date: Fri, 18 Sep 2020 18:05:10 +0800
Subject: [PATCH] WIP: get Z/A of current material from Geant4.
---
.../src/BetheBlochEquationDedxSimTool.cpp | 16 ++++++++++++----
1 file changed, 12 insertions(+), 4 deletions(-)
diff --git a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
index 5b64de69..6aecca2b 100644
--- a/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
+++ b/Simulation/DetSimDedx/src/BetheBlochEquationDedxSimTool.cpp
@@ -1,5 +1,6 @@
#include "BetheBlochEquationDedxSimTool.h"
#include "G4Step.hh"
+#include "G4SystemOfUnits.hh"
// https://folk.uib.no/ruv004/
DECLARE_COMPONENT(BetheBlochEquationDedxSimTool)
@@ -10,6 +11,14 @@ double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep)
G4Track* gTrack = aStep->GetTrack() ;
G4int z = gTrack->GetDefinition()->GetPDGCharge();
if (z == 0) return 0;
+
+ G4Material* material = gTrack->GetMaterial();
+ G4double material_density = material->GetDensity() / (CLHEP::g/CLHEP::cm3); // conert from G4 unit.
+ G4double material_Z = material->GetZ();
+ G4double material_A = material->GetA();
+
+ m_I = material_Z*10; // Approximate
+
G4double M = gTrack->GetDefinition()->GetPDGMass();//MeV
M = pow(10,6)*M; //to eV
G4double gammabeta=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma();
@@ -17,10 +26,10 @@ double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep)
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));
+ 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*m_material_density; // MeV / cm
+ return dedx*material_density; // MeV / cm
}
StatusCode BetheBlochEquationDedxSimTool::initialize()
@@ -28,12 +37,11 @@ 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
- m_I = m_material_Z*10; // Approximate
+
info() << "Initialize BetheBlochEquationDedxSimTool with following parameters" << endmsg;
info() << "-> m_me: " << m_me << endmsg;
info() << "-> m_K: " << m_K << endmsg;
- info() << "-> m_I: " << m_I << endmsg;
return StatusCode::SUCCESS;
}
--
GitLab