Scoring the information of all the secondaries from inelastic nuclear interactions

Advanced Features and User Routines
1

Dear FLUKA experts,

I am trying to simulate the interaction of carbon ions with hydrogen (or carbon or oxygen) by considering only inelastic nuclear interactions and to kill the event after one interaction for each projectile (namely, to interact only once).
I am trying to score the type, energy distribution, and angular distribution of all the secondaries (per nuclear event).

I tried to use THRESHOLd to switch off the elastic interactions, EMF=off to switch off the EM interactions, LOW-BIAS with 20MeV to switch off the low energy neutron transport, and MULSOPT with what(2)=what(3)=3 to switch off MCS. But I’m not sure whether these settings are enough (only inelastic nuclear interactions survive & interact only once).

In addition, I used USRDUMP with mgdraw.f to score the secondaries. I found “KPART” denotes the secondary IDs. I also found that “KHEAVY” and “ICRES” should be considered to include the fragments and residuals, respectively. I am a little confused of the definition of “secondaries” in USRDUMP. For example, in the reaction,

12C + 12C = p + n + gamma + 3He + 7Be + 7Li

I think all the particles on the right hand side are called secondaries. But my testing from USRDUMP showed that “p + n + gamma” came from “KPART”, while “3He + 7Be + 7Li” from “KHEAVY”.

Could someone help me in the settings of physics, transport and scoring?

Thanks in advance!

3

Hello Kaiwen,

Welcome to the FLUKA Forum.

To attain your goal, the missing relevant piece I did not see in your suggested course of action is to use a thin target as your geometry and to use the LAM-BIAS card to shorten the inelastic scattering length of your beam particles accordingly. Feel free to take a look at the relevant section of the manual or this lecture from the last FLUKA Beginner Online Training.

A possible/natural way to score secondaries from a nuclear inelastic interaction consists in issuing a USRYIELD card, requesting Type: EMERGING. You will need need several cards for the various angular domains and particle species you’re interested in. With a sufficiently thin target and a sufficiently biased (nuclear) inelastic scattering length for beam particles, the risk of scoring spurious contribution from reinteraction secondaries should be negligible (going overboard one can do a sensitivity study to make sure spectra do not change significantly). To be absolutely sure, though, you could look into the usrmed.f routine (needs to be requested with MAT-PROP) and set the statistical weight to zero for those particles you wish to effectively discard. The USRYIELD approach has the advantage that all the well-tested scoring, processing, and plotting machinery is at your fingertips with Flair.

Alternatively, in a more “do it yourself” way, you could proceed via ENTRY USDRAW in mgdraw.f as you report, keeping an eye on / making sure that you filter by JTRACK=-2, IBARCH(-2)=12, and ICHRGE(-2)=6 (your primary 12C), as well as by LTRACK=1 (beam particle) and by incoming energy ETRACK-AM(-2) not too far away from the beam kinetic energy (kinetic energy as is, not per nucleon). You can e.g. print all secondary info (including their statistical weight) out to a file and apply your own careful post-processing (it’s very easy to slip here and there, so be on your toes). With this filtering conditions prior to your WRITE statement you’re sure to be free from reinteraction secondaries; you can still look at usrmed.f to kill particles at your discretion if you wish. Do keep an eye open in case you stumble upon surviving biased primaries (see manual LAM-BIAS).

Regarding the definition of secondaries, there is no formal problem, it’s only an operational question. In FLUKA you get heavy fragments in include/fheavy.inc, the rest of secondaries in include/genstk.inc, and the residual (if any) in include/resnuc.inc, as you indeed witnessed in your attempt.

Incidentally, do not forget to add a PHYSICS card with SDUM=COALESCE and WHAT(1)=1, and another PHYSICS card with SDUM=EVAPORAT and WHAT(1)=3. Beware to link with ldpmqmd / run with flukadpm if your run ecounters ions with kinetic energy in excess of 100 MeV/n.

Cheers,

Cesc

PS1: If you call your sample thickness dz, a reasonable estimate for the factor f<1 by which to shorten the inelastic scattering length lambda_i (“between consecutive nuclear inelastic interactions” implicitly understood) can be obtained as follows. You extract lambda_i from the table of inelastic scattering lengths in the output file for beam particles in the relevant material. You then set f such that dz/ (f * lambda_i) ~ 1, so that every primary will (on average) undergo roughly one nuclear inelastic interaction. Note that this deliberate distortion of the interaction length is automatically accompanied by the corresponding modification of the statistical weight of the concerned particles and their secondaries, so that the expectation value of physical observables is respected.

