Monte Carlo modelling of pixel clusters in Timepix detectors using the MCNP code

The track structure of the signal measured by the semiconductor pixel detector Timepix3 was modelled in the Monte Carlo MCNP ® code. A detailed model at the pixel-level (256 × 256 pixels, 55 × 55 µ m 2 pixel size) was developed and used to generate and store clusters of adjacent hit pixels observed in the measured data because of particle energy deposition path, charge sharing, and drift processes. An analytical model of charge sharing effect and the detector energy resolution was applied to the simulated data. The method will help the user sort the measured clusters and distinguish radiation components of mixed fields by determining the response of Timepix3 detector to particular particle types, energies, and incidence angles that cannot be measured separately.


Introduction
In the planning and delivery of external radiation therapy with ionizing radiation (e.g., photons, electrons, protons), i.e. "radiotherapy", the main focus placed on accurate delivery of the dose to the target volume and minimization of the dose to the surrounding critical organs. Currently, another topic of great interest in radiotherapy is the dosimetry outside the primary particle beam, or out-of-field dosimetry [1].
Out-of-field dose is caused by largely scattered primary radiation and by secondary radiation generated in interactions of primary beam with material. Mixed radiation fields are composed of particles of different type and energy, also called "stray radiation" or "scattered radiation", and their composition heavily depends on the parameters of the primary beam, target material, and the distance from the beam axis. The stray radiation delivers extra dose to healthy tissue and critical organs outside the target volume. The dosimetry of mixed fields is challenged by the need of the identification of the field components.
The delivered amount of radiation to reach the therapeutic effect is expressed by absorbed dose, D T , (in Gray, Gy). The result has a deterministic effect visible in a short period of time after irradiation, e.g., death or degeneration of a large number of cells at a time that may lead to the recovery from cancer. D T is directly measurable by a dedicated instrument, e.g., ionization chamber. On the other hand, out-of-field dosimetry is covered by radiological protection dosimetry because the radiation dose is much lower, in the region of stochastic effects below the threshold of deterministic effects. The equivalent dose in an organ or tissue, H T (in Sievert, Sv), is used instead. This quantity is not directly measurable but is obtained from the D T by applying radiation weighting factors, w R , dependent on radiation type and energy, to reflect the relative biological effectiveness of the radiation [2]. This means that the determination of H T requires the knowledge of the field composition.
In personalized dosimetry, the capability of measurement and proper identification of field components plays an important role in the risk assessment for various target categories of patients.
The development of traceable and validated methods for characterization of stray radiation is one of several goals of the joint European research project "Metrology for advanced radiotherapy using particle beams with ultra-high pulse dose rates" (18HLT04, UHDpulse) [34] aiming to develop a metrological framework for traceable absorbed dose measurements in ultra-high-dose-per-pulse (UHD) particle beams. It runs from 2019 to 2023 within the European Metrology Programme for Innovation and Research (EMPIR) organized by the European association of national metrology institutes EURAMET.
The project partner company ADVACAM, Czech Republic, aims to design, manufacture, test and validate an optimized pixel Timepix3based detector for measurements of the absorbed dose from stray radiation fields generated by UHD beams [5] and its recalculation to the equivalent dose using the knowledge of the field composition. The development includes the work on methods to distinguish and characterize radiation components of mixed and time-dependent stray radiation fields occurring close to the UHD beam by the analysis of the measured signal. Due to the complexity of the task, Monte Carlo simulations were used in various steps of the work. This paper describes the Monte Carlo modeling of the real track structure of the signal measured by the Timepix3 detector. The method includes an adaptation of an analytical model of charge sharing effect, which facilitates precise modelling of the so-called clusters of the adjacent hit pixels observed in the measured data caused by the particle path itself, charge sharing to adjacent pixels, and drift processes in the sensor material. The method will allow to optimize the sorting of measured clusters to distinguish radiation components by determining the response of the Timepix3 detector to the particular particle types, energies, and incidence angles, which cannot be measured separately.

