LET calculation for incident neutrons

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).

  1. What should I do differently with regard to the physics to obtain non-zero results?

  2. 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?

  3. 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!

Neutrons don’t do dE/dx, and they don’t interact electromagnetically. They interact via elastic scattering off of nuclei, capture by a nucleus, or quasi-elastic scattering where they leave the target nucleus in an excited state.

Thank you for the clarification! In this case, how could I obtain LET values from neutron-induced secondary ions, as in the fig below?

Reference paper: C. Peng, Z. Lei, Z. Zhang, Y. En and Y. Huang, “Investigating Neutron-Induced Single Event Transient Characteristics by TCAD Simulations in 65 nm Technology and Below,” 2019 3rd International Conference on Radiation Effects of Electronic Devices (ICREED) , 2019, pp. 1-4, doi: 10.1109/ICREED49760.2019.9205173.