Mgdraw routine to determine the number of Compton events before absorption

BV96 · 11 May 2023 08:34

Dear Community

My name is David and I have the following goal:

Goal
I have to determine the mean number of Compton scattering events for a given photon primary before it gets absorbed by photoelectric absorption (or drops below the simulation threshold) in a given scintillator.

Approach
For that purpose, I took a look at the available scoring cards in FLUKA and came to the conclusion that I have to use the userdump card together with a custom mgdraw routine. So, I have written a first draft of a mgdraw routine using the corresponding USDRAW section for txt-readable output:

My idea is to log all Compton scattering events of the primary particle (photon) together with the photoelectric events. So I have used the ICODE.EQ.219 flag for Compton scattering and the ICODE.EQ.221 flag for photoelectric absorption.
To identify, which primary is involved in these events, I plan to use NCASE variable from the SODRAW entry.
To ensure that the Compton and photoelectric events are caused by the primary and not the secondaries, I tested to log the stack pointer (NPFLKA), the particle generation (Loflk) and the particle identity (ILOFLK). I would have expected that Loflk>1 or NPFLKA>0 would indicate that secondaries are involved…
As a source, I use the BEAMPos card with the flood option to get an isotropic photon flux within a sphere containing my scintillator.

Problem
With my first mgdraw routine draft , I obtained unexpected results, i.e. multiple photoelectric absorption events for the same photon primary or even Compton scattering events after a photoelectric absorption event (both shouldn’t be possible from a physical point of view).

Problem Solving Status
I have the suspicion that these problems are either caused by my rudimentary Fortran coding skills, i.e. bugs, or that I score also secondaries and can’t identify them properly with my approach (no. 3 above).

Questions

Is my general approach correct or do you have maybe an easier way to achieve my goal?
Is there a problem with my mgdraw routine?
How do I access/log the particle identity (i.e. a unique identifier for the primary photon) properly? I have the suspicion that my indexing for stack variables is wrong (e.g. Loflk(I))…
I’ve tested also a second way to check for Compton scattering and photoelectric events by using the LCMPTN and LPHOEL flags from EVTFLG. However, the result is unfortunately the same. So it seems my two methods are identical, aren’t they?

Attachments
For your convenience, I attach here all input, flair and output files:
input_file.flair (5.6 KB)
input_file.inp (4.2 KB)
mgdraw_ed.f (15.4 KB)
test1001.out (77.5 KB)
test1001_fort.txt (4.0 KB)

In addition, I add here my custom section from the mgdraw routine:

ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
      IF ( .NOT. LFCOPE ) THEN
         LFCOPE = .TRUE.
         IF ( KOMPUT .EQ. 2 ) THEN
            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
         ELSE
            FILNAM = CFDRAW
         END IF
         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
     &          'UNFORMATTED' )
      END IF
	  
	  IF(ICODE.EQ.219) THEN
*	  IF(LCMPTN) THEN
		 WRITE(22,*) 'NCASE ',NCASE,'NPFLKA ',NPFLKA,'Jtrack ',Jtrack,'gen ',Loflk(I),'id ',ILOFLK(I),NSTMAX,' Compton scattering'
*		 WRITE(22,*) 'NCASE ',NCASE,'NPFLKA ',NPFLKA,'Jtrack ',Jtrack,'gen ',Loflk(NPFLKA),NSTMAX,' Compton scattering'
	  END IF
	  
	  IF(ICODE.EQ.221) THEN
*	  IF(LPHOEL) THEN
	     WRITE(22,*) 'NCASE ',NCASE,'NPFLKA ',NPFLKA,'Jtrack ',Jtrack,'gen ',Loflk(I),'id ',ILOFLK(I),NSTMAX,' Photoelectric' 
*		 WRITE(22,*) 'NCASE ',NCASE,'NPFLKA ',NPFLKA,'Jtrack ',Jtrack,'gen ',Loflk(NPFLKA),NSTMAX,' Photoelectric' 
	  END IF
      RETURN

Additional Information
I performed all simulations with the FLUKA CERN version 4-3.2 together with the FLAIR interface version 3.2-4.1 on Ubuntu 18.04 LTS WSL.

Thank you very much for your support
Cheers,
David

gmazzola · 15 May 2023 16:34

Hello David,