Timepix3 pixel detector
For the purpose of the dosimetry in mixed radiation fields generated by UHD beams, a specific customized version of MiniPIX Timepix3 Flex detector was developed (Fig. 1). The readout chip contains a matrix of 256 × 256 pixels (total 65 536 independent channels, pixel size 55 µm) and an active sensor area of 14 mm × 14 mm for a total sensitive area of 1.98 cm 2 [6]. This detector has the readout electronics detached from the Timepix3 ASIC and a sensor assembly to minimize scattering and self-shielding by using a 5-cm-long flexible printed circuit board [7]. Thanks to the timestamp of each hit with the precision of 1.5 ns and the energy absorbed in each individual pixel, the pixelized structure of the sensor allows recording of particle tracks. This version of the detector utilized the Si sensor 100 and 500 µm thick. The sensor was fully depleted as the high bias voltage was set to 50 V and 200 V, respectively. However, the pixel cluster modelling showed below was also tested for other sensors. If so, it is mentioned in the respective results section.
Particle tracking with Timepix-family detector is widely used and continuously improved, allowing up to 3D track mapping [8]. This paper deals with the standard 2D tracking with the emphasis on the charge sharing effect.

Monte Carlo code
Monte Carlo (MC) simulations were performed with the code MCNP® in version 6.2 [9]. It is a general-purpose code for coupled ionizing radiation particle transport, developed in Los Alamos National Laboratory, USA. Specific areas of application include, but are not limited to, radiation protection and dosimetry [10], radiation shielding [11], medical physics [12], detector design and analysis [13], decommissioning [14], radiography, nuclear criticality safety, etc. Important standard features that make MCNP® very versatile and easy to use include a powerful general source, criticality source, and surface source, both geometry and output tally plotters, a rich collection of variance reduction techniques, a flexible tally structure, and an extensive collection of cross-section data. MCNP® contains numerous flexible tallies: surface current and flux, volume flux, point and ring detectors, particle heating, fission heating, pulse height tally for energy or charge deposition, mesh tallies, and radiography tallies [15].
The lower cut-off energy in the detection sensor was set as follows: 1 keV for electrons, 10 keV for photons, 300 keV for all other charged  particles. Higher energy cut-off was set in materials surrounding the sensor to include the signal from particles created outside the sensor as well. No energy cut-off was set for neutrons. Generating of δ-electrons was turned off. Other parameters of particle physics on PHYS cards [9] were kept default.
No variance reduction methods were used because the weight of particles generating the signal in the sensor had to remain equal to one. Otherwise, the energy deposition and the ratio of different track types would be biased.

Timepix detector model
A detailed MC model of the MiniPIX Timepix3 FLEX device used in the UHDpulse project was developed comprising all parts relevant for ionizing particles transport (Fig. 2). Dimensions and material specification corresponded to the real device. The model included the sensor chip of the given material (Si, CdTe, or GaAs) and the thickness equal to the given sensor. The readout chip (720 µm thick Si), glue (250 µm) and bonding (20 µm) of the sensor chip were also included. A graphite cooling plate below the readout chip, flex cables connecting the readout chip with electronics, and a 3D-printed holder and detector cover were modelled as well.

Particle tracks modelling
The goal was to model tracks of particles by means of scoring the deposited energy in each pixel of the sensor and to store all information for further data processing. To achieve this goal, a lattice (LAT card) was placed onto the whole sensor chip with 55 × 55 × Z µm 3 mesh, where Z is the actual sensor chip thickness (Fig. 3) that varies depending on the detector expected application (typically between 100 and 2000 µm).
The particle track file (PTRAC card) was used to write and store information about the energy deposited in each pixel. The PTRAC card in simulation for mixed radiation field generated by UHD proton beam was set up as follows: Events written into the file included any particle crossing surfaces 1 to 9 (specifically reserved for lattice mesh and sensor chip borders only) and any interaction, generation, and termination event in the sensor that resulted in energy deposition higher than 3 keV by any charged particle. Photons were included to correctly count also electrons generated by photons originating in the sensor, e.g., characteristic and bremsstrahlung photons. The VEL parameter (limit on particle velocity) allowed to omit writing mesh surface crossings of photons into the PTRAC file and thus reduce its size. For simulations with primary particles of a different type or different energy, the parameter TYPE was modified to match all particles that might occur in the detector. There is no information about particle track inside the pixel except for the information about particle termination or particle generation.
The parameters written in the PTRAC file relevant for data evaluation included particle type, event type, spatial coordinates of the event, direction vector, particle energy, and time.
Typically, 5 × 10 5 events were stored into the PTRAC file. The running time depended on the type of the source and the geometry.

