Flag secondary particles which are the same kind of the one initiating the collision

[thomas.f (Thomas)]

Dear Fluka Team,

I am trying to flag secondary particles based on a specific logic. My objective is for the flag to be TRUE when a secondary particle is the same type as the particle that initiated the collision.

From various forum discussions, I gathered that STUPRE/STUPRF, MDSTCK, and MGDRAW>USRDRAW could be utilized for this purpose. However, I’ve faced challenges in implementing this, primarily due to uncertainty about the correct parameters to employ (e.g., from EMFTSK or GENTSK or …).

Below are the details of my attempts:

1. STUPRE:

• Encountered an issue with the values stored in KPART(NPNW - NPSTRT + 1) or KPART(NPNW) that does not return the correct particle. I guess I do not consider the correct stack here.

      DO 500 NPNW = NPSTRT,NPEMF
         DO 51 ISPR = 1, MKBMX1
            ESPARK (ISPR,NPNW) = SPAUSR (ISPR)
 51      CONTINUE
         DO 52 ISPR = 1, MKBMX2
            IESPAK (ISPR,NPNW) = ISPUSR (ISPR)
 52      CONTINUE
         IP = NPNW - NPSTRT + 1
	     NP = NPEMF - NPSTRT + 1
         WRITE(99,*) 'STUPRE ',REACTION,NP,JTRACK,KPART(IP)
     &            ,TKI(IP),WEI(IP)
	 IF ( JTRACK .NE. KPART(IP)) THEN
            LOUEMF (NPNW) = TRUE
	 ELSE
            LOUEMF (NPNW) = FALSE
	 ENDIF
 500  CONTINUE

2. MDSTCK:

• Encountered a similar problem when accessing KPART(NPSECN).

      DO IP=1, NPSECN  !* SPECN is linked to genstk
        WRITE(99,*) 'MDSTCK ',REACTION,NPSECN,JTRACK,KPART(IP)
     &  	   ,TKI(IP),WEI(IP)
        IF ( JTRACK .NE. KPART(IP)) THEN
           LOUEMF (IP) = TRUE
        ELSE
           LOUEMF (IP) = FALSE
        ENDIF
      END DO

3. MGDRAW:

• While I seem to capture the right secondary, accessing its flag is problematic.
• LOUEMF might not be the right choice as MGDRAW is invoked post-STUPRE. LLOUSE doesn’t seem effective either.

      DO IP=1, NP
        WRITE(99,*) 'MGDRAW ',REACTION,NP,JTRACK,KPART(IP)
     &             ,TKI(IP),WEI(IP)
	IF ( JTRACK .NE. KPART(IP)) THEN
*	   LLOUSE(IP) = TRUE
	ELSE
*	   LLOUSE(IP) = FALSE
	ENDIF
      END DO

In comparing the WRITE function outputs across these routines, I found consistent results for Moller and Bremstrahlung interactions. However, discrepancies arose for other interactions. Interestingly, MGDRAW appears to yield correct results in such cases.

> ROUTINE INTERA                           NSEC      JTRACK        JSEC                      KINE    
>  MGDRAW Moller                              2           3           3   2.3284173093807126E-005  
>  MGDRAW Moller                              2           3           3   1.9751929487762943E-006  
>  STUPRE Moller                              2           3           3   2.3284173093807126E-005  
>  STUPRE Moller                              2           3           3   1.9751929487762943E-006  
>  MDSTCK Moller                              2           3           3   2.3284173093807126E-005  
>  MDSTCK Moller                              2           3           3   1.9751929487762943E-006  
>  
>  STUPRE PhotoElectric                       2           7           3   2.3284173093807126E-005  
>  STUPRE PhotoElectric                       2           7           3   1.9751929487762943E-006  
>  MDSTCK PhotoElectric                       2           7           3   2.3284173093807126E-005  
>  MDSTCK PhotoElectric                       2           7           3   1.9751929487762943E-006  
>  MGDRAW PhotoElectric                       2           7           3   1.2797920875950585E-004  
>  MGDRAW PhotoElectric                       2           7           3   1.5254000000000786E-006  
>  
>  STUPRE Rayleigh                            1           7           3   1.9361811555252762E-006     
>  MDSTCK Rayleigh                            1           7           3   1.9361811555252762E-006     
>  MGDRAW Rayleigh                            1           7           7   6.6306973014815596E-005