sorry for the delay in the reply but I was looking into the details of your request together with @dcalzola.

Let me start pointing some other threads in the forum that could help in describing the solution I will propose you, in particular this one and a second one.
Regarding your questions, I will try to answer them explaining the solution I propose to you.

Your general approach is correct but unfortunately not sufficient.
mgdraw.f is the right tool to print interaction by interaction the Compton scattering (and photoelectric effect) happening for the primary photon interacting in the detector region.
As you state, to reach your aim, it is important to correctly identify just the Compton scattering correlated to the primary photon excluding instead any contributions that can be generated by photons created by secondary electrons. This can be done in strupre.f.
stupre.f is a routine that gives the possibility to access information of the primary and secondary EMF particles placed in the EMFSTK (the stack for emf particle - LOFLK and ILOFLK that you are using are instead for the main stack of FLUKA) and it can be used to assign user-properties to this kind of particle. With this routine, a logical flag can be created to signal that the primary photon had a Compton scattering and to associate it to the outgoing photon created in the reaction. This process then can be followed for any Compton scatterings happening afterwards until the primary photon is absorbed (or leave the detector region!).
Together with the flag, in the stupre.f routine we can also define a counter for the Comptons.

I am attaching a combination of the two routines that works for your case:
mgdraw_ed.f (12.4 KB) → it prints the Compton events and photoelectric effect of the primary photon
stupre.f (6.7 KB) → it generates the flag to be associated to the outgoing photon from the primary photon’s Compton scattering and it counts the number of Comptons before the photoelectric effect

Finally let me just state that NCASE is a variable accessible from caslim.inc and that some of the primaries don’t even enter the detector region or they could enter without interacting or interacting for then leaving.

I hope this can help you.
Cheers,
Giuseppe

BV96 · 16 May 2023 10:54

Hi Giuseppe

First, thank you very much for your detailed reply. I really appreciat the time and effort you invested in answering my questions. That said, I tested your routines and after studying the results, I have still a couple of small questions:

Questions

I would expect that the “Compton done” column in the output file for a photoelectric event is always equal to the corresponding preceeding Compton event number. However, this is not always the case, e.g. in the last photoelectric event for the primary 2098 in the output file attached below. In this case, no Compton event is registered before. Do I miss something here?

         NCASE         Jtrack            gen            flag   Compton done   Energy gamma          ICODE               
           2096              7              3              1              2    0.18625E-03            221  Photoelectric
           2098              7              1              0              0    0.30000E-02            221  Photoelectric
           2098              7              2              1              2    0.32294E-04            221  Photoelectric
           2100              7              1              0              1    0.30000E-02            219        Compton

I would expect that for a unique primary NCASE, there is always only 1 photoelectric event possible (because we are only considering the primary photon and no secondary photons). However, I found some counter-examples in the output file, e.g. 2078, 2098, 2520, 3558, 3575 or 4371. How can you explain these events?

         NCASE         Jtrack            gen            flag   Compton done   Energy gamma          ICODE
           2072              7              3              1              2    0.25493E-03            221  Photoelectric
           2078              7              1              0              1    0.30000E-02            219        Compton
           2078              7              2              1              2    0.50821E-03            219        Compton
           2078              7              3              1              2    0.28306E-04            221  Photoelectric
           2078              7              3              1              3    0.19570E-03            219        Compton
           2078              7              4              1              3    0.17621E-03            221  Photoelectric
           2079              7              1              0              1    0.30000E-02            219        Compton

         NCASE         Jtrack            gen            flag   Compton done   Energy gamma          ICODE
           3555              7              1              0              1    0.30000E-02            219        Compton
           3558              7              1              0              0    0.30000E-02            221  Photoelectric
           3558              7              2              1              1    0.28610E-04            221  Photoelectric
           3575              7              1              0              0    0.30000E-02            221  Photoelectric
           3575              7              2              1              1    0.28306E-04            221  Photoelectric
           3580              7              1              0              1    0.30000E-02            219        Compton

         NCASE         Jtrack            gen            flag   Compton done   Energy gamma          ICODE
           4371              7              1              0              1    0.30000E-02            219        Compton
           4371              7              2              1              2    0.28581E-03            219        Compton
           4371              7              4              1              2    0.28610E-04            221  Photoelectric
           4371              7              3              1              2    0.16430E-03            221  Photoelectric

