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

thomas.f · 21 October 2023 02:01

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:

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

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

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 · 27 October 2023 11:09

Dear @thomas.f,

Thanks for your question.

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…
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.
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 · 28 October 2023 06:05

Dear Andre,

Thank you for your timely responses.

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.
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.
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 · 22 November 2023 15:26

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 · 23 November 2023 10:38

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.