MCLTDCR: A Monte Carlo code for generation of list mode TDCR files

Contents lists available at ScienceDirect Applied Radiation and Isotopes journal homepage: www.elsevier.com/locate/apradiso MCLTDCR: A Monte Carlo code for generation of list mode TDCR files K. Mitev a,∗, V. Todorov a, P. Cassette a, B. Sabot b aSofia University ‘‘St Kliment Ohridski’’, Faculty of Physics, 5 James Bourchier Blvd., Sofia, 1164, Bulgaria bUniversité Paris-Saclay, CEA, LIST, Laboratoire National Henri Becquerel (LNE-LNHB), Palaiseau, F-91120, France A R T I C L E I N F O Keywords: TDCR counting Monte Carlo simulation List-mode acquisition Correction for accidental coincidences A B S T R A C T This work describes a Monte Carlo code developed for generation of list-mode TDCR files. The structure of the code, along with its input and output formats, is described. The significance of adequate modeling of the time distributions between the signals of the photomultiplier tubes is highlighted, and the options for achieving this with the current code are presented. The advantages of having datasets with known basic reference, such as those provided by Monte Carlo simulation, are demonstrated through two examples. It is shown, first, how the results of the code can be used for benchmarking list-mode TDCR analysis software, and second, how the code can be used to study the effect of correcting accidental coincidences on TDCR counting results. 1. Introduction The triple-to-double coincidence ratio (TDCR) method within liquid scintillation counting (LSC) is widely used for the standardization of a large number of pure 𝛽− and electron-capture radionuclides, such as 3H, 14C, 32P, 33P, 35S, 45Ca, 55Fe, 59Ni, 99Tc, 147Pm, 89Sr, 90Sr, 90Y, 241Pu, and others (Broda et al., 2007). For this reason, many national metrology laboratories and other institutions maintain and develop TDCR counting systems for radionuclide standardization. Over the last decade, the use of digitizers for acquisition of data in list-mode format for TDCR applications has become increasingly popular. This is mainly because such acquisition allows for flexible post-processing of the recorded data with different settings, such as varying the base deadtime duration or the coincidence time window. For instance, list-mode TDCR data acquisition is employed at the BIPM (Coulon et al., 2023), ENEA, Italy (Mini et al., 2014), LNHB, France (Sabot et al., 2024), NIM, China (Zhou et al., 2022), PTB, Germany (Takács and Kossert, 2021), SU, Bulgaria (Dutsov et al., 2019), and other institutions. Typically, the list-mode files are processed using custom-built software to extract the information relevant to the application of the TDCR method, and most of the institutions mentioned have developed their own processing codes. For the development and testing of custom-built software for processing TDCR list-mode files, it is useful to have a dataset with known basic reference (e.g., counting rates, source activity) against which to compare the results obtained from the software under test. Obtaining such datasets experimentally is challenging, as the exact values of activity, counting rates, and other parameters are typically not known. ∗Corresponding author. E-mail address: [email protected] (K. Mitev). The Monte Carlo (MC) method, on the other hand, allows for adequate statistical modeling of physical experiments and provides accurate information about the variables involved in the simulated process. Therefore, alongside the development of the cdt_logic code (Dutsov, 2021), which is a software for determination of TDCR related quantities (e.g. counting rates and time distributions) from list-mode files, it was decided to also develop a MC code for simulation of TDCR counting experiments. The code was intended to simulate as realistically as possible the TDCR counting process and to generate list-mode files with known MC basic reference. The objective of this work is to present the basic features of the developed MC code and to describe its input and output. Over the years, the approach of using MC-generated list-mode files has proven useful not only for testing TDCR analysis software, but also for studying various physical phenomena involved in TDCR counting. Therefore, this paper also presents examples of several applications of the code for investigating the underlying physics of the TDCR method. 2. Methods and materials 2.1. Input for the MC simulation The MCLTDCR code is written in Fortran 77 programming language. It reads the input values needed for the simulation from an input file, the content of which is shown in Fig. 1. The program reads the duration of the simulated interval 𝑇, the decay rate of the MC source (𝐴𝑆𝑅𝐶 ) and its half-life (𝑇1 2 ). Note that 𝐴𝑆𝑅𝐶 is supposed to be estimated from https://doi.org/10.1016/j.apradiso.2025.112094 Received 17 June 2025; Accepted 4 August 2025 Applied Radiation and Isotopes 226 (2025) 112094 Available online 18 August 2025 0969-8043/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC license ( http://creativecommons.org/licenses/bync/4.0/ ). K. Mitev et al. Fig. 1. Description of the input for the MC simulation. the TDCR measurement we want to simulate at the moment of its beginning. Next the program reads the so-called noise counting rates of the PMTs, if any. Further, the program takes as input the values of the net (background corrected) counting rates of individual PMTs, double and triple coincidences as they were determined in a given (performed) TDCR measurement. These counting rates are: 𝑛𝐴, 𝑛𝐵, 𝑛𝐶, for the PMTs A, B and C; 𝑛𝐴𝐵, 𝑛𝐵𝐶 , 𝑛𝐴𝐶 , for the double coincidences in channels AB, BC and AC and 𝑛𝐴𝐵𝐶 or 𝑛𝑇 for the triple coincidences. After that, the program reads the selected function for the time distribution between the pulses in the PMTs. The available distributions are: Voigt, Voigt+exponentially modified Gaussian distribution and an exponentially modified Voigt distribution. The blocks after the function selection specify the parameters of the distributions. In addition, the program can read experimental time distributions between the PMTs and sample from them. The next block sets the random number generator seeds and the last block specifies the output file names and type. 2.2. Setting up the probabilities for the MC simulation The code uses the values of the net (background corrected) counting rates of individual PMTs, double and triple coincidences (𝑛𝐴, 𝑛𝐵, 𝑛𝐶, 𝑛𝐴𝐵, 𝑛𝐵𝐶 , 𝑛𝐴𝐶 , 𝑛𝑇) as they were determined in a given (performed) TDCR measurement for a given coincidence window (CW). It also uses the value of the MC source decay rate 𝐴𝑆𝑅𝐶 , estimated from the TDCR measurement at the moment of its beginning. The counting rate of the logical sum of double coincidences 𝑛𝐷 is calculated as: 𝑛𝐷=𝑛𝐴𝐵 +𝑛𝐵𝐶 +𝑛𝐴𝐶 − 2𝑛𝑇(1) The logical sum of the single events is then calculated as: 𝑛𝑆=𝑛𝐴+𝑛𝐵+𝑛𝐶−𝑛𝐷−𝑛𝑇(2) In order to perform a MC simulation, we need first to evaluate the counting rates and the probabilities of elementary events. In the context of TDCR counting, the elementary events (hereafter referred to as ‘‘pure events’’) can be considered as events that belong to only one of the categories ‘‘Single’’, ‘‘Double’’ or ‘‘Triple’’. The concept for calculating the counting rates of the pure events is sketched in Fig. 2. The counting rates of the pure ‘‘Single’’, ‘‘Double’’ or ‘‘Triple’’ events are calculated as: 𝐶𝑇 𝑅𝐼𝑃 𝑆 =𝑛𝑇 𝐶𝐷𝐵𝐿𝐸𝑆 =𝑛𝐷−𝑛𝑇 𝐶𝑆𝐼𝑁𝐺𝐿𝐸𝑆 =𝑛𝑆−𝑛𝐷 (3) Here 𝐶𝑇 𝑅𝐼𝑃 𝑆 , 𝐶𝐷𝐵𝐿𝐸𝑆 and 𝐶𝑆𝐼𝑁𝐺𝐿𝐸𝑆 are the counting rates of the events that contribute to only one of these categories. For instance, if an event contributes to 𝐶𝐷𝐵𝐿𝐸𝑆 it is a pure double event, meaning that it is not and cannot be a triple event. Similarly, we define also pure A, B, C, AB, BC and AC events: 𝐶𝑃 𝑈𝑅𝐴 =𝑛𝐴−𝑛𝐴𝐵 −𝑛𝐴𝐶 +𝑛𝑇 𝐶𝑃 𝑈𝑅𝐵 =𝑛𝐵−𝑛𝐴𝐵 −𝑛𝐵𝐶 +𝑛𝑇 𝐶𝑃 𝑈𝑅𝐶 =𝑛𝐶−𝑛𝐴𝐶 −𝑛𝐵𝐶 +𝑛𝑇 𝐶𝑃 𝑈𝑅𝐴𝐵 =𝑛𝐴𝐵 −𝑛𝑇 𝐶𝑃 𝑈𝑅𝐵𝐶 =𝑛𝐵𝐶 −𝑛𝑇 𝐶𝑃 𝑈𝑅𝐴𝐶 =𝑛𝐴𝐶 −𝑛𝑇 (4) It is also necessary to define undetected events, i.e. true decays that have not been detected by the TDCR system. The rate of these undetected events (𝐶𝑍𝐸𝑅𝑂𝐸𝑆 ) is estimated as: 𝐶𝑍𝐸𝑅𝑂𝐸𝑆 =𝐴𝑆𝑅𝐶 − (𝐶𝑇 𝑅𝐼𝑃 𝑆 +𝐶𝐷𝐵𝐿𝐸𝑆 +𝐶𝑆𝐼𝑁𝐺𝐿𝐸𝑆 )(5) Note that Eq. (5) implies that a decay can contribute to single, double, triple or undetected event. Apparently, no decay schemes are taken into account here and the quantity 𝐴𝑆𝑅𝐶 is actually the decay rate of the MC source, rather than the activity of the source measured in the TDCR counting experiment. The link between the two has to be established by the user considering the decay data of the particular radionuclide and what is detected in the TDCR measurement. With this data, we can now estimate the probabilities 𝑃 for pure single, double, triple or undetected events: 𝑃𝑆𝐼𝑁𝐺𝐿𝐸𝑆 =𝐶𝑆𝐼𝑁𝐺𝐿𝐸𝑆 𝐴𝑆𝑅𝐶 , 𝑃𝐷𝐵𝐿𝐸𝑆 =𝐶𝐷𝐵𝐿𝐸𝑆 𝐴𝑆𝑅𝐶 𝑃𝑇 𝑅𝐼𝑃 𝐿𝐸𝑆 =𝐶𝑇 𝑅𝐼𝑃 𝐿𝐸𝑆 𝐴𝑆𝑅𝐶 , 𝑃𝑍𝐸𝑅𝑂𝐸𝑆 =𝐶𝑍𝐸𝑅𝑂𝐸𝑆 𝐴𝑆𝑅𝐶 (6) Then, we can also define the conditional probabilities for pure A, B or C events (𝑃𝑖, 𝑖=𝐴, 𝐵, 𝐶), given the occurrence of a single event: 𝑃𝐴=𝐶𝑃 𝑈𝑅𝐴 𝐶𝑆𝐼𝑁𝐺𝐿𝐸𝑆 , 𝑃𝐵=𝐶𝑃 𝑈𝑅𝐵 𝐶𝑆𝐼𝑁𝐺𝐿𝐸𝑆 , 𝑃𝐶=𝐶𝑃 𝑈𝑅𝐶 𝐶𝑆𝐼𝑁𝐺𝐿𝐸𝑆 (7) Finally, we can define the conditional probabilities for pure AB, BC or AC events (𝑃𝑗, 𝑗=𝐴𝐵, 𝐵𝐶, 𝐴𝐶) given the occurrence of a pure double event: 𝑃𝐴𝐵 =𝐶𝑃 𝑈𝑅𝐴𝐵 𝐶𝐷𝐵𝐿𝐸𝑆 , 𝑃𝐵𝐶 =𝐶𝑃 𝑈𝑅𝐵𝐶 𝐶𝐷𝐵𝐿𝐸𝑆 , 𝑃𝐴𝐶 =𝐶𝑃 𝑈𝑅𝐴𝐶 𝐶𝐷𝐵𝐿𝐸𝑆 (8) The probabilities defined in Eqs. (6) and the conditional probabilities (Eqs. (7) and (8)) are used to simulate the TDCR counting experiment. It should be noted here that these probabilities, as well as the source Applied Radiation and Isotopes 226 (2025) 112094 2 K. Mitev et al. Fig. 2. Description of the probabilities used in the MC simulation. decay rate, are determined from the TDCR counting rates. These rates, in turn, are determined for a given coincidence window, meaning the probabilities in Eqs. (6), (7), and (8) are, in principle, related to the coincidence window used in the TDCR analysis. Note also that the occurrence of after-pulses in the PMTs is not considered in the simulation. As long as after-pulses typically follow detected events, it is assumed that they do not interfere with the analysis because they are ignored by the imposed extending-type dead-time in the TDCR analysis software. 2.3. Structure of the MC simulation The idea of the simulation is to simulate a TDCR measurement of a radioactive source containing one radionuclide with known half-life (𝑇1 2 ). Thus, we want to generate a time series of events that happen within the TDCR counter as a result of the decays happening in the radioactive source. Each of these events represents a radioactive decay followed by a particle emission within the sample. The occurrence of the events is considered an inhomogeneous Poisson process with nonstationary rate 𝐴𝑆𝑅𝐶 (𝑇) (Furbish, 2013). Then, it follows that the time interval 𝑡 between events occurring at moments 𝑇 and 𝑇+𝑡 is exponentially distributed with probability density function 𝑝(𝑡) (McQuighan, 2010): 𝑝(𝑡) = 𝐴𝑆𝑅𝐶 (𝑇)⋅𝑒−𝑡⋅𝐴𝑆𝑅𝐶 (𝑇)(9) Thus, the sampling of the time until the next event 𝑡𝑛𝑒𝑥𝑡 is performed by the inverse transform method: 𝑡𝑛𝑒𝑥𝑡 = − ln 𝑅 𝐴𝑆𝑅𝐶 (𝑇)(10) where 𝑅 is a number between 0 and 1 generated by the pseudo random number generator. The pseudo random number generator used in this code is the one from the PENELOPE code system (Salvat, 2018). A flow diagram of the simulation is shown in Fig. 3. The simulation begins with the initialization subroutine (ReadIn) which reads the input file, calculates the needed probabilities and sets the global time to zero. Then subroutine TEVENT samples the time to the next decay 𝑡𝑛𝑒𝑥𝑡 (by Eq. (10)), sets the global time 𝑇=𝑇+𝑡𝑛𝑒𝑥𝑡 and samples the type of the next event. If the event is a source decay, the subroutine ITYPEC samples what kind of decay event has occurred. The decay events are: undetected decay, single, double or triple coincidence and the type of decay event is sampled from the probabilities determined in Eqs. (6). If an undetected event occurs, it is written in the list mode file with its time stamp and associated information (see Fig. 4 Right). If a single event is deemed to occur, the subroutine IMPTs samples which PMT is triggered from the conditional probabilities calculated in Eqs. (7) and the event is written in the list mode file. If a double coincidence is deemed to occur, the triggered channel (AB, BC or AC) is sampled from the conditional probabilities calculated in Eqs. (8). In the case of double or triple coincidence events, one needs to simulate correctly the time interval (𝛥𝑡) between the pulses generated by the PMTs. This is done in the code by a function named TBP, acronym for Time Between Pulses, described below. In the case of a double coincidence, the first triggered PMT takes the timestamp of the decay event (𝑇) and the second PMT takes the timestamp 𝑇+𝛥𝑡, where 𝛥𝑡 is sampled by the function TBP. All triggered PMTs are then written as separate lines in the list mode file. The same strategy is adopted and for the triple coincidences. The same subroutine (WRITELINE) is used to write a line in the list mode file for all types of events: undetected, singles, etc. When an event is fully simulated, the simulation continues with the sampling of the time to the next event by the subroutine TEVENT. The cycle continues until reaching a user defined time 𝑇𝑚𝑎𝑥. When the simulation is complete, the code generates two output files. The first file contains detailed information about the simulated events (see Fig. 4, left), including the number of total events, decays, noise events, pure singles, pure doubles, and pure triples, along with their corresponding simulated counting rates. It also reports the number of simulated singles, doubles, and triples, their associated counting rates, the TDCR value from the simulation, and the ratios T/AB, T/BC, and T/AC. This output is intended for benchmarking TDCR analysis codes for list-mode files, as it provides the basic reference (for a given coincidence window used) against which the results of these analysis codes can be compared. The second output file is the generated listmode file itself. The information written in it is as follows (see Fig. 4, right): Column 1: Time of the event from start, s; Column 2: Event type: 0 = Noise, 1 = Decay; Column 3: Decay type: -1 = Noise, 0 = Undetected decay, 1 = S, 2 = D 3 = T; Column 4: PMT triggered 0 = none, 1 = A, 2 = B, 3 = C; Column 5: Number of simulated decay events; Column 6: Number of simulated noise events, Column 7: Number of all simulated events. Applied Radiation and Isotopes 226 (2025) 112094 3 K. Mitev et al. Fig. 3. Description of the simulation flow. Fig. 4. Example of the output files: (a) the file containing statistics for simulated events (SimStatis1.dat), and (b) the simulated list mode file (Pulses_TDC1.dat). Applied Radiation and Isotopes 226 (2025) 112094 4 K. Mitev et al. Fig. 5. Example of MC sampling from Voigt distribution. 2.4. Sampling the time distribution between pulses from the PMTs The time distribution between the signals of the triggered PMTs (termed also cross-correlation distribution) is of paramount importance in the TDCR counting process. It has been shown to contain information about the probability to detect a scintillation event (i.e. the detection efficiency) in a TDCR counting experiment (Mitev et al., 2021). The shape of this distribution depends on many factors, such as: the scintillation cocktail (i.e. the ratio between its prompt and delayed scintillation emission), the scintillation non-linearity of the light emission of the cocktail (which is due to different ionization quenching for different stopping powers, hence the dependence on the beta-spectrum of the measured radionuclide), the counting statistics, combinatorial term accounting for the probability of detecting in the PMTs and the time response function of the PMTs. A theoretical derivation of the shape of the cross-correlation distribution is given in Mitev et al. (2021). Due to the complex shape of the time distribution between pulses, we have adopted a simplified approach for its simulation in the present work. We chose to sample the time between pulses from simple distributions that were previously found to approximate the measured time distributions reasonably well. The predefined distributions include the Voigt distribution (VGT), the exponentially modified Gaussian distribution (EM_Gauss), a combination of Voigt and exponentially modified Gaussian (VGT+EM_Gauss), and the exponentially modified Voigt distribution (EM_VGT). If none of these distributions adequately match the experimental time distribution, we have also implemented an option to sample directly from the experimentally measured time-distribution histogram. The Gaussian distribution is defined as 𝐺(𝑥;𝑥0, 𝜎𝑔) = 1 √2𝜋𝜎𝑔2 exp −(𝑥−𝑥0)2 2𝜎𝑔2(11) where 𝑥0, 𝜎𝑔2 are the mean (i.e. position of the centroid) and the variance of the distribution. The sampling from the Gaussian distribution is performed by the Box–Muller method (Box and Muller, 1958). The Voigt distribution is a convolution between a centered Gaussian distribution and a Lorentzian distribution (𝐿): 𝑉(𝑥;𝜎𝑔, 𝛾𝐿) = ∫∞ −∞ 𝐺(𝑥′;𝜎𝑔)𝐿(𝑥−𝑥′;𝛾𝐿)𝑑𝑥′(12) with 𝐿(𝑥;𝑥0, 𝛾𝐿) = 𝛾𝐿 𝜋((𝑥−𝑥0)2+𝛾𝐿2) and 𝑥0 being the centroid and 𝛾𝐿 being the half-width at half-maximum of 𝐿(𝑥). The composition method proposed by Lee (1974) is used for sampling from the Voigt distribution. Example of the MC sampling from the Voigt distribution with our algorithm is shown in Fig. 5. The normalized residuals shown in Figs. 5–7 are calculated as: 𝑁𝑜𝑟𝑚.𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙𝑠 =𝑝(𝑥) − 𝑝𝑀𝐶 (𝑥) 𝜎𝑀𝐶 (𝑥)100 (13) where 𝑝(𝑥) is the value of the distribution (theoretical or experimental) at point 𝑥, 𝑝𝑀𝐶 (𝑥) is the MC estimate of the value of the distribution at the same point and 𝜎𝑀𝐶 (𝑥) is the statistical uncertainty of the MC estimate. The EM_Gauss and EM_VGT distributions are defined as: 𝐸𝑀_𝐺𝑎𝑢𝑠𝑠(𝑥;𝜇𝑔, 𝜎𝑔, 𝜆) = ∫𝑥 −∞ 𝐺(𝑥′;𝜇𝑔, 𝜎𝑔)⋅𝜆exp−𝜆(𝑥−𝑥′)𝑑𝑥′(14) 𝐸𝑀_𝑉 𝑜𝑖𝑔𝑡(𝑥;𝜎𝑔, 𝜆) = ∫𝑥 −∞ 𝑉(𝑥′;𝜎𝑔, 𝛾𝐿)⋅𝜆exp−𝜆(𝑥−𝑥′)𝑑𝑥′(15) The samplings from the exponentially modified distributions are performed by the composition method, similarly to the approach proposed by Lee (1974). The usage of the composition method is preferred here because it avoids the explicit calculation of the probability density function of the Voigt or EM_Gauss distributions. The code can sample also from superimposed Voigt and EM Gaussian distributions (denoted VGT+EM_Gauss). In this case, it is assumed that the random variable to be sampled (𝑍) is a sum of two other independent random variables 𝑋 (Voigt distributed) and 𝑌 (Gauss distributed). Explicitly, 𝑍=𝑎𝑋 + 𝑏𝑌 , where 𝑎 and 𝑏 are user-selectable areas of the Voigt and Gauss distributions. An example of such sampling is shown in Fig. 6. An adequate simulation of the time distribution between the signals of the triggered PMTs is important for the generation of list mode files. Sometimes, the above functions are sufficient to model the experimental time distribution. For instance, it was found in Dutsov et al. (2019) that a sum of Voigt+EM_Gauss distributions may approximate well the measured time distribution. However, this may not always be the case and in order to handle these exceptions, the code has also the functionality to sample directly from histograms of the measured time distributions. The sampling is performed with the subroutine IRND of the PENELOPE-2018 code system (Salvat, 2018). This subroutine is used here because it provides fast and accurate sampling from discrete distributions using the Walker’s aliasing method. An example of sampling from an experimental distribution is shown in Fig. 7. It should be emphasized that the MCLTDCR code was developed to simulate the time distribution using empirical functions that fit experimental data. As such, the code is not intended to simulate the fundamental properties of the scintillation process—such as modeling the prompt and delayed components of scintillation emission from first principles. A separate code, developed by Dutsov (2021), takes Applied Radiation and Isotopes 226 (2025) 112094 5 K. Mitev et al. Fig. 6. Example of MC sampling from Voigt+EM_Gauss distributions. Fig. 7. Example of MC sampling from an experimental time distribution histogram with the function IRND from the Penelope-2018 code system (Salvat, 2018). a different approach: it aims to simulate the timing of events using a minimal set of fundamental parameters and physical models. This code is described by Dutsov (2021) and Dutsov et al. (2020a), and it generates the time distribution between pulses based on physical models of prompt and delayed fluorescence emissions in scintillation cocktails. Both codes have been applied in various studies. The MCLTDCR code was used in research on different TDCR counting algorithms (Dutsov et al., 2019), as well as during the development of corrections for accidental coincidences (Dutsov et al., 2020b). Dutsov’s code was employed in studies on the selection of optimal coincidence resolving time (Dutsov et al., 2020a), investigations of the time domain and the evaluation of the cross-correlation function (Mitev et al., 2021), in research into the significance of correcting for accidental coincidences in TDCR counting (Dutsov et al., 2022) and in the evaluation of time properties of scintillator in the Compton-TDCR experiment (Sabot et al., 2024). Overall, both codes have proven their usefulness over time in advancing the understanding of the specific characteristics of TDCR counting. 3. Results In the following we will present some results obtained with the MCLTDCR code and will discuss their usage in TDCR counting studies. Applied Radiation and Isotopes 226 (2025) 112094 6 K. Mitev et al. Fig. 8. Relative difference between simulated and calculated values as a function of the coincidence window (CW), for realistic (Fig. 8(a)) and narrow (Fig. 8(b)) time distributions. The solid and open symbols on the figure indicate data that is obtained with or without correction for accidental coincidences, respectively. Table 1 Input values used in the simulations for benchmarking the cdt_logic code. Parameter Input name 3H-like time distribution Narrow time distribution Centroid Voigt peak, ns Position of Voigt peak 0.0 0.0 Area Voigt peak Area of the Voigt peak 1.0 1.0 Gaussian FWHM, ns Gaussian FWHM 2.82E−9 2.82E−15 Lorentzian HWHM, ns Lorentzian HWHM 2.12E−9 2.12E−15 Centroid EM_Gauss, ns Centroid AGS 0.0 0.0 Sigma EM_Gauss, ns Sigma AGS 2.82E−9 2.82E−15 Area EM_Gauss Area AGA 0.141 0.141 Left exponential parameter Left exp. parameter 0.0E0 0.0E0 Right exponential parameter Right exp. parameter 862E−9 862E-15 3.1. Benchmarking list-mode TDCR analysis software A straightforward application of the MCLTDCR code is to evaluate the performance of software used for analyzing list-mode TDCR data files. The code used at Sofia University is developed by Ch. Dutsov and is called cdt_logic (the source code and its description are provided in Dutsov, 2021). To benchmark the output of this analysis code against Monte Carlo basic reference, we conducted two simulations using different time distributions. The input parameters for both simulations (listed in Table 1) are identical, except for the width of the time distribution. The first simulation uses a realistic time distribution that mimics 3H acquisition (the input file is shown in Fig. 1). The second simulation uses a non-physical, extremely narrow time distribution with widths on the order of 10−15 seconds, chosen for demonstration purposes. For the purpose of the analysis, the relative difference (𝛥) is defined as: 𝑅𝑒𝑙.𝐷𝑖𝑓𝑓 =𝛥=𝑋𝑠𝑖𝑚𝑢𝑙𝑎𝑡𝑒𝑑 −𝑋𝑐𝑎𝑙𝑐𝑢𝑙𝑎𝑡𝑒𝑑 𝑋𝑠𝑖𝑚𝑢𝑙𝑎𝑡𝑒𝑑 100 (16) where 𝑋𝑠𝑖𝑚𝑢𝑙𝑎𝑡𝑒𝑑 is the value of a given variable 𝑋 obtained from the Monte Carlo (MC) simulation (i.e., the basic reference), and 𝑋𝑐𝑎𝑙𝑐𝑢𝑙𝑎𝑡𝑒𝑑 is the value of the same variable determined from the analysis of the list-mode file using the cdt_logic code. The two list-mode files generated with MCLTDCR with the inputs shown in Table 1, were analyzed with the cdt_logic code, applying a 10 μs base dead-time duration and 40 ns coincidence resolving time (coincidence window, CW). We chose a 40 ns coincidence resolving time because it is the value set in the MAC3 TDCR acquisition module (Bouchard and Cassette, 2000) and has been widely used previously. The results of the comparison are shown in Table 2. Note that the activity of the source for the 3H-like distribution is determined using a value of 𝑘𝐵 = 0.010 cm/MeV, estimated from the results shown in Fig. 10(a). The activity of the source for the narrow distribution is determined using an optimal value of 𝑘𝐵 = 0.012 cm/MeV, estimated from the results shown in Fig. 9(a). The results in Table 2 show that for the narrow distribution the estimated values with the cdt_logic code practically coincide with the MC input. This is an indication that the results obtained with the cdt_logic code are correct. For the realistic 3H-like distribution however, the estimated counting rates are systematically smaller than the MC reference values. This is because with 40 ns coincidence window there is a chance to lose some events with long intervals between the detection of pulses in the PMTs. Nevertheless, the TDCR method in this case somehow manages to compensate the missed events and gives a result for the calculated activity that is quite coherent with the MC reference. 3.2. Applying the code to study the properties of the TDCR method To take the comparison a step further, each of the two output files generated using the input parameters from Table 1 was analyzed using different coincidence windows, both with and without applying the analytical correction for accidental coincidences proposed in Dutsov et al. (2020b). The results are presented in Fig. 8. It is immediately evident from the figure that correct application of the TDCR method requires a correction for accidental coincidences. Interestingly, for the simulation with the narrow time distribution, the correction works very well—after correction, all calculated quantities closely match the Monte Carlo reference values (Fig. 8(b)). This is not exactly the case for the simulation using the realistic 3H-like time distribution. In this case, the correction performs better with larger coincidence windows and less effectively with smaller ones (Fig. 8(a)). This is likely due to the aforementioned missed events when using narrow coincidence windows. Nevertheless, as shown in Table 2, the TDCR method appears to compensate for the effect of missed events when it comes to the final activity calculation. Applied Radiation and Isotopes 226 (2025) 112094 7 K. Mitev et al. Table 2 Comparison between MC simulated and estimated with cdt_logic TDCR values for the inputs from Table 1. Corrections for accidental coincidences are not applied. Activity estimates use CW = 40 ns, with 𝑘𝐵 = 0.012 cm/MeV for the narrow distribution and 𝑘𝐵 = 0.010 cm/MeV for the 3H-like distribution. Variable 3H-like time distribution Narrow time distribution MC simulated Estimated 𝛥MC simulated Estimated 𝛥 A, cps 10680 10612 0.6% 10680 10682 −0.01% B, cps 11914 11856 0.5% 11914 11915 −0.01% C, cps 10619 10 547 0.7% 10619 10620 −0.01% AB, cps 6365 6262 1.6% 6365 6369 −0.05% BC, cps 6221 6116 1.7% 6221 6225 −0.06% AC, cps 5677 5568 1.9% 5677 5681 −0.07% T, cps 4370 4248 2.8% 4370 4374 −0.09% D, cps 9523 9450 0.8% 9523 9527 −0.04% T/AB 0.687 0.678 1.2% 0.687 0.687 −0.03% T/BC 0.703 0.695 1.1% 0.703 0.703 −0.02% T/AC 0.770 0.763 0.9% 0.770 0.770 −0.02% T/D 0.459 0.450 2.1% 0.459 0.459 −0.05% Activity, Bq 19022 18998 0.1% 19022 19014 0.04% Fig. 9. Relative difference between simulated and calculated activity as a function of the 𝑘𝐵 value for the narrow time distribution and different coincidence windows without (Fig. 9(a)) and with (Fig. 9(b)) correction for accidental coincidences. An interesting and subtler example of the code’s application is in demonstrating the importance of correcting for accidental coincidences when selecting the 𝑘𝐵 value in TDCR counting. The list mode files from the previous example are analyzed using the cdt_logic code, and the TDCR17 code (Cassette, 2017) is then applied to determine the source activity for various 𝑘𝐵 values. The results obtained for the narrow distribution are shown in Fig. 9. They indicate that when raw data, uncorrected for accidentals, is analyzed there is a spread in the estimated activity for each 𝑘𝐵 value. This impedes the choice of the correct 𝑘𝐵 value. In contrast, when the correction for accidentals is applied, the spread is eliminated and the correct 𝑘𝐵 value (0.012 cm/MeV in this case) is easily identified. Similar tendency, although not so strongly manifested, is seen in the results with the 3H-like time distribution (Fig. 10). The analysis of the raw data (i.e. uncorrected for accidentals) exhibits the same spread of deviations from the MC reference for each 𝑘𝐵 value (Fig. 10(a)). Applying the accidental coincidences correction slightly reduces the spread, especially for larger coincidence windows, but a significant amount still remains in the results (Fig. 10(b)). Interestingly, the accidentals correction does not significantly alter the 40 ns and 100 ns data points. However, it clearly tightens the grouping of the longer CW points. A possible explanation for this lies in Fig. 9(a). In this figure, the T/D values corrected for accidental coincidences are significantly biased for small CW values. Conversely, even with a higher number of accidentals for large CW values, the corrected T/D values for large CW do not show bias when compared to the MC results. This likely accounts for the tight grouping of points in Fig. 10(b) when CW ≥ 200 ns. 4. Conclusions This work describes the development of a Monte Carlo code designed for the generation of list-mode TDCR files. The code’s structure, its input and output formats, and the types of time distributions between PMT pulses that it can model are thoroughly outlined. Illustrative examples of the code’s application are subsequently presented and discussed. We demonstrate its utility for benchmarking TDCR analysis software by comparing its outputs to Monte Carlo reference data, and for investigating the influence of accidental coincidences correction. Additionally, this Monte Carlo code is employed to elucidate the potential impact of accidental coincidences correction on the determination of the 𝑘𝐵 value in TDCR counting experiments. Based on our experience, three directions for future developments appear interesting. The first is the implementation of a method for sampling from arbitrary functions defined within a function program. This enhancement is essential for more accurate modeling of experimentally measured time distributions between PMT pulses. The second proposed development involves the inclusion of a fourth channel in the simulation, enabling the generation of list-mode files compatible with four-channel acquisition systems, such as TDCR-gamma and ComptonTDCR counters. The third development is related to implementation of algorithms to handle time-correlated decays, which may occur for instance when counting two radionuclides in secular equilibrium and the progeny has a very short half-life. These developments are planned for the future. Applied Radiation and Isotopes 226 (2025) 112094 8 K. Mitev et al. Fig. 10. Relative difference between simulated and calculated activity as a function of the 𝑘𝐵 value for realistic 3H-like time distribution and different coincidence windows without (Fig. 10(a)) and with (Fig. 10(b)) correction for accidentals. CRediT authorship contribution statement K. Mitev: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft. V. Todorov: Data curation, Formal analysis, Investigation, Software, Validation, Writing – review & editing. P. Cassette: Software, Validation, Writing – review & editing. B. Sabot: Methodology, Resources, Validation, Writing – review & editing. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: K. Mitev reports financial support was provided by National Recovery and Resilience Plan of the Republic of Bulgaria. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper Acknowledgments This work is supported by the European Union NextGenerationEU, through the National Recovery and Resilience Plan of the Republic of Bulgaria, project No BG-RRP-2.004-0008-C01. This work is supported by the project RadonNET (23IND07), which received funding from the European Partnership on Metrology (ID: 10. 13039/100019599), co-financed from the European Union’s Horizon Europe Research and Innovation Programme and by the Participating States. Data availability Data will be made available on request. References Bouchard, J., Cassette, P., 2000. MAC3: an electronic module for the processing of pulses delivered by a three photomultiplier liquid scintillation counting system. Appl. Radiat. Isot. 52 (3), 669–672. http://dx.doi.org/10.1016/s0969-8043(99) 00228-6. Box, G.E.P., Muller, M.E., 1958. A Note on the Generation of Random Normal Deviates. Ann. Math. Stat. 29 (2), 610–611. http://dx.doi.org/10.1214/aoms/1177706645. Broda, R., Cassette, P., Kossert, K., 2007. Radionuclide metrology using liquid scintillation counting. Metrologia 44 (4), S36. http://dx.doi.org/10.1088/0026-1394/44/ 4/S06. Cassette, P., 2017. TDCR17– Detection efficiency calculation for pure-beta radionuclides, LNHB program with short tutorial, 2017 version. http://www.lnhb.fr/home/ conferences-publications/icrm_lsc_wg/icrm_lsc_software/. Coulon, R., Gressier, V., Michotte, C., 2023. Test of a digitizer to process the pulse signal from the 3 photomultiplier tubes of a TDCR liquid scintillation counter. Appl. Radiat. Isot. 192, 110598. http://dx.doi.org/10.1016/j.apradiso.2022.110598, URL: https://www.sciencedirect.com/science/article/pii/S0969804322004833. Dutsov, C., 2021. Studies on the Application of the Triple-to-Double Coincidences Ratio Method for Primary Standardisation Using Liquid Scintillation Counting (Ph.D. thesis). Sofia University ‘‘St. Kliment Ohridski’’, Sofia, URL: https://physica.dev/ files/phd_thesis.pdf. Dutsov, C., Cassette, P., Mitev, K., Sabot, B., 2020a. In quest of the optimal coincidence resolving time in TDCR LSC. Nucl. Instrum. Methods Phys. Res. Sect. A 164846. http://dx.doi.org/10.1016/j.nima.2020.164846. Dutsov, C., Cassette, P., Sabot, B., Mitev, K., 2020b. Evaluation of the accidental coincidence counting rates in TDCR counting. Nucl. Instrum. Methods Phys. Res. Sect. A 977, 164292. http://dx.doi.org/10.1016/j.nima.2020.164292. Dutsov, C., Mitev, K., Cassette, P., Jordanov, V., 2019. Study of two different coincidence counting algorithms in TDCR measurements. Appl. Radiat. Isot. 154, 108895. http://dx.doi.org/10.1016/j.apradiso.2019.108895. Dutsov, C., Sabot, B., Cassette, P., Mitev, K., 2022. Significance of the corrections for accidental coincidences in liquid scintillation counting measurements. J. Radioanal. Nucl. Chem. 331 (8), 3303–3311. http://dx.doi.org/10.1007/s10967-022-08316-y. Furbish, D.J., 2013. Cool probabilistic things we typically don’t teach our students about radioactive decay, but should. URL: https://my.dev.vanderbilt.edu/ davidjonfurbish/wp-content/uploads/sites/978/2013/06/Radioactive-Decay.pdf. Lee, J.S., 1974. Monte Carlo Simulation of Voigt Distribution in Photon Diffusion Problems. Astrophys. J. 187, 159–162. http://dx.doi.org/10.1086/152603. McQuighan, P., 2010. Simulating the Poisson process. URL: https://www.math. uchicago.edu/~may/VIGRE/VIGRE2010/REUPapers/Mcquighan.pdf. Mini, G., Pepe, F., Tintori, C., Capogni, M., 2014. A full digital approach to the TDCR method. Appl. Radiat. Isot. 87, 166–170. http://dx.doi.org/10.1016/j.apradiso. 2013.11.103. Mitev, K., Dutsov, C., Cassette, P., Sabot, B., 2021. Time-domain based evaluation of detection efficiency in liquid scintillation counting. Sci. Rep. 11 (1), 12424. http://dx.doi.org/10.1038/s41598-021-91873-1. Sabot, B., Dutsov, C., Cassette, P., Mitev, K., Hamel, M., Bertrand, G.H.V., Lebbou, K., Dujardin, C., 2024. A compact detector system for simultaneous measurements of the light yield non-linearity and timing properties of scintillators. Sci. Rep. 14 (1), 6960. http://dx.doi.org/10.1038/s41598-024-57186-9. Salvat, F., 2018. PENELOPE-2018: A code system for Monte Carlo simulation of electron and photon transport. Issy-les-Moulineaux, France: OECD/NEA Data Bank, http://www.nea.fr/lists/penelope.html. Takács, M.P., Kossert, K., 2021. Half-life determination of short-lived nuclear levels in 237Np (59.54 keV), 233Pa (86.47 keV) and 227Ac (27.37 keV). Appl. Radiat. Isot. 176, 109858. http://dx.doi.org/10.1016/j.apradiso.2021.109858, URL: https: //www.sciencedirect.com/science/article/pii/S0969804321002591. Zhou, Q., Liu, H., Yang, Z., Zeng, W., Liang, J., 2022. Development of a portable TDCR system at NIM, China. Appl. Radiat. Isot. 187, 110315. http://dx.doi.org/10.1016/ j.apradiso.2022.110315, URL: https://www.sciencedirect.com/science/article/pii/ S096980432200207X. Applied Radiation and Isotopes 226 (2025) 112094 9