Positron annihilation in various layers

Advanced Features and User Routines
1

Dear all
I have a setup where several micron thick layers of materials surround a Na22 source. I am trying to see how many positrons are stopped in each layer. For that, I have using USDRAW in mgdraw.f code. THe problem is i am not getting any output in the collision files. can some one please take a look and let me know what mistake I am doing? The input and mgdraw files are attached.
thanks & regards
saurabh
silicon_poly_samwich.inp (3.6 KB)
mgdraw_pos_src_corr1.f (14.8 KB)

3

Dear @sxm3816,

You are on the right track, with the following caveats:

  1. In your userroutine, you are selecting the ‘ICODE .EQ. 110’, corresponding to ‘radioactive decay products’. Note that your source is in the center of your geometry and only there you have the decay itself! Afterwards you also select ‘MREG .EQ. 3’ (the ‘kpart’ region), but here the particles that arrive (photons and positrons) are not decaying themselves - they are decay products of your ^{22}Na source but do not decay in the ‘kpart’ region, hence why you have nothing in the output collision files.
  2. I adjusted your userroutine, please find it here.
    mgdraw_pos_src_corr1_adjusted.f (14.9 KB)
  3. Note that if you want to print out ‘TKI’ and not ‘TKI(KP)’, you will print out the entire array of kinetic energy for the entire stack (including zeroes)!

Here, I score every particle that enters ‘MREG .EQ. 3’, and from the results it seems that there are mostly doubly 511 keV photons detected, meaning that the positrons get already annihilated before entering the ‘kpart’ region. Some positrons do survive.

I hope this helps.

Best,
Daniel

4

Dear Daniel
thank you for the wonderful solution. This part worked nicely. There is something more that I was interested in. I did not make that clear in the question. sorry about that.
what I was looking for is the number of positrons annihilating in the various layers. For this, I tried few variation on what you suggested.

  1. The annihilation would be mostly because the positron is below the transport threshold. I am not sure about the in-flight annihilation probability. To record the positron’s position at the moment its KE drops below say 2 keV, I modified the code to have TKI (kp) < 2e-6. With this code, the total number of positrons I got after summing the hits in all the layers was about 10% of the total number of primaries. It should be ~ 90% of total number of primaries.
  2. Every time positron annihilates, one will get the 511 keV gammas. So, I set the USDRAW to look at these events. The mgdraw file is attached. The relevant part of the code is-
    IF (ICODE .EQ. 215 .AND. KPART (KP) .EQ. 7) THEN
    IF (MREG .EQ. 3 ) THEN
    WRITE(93,‘(1P,4G25.15)’) ZSCO, TKI (KP)
    With this code, I get about 82% of the primaries (accounting for emission of 2 annihilation photons).
    One thing could be happening is that some positrons might be escaping the layers. I made them fat and looked at the relevant USRBIN data. There was no positron escape.

Can you please suggest what mistake I might be doing?
thanks & regards
saurabh mukherjee
mgdraw_pos_src_corr1_adjusted.f (16.1 KB)
silicon_poly_samwich.inp (3.8 KB)

5

I think I solved it. I changed the code so that the code writes the 511 keV gammas originating in various regions. From this I calculate the positrons dying in various layers. The relevant part of the code is
IF ( KPART (KP) .EQ. 7 .AND. ICODE .EQ. 214 .OR. ICODE .EQ. 215) THEN
IF (MREG .EQ. 3 ) THEN
WRITE(93,‘(1P,4G25.15)’) ZSCO, TKI (KP)

BTW, is there a more direct way of doing the same thing?
thanks & regards
saurabh
mgdraw_pos_src_corr1_adjusted.f (16.2 KB)

6

Dear @sxm3816,

  1. Note that you are using PRECISION defaults (with a minimum threshold energy E_{thres}=100 keV), so if you want to further reduce the transport/production threshold you can use the EMFCUT card, allowing you to have E_{thres}=1 keV for electrons and positrons. If you are interested in the location of the positron annihilation, I recommend lowering as much as possible the positron threshold.
  2. Indeed, in your second message, you are scoring the locations where your positrons annihilate (thereby producing the two gammas), so it should be correct.
  3. At this stage, using mgdraw presents itself as the most direct way of achieving what you are interested in.

Best,
Daniel

7

Dear Daniel
thank you very much!
sincerely
_S