Stupre.f and stuprf.f routines / FLUKA stacks

acelenta · 28 March 2023 08:12

Dear colleagues, I am trying to implement the same modification to the FLUKA code as reported in this topic, for inelastic interactions.

In my setup, I have an electron beam impinging on a very thin target, and I want to bias photonuclear and electronuclear interactions just in that material via LAM-BIAS. However, I do not want to have the projectile being “cloned” in the final state. LAM-BIAS, by default, does that either via Russian Roulette or by always cloning the projectile in the final state, but with reduced weight.

I wanted to replicate the solution proposed in the aforementioned topic for \gamma N \rightarrow \mu^+\mu^- N, and thus I started to look at stupre.f and stuprf.f. According to the manual, stupre.f should be called when secondary electrons/photons are loaded in the stack, while stupr.f deals with hadrons, muons, …

I tried to understand more about this. First, I wrote a simple code for mgdraw.f that prints out all secondaries. This works properly, when the code is executed the secondaries from any interaction are printed.

ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
      IF (JTRACK.EQ.3) THEN
        CALL GEOR2N(MREG,REGNAME,IREGERR)
      	PRINT *,"MGDRAW_1",ICODE,PTRACK,cztrck,REGNAME,NCASE
      	do kp=1,np
      	   PRINT *,"MGDRAW_1_SEC ",kp,kpart(kp),plr(kp),czr(kp)
      	end do
      END IF

I also modified stupre.f and stuprf.f as follows:

stupre.f

[...]
 52      CONTINUE
         LOUEMF (NPNW) = LLOUSE
c THE FOLLOWING IS THE MODIFIED LINE
	 PRINT *,"STUPRE.f called, TRACK:",JTRACK,PTRACK,ZFLK 
            
 500  CONTINUE
[...]

stuprf.f

[...]
  200 CONTINUE
  
c THE FOLLOWING ARE THE THREE MODIFIED LINES
			PRINT *,"STUPRF.f called IJ-NPSECN: ",IJ,NPSECN
			PRINT *,"STUPRF.f TRACK: ",JTRACK,PTRACK,cztrck,Wtrack
			PRINT *,"STUPRF.f GENSTK : ",KPART(NPSECN),PLR(NPSECN),czr(npsecn),wei(NPSECN)		
[...]

Naively, I would have expected that, following an inelastic interaction, both stupre.f and stupr.f are called, if in the list of secondaries there are both photons/e+/e- and other particles.

For example, the following is the result of an interaction, as printed by mgdraw.f. There are 7 secondaries.

The first is actually the primary electron, that is “cloned” in the list of secondaries (same momentum, same cz).
The second is the final-state electron that is present in the interaction (interesting, different momentum but same cz).
Then, one neutron and 4 photons.

MGDRAW_1         101   99.999184647158955       0.99999998433983484      MMPCB40           14
 MGDRAW_1_SEC            1           3   99.999184647158955       0.99999998433983472     
 MGDRAW_1_SEC            2           3   99.980113297302054       0.99999998433983472     
 MGDRAW_1_SEC            3           8   2.4204667357567395E-002 -0.35656348784275327     
 MGDRAW_1_SEC            4           7   3.7739701897341865E-003  0.89226871675986441     
 MGDRAW_1_SEC            5           7   1.2530090283270269E-003  0.13674998267729363     
 MGDRAW_1_SEC            6           7   2.2254216640014818E-003 -0.87377399096657848     
 MGDRAW_1_SEC            7           7   6.4455971471241453E-004 -0.42476474927701219

I would have expected that there were 6 calls to stupre.f and 1 call to stuprf.f (for the neutron).

However, this is the result I see.

 STUPRF.f called IJ-NPSECN:            3           1
 STUPRF.f TRACK:            3   99.999184647158955       0.99999998433983484        1.0000000000000000     
 STUPRF.f GENSTK :            3   99.999184647158955       0.99999998433983472       0.99990000000000001     
 STUPRF.f called IJ-NPSECN:            3           2
 STUPRF.f TRACK:            3   99.999184647158955       0.99999998433983484        1.0000000000000000     
 STUPRF.f GENSTK :            3   99.980113297302054       0.99999998433983472        1.0000000000000000E-004
 STUPRF.f called IJ-NPSECN:            3           3
 STUPRF.f TRACK:            3   99.999184647158955       0.99999998433983484        1.0000000000000000     
 STUPRF.f GENSTK :            8   2.4204667357567395E-002 -0.35656348784275327        1.0000000000000000E-004
 STUPRF.f called IJ-NPSECN:            3           4
 STUPRF.f TRACK:            3   99.999184647158955       0.99999998433983484        1.0000000000000000     
 STUPRF.f GENSTK :            7   3.7739701897341865E-003  0.89226871675986441        1.0000000000000000E-004
 STUPRF.f called IJ-NPSECN:            3           5
 STUPRF.f TRACK:            3   99.999184647158955       0.99999998433983484        1.0000000000000000     
 STUPRF.f GENSTK :            7   1.2530090283270269E-003  0.13674998267729363        1.0000000000000000E-004
 STUPRF.f called IJ-NPSECN:            3           6
 STUPRF.f TRACK:            3   99.999184647158955       0.99999998433983484        1.0000000000000000     
 STUPRF.f GENSTK :            7   2.2254216640014818E-003 -0.87377399096657848        1.0000000000000000E-004
 STUPRF.f called IJ-NPSECN:            3           7
 STUPRF.f TRACK:            3   99.999184647158955       0.99999998433983484        1.0000000000000000     
 STUPRF.f GENSTK :            7   6.4455971471241453E-004 -0.42476474927701219        1.0000000000000000E-004
 STUPRE.f called, TRACK:           7   1.2530090283270269E-003   0.0000000000000000     
 STUPRE.f called, TRACK:           7   1.2530090283270269E-003   0.0000000000000000

