Imitate tcquench with userweig card

Dear FLUKA community

As already discussed in a previous post (Non-Linear Electron Scintillation Efficiency in NaI(Tl) Crystals), I’m trying to implement my own weighting function to account for the non-proportional scintillation efficiency in NaI(Tl) crystals using a modified version of Birks law (which is unfortunately not available in the tcquench card). For that purpose, in a first step,I would like to validate my approach by imitating the tcquench card by the userweig card and the corresponding user routine comscw.f.

I’ve already implemented a simple test case with four NaI crystal blocks for both cases, i.e. one simulation with the eventbin + tcquench card and another simulation with the detect + userweig card (and of course the corresponding user routines usreou.f and comscw.f). Unfortuantely, I didn’t get exactly the same results for the both cases (see two figures below). Therefore, I’m wondering, if I have made an error somewhere. I have the following specific questions for you regarding the comscw.f user routine and the implementation of Birks law in FLUKA:

  1. I’ve implemented Birks law considering only the continuous energy deposition events by electrons in my NaI crystals (Ntrack = 1, Mtrack = 1,Jtrack=3),.i.e.:

IF (( Mtrack .EQ. ONEONE ) .AND. ( Ntrack .EQ. ONEONE ) .AND.
& ( Jtrack .EQ. THRTHR )) THEN

COMSCW = ONEONE / ( ONEONE + kB * Dtrack(1) / Ttrack(1))

ENDIF

More details can be found directly in the attached comscw.f routine. Note, that I’m using a Birks coefficient of 10^-2 [g MeV^-1 cm^-2] for all simulations (chosen arbitrarily for test purposes).

Is this approach correct? Does FLUKA consider also local deposition events (Ntrack = 0, Mtrack = 0) for the implementation of Birks law?

  1. I’ve recognized that the Ttrack and Dtrack variables contain only non-zero variables in the first column of the arrays (for Ntrack = 1, Mtrack = 1). Therefore, I’m only considering Ttrack(1) and Dtrack(1) for my computations. Is this correct?

  2. I’m assuming that the units of Ttrack and Dtrack are default, i.e. cm and GeV.

I have to add that I’ve validated both standard pulse height cards (i.e. eventbin and detect+usreou.f card) and they give exactly the same results. So, the error has to be somewhere in the comscw.f routine.

For your convenience, I’ve attached the input files as well as the corresponding user routines to this post.

Thank you so much for your support

Best,
David

detect_userweig.flair (6.5 KB)
detect_userweig.inp (6.6 KB)
comscw_ed.f (6.7 KB)
usreou.f (2.8 KB)
detect_userweig.pdf (77.5 KB)

eventbin_tcquench.flair (6.8 KB)
eventbin_tcquench.flair (6.8 KB)
eventbin_tcquench.pdf (76.7 KB)

David @BV96, are you actually using the comscw version you uploaded? If so, there is a mistake in your Birks coefficient, which is assigned to the integer variable kB rather than the intended bK.
Note that in FLUKA and its user routines the systematic inclusion of the dblprc.inc file carries automatically the statement

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

making your REAL definition useless (actually spoiling the double precision requirement) and kB truncated to the lower integer value.

As a side note, not affecting the outcome and answering your question 2, NTRACK is normally 1 (not ONEONE, which is a double precision parameter, opposite to integer) in the absence of magnetic/electric field.

I’ll come back on other points (e.g. local deposition events) after your feedback on the above correction.

Dear Francesco

Thank you very much for your reply. I really appreciate that. Regarding your first inputs, I’ve modified the comscw.f routine as suggested by you, i.e. I’ve deleted the REAL definition but for now kept the ONEONE value for the NTRACK variable (could also integer values like 1 and 3 be adopted for the MTRACK and JTRACK variables?).

Ater that, I’ve compiled and run the detect+comscw simulation again. You can find the modified comscw.f routine as well as the input scripts in the attachment to this post. Unfortunately, there is still a significant difference between the detect+comscw and the eventbin+tcquench simulation (cf. pdf figures)…