Then a more physics-related question: Does your code account for “below EM threshold” events (Jtrack=211), i.e. a primary photon, which after some Compton scattering events might end up with an energy below the EM threshold? Because these events would most likely get absorbed by photoelectric events in the scintillator (if we neglect the escape of these photons from the scintillator), one should count these events also as “photoelectric absorption”, altough in the MC code, no photoelectric event has happened. What do you think?
Then I have a more technical question: I checked you code and found several locations, where you have used the ISPUSR(1,2) variable. To my knowledge, this is a user defined spare flag for the current particle in TRACKR. So far so good. However, I didn’t find the place, where you actually define/assign variables to these flags. Where do they get assigned and how are they related to IESPAK(1,2)? I have the suspicion that IESPAK(1)==ISPUSR(1) and IESPAK(2)==ISPUSR(2)… A follow-up question would then be, why you use at several occasions the logic “ISPUSR(2).EQ.1 .OR. LTRACK.EQ.1” in the mgdraw routine? Don’t give these two variables the same information (under the assumption IESPAK(2)==ISPUSR(2))?

Thanks again for your great support. I’m looking forward to your answer.
Cheers,
David

Attachment
input.flair (5.6 KB)
input_file.inp (4.2 KB)
mgdraw_ed.f (12.4 KB)
stupre.f (6.7 KB)
test1001.err (22.4 KB)
test1001.out (82.9 KB)
test1001_fort.txt (254.1 KB)

gmazzola · 17 May 2023 16:02

Hello David,

Thank you for your feedback. Indeed the routine could have been optimized for your aim.

I will answer your questions before proposing you an updated version of the stupre.f routine. In particular let me start replying to your last question that can be useful to better understand how the routines system works.

You are right in your suspicion and in fact you can find the correlation between IESPACK and ISPUSR in line 95 and line 117 of the new stupre.f I am attaching in this reply. These variables are used in order to associate the flag for the primary photon and the counter of Comptons to the transported particle.
The reason why the flag is needed is because the primary photon has LTRACK (the generation number of a particle) equal to 1 only until a first Compton (or photoelettric) reaction takes place. The photon outgoing the Compton scattering is treated by FLUKA as a particle with LTRACK increased of 1.
This means that LTRACK can not be used as an unique identifier for the primary photon along the different Compton scatterings and, instead, it makes necessary to create a flag in order to follow the primary photon along the reaction chain.
This explains also why is needed the logic “ISPUSR(2).EQ.1 .OR. LTRACK.EQ.1”. The first part is to consider the photon outgoing the Compton scattering of the primary photon (so it is a flagged photon since LTRACK is increased) while LTRACK.EQ.1 is to take into account also the primary particle (that is unflagged, so it has ISPUSR(2).EQ.0 - this can be seen also in the printing from mgdraw.f). I hope this is helpful to understand the logic used for the routines.

Now let me comment the strange results printed by mgdraw. All of them can be explained by the fact that a Compton scattering generates a scattered photon, an electron and - eventually- fluorescence and/or an Auger electron (same withtout the scattered photon for the photoelectric - slide 21 here). The logic of the routine as implemented before was not able to distinguish between the scattered photon and a fluorescence photon that then was being flagged and its reactions printed:

In this case the second photoelectric effect is due to a photon emitted by fluorescence (it can be suspected also by the gen number - LTRACK - that is increasing) and the Compton counter was not reset to zero.
For 2078 → the first photoelectric effect is due to fluorescence generated in the Compton before.
For 3558 and 3575 → same as 1.
For 4371 → In this case the fluorescence photon was even having a Rayleigh scattering before being absorbed (this can be seen since the first photoelectric is due to a photon with gen = 4 and not 3 as the second photoelectric - aka the primary photon)

In order to consider also the possible creation of fluorescence in the reactions, I attach an updated stupre:
stupre.f (7.1 KB) → now we create the flag only once Compton scattering has happened and we associate it to the only photon present between the first two particles stacked (they are in any case the scattered photon and the electron)

Finally a comment on the last question:
3. With default PRECISIOn the particle transport threshold is 100 keV and the transprort limit in FLUKA is set to 1 keV as you are doing in the input. I would say that it is very unlikely that a particle with such energy is created before the photon is being absorbed.