Here, stuprf.f is called 7 times, for all the secondaries of the track that was printed in mgdraw.f. Then, I also see stupre.f called, but for a different track.

I am quite puzzled by this, since, as reported by the manual, stupre.f should be called when electrons/photons are loaded in the stack. Here, it seems stuprf.f is called instead.

Thanks,
Bests,
Andrea

ghugo · 31 March 2023 09:31

Dear Andrea,

You might want to update your implementation in STUPRE. (1) You have put a printout about the primary only within the secondaries loop, hence the printouts from STUPRE can be a bit confusing. (2) In addition, you might have wanted to print CZTRCK (from trackr.inc) instead of ZFLK (from flkstk.inc).
For each “pure” electromagnetic interaction, ie the electromagnetic interactions with atoms and electrons, STUPRE is called.
It is called once once per interaction, outside of the secondaries loop (hence why there is such a loop inside STUPRE).
After the call to STUPRE, the EMFSTK (including its user variables) is flushed into the TRACKR stack.
For photonuclear and electronuclear interactions, STUPRF is called.
It is called for each secondary.
After the calls to STUPRF, the FLKSTK (including its user variables) is flushed into the TRACKR stack.
Hence in your case, since you have an electronuclear reaction: you should have 7 calls to STUPRF, and no call to STUPRE.
You should hence place your handling inside STUPRF, not STUPRE. You could hence use LOUSE(NPFLKA), from FLKSTK, for tagging the primary electron. FLKSTK will in turn be flushed into the TRACKR stack at the beginning of tracking.
I agree the manual can be a bit confusing on this, we should update it.
It should notably specify when STUPRE and STUPRF are called, ie for which interactions, and how many times (inside/outside secondaries loop). It should stress that it is the interaction defining which between STUPRE and STUPRF is called, not the content of the stacks (notably, e+/e- can be secondaries observed in STUPRF). Finally, it should explicit which stack is being copied to TRACKR after STUPRE / STUPRF call.

Thanks a lot for noticing this,
Hope this helps,

Gabrielle

NB: STUPRE/ STUPRF versus USDRAW calls ordering:

“pure” electromagnetic interactions: for γ projectiles, STUPRE is called before USDRAW, while for e+/e- projectiles, STUPRE is called after USDRAW.
photonuclear and electronuclear interactions: STUPRF is called after USDRAW.

acelenta · 31 March 2023 12:51

Dear ghugho,
thanks, your reply is very clear! I think it would be really helpful to fix these details in the manual. For example, it is now clear to me that a particle is inserted in a given stack not depending on its particle ID, rather due to the process that produced it. Electron-nuclear reactions, for example, will insert an electron (the scattered one) in the GENSTK stack, and not in the EMFSTCK.

I have a further question concerning this.

For biased photo-nuclear reactions, based on your answer I expect STUPRF to be called N+1 times, where N is the number of final state particles, and the factor +1 accounts for the biased photon, that is replicated in the final state. May I ask you if this is correct? I was not able to find this call for the biased photon. Instead, I was seeing this for electron-nuclear reactions, where I see two calls referring go an electron in the GENSTCK, the first being the replicated electron, the second being the scattered one.

Thanks,
Andrea

ghugo · 3 April 2023 15:23

Dear Andrea,

Sure, will transmit to the colleagues working on the manual!

Regarding electronuclear interactions, I confirm that with WHAT(2) > 0, one should always see the primary electron in the final state.

Regarding photonuclear interactions, the biasing is achieved slightly differently (no “trick” of adding the primary to the final state), and even with WHAT(2) > 0, the primary photon is guaranteed to never be in the final state.

Hope this helps,
Cheers,
Gabrielle