Do you have any idea why? Might it be an error in the unit conversion or some neglected physical effect, which is considered in the tcquench implementation in FLUKA (as already described, maybe the local energy deposition events or delta rays?).

Thank you very much for your support
Best,
David

comscw_ed.f (6.7 KB)
detect_userweig_new.inp (6.6 KB)
detect_userweig_new.flair (6.4 KB)
detect_userweig_new.pdf (77.0 KB)

Hi David, you did not correct at all the mistake: you keep assigning the Birks coefficient to the kB variable, which by default is an integer one (since its name starts with k), so its value gets transformed into an integer. Please replace everywhere kB with bK, as you initially meant to call the variable (you had defined bK as real, which is useless, but then used kB).

Dear Francesco

Now I think I got it. I was not aware that the implicit definition of variables is sensitive to the order in the characters. I’ve exchanged all kB with bK in a new simulation. The simulation results get closer but there is still a significant difference. Any ideas why?

Best,
David

comscw_ed.f (6.7 KB)
detect_userweig_new2.flair (6.4 KB)
detect_userweig2.inp (6.6 KB)
detect_userweig_new2.pdf (77.8 KB)

The reason for the remaining discrepancy concerns indeed the local energy deposition events by electrons below threshold or photons, which are also quenched when using TCQUENCH.
This is a comscw version I modified to account for them as well, with some explanatory comments. Now it should perfectly match the TCQUENCH effect in EVENTBIN, and is prone to further customization of your own.

Note that I’ve also corrected your stopping power calculation by using CTRACK (actual curved path) instead of TTRACK.

Moreover, note that your EMFCUT card is ineffective, since you did not specify the geometry regions where it’s expected to apply.

Dear Francesco

Thank you so much for your help and your support. I really appreciate that! I’ve tried your modified comscw routine. I have only one small problem: Sometimes, the simulation stops because an “unexpected track” was found, i.e. the “else” case in your routine got triggered.

Therefore, I’ve tried to add also positrons to the local deposition events as well as to the continous deposition events (because I’ve detected some events with JTrack = 4 and we actually expect to see some positrons by pair production). I don’t know if this is physical at all, because I expect positrons to get annihilated very fast. This modification made the routine more stable (i.e. less stops), but sometimes the simulation still doesn’t work because of the comscw routine. Do you have any idea why (other particles deposit energy in the crystal?)? Did my input scripts work with your routine in your case?

Side question: What exactly are you doing with this line in your code:

     ELKE = LOG ( MAX ( 0.6D+00 * ( EEPMIN (MEDEMF (MREG,IPRODC))
 &                                  - AMEMEV ), ENKEMF ) )

I didn’t get that. I guess you are transforming the kinetic energy of the electron in some way to be adaptable to the look-up table…

Best,
David

comscw_Cerutti_ed.f (8.2 KB)

Sorry, I was somewhat sloppy and did not accumulate enough statistics as to see positrons, which are perfectly legitimate, as you say, and are (physically) tracked until annihilation.
This adjusted comscw version differs from yours, since for positrons the stopping power tabulation is another, and happily survived one million histories (I used your input with the aforementioned EMFCUT correction). Should any additional track type (i.e. JTRACK value not anticipated here) trigger again an abort, thanks in advance for reporting.

Note that the second index of the stopping power tabulation refers to the material, being the dependence on the latter not limited to the average ionization potential.

Stopping power tabulations feature a logarithmic energy grid, with a lower limit to be taken into account.

Dear Francesco

Thank you very much for your help. I will test the modified comscw.f including positrons in the next days and will come back, if I detect some more issues. I expect that for some more exotic applications (e.g. cosmic ray simulations), additional charged particles have to be considered (e.g. muons, protons,…). That’s not in the scope of my work, but I plan to update the modified comscw routine with the most common charged particles to have a nice template also for other FLUKA users.

Based on your template that should be fairly easy, i.e. I guess I have only to add some more lines with the corresponding stopping power functions for those particles. Could you give me a hint, for which module (*.inc file) I have to look to find the stopping power tables/functions for other charged particles than electrons and positrons?

