Inelastic scattering length for user defined source

Dear all,

Is it possible to score the inelastic scattering length of the beam particles with a source defined energy spectrum? (I couldn’t find the name of the variable in order to extract it by using mgdraw.)
In the .out file I can find this quantity only for the incident particles at beam energy, but what is the value of the energy in this case? Is it the maximum energy defined with the BEAM card?

Best regards,
Mihaela

Dear Mihaela,

Let’s start with the latter question. The inelastic scattering lengths reported in the .out file are calculated with the energy defined in the WHAT(1) of the BEAM card.

However, you can obtain the inelastic scattering length at practically any energy below that value. To do so, you can call the SGTINL function directy in the user defined source.f routine.

The function definition is:

SGTINL ( IJ, EKE, ELKE, MMAT, LENOLD )

It requires several variables:

  1. IJ = particle type (paprop numeration)
  2. EKE = particle kinetic energy (GeV)
  3. ELKE = log (eke)
  4. MMAT = material number of the medium
  5. LENOLD = if .true. we are at the energy of the previous call. For your purposes, it must be always false

The output of this function is the inelastic macroscopic cross section in cm^-1. To get the scattering length you just need to divide 1 by this number.
Notice that the LENOLD logical variable needs to be reinitialized to false each time you perform this calculation, since the function changes its value to true each time a new energy tabulation index is computed.

Let me know if you have any further doubt!
Cheers,
Daniele

Dear Daniele,

Thank you for your response.

This works even better for me since I am willing to extract the cross sections in the end.

Unfortunately, I don’t understand how to use the SGTINL function. When I call it in the source_newgen.f routine and I try to compile it I receive 5 errors regarding the variables of the function:
image

In which file is this function defined? I saw that the only place where I can find its name is in sgtbcm.inc.
I tried to find information regarding this in the manual or on the forum but I couldn’t find anything.

Any further help would be really helpful.

Mihaela

This function is defined in the source code and it is, therefore, not accessible.
It seems that you haven’t defined the variables that you are providing to the function. Additionally, you are using a CALL instruction, which is applicable only to routines, while this is a function.

Would you mind sharing your source.f user routine?

You can find the source file attached.
source_newgen_SGTINL.f (7.6 KB)
.
Also, I have another question: What do you mean by IJ = particle type (paprop numeration)? Is it the same as IJBEAM?
If I replace IJ with -6 (alpha particle) I get this error:

Thank you,
Mihaela

Thank you for your file. I have attached a working copy in the end. A few remarks:

  1. You have declared MMAT as a double precision, while the material index is actually an integer. This discrepancy is the source of the error that you attached.
  2. MMAT has to be the material number that you need. Fluka has 25 pre-defined materials, you put 40. This is perfectly legitimate. However, for me it is impossible to guess which material you need without the inputfile. In the example I attached I just considered carbon as target material (i.e. MMAT = 6)
  3. SGTINL is a function already defined in the code. It should be declared as

double precision, external :: SGTINL

  1. To invoke the function, one should explicitly provide the required variables. Let’s call SIGIN the inelastic scattering length, in this case:

SIGIN = SGTINL( particle_code, EKE, ELKE, MMAT, LENOLD )

I attach here the example I run with a dummy geometry. The computed values of the inelastic cross sections seem reasonable to be.
Let me know if it works!
inScattLen.inp (1022 Bytes)
source_newgen_SGTINL.f (7.9 KB)

1 Like

Dear Daniele,

Thank you very much for your help! It is very clear to me now.

However, just to be sure, I would like to check the obtained results.

In my simulation, I used a uniform energy distribution for the alpha particles incident on a 40Ar target. I compared the cross sections obtained with Fluka with the ones from TENDL-2019 (extracted from ENDF) and the cross sections calculated with Empire 3.2.2 nuclear reaction code. You can see the results in the attached image.
Is it normal to obtain such a difference in the energy where the inelastic channel opens?

Regards,
Mihaela

Hello Mihaela,

In all the three cases the nuclear data are evaluated. Since different models and parametrization are applied, a certain difference might be expected.

Furthermore, from your plot, the cross section calculated in Fluka seems to be always between the other two cases (for energies higher than 6 MeV).

