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?

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:

IJ = particle type (paprop numeration)

EKE = particle kinetic energy (GeV)

ELKE = log (eke)

MMAT = material number of the medium

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

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:

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.

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 for your file. I have attached a working copy in the end. A few remarks:

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.

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)

SGTINL is a function already defined in the code. It should be declared as

double precision, external :: SGTINL

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

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)

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?

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?

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:

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

Missing library (but the compiler should complain).

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.

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?

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:

With the same mgdraw.f, you score all the inelastic scattering events in a dedicated file (the only variable to save is the energy).

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.

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.

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.