'stp_lib' Introduction

stopping power or dE/dx is to describe how a particle loses its energy when it is traveling in a medium. For example, it can be used to estimate how much energy a particle will deposit when it passes through your detector. In nuclear physics, the unit of dE/dx is $\frac{MeV}{mg/cm^{2}}$, since the beam normally in unit of $MeV$, and the target thickness is in unit of $mg/cm^{2}$.

There are several programs can simulate dE/dx for all kinds of particles in a medium. For example, SRIM ( official website) and LISE++ ( official website). However, sometimes, I need to generate dE/dx data for internal use for my program, and reading stopping power from external files is not genuinely convenient. So I search the library that can help me generate dE/dx file, but I don't have any good luck to find it.

I come across that VIKAR program ( Virtual Instrumentation for Kinematics and Reactions, by Dr. S.D. Pain ) has its dE/dx generator, and I have the source code. Not only I would like to study how VIKAR's dE/dx data compares to those from LISE++, but also I would like to extract VIKAR's dE/dx generator. The VIKAR is written in Fortran, and some parts of code is in Fortran 77. The main part for dE/dx calculation is the subroutine ncdedx. I rewrite and simplify the ncdedx subroutine, and translate the code into a C++ class, which provides easy-to-use methods to get the dE/dx and stopping range data.

Run the code:

I prepare two versions. One can be used in ROOT (i.e. using TMath library), and one is just pure c++ for your convenient. You can download the source from my GitHub repository (link) or here

In the following, I will just use the ROOT version to do a quick demo. The case is for $^{86}Kr$ particle to travel in a $CH_{2}$ medium. The class constructor will need 6 pieces of information; three from the target (CH2), and three from the beam. The CH2 target is composed of one carton and two hydrogen. Hence, the average $A = 14/3$, average $Z = 8/3$.

The class name is "stp_class".
to calculate the dEdx data, we invoke process() method.
to get the dEdx, we invoke get_dEdx() method.
to get the range, we invoke get_range() method.
to set the beam energy, we invoke set_beam_energy( eng )

    // target = CH2
    // C ==> A = 12, Z = 6
    // H ==> A =  1, Z = 1
    // A_avg = ( 12*1 + 1*2) / 3
    // Z_avg = (  6*1 + 1*2) / 3
    // beam = 86Kr, A = 86, Z = 36
    double target_density = 0.93; // [g/cm^3]
    double target_A_avg = 14./3;
    double target_Z_avg = 8./3;
    double beam_A = 86;
    double beam_Z = 36;
    double beam_eng = 90; // [MeV]
    stp_class*  demo_obj;
    demo_obj = new stp_class( target_density, 
                              beam_eng );

    TGraph* gr = new TGraph();
    int pointN = 0; 
    for( double eng = 1; eng < 5000; eng+=1) {
        demo_obj->set_beam_energy( eng );
        double dEdx = demo_obj->get_dEdx();
        gr->SetPoint( pointN, eng, dEdx );  pointN+=1; }

compare the results.

The following are the results of stp_lib which is derived from the ncdedx subroutine of VIKAR program. The curves labeled as "He_base", "H_base", "ATIMA-LS", and "ATIMA-no LS" are from the LISE++ calculations.

The reference:
[He-base] F.Hubert et al, AD&ND Tables 46 (1990) 1
[H -base] J.F.Ziegler et al, Pergamon Press, NY (low energy)
ATIMA 1.2 LS-theory (recommended for high energy)
ATIMA 1.2 without LS-correction

It is clear the dE/dx from VIKAR rises and reach its max too fast ( see the energy under 50 MeV), and it also drops too fast for the higher energy part ( see the energy above 3000 MeV). The good agreement comes at the range between around 1000 to 2000 MeV. If you would like to have precise dE/dx data, I would recommend to use the results from LISE++.