Let me know if you have any further doubt!

Cheers,
Daniele

Hello again :grin:!

Continuing this subject, I created a mgdraw.f file in order to see what is the probability for the (alpha,n) reaction. The problem is that all the cross sections are 0.

Can you please help me understand what am I doing wrong?

Here is the code:

  ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
  MMAT = 20
  LENOLD = .false.

  particle_code = -6
  • this is a alpha-induced reaction

    IF(ICODE.EQ.101.AND.JTRACK.EQ.-6) THEN
    DO I=1,NPHEAV
    
  • if Ca-43 is produced, print the cross section

    IF((IBHEAV(KHEAVY(I))).EQ.43.AND.(ICHEAV(KHEAVY(I))).EQ.20) THEN
    
  • kinetic energy = total energy - rest energy

    EKE = ETRACK - 3.727379378
    ELKE = log(EKE)
    
    SIGIN = SGTINL( particle_code, EKE, ELKE, MMAT, LENOLD )
    
  • cross section in b

    inxs = SIGIN*40/(1.39*6.023*1e23)*1e24
    
    WRITE(72,*) 'Neutron production', EKE*1000, particle_code, inxs
    
    END IF
    END DO
    
    END IF
    RETURN
    

Hello Mihaela,

If possible, could you attach your mgdraw.f file?
The only reasons (as far as I can foresee) because of the code could not possibly work are the following:

  1. Wrong type definition for any of the variables used when calling the function (i.e. double instead of integer for the MMAT and particle_code).
  2. Missing library (but the compiler should complain).
  3. Notice that, when you compute the inelastic microscopic cross section, you are dividing and multiplying by huge numbers. I would avoid that, since we are in the domain of the machine error. The optimization of the compiler should save you, but I would rewrite it.

Also, notice that, despite you enter in the if condition only in case of (alpha, n) reactions, the computed cross section is the one which contains all inelastic scattering processes.

Cheers,
Daniele

Dear Daniele,

The problem was that I forgot to define particle_code, but the program was printing the right value for it, so I thought that it is already correctly defined.

Regarding this statement,

I thought that this condition would filter the results.

IF((IBHEAV(KHEAVY(I))).EQ.43.AND.(ICHEAV(KHEAVY(I))).EQ.20) THEN

Since this is not the case, what would be the best way to extract individual cross sections?

Thank you,
Mihaela

Hello Mihaela,

I thought that this condition would filter the results.

What the code is doing is checking whether the (alpha, n) reaction has occurred and, in that case, retrieving the inelastic cross section (SGTINL). However, the calculated cross section depends only on the projectile, its energy, and the target.

In case you want to extract individual cross sections I would suggest the following:

  1. With the same mgdraw.f, you score all the inelastic scattering events in a dedicated file (the only variable to save is the energy).
  2. For each (alpha, x) reaction you are interested in, you select the reaction products in an if condition. The one you already selected should work for (alpha, n). Then, you save the projectile energy in a different dedicated file.
  3. You post-process the saved data. Particularly, you need to define energy binning and score the ratio between the number of (alpha, n) reactions and the total number of inelastic scattering ones.
  4. To obtain the (alpha, n) cross section, you need to multiply this ratio with the inelastic scattering cross section calculated for the same energy.

Translated in code, this would be:

IF(ICODE.EQ.101.AND.JTRACK.EQ.-6) THEN
   WRITE(72,*) EKE
   DO I=1,NPHEAV
      IF((IBHEAV(KHEAVY(I))).EQ.43.AND.(ICHEAV(KHEAVY(I))).EQ.20) THEN
         WRITE(73,*) EKE
      ENDIF
   END DO
ENDIF

Please, be aware that this process might be computationally intensive. I would further suggest modifying the other thresholds (EMFCUT, DELTARAY, etc.) to reduce the computational times in calculating EM interactions and delta ray production.

Finally, I just want to remind you that these obtained cross sections are subject to statistical uncertainty.

Let me know in case of troubles.
Cheers,
Daniele

1 Like

Dear Daniele,

Thank you very much for your help!

Best regards,
Mihaela