Non-Linear Electron Scintillation Efficiency in NaI(Tl) Crystals

Dear FLUKA community

As already described in previous posts, I’m currently working on simulating the gamma-ray spectral detector response for NaI(Tl) detectors using an event-by-event approach. For that purpose, I’m adopting the detect card together with the usreou.f user routine (cf. https://fluka-forum.web.cern.ch/t/prompt-decay-radiation-from-isotopes/1479; Energy deposition: DETECT or MGDRAW). So far, everything is working perfectly and I got already good results in benchmark comparisons with experimental data.

There is only one small issue, i.e. a systematic bias in the valley between the Compton edge and the full energy peak/photopeak in the simulated spectra. This bias is expected and can be linked to the neglected nonlinear electron scintillation efficiency [1-4] in my approach, i.e. I’m scoring the full energy deposited by gamma-rays and electrons per event and to my knowledge, FLUKA doesn’t account for scintillation efficiency by default . Therefore, I’m wondering, if it is possible to incorporate the nonlinear electron scintillation efficiency in a FLUKA simulation? At the moment, I have the following two ideas:

  1. Offline approach: Try to obtain the deposited energy for each electron interaction in the scintillation crystal during an event by a scoring user routine (maybe also usreou.f?). Afterwards, postprocess this electron data by the proper scintillation efficiency function, e.g. the one published by [4]. This offline approach has two disadvantages. First, I guess it will produce huge files. Second, I’m not sure if it is possible to obtain the electron data at all.

  2. Online approach: Try to change the transport physics in the FLUKA simulation by adopting again a user routine, i.e. incorporate the nonlinear electron scintillation efficiency function in the code itself. Unfortunately, I have no idea if this is possible and if yes how to do that.

What do you think? Do you have any other idea or did any of you already implemented something similar? Or is there a standard FLUKA card, which can model scintillation efficiency functions?

Thanks for your support

Best

David

[1] doi10.1016/0029-554X(81)91225-8

[2] doi10.1016/S0969-8043(02)00140-9

[3] doi10.1016/J.APRADISO.2011.12.006

[4] doi10.1016/J.APRADISO.2011.11.003

Dear David @BV96 ,
there is indeed the dedicated TCQUENCH card, which I already advertised for a different purpose at the the end of the first thread you recalled. It implements the more convenient online approach, as in the manual. Nevertheless, it’s coupled to the EVENTBIN (and USRBIN) scoring, and not to DETECT.

Dear Francesco

Thank you very much for your reply. As suggested by you, I’ve checked the TCQUENCH card. According to the Fluka manual, the TCQUENCH card provides calibration constants, time cutoff parameters or Birk’s law parameters. I guess you are referring to Birk’s law coefficents. I was checking this topic in some textbooks and I’m actually not sure, if Birk’s law is equivalent to my problem. To clarify, I repeat here my goal/problem in more details:

According to some publications already cited in the previous post, there is a non-linear scintillation efficiency present for electrons in inorganic scintillators (Nai(Tl), BGO,…), i.e. the number of scintillation photons is not proportional/linear to the deposited energy of the electrons. Therefore, to correct for that issue, one has to multiply each deposited energy by electrons (events, where electron energy goes below the energy threshold (local) as well as continuous energy loss) by the corresponding scintillation efficiency:

E_detector(E_e_deposited) = ScintillationEfficiency(E_e_deposited) * E_e_deposited

For the ScintillationEfficiency(E_e_deposited), there exist empirical functions (cf. [4] from above). I’m also attaching a figure from this paper, which shows the general shape of the ScintillationEfficiency(E_e_deposited). To obtain the total energy deposited per history/primary (E_event), one has to sum all electron interactions in the detector crystal volume plus some other energy deposition events by photons:

E_event= sum (E_detector) + …

This E_event is provided by the detect (+usreou.f routine) or eventbin card, but of course without considering the effect described above, i.e. ScintillationEfficiency(E_e_deposited) =1.

So, I have the following questions for you:

  1. Is the TCQUENCH card able to account for the described effect? I’m not entirely sure, if Birk’s law is the same as the effect described above. Birk’s law seems to be related to quenching, which seems to be only relevant for ions at high energies. In addition, Birk’s law seems to cover only continuous energy loss and no local energy deposition events by electrons. Or am I wrong on this? To be honest, I do not see a way to enter an energy dependent scintillation efficiency in the TCQUENCH card.
  2. Is there maybe another way than the TCQUENCH card. In the meantime, I’ve found this discussion on the forum: Charge trapping using comscw.f routine. It’s covering a similar topic, but unfortunately, no final solution was provided.

Thank you very much for your support. I really appreciate that!
Cheers,
David

I admit that I did not manage to find your ref 4 and I was indeed referring to Birks law quenching of the energy deposited in a charged particle step, by means of the factor (1 + BS + CS**2) dependent on the particle stopping power S = dE/dx (that is energy dependent in turn), where B and C are the Birks parameters to be input [J.B. Birks, The theory and practice of scintillation counting, Pergamon Press, Oxford 1964].
One can also implement any efficiency function by means of the comscw user routine, to be activated by the addition of the USERWEIG card in the input (see the manual). The provided template just contains the two instructions:

      LSCZER = .FALSE.
      COMSCW = ONEONE

that mean ‘do NOT discard the present scoring contribution’ and ‘multiply the present scoring contribution by 1 (i.e. do not alter it)’, respectively. You can implement there your scintillation efficiency function and assign the COMSCW variable the corresponding value, which will be used right after by FLUKA to multiply the deposited energy, as scored either by EVENTBIN and DETECT. In order to implement the aforementioned function, several relevant variables are available, such as JTRACK (or J0TRCK, read the include/trackr.inc header) - integer giving the particle type - and ETRACK - double precision real number giving the particle total energy (mass + kinetic energy) -.

Dear Francesco

After some more reading, I have to correct my previous statement. Birk’s law can actually be used to simulate also the non-linear scintillation effect in inorganic scintillators. There is only one problem. The relative light yield response curve for NaI(Tl) possesses a distinct peak (see figure below). Therefore, most of the studies apply a modified Birk’s law, i.e. Birk’s law +Onsager exponential term (cf. Payne et al., 2009, Nonproportionality of Scintillator Detectors:Theory and Experiment). The modified version of Birks law, which is implemented in FLUKA, is unfortunately not capable to describe the shape of the response curve for halide scintillators (the second term was introduced by Crau and Smith to get better results for organic scintillators). Therefore, I wonder how I can implement my own modified version of Birks law in FLUKA.

You suggested to use USERWEIG card for that purpose. I’m wondering, is there a USERWEIG template for the default Birks law quenching correction used by FLUKA? Unfortunately, I’m a beginner in FORTRAN coding and I’m not that familiar with the user routines. Therefore, I try to explain in the following lines, how I would implement this custom weighting function with the USERWEIG routine:

  1. First, I have to access the energy deposited in a charged particle step, let’s say “dE”. This, I can do by the variable “Dtrack” from the trackr module. But do I have to filter this array to get some specific deposition events (only charged particles?) or the events from a specific region, i.e. the detector region? I guess, this variable will contain all particle tracks from the entire geometry and event, because the USERWEIG is called at the end of an event, won’t it?
  2. Then I have to weight this deposited energies by the term provided by the (modified) Birks law to get the corrected deposited energy, i.e. dE’ = f(dE).
  3. In the third step, the total energy deposited in all charged particle steps per event have to be computed, i.e. dE_sum = sum(dE) and dE’_sum = sum(dE’).
  4. As a last step, the weighting factor has to be computed to correct the total deposited energy per event E_tot ( can be obtained by the variable RULL, I guess) :

COMSCW = E’_tot / E_tot = (E_tot - dE_sum + dE’_sum) ./ E_tot

In this way, the corrected total energy deposition per event E’_tot will be computed by FLUKA as follows:

E’_tot = COMSCW * E_tot

Is this general approach correct? How does the FLUKLA compute this correction factor in the default case? I’m wondering, if and how FLUKA considers the local deposition events (deposition events, which are caused by the energy threshold)? Because, in my understanding, the trackr module doesn’t provide information about local energy deposition events, does it?

I know, it’s a lot of questions, but, as I said, unfortunately I’m a novice in user routines. Maybe, we can do a video call to answer all my questions. What do you think?

Thank you so much for your support
Best,
David

Payne et al., 2009, Nonproportionality of Scintillator Detectors:Theory and Experiment

Dear FLUKA community

In the meantime, I was able to implement a first version of a comscw.f user routine to imitate the tcquench card. The idea is that I validate first this approach, before implementing my own custom electron response function. Due to the shift of topic during the discussion with Francesco, I have shifted this discussion to the user-routine section: Imitate tcquench with userweig card.

Cheers
David