I would greatly appreciate your help about this and advices to flag the secondary as a function of its type compared to the initial particle.
In addition, I do not know which of those routines is the most suitable for my purpose and maybe you could guide me about it.

Thanks for your help.
Best regards.

[adonadon (Andre Donadon Servelle)]

Dear @thomas.f,

Thanks for your question.

1. Why do you need to flag the particles in this way? In case you just need to achieve this and dump the information in a file, this is already done in your approach with MGDRAW. I guess you know that, after being transported, the particles’ information won’t accesible anymore…

2. In any case, given your output where only EM particles where utilised, and the fact that you tried to use the routine stupre.f, I guess you are only interested in these particles. Considering this, I would continue your stupre.f option. As MDSTCK is called “after a nuclear interaction in which at least one secondary particle has been produced”, I would not use this option for your problem.

3. I still do not see the reason why it is not working with the stupre routine. Could you provide me the files (.inp and .f) to observe more in detail where the problem is coming from? Also the files will be helpful in case that what I supposed in (1) and (2) are not correct.

Thanks,
Best,

André

[thomas.f (Thomas)]

Dear Andre,

Thank you for your timely responses.

1. I’m currently working on a biasing technique using the USIMBS user routine. In this technique, the particle weights rely on the preceding particle track. However, a challenge arises when an interaction spawns a different type of particle, introducing a discontinuity that I need to rectify. My approach to address this is by flagging particles when this discontinuity occurs. Doing so will allow me to detect these processes within the USIMBS routine and achieve accurate renormalization.

2. While my current test emphasizes em particles for the sake of simplicity, I intend to broaden the scope to encompass all particles in due course. However, for other particles, I recognize that I can leverage STUPRF. The distinction being, STUPRF is invoked for each secondary, whereas STUPRE is called at every interaction.

3. Upon reviewing my code, I’ve observed that the way I’ve implemented the STUPRE routine results in inconsistent outcomes. A case in point is the Rayleigh scattering interaction: in both STUPRE and MDSTCK, the final state comprises an electron and lacks a photon, whereas in MGDRAW, the final state has a photon. I also noticed a difference in the secondary energies between [STUPRE&MDSTCK] and MGDRAW.

Below, the first integer is the number of secondaries, the second one is the code of the initial particle and the third one is the code of the final state particle. The last number is the kinetic energy. You will notice that only the MGDRAW routine provides the correct final state.

>  STUPRE Rayleigh                            1           7           3   1.9361811555252762E-006     
>  MDSTCK Rayleigh                            1           7           3   1.9361811555252762E-006     
>  MGDRAW Rayleigh                            1           7           7   6.6306973014815596E-005

As per your request, I’ve attached a simplified example of the .inp and .f files to this email. Despite its simplicity, it sufficiently highlights the discrepancies observed with MGDRAW, STUPRE and MDSTCK.

benchmark.inp (1.9 KB)
mdstck.f (2.6 KB)
mgdraw.f (12.5 KB)
stupre.f (4.9 KB)

[adonadon (Andre Donadon Servelle)]

Dear Thomas,

Sorry for the delay in replying!

The problem in your routines lies in how you are dumping the information related to the secondary particles.

In STUPRE and MDSTCK you are utilising the KPART variable, which actually represents the particles in your in genstck. What you have to do is to use the variable ICHEMF:

             Ichemf(i) = charge (-1 for e-, 1 for e+, 0 for photon).

What you observe in MGDRAW is correct because the information from the common ‘emfstk.inc’ is copied into the common ‘genstk.inc’ before calling USDRAW.

Thus, you can just use the MGDRAW with the routine that you already have, or modify the others changing KPART to ICHEMF.

Finally, it is important to note that the energies in ‘emfstk.inc’ are given in MeV (in case you need this information).

I hope this helps,
Have a nice day,

André

[thomas.f (Thomas)]

Thank you Andre for your response.
However, I do not understand very well.

My issue in MDSTCK and STUPRE is that the particles and energies of the secondaries are not consistent with those provided in MGDRAW. Using ICHEMF to get the charge could eventually help me to know whether the secondary is an electron, positron or photon. Does that mean there is no variable in this routine to get the particle code? Also, I still do not understand why the energies are not the same as those for MGDRAW. The difference does not seem to be a GeV/MeV conversion.

You suggest me to use MGDRAW, which I can eventually do. In this case however, the LLOUSE variable (from TRACKR) seems to not be recognized. Same for LOUEMF. Which variables should I be using to flag the secondaries in MGDRAW/USRDRAW?