MC data analysis
The PTRAC files generated by the MCNP® were processed in a purpose-written script in Matlab® software (The MathWorks®, USA). The individual steps of the analysis included: • spatial coordinate transformation to the default coordinates of the software package Pixet Pro [23]; • calculation of a) energy deposited in each pixel and b) effective centre of the track in each pixel; • application of the model of charge-sharing effect; • application of the detector energy resolution.
The deposited energy in each pixel was calculated as the difference between the particle energy when entering the pixel (or particle initial energy if created there) and particle energy when leaving the pixel (if it crossed the sensor surface).
The effective centre of the track in each pixel was determined as the average X, Y, Z coordinates of particle location when entering the pixel (or location of particle origin if created there) and particle location when leaving the pixel (or location of particle end of path if terminated there).

Model of charge sharing effect
Charge sharing effect in pixel detectors indicates the sharing of initial charge cloud generated by ionizing particle between two or more adjacent pixels caused by the lateral spread of the cloud because of electrical repulsion when traveling towards the read-out electronics accelerated by the electric field. The charge sharing in pixel detectors has been studied by many authors, e.g. [24][25][26]. The effect leads to the degradation of spatial and energy resolution of the detector [27] that is of high importance in medical applications where it causes reduction of dose efficiency [28] and image quality [29]. The charge sharing also affects the pixel detector energy calibration [30].
The analytical model of the expansion of the charge cloud generated in the sensor by ionizing radiation and the cloud's sharing between the adjacent pixels was developed (Fig. 4). It assumed the point generation of the initial charge. The model took into account electrical properties of the sensor material, bias voltage, depletion voltage, and the depth of particle interaction in the sensor.
The main principle relied on electric repulsion of the same charge within the charge cloud. The model assumed the initial charge distribution described by a Gaussian function with a certain known dispersion (potentially energy dependent). The model showed that the Coulomb repulsion dominated the charge sharing process especially for energies larger than units of keV. The Coulomb repulsion dominates shortly (in a fraction of nanosecond) after the initial ionization and positive-negative charge separation. The charge distribution is shaped as homogeneous sphere with radius r: where R is material constant incorporating charge mobility and ionization energy for electron-hole production, E is energy of primary particle spent for ionization, t is charge collection time, and C is the initial charge size.
The charge diffusion is relevant in low energy range. The electrostatic induction in combination with charge sensitive preamplifier properties didn't contribute to the signal registered by pixels for deposited energy lower than about 1 MeV.
The model was applied on the simulated data. For each pixel with any energy deposition, the radius of the final charge cloud was calculated from the effective centre of the track. Then, the energy deposited at   6. Comparison between measured (above) and simulated (below; with charge sharing applied) pixel clusters from 241 Am radionuclide (main energies: gamma − 60 keV, X-ray 12-22 keV). Data for a Timepix3 detector with 1 mm Si sensor. a point in the simulation was spread over the final radius and recalculated into the pixels lying within this radius. The energy shared by each pixel due to the charge sharing effect was determined as the original total deposited energy multiplied by the relative volume of the charge cloud sphere in the given pixel. The relative volume was determined by random sampling of points inside the pixel and counting the hits inside the sphere, which proved to be the fastest method.
The energy resolution of the detector was applied on the simulated data after the charge sharing model had been employed. It describes the accuracy of the absorbed energy measurement and comprises the statistical noise, electronic noise, and operational drift components. The energy-dependent sensor-wise energy resolution function was measured out of the scope of this work for several Timepix3 detectors using selected suitable radionuclides emitting photons (X-ray fluorescence and γ) in the energy range from 14 keV to 662 keV. It was determined from the full-width-at-half-maximum of full-energy peaks in the measured energy spectra of clusters in the whole sensor area. The energy resolution was applied as follows: First, for each cluster, a random difference from the true total deposited energy was determined by random number generation from Gaussian distribution using sigma (variation from the mean) taken from the energy-dependent energy resolution function for the given deposited energy and for the given sensor (material, thickness, settings). Then, the difference was split among all pixels of the given cluster linearly with respect to the energy in each pixel.

