(p.n) reaction position-MGDRAW

Dear FLUKA experts,
My simulation is a proton beam incident on a metal target(there exist water channels inside target).
I want to get the information about the position of (p,n) reaction, and I have read manual and the forum.
it seems that I need a customized mgdraw (usdraw entry) user routine.
I have never used mgdraw.f.
for this purpose, could you please help me have a look if the code is correctly modified?
here you can find my input file and mgdraw.f.10_1.inp (10.1 KB) mgdraw_pn.f (3.4 KB)
any help will be appreciated.


Dear FLUKA experts,

I used mgdraw (usdraw entry) user routine to record the information of (p,n) reaction.
and I already got USERDUMP output files with five cycles running.
the datas in five output files are different.
I don’t know how to access the datas and how to get the relevant plots.
and I tried to process it on flair as USRBIN output file, but it didn’t work.

could you please tell me how to analysis the datas and get the plots.

Best regards,

Qi Ding

Dear Qi,

Note that your mgdraw_pn.f is essentially empty, so it cannot do much. The idea is to populate the relevant entries with custom user code as needed.

Since nuclear inelastic interactions of protons in Ta and water at the 70 MeV you request are in the order of 10 and 80 cm respectively (search output file for “Inelastic” and you’ll see a table) and your geometry is thinner than 1 cm (adapt as needed if I misread), you may want to bias the interaction length between consecutive nuclear inelastic processes in Ta and water, so as to effectively steer CPU time into sampling these otherwise seldom interactions:

LAM-BIAS                 -0.01  TANTALUM    PROTON    PROTON        1.INEPRI
LAM-BIAS                -0.001     WATER    PROTON    PROTON        1.INEPRI

This will shorten the inelastic scattering length of protons in Ta and H2O by .01 and .001 respectively. The negative sign in front means that the primary proton will survive/die playing Russian roulette (see manual). The INEPRI means the biasing will apply to the primary particles (protons) only. Reduce the number of primaries to, say, 100 and scale up as needed once you see how much CPU time per primary is needed.

Try putting something along the following lines in ENTRY USDRAW in your mgdraw_pn.f (logics a bit exaggerated for the sake of clarity, adapt as needed):

*     If there was a nuclear inelastic interaction:
      IF (ICODE.EQ.101) THEN
*        If the projectile is a proton:  
         IF (JTRACK.EQ.1) THEN
*           If the first secondary is of the same species and has the
*           same energy as the primary proton, it means the biased
*           primary proton survived Russian roulette and is in the
*           secondaries list, so we examine secondary number 2.
            IF (KPART(1).EQ.1 .AND.
     &         ABS(ETRACK-AM(1)-TKI(1)).LT.1.0D-13*TKI(1)) THEN
*           Otherwise, we examine secondary number 1.
*           If there was only one secondary (possibility of surviving
*           biased primary accounted for by the logics above):
            IF (NP.EQ.ICHECK) THEN
*              If said secondary is a neutron:
               IF (KPART(ICHECK).EQ.8) THEN
*                 Write the (p,n) interaction point to a file *_fort.99
                  WRITE(99,*) XSCO,YSCO,ZSCO
               END IF
            END IF
         END IF
      END IF

You are then supposed to compile your mgdraw_pn.f and bundle it into a custom binary either comfortably in Flair or if you’re a shell fan via e.g.:

fff mgdraw_pn.f
lfluka -m fluka -o myflu mgdraw_pn.o

And then use the resulting myflu binary for your run.

On the physics side: you are shooting a 70-MeV proton beam onto Ta and water. At this energy, the explicit (p,n) channel, i.e. emitting a single neutron, will be suppressed compared to other channels producing many more secondaries, so you will stumble upon a naked (p,n) event once in a blue moon. You can examine this by completing the block of code above with:

*        If the projectile is a proton:  
         IF (JTRACK.EQ.1) THEN
*           Write out number of secondaries and their respective species 
            WRITE(98,*) NP,KPART(1:NP)

so as to have a general idea of what comes out (see the resulting *_fort.98 file).

Adapt the logics above accordingly in case you meant (p,xn).

Note that with USDRAW you are dumping customized output in whichever format you wish, so it’s not immediate to plot with Flair. You will have to post-process your files manually (i.e. make a histogram of positions or whatever you need in the end).


1 Like

Dear Cesc,

Thanks for your detailed response.
regarding my simulation I need record the positions where neutrons are produced by (p,n) reaction.
for this purpose, I modify mgdraw (usdraw entry) user routine accordingly.
could you please tell me if it’s right modifications in my mgdraw.f filemgdraw.f (3.6 KB) to record the positions of (p,n) reaction?

best regards,

Qi Ding

Dear Qi,

The mgdraw.f you submitted goes in the right direction. Here are a few comments:

  • The INCLUDE header structure implies your mgdraw routine is not compatible with FLUKA-4 as distributed by CERN under https://fluka.cern and supported in this forum. Adapt the INCLUDE headers accordingly (compare with src/user/mgdraw.f). You may comfortably use the provided python script bin/fluka_src_converter.py

  • Note that as you have it now you print the position (along with kinetic energy, direction unit-vector, and statistical weight) for neutrons produced during a nuclear inelastic interaction regardless of the species of the primary particle. You may want to filter by proton, and if you really want just the p,n reaction please refer to the detailed reply above, where the source code is highly commented to help you navigate it.


Dear Cesc,

thanks for your detailed response.
I’ll try to adjust the magraw.f file according to your comments.

Best regards,

Qi Ding