Hello,
I am looking to obtain the LET statistics for incident neutrons (2 MeV, 2000 events) on a Si cube (1 um) I am working on the modified extended Hadronic example (Hadr04) such that the physics list includes NeutronHP physics, Neutron transport and thermal scattering. I get the following result for LET (elaborate output at output.txt):
unrestricted dE/dx : | 0 eV/cm | 0 eV/cm | 0 eV/cm | 0 eV/cm | 0 eV/cm |
(MeV/g/cm2) : | 0 eVcm2/g| 0 eVcm2/g| 0 eVcm2/g| 0 eVcm2/g| |
I obtain the Stopping Power (and hence, LET) as follows. Note that LET is obtained from extended Electromagnetics example TestEm0. The processes in consideration are are fHadronic (from G4ProcessType.hh).
-
What should I do differently with regard to the physics to obtain non-zero results?
-
The below code for LET is implemented in EndofRun fuction of Run.cc. Is this the right approach to obtain the stopping power for incident neutrons?
-
I get a segmentation fault when using fElectromagnetic instead of fHadronic as process type. Why is the Em process no longer applicable here?
Code snipped for Stopping Power/LET calculation in Run.cc:
G4EmCalculator emCal; //instantiate EmCalculator
G4ProcessVector *plist = particle->GetProcessManager()->GetProcessList();
G4String pName;
G4double cut;
std::vector<G4String> emName;
std::vector<G4double> enerCut;
size_t length = plist->size();
for (size_t j = 0; j < length; j++)
{
pName = (*plist)[j]->GetProcessName();
if (((*plist)[j]->GetProcessType() == fHadronic) )
{
emName.push_back(pName);
}
// define material
// const G4Material *material = fDetector->GetMaterial();
std::vector<G4double> dedx1;
std::vector<G4double> dedx2;
G4double dedx, dedxtot = 0.;
size_t nproc = emName.size();
for (size_t j = 0; j < nproc; j++)
{
dedx = emCal.ComputeDEDX(energy, particle, emName[j], material, energy);
dedxtot += dedx;
dedx1.push_back(dedx);
dedx2.push_back(dedx / density);
dedx1[j] = dedx;
dedx2[j] = dedx / density;
}
dedx1[nproc] = dedxtot;
dedx2[nproc] = dedxtot / density;
// print stopping power
G4cout << "\n \n unrestricted dE/dx : ";
for (size_t j = 0; j <= nproc; j++)
{
G4cout << "\t" << std::setw(13) << G4BestUnit(dedx1[j], "Energy/Length");
}
G4cout << "\n (MeV/g/cm2) : ";
for (size_t j = 0; j <= nproc; j++)
{
G4cout << "\t" << std::setw(13) << G4BestUnit(dedx2[j], "Energy*Surface/Mass");
}
I also tried to implement the above code for LET in RunAction.cc, in accordance to the extended Electromagnetic example TestEm0, but get the same result. Any guidance would be most appreciated!