Best
David

Well, actually this is not so straightforward and such an extension would immediately call, as far as local energy deposition events are concerned, for dealing with recoiling nuclei from nuclear interactions (falling below transport threshold). With regard to stopping power retrieval, a more transparent way to consider is invoking in the user code the

      DOUBLE PRECISION FUNCTION GETLET ( IJ, EKIN, PLA, TDELTA, MATLET )

where

*         Input variables:                                             *
*                                                                      *
*               Ij = particle index (Paprop)                           *
*             Ekin = particle kinetic energy (GeV)                     *
*              Pla = particle momentum (GeV/c)                         *
*           Tdelta = maximum secondary electron energy (GeV)           *
*                   (unrestricted if =< 0)                             *
*           Matlet = material index for which LET is requested         *
*                                                                      *
*         Output variables:                                            *
*                                                                      *
*           Getlet = (un)restricted LET (keV/(um g/cm3))               *

with attention to be paid to the delta ray (i.e. secondary electron) threshold.

Dear Francesco

Based on your input, I’ve decided not to extend the routine to more exotic particles. During some simulations, I’ve encountered a small issue, i.e. for higher EMF transport thresholds, local deposition cases of photons (JTrack=7) are possible and have to be considered therefore in the subroutine. Consequently, I’ve included the photons for the local deposition events (see routine below). However, I have still three small questions for you:

  1. I wonder, if it would be justified to alter only the electrons and positron events with the userweig routine (i.e. keeping local photon deposition events unaltered). My concern is that local energy deposition events by photons with a high energy, i.e. not caused by the EMF threshold, aren’t well described by electron response curves.

  2. Is it possible to submit somehow custom numerical input parameters (like for example the Birks coefficient) from the input file to the (compiled) userweig routine (maybe with another routine)?

  3. Is the initialization and definition of variables/parameters with the PARAMETER statement (as used in the routine below) a good idea in a FLUKA routine?

Thank you so much for your support
Best,
David

comscw_BirkDefault.f (12.8 KB)

Dear David,
nice and clean work! [As a quite secondary correction, you may want to replace my own postal address from LinkedIn with CERN’s one].

  1. I’d be inclined to rather lower the thresholds to increase accuracy.
  2. Very good point. For the source user routine, not relevant here, this can be directly done by means of the respective SOURCE card. This does not apply to USERWEIG. Nevertheless, either the USRGCALL and USRICALL cards allow one to access input parameters in the usrglo and usrini routines, respectively (see the manual), and, after assigning them there to own custom variables, share the latter ones with other routines, such as comscw, through COMMON instructions (as you can find in the .inc files located in the include directory) to be included in all relevant routines. Evidently, the advantage of such an elegant practice is the possibility to change the parameter values in the input file with no need to recompile and regenerate the executable.
  3. I do not see any explicit PARAMETER statement in your routine. These can certainly be added on the top to set unchangeable elements (unless the above practice fits better). For variables, which opposite to the PARAMETER type can change their value runtime, the DATA statement allows initialization.

Dear Francesco

Thank you so much for your reply and your support. Regarding the input variable task, I’ve followed your approach and tried to implement a USRGCALL card together with the corresponding usrglo and comscw user routines. Unfortunately, due to limited/inexistent FORTRAN programming skills, I was not successful. I’ve attached the routines with my trials (some of which are commented). Could you maybe give me a hint, where I made errors (I guess my variable declarations in the usrglo routine are not correct).

Best
David

P.S. I didn’t find your CERN credentials online. Could you maybe send them to me. Alternatively, you could also fill in your data directly into the routine and send it back to me.

usrglo_ed.f (1.6 KB)
comscw_BirkDefault_usrglo.f (12.9 KB)
test_4box_detectweight_usrglo.flair (7.6 KB)
test_4box_detectweight_usrglo.inp (6.7 KB)

Here we go:
usrglo_Birks.f (1.3 KB)
comscw_Birks.f (12.8 KB)

Dear Francesco

Thank you so much. I will come back, if I find some more issues/problems, which exceed my capabilities.

Best
David