Let me know if you have any feedback.

Cheers,
Giuseppe

BV96 · 18 May 2023 09:21

Dear Guiseppe

Perfect, thank you so much for your additional work. I tested the updated stupre routine and the result is now exactly what I was expecting. Your analysis of the unexpected results from the previous test is also intuitive and well explained. Last but not least, I did some additional tests and have now enough confidence to share your thoughts about the low probability of “below EM threshold” for primary photons.

Thanks again for your great support
Cheers,
David

BV96 · 24 May 2023 14:17

Dear Community

My name is David and due to newly encountered rare errors with the user routines described in this post here, I would like to reopen the discussion.

Problem
During my simulations, I encountered again some rare unphysical events, i.e. multipe photo-electric absorptions (ICODE=221) for the same primary photon. The problematic NCASE numbers are as follows: NCASE =[15716,21259,34326,36724,51853,61132,65134,76308,84728,90142,99440].

Here are two examples from the output file (*_fort.txt file attached below) for the events 15716 and 99440, respectively:

          NCASE         Jtrack            gen            flag   Compton done   Energy gamma          ICODE               
          15716              7              1              0              1    0.20000E-02            219        Compton
          15716              7              2              1              1    0.28306E-04            221  Photoelectric
          15716              7              2              1              2    0.28325E-03            219        Compton
          15716              7              3              1              2    0.26035E-03            221  Photoelectric

          99440              7              1              0              1    0.20000E-02            219        Compton
          99440              7              3              1              1    0.28610E-04            221  Photoelectric
          99440              7              2              1              2    0.97397E-03            219        Compton
          99440              7              3              1              3    0.90240E-03            219        Compton
          99440              7              4              1              3    0.52765E-03            221  Photoelectric

Questions

Could this unexpected behavior be caused again by unintentionally scored secondary photons (e.g. due to pair production) or by the already discussed “below threshold” events?
How can I adapt/correct the user routines (stupre & mgdraw attached below) to prevent these errors?

Attachment
test1001.out (83.2 KB)
test1001_fort.txt (2.8 MB)
input.flair (6.0 KB)
input.inp (4.4 KB)
stupre.f (7.1 KB)
mgdraw_ed.f (12.4 KB)

Thank you very much for your support
Cheers,
David

P.S. Thank you Dávid Horváth and Roberto Versaci for reopening and moving the post.

gmazzola · 25 May 2023 12:04

Hello David,

Thank you for your feedback.
I double checked the routine and indeed the problem was still related to flagging of some fluorescence photon created after Compton scattering of the primary photon.
The (eventual) fluorescence photon and/or Auger electron were not correctly flagged due to the logic of the routine and in some rare cases the photon was inheriting the positive flag causing then the printing of its photoelectric effect from mgdraw.

You can find here a stupre.f that now solves this problem (you can actually see that the difference in the logic is minimal):
stupre.f (7.1 KB)

Thanks,
Giuseppe

BV96 · 25 May 2023 12:33

Hi Guiseppe

Thanks for your swift reply and support. I really appreciate your help. After testing the revised routine, I can confirm that I don’t get any more any unexpected events. However, I wonder, why the same logic adaption is not required for line 93 in the stupre.f routine, i.e. why we have to keep “IF ( NPSTRT .LT. NPEMF) THEN” but change “IF ( NPSTRT+2 .LT. NPEMF) THEN” to “IF ( NPSTRT+2 .LE. NPEMF) THEN”?

Cheers,
David

gmazzola · 25 May 2023 13:19

Hello David,

the first if condition basically takes into account that after the Compton scattering one photon and one electron are always uploaded in the EMFSTK since this means then that NPSTRT will be always lower than NPEMF (for example starting with an empty stack, after a Compton that creates a photon and an electron, NPEMF will be 2 and NPSTRT is 1).

The second if condition takes into account all the possible fluorescence/Auger created (and always stacked after the photon and electron). With the LE we are now taking into account that one single fluorescence (or Auger) could have been created in the interaction. Using LT before, we were considering situation in which both the Auger and fluorescence were created in the same Compton.

Hope this help,
Giuseppe

BV96 · 25 May 2023 13:33

Hi Guiseppe

I see. Thanks for the additional explanation and again for your great support.

Cheers,
David