Comparison of measured and simulated clusters in various fields
The method for pixel clusters modelling was initially tested on imaginary point isotropic sources of monoenergetic 50 MeV protons and 2.9 MeV electrons originating at the sensor normal at 1 cm distance from the sensor. Fig. 5 shows the visualization of particle tracks. Clusters created by protons are nearly straight and heading from the source. The cluster length increases with the increasing angle of the proton path from the sensor normal, intersecting more pixels. The electron tracks, on the contrary, are largely curved because of straggling at numerous interactions and therefore the direction to the source is not obvious. The clusters are not broken and their shape, neglecting charge sharing effects, corresponds with our experience.
The comparisons of simulated clusters and real measured clusters created in various particle fields are presented in Figs. 6-11. Fig. 6 shows the integral frames with clusters generated by photons from 241 Am quasi-point source in 1 mm thick Si sensor. The source was positioned in air at the distance of 10 cm from the detector on the sensor normal. Because of low energy of photons (<60 keV) and small radius of charge sharing, the clusters are typically composed of 1 to 3 pixels. Fig. 7 shows the comparison of the measured and simulated detector pulse-height spectra, where each pulse represents the sum of energy deposited in all pixels of one cluster. Energy resolution of the given detector was applied on the simulated data. A good agreement was achieved. One of the main sources of a minor discrepancy was caused by the energy calibration of the detector. It can be observed in the position of the main two peaks (~14 keV and 60 keV), which should be in the ideal case in accordance with the simulation.
A similar result is depicted in Fig. 8 for 137 Cs source. The experimental setup was the same as in the case of 241 Am source. Some clusters are longer because of higher maximum photon energy emitted by 137 Cs (662 keV) compared to 241 Am. In addition, the 137 Cs also emits beta radiation with the maximum energy of 1175 keV (5.6 % yield). Higher   Fig. 7. Relative comparison between measured (red) and simulated (blue) energy spectrum of clusters for 241 Am photon source. Charge sharing and energy resolution were applied on simulated data. Data for a Timepix3 detector with 1 mm thick Si sensor. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 8. Comparison between measured (above) and simulated (below; with charge sharing applied) pixel clusters from 137 Cs radionuclide (main energies: gamma − 662 keV, X-ray 32-37 keV). The beta radiation was included in simulated data. Data for a Timepix3 detector with 1 mm thick Si sensor. energy beta particles also reached the detector because there was no absorber: there was no detector input window and the 1.5 mm thick polystyrene cover of the point source did not completely absorb the electrons. Beta radiation was included in the simulation. Fig. 9 compares the measured and the simulated maximum deposited energy in a pixel in cluster. Fig. 10 presents the comparison between the measured and the simulated clusters resulting from thermal neutrons. Neutron detection was possible using the Timepix3 detector equipped with suitable neutron convertors [31]. For thermal neutron detection, a convertor made from Lithium Fluoride enriched by 6 Li nuclide was used in a form of a thin (~10 µm) layer placed above a part of a sensor (~300 µm distance). In Fig. 10  Thermal neutrons interact with 6 Li creating triton and alpha particle ions. The ions are emitted in the direction heading away from each other and one of them usually reaches the detector and generates a thick circular cluster. Fig. 10 shows such tracks recognised as ions created by the thermal neutron interactions in the converter.
Thermal neutron detection efficiency strongly depends on the 6 LiF layer properties (thickness, distance, enrichment, homogeneity, glue concentration, etc.) and each piece of the converter is unique and has to be calibrated individually together with the given piece of the detector/ sensor. That is why the parameters of MC model of the converter (thickness, distance, area) have to be always tuned to match the measured detection efficiency (number of clusters per cm 2 per unit fluence of thermal neutrons) of the given device/converter. The presented results were obtained for 500 µm thick Si sensor with a converter providing the thermal neutron detection efficiency of about 1.6 % that is similar to the data presented in [31]. Fig. 11 shows the comparison between the measured and the simulated clusters created by 220 MeV proton beam delivered to a water phantom at a therapeutic proton therapy facility. Specifically, a complex mixed radiation field was measured and simulated on the beam axis about 4 cm behind the Bragg peak. The frames show a large number of different cluster types generated by various particlesthick circular and long tracks from recoiled protons, thin long tracks from high energy electrons, and short tracks from low energy electrons.
There was a noticeable difference in the shape of large circular clusters created by low energy (E <~10 MeV) recoiled protons. Small amount of energy was registered also by pixels farther from the interaction point creating a corona effect caused by local distortion of electric field caused by a huge amount of charge generated at one place. Such effect was not observed in simulated clusters because this phenomenon was not described by the charge sharing model.
To provide accurate simulations of the clusters made by scattered particles to characterise the mixed radiation fields created in therapeutic proton and electron beams is the main goal of this work. In the case presented in Fig. 11, the radiation field was composed by primary protons, recoiled nuclei, secondary photons and neutrons generated in proton nuclear interactions, and tertiary charged particles generated by secondary photons and neutrons. Fig. 9. Relative comparison between measured (red) and simulated (blue) spectrum of maximum deposited energy in a pixel in cluster for 137 Cs photon source. The charge sharing and energy resolution were applied on simulated data. Beta radiation was included in simulations. Results for Timepix3 detector with 1 mm thick Si sensor. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 10. Comparison between measured (above) and simulated (below; with charge sharing applied) clusters generated by tritons and alpha particles originating from 6 Li disintegration from interaction of thermal neutrons. The 6 LiF layer was placed over the top part of the sensor at the distance of about 300 µm. Different areas were covered in measurement and simulation. Data for a Timepix3 detector with 500 µm thick Si sensor.