PS2: when using the LAM-BIAS card, make sure you select SDUM=INEPRI to just bias the primaries (beam particles).

4

Dear Francesc,

Thank you so much for your quick and detailed response!

I followed your instructions to set the options in flair and mgdraw.f. And I’ve got another two questions.

  1. for the value of “f” in LAM-BIAS, since I need to calculate various kinetic energies for the incident particle (Ek=0~450MeV with an interval of 50MeV), it might be difficult to estimate the value “f” for each energy. As an alternative, I tired to set “f” to a very small value in order to make the number of inelastic interactions greater than the number of incident particles (all the primary particles were interacted inelastically with the target). (In the ***001.out file: particle number=1000, number of inelastic interactions= 3862) But I am confused about the underlying physics. Does it mean not only the primary particles have inelastic interactions, but also the inelastic interactions from the secondaries appear?
    Then in the mgdraw.f file, I set

ICODE.EQ.101 .AND. JTRACK.EQ.-2 .AND. IBARCH(-2).EQ.12
* .AND. ICHRGE(-2).EQ.6 .AND. LTRACK.EQ.1

in order to score the inelastic yields from primary particles only. But the results showed that still 3862 inelastic interactions appear, whose number is the same as the number of inelastic interactions in the ***001.out file. So what did the 1000 primary particles happen?

  1. in the generated file (.72), KPART(1) is always “***”, I don’t know why…
    image
    image1320×218 21.6 KB

I have uploaded the .inp file and the .f file. I would be so appreciated if you could have a look at them.
Thank you very much!

Best,
Kaiwen
mgdraw_carbon_ine.f (10.7 KB)
carbon_ine.inp (1.4 KB)

5

Dear Kaiwen,

I was not clear in my comment above

Do keep an eye open in case you stumble upon surviving biased primaries (see manual LAM-BIAS)

Looking at the section of the manual on the LAM-BIAS card, specifically at WHAT(2), you’ll see that if you enter the “f” factor as a positive number, the biased particle (in your case the primary 12C since you’re requesting INEPRI) always survives, but with a reduced statistical weight (secondaries from subsequent interactions of the surviving biased primary see their statistical weight reduced accordingly). The surviving biased particle appears as the first secondary, in your case with KPART(1)=-601240 (encoding Z=6 and A=12 in an integer, nevermind the dangling 40, it’s irrelevant here). This is why your WRITE statement format failed (cannot fit such an integer in just 3 allocated spaces for output).

FYI, if you enter WHAT(2) as a negative number (-f if you will) then Russian roulette is played against the probability f, killing the primary at some point. Unfortunately, in your case with f=1e-6 the primary will mostly survive and keep on reducing its weight, so you’ll probably see no effect puting WHAT(2) positive or negative, other than the random-number history changing.

Being practical, you can simply adapt your mgdraw.f such that the reported number of secondaries is NP-NP0 and that your WRITE statement has (KPART(I),I=NP0+1,NP). Nothing to touch for the heavies and the possible residual. Variable NP0 is such that NP0+1 is the first actual secondary.

It may be interesting to consider the USRYIELD path above (sparing yourself these troubles).

Season cheers,

Cesc

6

Dear Francesc,

Thank you for your explanation!

Let me try to understand the physics process. When I set f=1e-6 and SUDM=INEPRI in LAM-BIAS, the following things happen

====================================================

  1. only inelastic nuclear interactions appear

  2. for each primary 12C (or most of the primary 12C), it survives by force after inelastic nuclear interaction and appears as one of the generated secondary particles. The kinetic energy and momentum of the 12C do not change, while the weight is reduced.

  3. then the survived 12C interacts again with the target inelastically and generates new secondaries, including itself. And the weight is reduced again.

  4. the above processes happen again and again until the 12C goes out of the target.

====================================================

But actually since I set f=1e-6 (a very small value), the weight is always 1 until the 12C goes out of the target. As I checked and printed WEI(NP0) after each interaction.

And if I set the value of “f” to an even smaller value (e.g. f=1e-8), the number of inelastic nuclear interactions will increase.

Is my understanding correct?

Best,

Kaiwen

7

Hello Kaiwen,

  1. only inelastic nuclear interactions appear

Of course, if you filter by ICODE.EQ.101 you only see the nuclear inelastic interactions, but keep in mind that the biased 12C is still subject to an electronic stopping power along the steps of its trajectory (very short as a result of the present biasing), which is e.g. why the kinetic energy at the interaction point is (ever so slightly) below the beam energy.

  1. for each primary 12C (or most of the primary 12C), it survives by force after inelastic nuclear interaction and appears as one of the generated secondary particles.

