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):
ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
*
* 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
ICHECK=2
*
* Otherwise, we examine secondary number 1.
ELSE
ICHECK=1
ENDIF
*
* 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
RETURN
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).
Cesc