The charge sharing model summary
The analytical model of charge sharing effect has recently been developed for point deposition of the charge, assuming the spherical expansion of the charge cloud based on the electric repulsion. That is applicable especially for low energy electrons (E <~50 keV). The model simplifies charge sharing a) for long tracks where the charge is generated continuously along the path of the charged particle, for example, high energy protons (E >~50 MeV) or electrons (E >~200 keV), and b) for very high energy depositions at one place, for example low energy protons (E <~10 MeV). The model can be improved also by including other physical effects, such as the charge diffusion, recombination, or the electrostatic induction.
However, the presented results are sufficient for the first approximation of the charge sharing effect in the Timepix-family semiconductor pixel detectors. The model will be further studied and improved in the next work.

Conclusions
In the MCNP radiation transport code, Monte Carlo model of MiniPIX Timepix3-FLEX pixel detector was developed. The model was voxelized into 256 × 256 pixels and allows simulating per-pixel deposited energy and storing the data for further analysis. A model of charge sharing between adjacent pixels was adapted for the use on simulated data to realistically simulate the shape of pixel clusters depending on the particle type, deposited energy, interaction depth, sensor parameters, and measurement settings. The method helps to optimize the sorting of measured clusters to distinguish radiation components especially in conditions where it is not possible to determine it experimentally. Ultimately, the developed method may improve the measurement of equivalent dose in mixed particle fields by Timepix3 detectors.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.