Gamma cascade energy from neutron capture on gadolinium

I’m investigating neutron captures on gadolinium-doped water and I am finding that the energy of the gamma cascade produced when a neutron captures on a Gd nuclei does not match with the expected energy of ~8 MeV.

I am using USDRAW to score the captures with the following code:

         ! Low energy neutron interaction
         IF (ICODE .EQ. 300) THEN
            ! Loop over secondaries
            IF(JTRACK .EQ. 8) THEN
               ! Reset gamma cascade energy counter
               GAMMA_NEUTRON_ENE = 0
               DO I = 1, NP
                  IF (KPART(I) .EQ. 7) THEN
                     GAMMA_NEUTRON_ENE = GAMMA_NEUTRON_ENE + TKI (I)
                  END IF
               END DO
               IF (ICRES .EQ. 64) THEN
                  ! Look at secondaries -- are they all gammas?
                  DO I = 1, NP
                     PRINT *, "Secondary with Gd residual ", KPART(I), " ene: ", TKI(I) * 1000
                  END DO
                  ! Print out secondary information -- is it around 8 MeV?
                  PRINT *, "Gamma cascade energy ", GAMMA_NEUTRON_ENE * 1000
                  PRINT *, "Number of photons from gad capture: ", NP
                  N_GADCAPS = N_GADCAPS + 1
                  PRINT *, "Number of gad captures: ", N_GADCAPS
               END IF
            END IF
         END IF

The printouts for the gamma cascade energy (Gamma energy below) range from ~4 MeV all the way to 15 MeV:

 Secondary with Gd residual            7  ene:    1.6283083033559327
 Secondary with Gd residual            7  ene:    1.0976518308104020
 Secondary with Gd residual            7  ene:   0.79538810537568339
 Gamma energy    3.52134824
 Number of photons from gad capture:            3
 Number of gad captures:           39
.
.
.
 Secondary with Gd residual            7  ene:    4.4279597575398153
 Secondary with Gd residual            7  ene:    1.4908416220897729
 Secondary with Gd residual            7  ene:    6.0881862942761211
 Secondary with Gd residual            7  ene:    3.6074345911704366
 Gamma energy    15.6144218
 Number of photons from gad capture:            4
 Number of gad captures:          575

I found this reply on the INFN forums for someone asking about a similar issue, though it makes little sense to me:

Neutron interactions below 20 MeV are treated on the basis of neutron cross section evaluations in a group-wise way, except for a few isotopes. This means that also photon emission is treated group-wise, I mean that the gamma energy is “grouped” in intervals, not lines. More than this, there is no correlation among emitted gamma lines.

I think I am scoring the captures correctly so my thinking is that the FLUKA handling of the gamma cascade after a Gd neutron capture.

Any information anyone has on this, whether I am scoring incorrectly or whether the cascades are in fact wrong, is appreciated

Hello Jack,

Your attempted filtering by ICRES was unfortunately ineffective: information on residual nuclei from point-wise non-elastic neutron interactions is stored in the fheavy.inc common (this will be better highlighted in our manual and training material shortly).

The filtering by ICODE.eq.300 (low-energy neutron interaction) is alright, the one over JTRACK.eq.8 (particle originating the interaction was a neutron) should be redundant (but it doesn’t harm). The gamma energy accumulation looks alright.

In order to intercept only neutron capture events on Gd, I’d suggest filtering by ICHTAR.eq.64 (the target was some isotope of Gd) and then loop over the heavies stack (see fheavy.inc) from 1 to NPHEAV, make sure there is only one particle there (it should be the residual) – this should be the signature of a capture event. Then fish out A and Z via IBHEAV(KHEAVY(1)) and ICHEAV(KHEAVY(1)), make sure that they match IBTAR+1 and ICHTAR respectively (then you have the explicit channel in your hands), and do your analysis under this scenario.

Doing this I see events with total energy of gammas of 7.8 MeV, 8.5 MeV, etc as a result of capture on the different isotopes of Gd in its natural composition. One does not stumble upon events with 15 MeV accumulated gamma energy.

Take it for a spin and let us know if issues persist or if you need further hints.

Cheers,

Cesc

1 Like

Hi Cesc,

Thanks for your reply. That has been really helpful and I now see the same as you.

It could be worth noting for future reference that I was also running an older version of FLUKA. I had 4-3.2. This was before the new neutron libraries started to be used by FLUKA and before point-wise interactions were activated for neutrons below 20 MeV. (Edit: Neutron libraries were in FLUKA in this version, though they were not enabled in the most relevant DEFAULTS lists. See Cesc’s reply below.) This could have contributed to what I was seeing as well as my failed filtering. Your method also does not work on the older version of FLUKA and I had to update to the latest release (4-4.1) to get information from ICHTAR, IBTAR, and the FHEAVY commons.

Just in case anyone else wants to implement something similar, the code I use is here:

         !--------------------------------------------------------------
         !===================== NEUTRON CAPTURES =======================
         !--------------------------------------------------------------
         ! Check if the interaction is a low energy neutron interaction
         IF (ICODE .EQ. 300 .AND. JTRACK .EQ. 8) THEN
            ! Check if the incoming particle was a neutron
               IF ((ICHTAR .EQ. 64 .OR. ICHTAR .EQ. 1) .AND.
     &                IBHEAV(KHEAVY(1)) .EQ. IBTAR + 1 .AND. 
     &                   ICHEAV(KHEAVY(1)) .EQ. ICHTAR) THEN
                  ! If we have a residual nuclei with 64 protons
                  !  then we have a Gd capture
                  !  Add the capture position to the neutron write-out
                  !   lists
                  ENE_NEUT (NEUT_CNT + 1) = ETRACK
                  X_NEUT (NEUT_CNT + 1) = XSCO
                  Y_NEUT (NEUT_CNT + 1) = YSCO
                  Z_NEUT (NEUT_CNT + 1) = ZSCO
                  X_DIR_NEUT (NEUT_CNT + 1) = CXTRCK
                  Y_DIR_NEUT (NEUT_CNT + 1) = CYTRCK
                  Z_DIR_NEUT (NEUT_CNT + 1) = CZTRCK
                  T_NEUT (NEUT_CNT + 1) = ATRACK
                  NEUT_CNT = NEUT_CNT + 1
               END IF
         END IF

Like you say, the JTRACK .EQ. 8 check is superfluous and is only there for my sanity. I’m also interested in capture on hydrogen, hence why I also have IF ICHTAR .EQ. 1 in the check.

Thanks for all your help,
Jack

Hello Jack,

Lovely! Great to see that this did the trick.

It could be worth noting for future reference that I was also running an older version of FLUKA. I had 4-3.2.

Jupp, make sure to stay up to date on the latest versions.

NB: the new point-wise low-energy neutron interaction treatment was already included in v4-3.0, but indeed it became default (for the most relevant DEFAULTS) a bit later in v4-4.0.

Sunny cheers,

Cesc

1 Like