Thanks for your help.

[ghugo (Gabrielle Hugo)]

Dear Thomas,

Thank you for your questions.

All issues ended up connected, and stemming from a deeper fact (mentioned as a side note in Stupre.f and stuprf.f routines / FLUKA stacks - #3 by ghugo):

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

Even more confused? Let’s go.

• MDSTACK / STUPRE versus MGDRAW results discrepancies

My issue in MDSTCK and STUPRE is that the particles and energies of the secondaries are not consistent with those provided in MGDRAW. Also, I still do not understand why the energies are not the same as those for MGDRAW.

In MDSTACK / STUPRE, you are printing secondaries information from GENSTK. You should use EMFSTK instead.
Using the following code (relying on EMFSTK) in MDSTACK / STUPRE (for electron or photon interactions indifferently) fixes your secondaries printouts:

do IT=NPSTRT, NPEMF
           WRITE(99,*) "MDSTCK IT =", IT, 
     &          ", ETEMF(IT) = ",ETEMF(IT), 
     &          ", PMEMF(IT) = ", PMEMF(IT), 
     &          ", UEMF(IT) = ", UEMF(IT), 
     &          ", VEMF(IT) = ", VEMF(IT), 
     &          ", WEMF(IT) = ", WEMF(IT), 
     &          ", WTEMF(IT) = ", WTEMF(IT)
end do
*     Off-topic: note the use of particle direction, better suited than energy
*     to uniquely identify a particle (e.g. Rayleigh scattering).

Why? What is wrong with GENSTK?

GENSTK is always populated correctly at the time MGDRAW is called, to be able to study any interaction.

** “pure” EM interactions with electron primary: MDSTACK / STUPRE are called after MGDRAW. Hence, when they are called, the secondaries are already properly loaded into GENSTK.
This is why your print of secondaries information for Moller or Bremstrahlung interactions from MDSTACK and STUPRE is correct, and matches your print from MGDRAW.

** “pure” EM interactions with photon primary: MDSTACK / STUPRE are called before MGDRAW. At that time, the secondaries are not already loaded into GENSTK.
This is why, for example for Rayleigh scattering or photoelectric interactions, your print of secondaries information from MDSTACK and STUPRE is not correct, while your print from MGDRAW is ok. The secondaries you are printing from MDSTACK / STUPRE with KPART(IP),TKI(IP),WEI(IP) (stored in GENSTK) in your message, are actually the secondaries from the previous electron interaction.

• Let’s now discuss what you would actually like to do for your study. Ie, being able to know, within USIMBS, whether a secondary is same type as its parent.

In MGDRAW, LLOUSE variable (from TRACKR) seems to not be recognized

The use of TRACKR variables (for instance LLOUSE) is licit within MGDRAW.
However, you are using LLOUSE as an array (LLOUSE(IP) = TRUE), while LLOUSE is a scalar.
This points to the next issue: you are trying to store information per secondary, while here we are still at the primary hierarchy level.

Which variables should I be using to flag the secondaries in MGDRAW/USRDRAW?

Instead, you can simply store the parent particle ID, for example by adding:

ISPUSR(1) = JTRACK

This information will be automatically propagated to each secondary before they are tracked.
At tracking time (e.g. in USIMBS), you will then be able to simply compare the tracked particle type (which will be in JTRACK) with its parent type (ISPUSR(1)). Hence, you will be able to know whether each secondary is same particle type as its parent or not.

• Wait, are we not done yet? Why an extra point?

This relates to how the information is automatically propagated with ISPUSR(1) from the primary hierarchy level to all of its secondaries. This propagation happens to be done in STUPRE (for “pure” EM interactions, otherwise it is in STUPRF):

IESPAK (1,NPNW) = ISPUSR (1)

Now, hold on, this is where we close the circle.
For photon “pure” EM interactions, STUPRE is called before MGDRAW.
Hence, if you only set ISPUSR(1) = JTRACK within MGDRAW, you will set it too late (the propagation will already have taken place in STUPRE). You will thus not be able to access the correct parent for any photon “pure” EM interaction.
Instead, all the other interactions will be OK, because for them, STUPRE or STUPRF is called after MGDRAW (hence the propagation happens after ISPUSR(1) is set).
To solve this, you thus need to also add ISPUSR(1) = JTRACK at the top of STUPRE.

The good news is that eventually, these are just 2 lines of code to add

Hope this helps,
Cheers,
Gabrielle