Indeed, if you enter WHAT(2)=1e-6 (positive) as you declare above, the biased particle always survives with a reduced weight, as clearly written in the manual. It appears as the first secondary (if NP0>0).

The kinetic energy and momentum of the 12C do not change, while the weight is reduced.

This is true in the immediate aftermath of the biased nuclear inelastic interaction. But when the surviving 12C takes the next (albeit very short) step, it still undergoes dE/dx, and it sees its energy decreased a bit. I did one more / final quick run of your case and picked an event (simply filtering by NCASE) where I spotted the following consecutive 12C kinetic energies in GeV:

9.9999777935689885E-002
9.9999258660380175E-002
9.9999148749514280E-002
9.9999113202211110E-002
9.9999030852911577E-002

As you can see the energy does decrease a little bit.

  1. then the survived 12C interacts again with the target inelastically and generates new secondaries, including itself. And the weight is reduced again.

Yes.

  1. the above processes happen again and again until the 12C goes out of the target.

Yes.

But actually since I set f=1e-6 (a very small value), the weight is always 1 until the 12C goes out of the target. As I checked and printed WEI(NP0) after each interaction.

I ran your input with a slightly modified mgdraw.f:

mgdraw_carbon_ine.f (11.0 KB)
run.inp (1.4 KB)

For an event where the 12C happened to undergo 5 nuclear inelastic interactions (of course each with a varying number of secondaries), the weight of the produced secondaries (looking just in genstk) is the following (where the first element is for the surviving biased 12C).

GENSTK WEI:   9.99999000E-001  1.00000000E-006  1.00000000E-006
GENSTK WEI:   9.99998000E-001  9.99999000E-007  9.99999000E-007  9.99999000E-007  9.99999000E-007
GENSTK WEI:   9.99997000E-001  9.99998000E-007  9.99998000E-007  9.99998000E-007  9.99998000E-007  9.99998000E-007  9.99998000E-007  9.99998000E-007  9.99998000E-007
GENSTK WEI:   9.99996000E-001  9.99997000E-007  9.99997000E-007
GENSTK WEI:   9.99995000E-001  9.99996000E-007  9.99996000E-007  9.99996000E-007  9.99996000E-007  9.99996000E-007  9.99996000E-007  9.99996000E-007  9.99996000E-007

As you can see, the weight of the surviving biased 12C is gradually reduced, as is that of the secondaries. It is then a matter of proper bookkeeping of these weights.

Again: the use of tried and true built-in scorings like USRYIELD as mentioned in the previous post remains the recommended way to go, as it spares you many pitfalls along the way.

With kind regards,

Cesc

8

Dear Francesc,

Thank you for your explanation and results!
I checked the following things after our last discussion.

  1. For the 12C transport, is it possible to switch off the electronic force to make dE/dx=0, so that the kinetic energy of 12C can always keep unchanged?

  2. After each inelastic interaction, although the survived 12C is stored as a secondary particle, its LTRACK still equals to 1, is it right?

image
3) I have outputted the value of LTRACK out of those “IF” conditions, as shown in the figure. And I found there are 9 particles with LTRACK=2. Where do these particles come from? On the other hand, when I outputted all the secondary particles from fheavy.inc (IBHEAV & KHEAVY), I found there are 486 12Cs in the secondaries. Where did these 12Cs go? They are not treated with LTRACK=2? Or are they just stopped?

The reason that I prefer to use USRDUMP instead of those built-in cards (such as USRYIELD) is, the former can perform customized scorings and output so many variables in detail. It helps to understand the underlying physics and check whether the results satisfy my requirements, although it is complicated. Thank you very much for your kindly suggestion and Merry Christmas!

Best,
Kaiwen

9

Hi Kaiwen,

Just a quick pre-christmas reply to close the thread.

I am not sure one can easily switch off dE/dx. But do note that the changes in the kinetic energy are essentially irrelevant: we’re speaking of minor variations in the 6th significant figure (!).

After each inelastic interaction, although the survived 12C is stored as a secondary particle, its LTRACK still equals to 1, is it right?

Jupp.

Regarding the LTRACK=2, is can well be that once in a blue moon (9 times out of thousands of interactions) there may be a secondary which manages to reinteract.

Assuming you stayed on 12C on oxygen (natural composition), from the various 12C+16O or 12C+17O or 12C+18O interactions I understand you may get 12C among the secondaries. This is a legitimate secondary particle or “fragment”, not to be confused with a surviving biased primary.

Enjoy the USDRAW adventure and happy holidays,

Cesc

10

Dear Francesc,

Thank you very much for your patience and happy holidays!

Best,
Kaiwen