Skip to content
Snippets Groups Projects
Commit e185ef0a authored by lintao@ihep.ac.cn's avatar lintao@ihep.ac.cn
Browse files

WIP: get Z/A of current material from Geant4.

parent af161167
No related branches found
No related tags found
No related merge requests found
#include "BetheBlochEquationDedxSimTool.h" #include "BetheBlochEquationDedxSimTool.h"
#include "G4Step.hh" #include "G4Step.hh"
#include "G4SystemOfUnits.hh"
// https://folk.uib.no/ruv004/ // https://folk.uib.no/ruv004/
DECLARE_COMPONENT(BetheBlochEquationDedxSimTool) DECLARE_COMPONENT(BetheBlochEquationDedxSimTool)
...@@ -10,6 +11,14 @@ double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep) ...@@ -10,6 +11,14 @@ double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep)
G4Track* gTrack = aStep->GetTrack() ; G4Track* gTrack = aStep->GetTrack() ;
G4int z = gTrack->GetDefinition()->GetPDGCharge(); G4int z = gTrack->GetDefinition()->GetPDGCharge();
if (z == 0) return 0; 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 G4double M = gTrack->GetDefinition()->GetPDGMass();//MeV
M = pow(10,6)*M; //to eV M = pow(10,6)*M; //to eV
G4double gammabeta=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma(); G4double gammabeta=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma();
...@@ -17,10 +26,10 @@ double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep) ...@@ -17,10 +26,10 @@ double BetheBlochEquationDedxSimTool::dedx(const G4Step* aStep)
float beta = gammabeta/sqrt(1.0+pow(gammabeta,2)); float beta = gammabeta/sqrt(1.0+pow(gammabeta,2));
float gamma = gammabeta/beta; float gamma = gammabeta/beta;
float Tmax = 2*m_me*pow(gammabeta,2)/(1+(2*gamma*m_me/M)+pow(m_me/M,2)); 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*m_scale;// the material density can be absorbed in scale
dedx = dedx*(1+((*m_distribution)(m_generator))); dedx = dedx*(1+((*m_distribution)(m_generator)));
return dedx*m_material_density; // MeV / cm return dedx*material_density; // MeV / cm
} }
StatusCode BetheBlochEquationDedxSimTool::initialize() StatusCode BetheBlochEquationDedxSimTool::initialize()
...@@ -28,12 +37,11 @@ StatusCode BetheBlochEquationDedxSimTool::initialize() ...@@ -28,12 +37,11 @@ StatusCode BetheBlochEquationDedxSimTool::initialize()
m_distribution = new std::normal_distribution<double>(0, m_resolution); m_distribution = new std::normal_distribution<double>(0, m_resolution);
m_me = 0.511*pow(10,6);//0.511 MeV to eV m_me = 0.511*pow(10,6);//0.511 MeV to eV
m_K = 0.307075;//const m_K = 0.307075;//const
m_I = m_material_Z*10; // Approximate
info() << "Initialize BetheBlochEquationDedxSimTool with following parameters" << endmsg; info() << "Initialize BetheBlochEquationDedxSimTool with following parameters" << endmsg;
info() << "-> m_me: " << m_me << endmsg; info() << "-> m_me: " << m_me << endmsg;
info() << "-> m_K: " << m_K << endmsg; info() << "-> m_K: " << m_K << endmsg;
info() << "-> m_I: " << m_I << endmsg;
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment