Track secondary neutron in (n,2n) reaction

Advanced Features and User Routines
1

Dear Fluka experts,

I want to use a polyethlene moderated 3He tube to detect emitted neutrons from a Au sample. Attached figure shows my detector’s configuration and you can also find input file and mdstck.f and mgdraw.f in the attachment files.

I used mdstck.f to record the interaction information of neutrons in the target region, and flag the neutron elastic scattering, neutron inelastic scattering, and (n, 2n) by 0,1,2 respectively in first column of the generated file ctry_001_3He_detector_mdstck.dat. I also used mgdraw.f at the detector position to record the detected neutron information and generate the ctry_001_3He_detector_mgdraw.dat file.

I don’t know how to track the secondary neutron produced by the (n,2n) reaction and if FLUKA exists some parameter like trackID in G4? Could you tell me how can I flag neutron reaction mark 0,1,2 in the ctry_001_3He_detector_mgdraw.dat file? In addition, I found NPFLKA doesn’t have the same value in sturpf.f and mdstck.f.

Best regards,

Jianqi

fig1_detector_model
fig1_detector_model828×618 41.4 KB

file1_ctry4.inp (4.2 KB)
mgdraw.f (15.2 KB)
mdstck.f (4.7 KB)

3

Hi Jianqui,

I will have a look into your problem and come back with an answer in a few days.
Thanks for providing the files!

Cheers,
Jerzy

5

Thank you, Looking forward to your answer.

6

Hi Jianqi,

First of all, sorry for the late reply, it took me a while to go through your request.
From your code I can see that you are quite an advanced FLUKA user already. I am sure that you would figure out the solution yourself, but I leave you some advice below.

Let me start from the end, I think that the difference in NPFLKA comes from the fact that between the place where these two functions are invoked in the main fluka code, some kind of main stack management is happening.

Even though among the files that you attached there is not sturpf.f, from what I understand, you are well aware of this user routine. If you want to have a variable which acts similarly to G4 trackID, you should in fact use one of the available user variable and assign them inside the sturpf.f. You could try to use a similar logic that you have in your mdstck.f function to flag the neutrons and then print out the flags when the neutrons get detected in your USDRAW subroutine of ‘mgdraw.f’. I am attaching a version of stupfre.f that should do the same thing as your mdstck.f and also flag the interaction type origin to all the produced particles according to your logic.
stuprf.f (3.1 KB)
this routine can be used together with your mgdraw.f where in your USDRAW you can print just one more variable LLOUSE

                     WRITE(10,*) NCASE,ATRACK,JTRACK,XSCO,YSCO,ZSCO,
     &                           NAM,NP,KPART(1:NP),TKI(1:NP), LLOUSE

to identify the origin of your detected neutrons.
When you use the above stupr.f, DO NOT include your mdstck.f to your executable.

There are a few small things in your code that I think don’t act as you would like them to.
Firstly, with such code

                     WRITE(7,*)  2, NCASE, XTRACK(NTRACK), YTRACK(NTRACK),
     &                           ZTRACK(NTRACK),CXR(NTRACK),CYR(NTRACK),
     &                           CZR(NTRACK),KPART(IPNEW), TKI(IPNEW)

the cosines are always taken from the parent particle and they end up to be all the same for all the daughters (interaction products). Instead of using CZR(NTRACK) you should rather take CZR(IPNEW), see the genstk variables in fluka_dir/include/genstk.inc.

Secondly, when you put conditions for the elastic scattering,

            ELSE IF ( NCASE .EQ. TEMP_INL(IPNCASE).AND. NP.EQ.1) THEN
*     &              .AND.JTRACK .EQ.8) THEN
              EKIN = ETRACK-AM(JTRACK)
*              PLOSS= DPTRCK(1,NTRACK)**2+
              TXX= (PTRACK*CXTRCK-PLR(NTRACK)*CXR(NTRACK))/PTRACK
              TYY= (PTRACK*CYTRCK-PLR(NTRACK)*CYR(NTRACK))/PTRACK
              TZZ= (PTRACK*CZTRCK-PLR(NTRACK)*CZR(NTRACK))/PTRACK
              
              PRINT *,"DTRACK=", DTRACK
              PRINT *,"DPTRCK=", DPTRCK

              DO IPNEW = 1, NPSECN
*                  IF(KPART(IPNEW) .EQ. 8) THEN     
                    WRITE(7,*)  0, NCASE, XTRACK(NTRACK), YTRACK(NTRACK),
     &                           ZTRACK(NTRACK),TXX, TYY,
     &                           TZZ, 8 , EKIN-TKI(IPNEW)
*                  END IF
              END DO
            END IF

you are printing the information about the only particle on the stack and you assume it’s a neutron, but in reality it is always a gold residual nucleus -7919707. You should maybe rethink the way you want to tag the interaction types that you are interested in.

I hope that helps, let me know.

Cheers,
Jerzy

7

Hello Jerzy,

Thank you very much for your very informative and inspiring reply. Tracking neutron from (n, 2n) reaction in FLUKA haunts me for a long time. Your answer helps me a lot.

In my post, I didn’t add stuprf.f, because I didn’t know how to use it properly at that time. One week ago, I found one way to track neutron flag through mdstck.f(attached_file_1) and stuprf.f(attached_file_2). I used stuprf.f to record the parent NCASE and energy info, and mdstck.f to record the neutron reaction info in the target and 3He region, so I can track the neutron reaction flag through tracking the energy of parent particle. Here I didn’t use one 3He tube but 110 3He tubes for real model, so I also add 3He tube position info in another file helium3_tube_110_pos_new.dat(attached_file_3). The principle behind it is same. Compared with your method, it appears dumb, but it can still solve my problem by an indirect way.

Let me first answer your question about tagging the neutron elastic scattering.

 ELSE IF ( NCASE .EQ. TEMP_INL(IPNCASE).AND. NP.EQ.1) THEN
*     &              .AND.JTRACK .EQ.8) THEN
              EKIN = ETRACK-AM(JTRACK)
*              PLOSS= DPTRCK(1,NTRACK)**2+
              TXX= (PTRACK*CXTRCK-PLR(NTRACK)*CXR(NTRACK))/PTRACK
              TYY= (PTRACK*CYTRCK-PLR(NTRACK)*CYR(NTRACK))/PTRACK
              TZZ= (PTRACK*CZTRCK-PLR(NTRACK)*CZR(NTRACK))/PTRACK
              
              PRINT *,"DTRACK=", DTRACK
              PRINT *,"DPTRCK=", DPTRCK

              DO IPNEW = 1, NPSECN
*                  IF(KPART(IPNEW) .EQ. 8) THEN     
                    WRITE(7,*)  0, NCASE, XTRACK(NTRACK), YTRACK(NTRACK),
     &                           ZTRACK(NTRACK),TXX, TYY,
     &                           TZZ, 8 , EKIN-TKI(IPNEW)
*                  END IF
              END DO
            END IF

You are right, CXR(NTRACK), CYR(NTRACK), and CZR(NTRACK) they are recording the emitting direction info of gold residual nucleus -7919707. In order to obtain the info of neutron direction, I used the conservation law of momentum to calculate them. In the meanwhile, I used the conservation law of energy to calculate the neutron energy. Here I have to say, the calculation process of emitting angle along different axis is not accurate, because the momentum of incident neutron is unknown(I don’t know how to get it), I was using the momentum of outgoing neutron(PTRACK) to substitute it, since the module of neutron momentum doesn’t change too much for neutron elastic reaction.

I have already tried your script, It works very well. Very impressive!!! I never thought it can be done like this way. But I still have one question. For some neutrons, they will travel the target directly without any reaction. So these neutrons info will not be recorded by stuprf.f, but they still have the possibility of being reflected by wall and then be detected by 3He, I don’t know how to tag this type of neutron, do you have any ideas for this condition?

Cheers,

Jianqi

file1_mdstck.f (7.2 KB)
file2_stuprf.f (2.2 KB)
file3_helium3_tube_110_pos_new.dat (2.1 KB)

8

Hello Jerzy,

Recently, I found one way that maybe capable of tracking neutron reaction based on your instructive and enlightening advice and routine lecture from this link at page 17 https://indico.cern.ch/event/1200922/contributions/5411837/attachments/2661810/4611453/08_User_routines_2023_Advanced_ANL.pdf.

Attached files marked with suffix “new” are my new codes. Here I was using ISPUSR(1) instead of LLOUSE in stuprf.f to record the neutron reaction info, and then using following loop method given by the above lecture to record the reaction parent info.

     DO 200 ISPR = 1, MKBMX2-1
                 ISPARK (ISPR+1,NPFLKA) = ISPUSR (ISPR)
 200     CONTINUE

After that, I was using mgdraw.f card to record the detection info with the neutron reaction info like this.

                     WRITE(10,*) NCASE,ATRACK,JTRACK,XSCO,YSCO,ZSCO,
     &               NAM,NP,KPART(1:NP),TKI(1:NP),ISPUSR(1),
     &               ISPUSR(2), ISPUSR(3), ISPUSR(4), ISPUSR(5)

So In principle, I can track the reaction by this method. But what confused me is some result given by this method doesn’t coincide with the previous method by tracking the parent neutron energy(marked with suffix old). I don’t know the reason behind it.

For example, when ncase is 538, the result (ctry4001_3He_detector_mgdraw_new.dat) given by the new method is

         538   2.2177272973285965E-005           8   14.320330604512852       0.68861043573417624        3.1649887574694944      c421               2           1           7   5.7251164476090910E-004   3.0181749290870368E-008           0           0           0           0           0

It seems no reaction happened inside sample.However, if we are using the old method of tracking parent neutron energy, the result given by the old method(ctry4001_3He_detector_mdstck_old.dat) is

           3           0   0.0000000000000000              538           8   1.4000000000000012E-002           1       36650           0           0  -2.8656763135969484E-002  0.10917639650257346        1.4472453757208876                8   5.1329621813515879E-003   3.8044877155283707E-010           0
           3           0   0.0000000000000000              538           8   1.4000000000000012E-002           1       36650           0           0  -2.8656763135969484E-002  0.10917639650257346        1.4472453757208876                8   7.0144552037515950E-004   3.8044877155283707E-010           0
           1         538   5.1329621813515879E-003         538           8   5.1329621813515879E-003           2       36652       53995           2  0.20805753887368736       0.25861619098946625        1.5068571066381633                8   5.1321995024180174E-003   4.7215761207757894E-010           0
           8         538   7.0144552037515950E-004         538           8   1.7421530474826023E-011          17       36699       54256          16   14.320330604512852       0.68861043573417624        3.1649887574694944               -4   1.9124064314324229E-004   2.2177272973285965E-005          21

It is easy to see that the (n,2n) reaction happened inside the sample.

So I am wondering if the new method works? Do you have some suggestions on it?

In addition, for the old method of tracking neutron parent energy, I used parameters like LOFLK, NUMPAR, NEVENT, and NPFLKA in stuprf.f, but I don’t know exact meaning of these parameters,
their interpretation inside flkstk is too simple and I don’t get the meaning of “particle generation” for LOFLK, and “particle number” for Numpar? “number of the event which created the particle” for Nevent? Would you please tell me more about them? When the ncase is 1, the result is given like this, I don’t know if these parameters can give us more hint to track the neutron?

   #reaction_type   ncase_parent energy_parent             ncase       jtrack     eptrack                     loflk      numpar    nevent        npflka   xtrack                     ytrack                     ztrack                  kpart(ip)            tki(ip)                 atrack 

           1           0   0.0000000000000000                1           8   1.4000000000000012E-002           1           1           0           0 -0.43092745688245293       0.61296445287527879       0.56656761436543246                8   1.3997189481580161E-002   2.0838377316117522E-010           0
           3           0   0.0000000000000000                1           8   1.3997189049446468E-002           2           2           1           1 -0.29207346433233183       0.81322272233927917        1.7838052862874818                8   1.5945015528584603E-003   4.5094797480441392E-010           0
           3           0   0.0000000000000000                1           8   1.3997189049446468E-002           2           2           1           1 -0.29207346433233183       0.81322272233927917        1.7838052862874818                8   1.2682917708824243E-003   4.5094797480441392E-010           0
           2           1   1.5945015528584603E-003           1           8   1.5945015528584603E-003           3           8           2           7 -0.23022193428593715       0.83887063882949031        1.8182507750976640                8   1.3157130034544551E-003   4.9411538951187750E-010           0
           8           1   1.3157130034544551E-003           1           8   3.4041937200804006E-011          59          38          58          35   10.023958647799379        3.6233448669981558        10.255316895700618               -4   1.9123024299272373E-004   6.0345241694901193E-005           2
           8           1   1.2682917708824243E-003           1           8   1.3708262497929762E-010         139          58         204          22  -5.8292830237726037        18.802277055506185        3.4813708765256806               -4   1.9120411559558903E-004   1.4806481390413736E-004          47

Best Regards,

Jianqi

stuprf_new.f (3.2 KB)
mgdraw_new.f (18.4 KB)
ctry4001_3He_detector_stuprf_new.dat (546.4 KB)
ctry4001_3He_detector_mgdraw_new.dat (46.1 KB)

stuprf_old.f (2.2 KB)
mgdraw_old.f (18.3 KB)
mdstck_old.f (6.9 KB)
file3_helium3_tube_110_pos_new.dat (2.1 KB)
ctry4001_3He_detector_mgdraw_old.dat (36.4 KB)
ctry4001_3He_detector_mdstck_old.dat (328.2 KB)

ctry4.inp (94.6 KB)

9

Sorry, I forgot to mark that the neutron reaction type in the new code has been changed. Now elastic, inelastic and (n, 2n) reaction is marked as 1, 2, 3, respectively rather than 0,1,2 in order to avoid confusion when the value of ISPUSR is equal to 0.

Best regards,

Jianqi

10

Hi Jianqi,

I can see that you are advancing with your search for a solution, seems like you are almost there. Sorry for the lack of reply, but I was on holidays and I will be able to have a look again at your questions next week.

Best regards,
Jerzy

11

Hi Jerzy,

Thank you for your encouragement. I am still working on it.

No worries, have a good vacation, Bon voyage!

Best regards,
Jianqi

12

Hi Jianqi,

I am slowly trying to unfold all the things that you were trying to do. Let’s maybe simplify the problem.

First of all, trying to save ncase_parent and ncase does not make much sense, because a given NCASE corresponds to a single primary or in other words to a single event.

I would say that saving LOFLK, NUMPAR, NEVENT, or NPFLKA will not give us much useful information. The latter is actually a pointer which points to the current particle on the stack. LOFLK is a stack variable which most likely corresponds to LTRACK of a given track, but I guess it’s better that we stick to LTRACK if you want to save it. With LTRACK == 1 you can filter the primary particles and this fact brings me to a few questions:

  1. Which particles/neutrons you actually want to track? Only the ones which originate from the primary interaction?
  2. What is the information that want to save? I understood that your goal was to flag the primary interaction inside the target and then when you save the detected neutrons you are able to tell what was the first interaction type. Is it enough for you? I am asking because I can see that in the _new version you are saving the interaction type of the last 5 generations, but why is that you are interested in this information?

Another thing is that maybe there is a way to simplify the tagging of neutron elastic interaction. I must say I don’t quite understand what is the interaction type that you flag as 8 in your mdstck_old.f.
I have to think about it, but there should be a way to use the mgdraw main routine to check whether a particle is still a primary and if it’s direction was modified with respect to the previous step. This would easily work for primaries, but if you also want to track elastic interactions of secondary neutrons, probably this could still be done with the usage of LTRACK or by assigning a unique track ID when the particle is loaded into the stack (in stuprf) or when it’s generated (in USDRAW). I tried to do it but for the moment it doesn’t work as expected, I need a bit more time to check how the things work.
Stay tuned!

Cheers,
Jerzy

13

Hi Jianqi,

Please try to use the following mgdraw routine:
mgdraw_dir_elastic.f (15.2 KB)
You can modify there your detector or target conditions as you like, I was still using your old FLUKA files for testing.

This approach uses only USDRAW to tag the very first interaction of the primary neutrons. Then, in the detector file you can identify the interaction type and the NCASE.
In the case of elastic interaction,

*        If the projectile is a primary neutron:  
         IF (JTRACK.EQ.8 .AND. LTRACK.EQ.1) THEN
*		 If interaction inside the target				
				  IF(NAM.EQ."target") THEN
*				  If elastic interaction
					IF(NP .EQ. 1) THEN
						ISPUSR(1) = 1
						DO IPNEW = 1, NP
						  WRITE(20, *) ISPUSR(1), NCASE, KPART(IPNEW), TKI(IPNEW),  CXR(IPNEW), CYR(IPNEW), CZR(IPNEW)
						END DO

The only particle on the stack turns out to be the scattered neutron, so you can directly retrieve all the information about it using the corresponding genstk variables.
For inelastic interactions:

*		If interacton inelastic				
					ELSE IF(NP .GT. 1) THEN
						DO IPNEW = 1, NP
						  IF(KPART(IPNEW) .EQ. 8) THEN
							NNEUT = NNEUT + 1
						  END IF
						END DO

						IF(NNEUT .EQ. 2) THEN
*		If Nto2N
							ISPUSR(1) = 3
*		Else inelastic						
						ELSE 
							ISPUSR(1) = 2
						END IF
						DO IPNEW = 1, NP
						  WRITE(20, *) ISPUSR(1), NCASE, KPART(IPNEW), TKI(IPNEW),  CXR(IPNEW), CYR(IPNEW), CZR(IPNEW)
						END DO
					END IF
				  END IF
		 END IF

the script checks the number of neutrons pushed to the stack by the registered interaction to set a flag either inelastic or nTo2n.
If you want to have the history of your interaction chain, you could still set up a stuprf routine where the only modification would be the one that you already mentioned:

DO 200 ISPR = 1, MKBMX2-1
   ISPARK (ISPR+1,NPFLKA) = ISPUSR (ISPR)
200   CONTINUE

In this case however, you need to lift the requirement for the primary particle LTRACK .EQ. 1 and modify the mgdraw code so that it only sets the ISPUSR(1) in the case of secondary neutrons without actually writing them to file (unless you them to be written).

If you still want the information about the gold nuclei residuals, I think it can be done with the old method by extending the stuprf.f routine for example. Possibly the residuals can also be accessed via the proper stack in USRDRAW.

Hope this helps,
Cheers
Jerzy

14

Hello Jerzy,

Thank you very much for patiently helping me solve this problem.

Let me answer your questions first.

Blockquote 1. Which particles/neutrons you actually want to track? Only the ones which originate from the primary interaction?

  1. I want to track neutron but only for the primary interaction, I want to track the whole reaction process of neutron inside the sample and also need to know which emitted neutron is captured by detector.

Blockquote 2. What is the information that want to save? I understood that your goal was to flag the primary interaction inside the target and then when you save the detected neutrons you are able to tell what was the first interaction type. Is it enough for you? I am asking because I can see that in the _new version you are saving the interaction type of the last 5 generations, but why is that you are interested in this information?

  1. The information I wanted to record is the reaction process of detected neutron inside the sample. So I pay attention not only the whole process of neutron propagating inside the sample but also need to know which emitted neutron is finally detected by detector.

Blockquote Another thing is that maybe there is a way to simplify the tagging of neutron elastic interaction. I must say I don’t quite understand what is the interaction type that you flag as 8 in your mdstck_old.f.

* conversation of energy
              EKIN = ETRACK-AM(JTRACK)
* conversation of momentum
              TXX= (PTRACK*CXTRCK-PLR(NTRACK)*CXR(NTRACK))/PTRACK 
              TYY= (PTRACK*CYTRCK-PLR(NTRACK)*CYR(NTRACK))/PTRACK
              TZZ= (PTRACK*CZTRCK-PLR(NTRACK)*CZR(NTRACK))/PTRACK

              DO IPNEW = 1, NPSECN
                    WRITE(7,*)  1, ISPUSR(1), SPAUSR(1),
     &                NCASE, JTRACK, EPTRCK,
     &                LOFLK(NPFLKA),NUMPAR(NPFLKA), NEVENT(NPFLKA),
     &                NPFLKA, XTRACK(NTRACK), YTRACK(NTRACK),
     &                ZTRACK(NTRACK),8,EKIN-TKI(IPNEW),ATRACK,0
               END DO
            END IF
         END DO

I guess you mean this part of the code. Here I am also using 8 to denote neutron, just because for elastic scattering reaction, there is no neutron produced in residual nuclei, but the information I want to record is the scattered neutron, so I obtained the neutron energy and direction information through the conversation of energy and momentum, respectively.

Best regards,

Jianqi

15

Very clever method, thank you very much!

I will try it.

Cheers,

Jianqi