A Lightweight Algorithm to Model Radiation Damage Effects in Monte Carlo Events for High-Luminosity Large Hadron Collider Experiments

Radiation damage significantly impacts the performance of silicon tracking detectors in Large Hadron Collider (LHC) experiments such as ATLAS and CMS, with signal reduction being the most critical effect; adjusting sensor bias voltage and detection thresholds can help mitigate these effects, generating simulated data that accurately mirror the performance evolution with the accumulation of luminosity, hence fluence, is crucial. The ATLAS and CMS collaborations have developed and implemented algorithms to correct simulated Monte Carlo (MC) events for radiation damage effects, achieving impressive agreement between collision data and simulated events. In preparation for the high-luminosity phase (HL-LHC), the demand for a faster ATLAS MC production algorithm becomes imperative due to escalating collision, events, tracks, and particle hit rates, imposing stringent constraints on available computing resources. This article outlines the philosophy behind the new algorithm, its implementation strategy, and the essential components involved. The results from closure tests indicate that the events simulated using the new algorithm agree with fully simulated events at the level of few %. The first tests on computing performance show that the new algorithm is as fast as it is when no radiation damage corrections are applied.


Introduction
Silicon radiation detectors are at the core of experiments conducted at high-energy colliders, such as the CERN Large Hadron Collider (LHC) [1].Multiple layers of pixel and strip detectors are utilized to accurately measure the momentum and vertices of charged particles.Examples of such detector systems are part of the CMS [2] Silicon Tracker [3] and the ATLAS [4] Inner Detector [5].
Hybrid pixel detectors, positioned close to the beam's interaction point, offer fine segmentation of a few tens of micrometers, ensuring spatial resolution on a single point down to 10 µm.Precision measurements and the discovery potential heavily rely on the performance of tracking devices, particularly for energy/momentum resolution and flavor identification.Notably, both ATLAS [6,7] and CMS [8] pixel detectors are equipped with thin planar n + -on-n silicon sensors (ATLAS Innermost B-Layer includes sensors made in 3D n + -on-p technology in the forward region).
The LHC is scheduled for an upgrade to a high-luminosity machine (HL-LHC) [9], with the goal of integrating a dataset of 4000 fb −1 ; this is twenty times larger than its current capacity, allowing for the full exploitation of the physics reach of the LHC machine.To integrate such a large dataset within a reasonable time, the instantaneous luminosity will need to increase by a factor of 5-7 beyond the nominal value, reaching 7 ×10 34 cm −2 s −1 .This increase in instantaneous luminosity will lead to a corresponding rise in particle flux, track densities, and hit rates.Consequently, the current inner tracking systems will need to be replaced to cope with this challenging environment, as radiation levels affecting both sensors and electronics will also increase by a similar amount.
The effects of radiation damage are already measurable at the fluences of LHC Run 2 and Run 3.For instance, the charge collection efficiency (CCE) of the innermost ATLAS Pixel layer has now reduced to approximately ∼70% at a fluence of about 1 × 10 15 n eq cm −2 [10].
The ATLAS, CMS, and LHCb collaborations have developed algorithms to incorporate radiation damage effects into their Monte Carlo (MC)-simulated events [11].Thanks to these algorithms, it is possible to track the evolution of pixel detector performance as radiation damage accumulates.Despite the different approaches adopted by the three collaborations, the level of agreement between data and Monte Carlo (MC) simulations is remarkable, achieving agreement at the O(%) level [10].These algorithms enable the prediction of operational parameters such as depletion voltage and make possible the generation of samples for future detector conditions; these are essential for training reconstruction algorithms.
Despite the importance and success of these algorithms, some may not be viable for the HL-LHC phase due to computing resource constraints.This is particularly true for the ATLAS Pixel detectors.Therefore, there is a pressing need for faster yet equally precise algorithms to meet the demands of the HL-LHC data-taking conditions.
In this paper, a new lightweight algorithm is presented which is able to address the challenges posed by the HL-LHC on computing resources for the inclusion of radiation damage effects in MC simulated events.The paper is organized as follows.
In Section 2, the ATLAS algorithm for simulating radiation damage effects to silicon pixel sensors is presented along with some key results.Section 3 introduces the new algorithm, which relies on Look-Up Tables (LUTs).This section also discusses the various types of simulations-TCAD (Technology-Computer-Aided Design) and Geant4 [12][13][14] based-that serve as inputs to the new algorithm.
In order to validate the new algorithm, a series of closure tests were conducted, comparing events simulated using the LUTs algorithm to fully simulated events.The results of these closures tests are reported in Section 4. Finally, the paper concludes (Section 5) with a discussion on the results and an outlook for the project.

The ATLAS IBL and Pixel Radiation Damage Digitizer
The ATLAS Collaboration developed and implemented an algorithm to include the effects of radiation damage to silicon pixel sensors in the Monte Carlo-simulated events.In the following (Section 2.1), the ATLAS Pixel and IBL detectors will be briefly discussed.The radiation damage digitizer for these detectors will be presented in Section 2.2.

The ATLAS IBL and Pixel Detector
The ATLAS Pixel detector [7] has been in operation since the beginning of LHC Run1.It comprises three barrel layers and three disks on each side.The barrel layer closest to the interaction point (IP) is called the B-Layer and is positioned at a radius R = 50.5 mm from the IP.All layers are composed of hybrid pixel modules, with sensors that are planar n + -on-n, 250 µm thick, and have a pitch of 50 × 400 µm 2 (bending and transverse planes, respectively).
For Run 2 of the LHC, which started in 2015, a new pixel layer was added: the Insertable B-Layer (IBL) [6], positioned 33 mm from the IP.The IBL is also composed of hybrid pixel modules.Compared to those in the Pixel detector, IBL modules feature thinner sensors, 200 µm instead of 250 µm, and a finer pitch of 50 × 250 µm 2 .Planar n + -on-n technology is used for IBL sensors, except in the forward region, where novel 3D n + -on-p sensors are employed.
By the end of 2023, the IBL (B-Layer) planar sensors had integrated a fluence of approximately 1.2 (1.0) ×10 15 n eq cm −2 .These fluence levels led to a signal loss of about 30% for both detectors.

The Radiation Damage Digitizer for Pixel and IBL Detectors
Signal reduction is the most important radiation damage effect for the performance of silicon pixel tracking detectors in ATLAS.It is important to have simulated data that reproduce the evolution of performance with the increase in radiation fluence.ATLAS collaboration developed and implemented an algorithm that reproduces signal loss and changes in Lorentz angle due to radiation damage [15].This algorithm is now the default for Run 3-simulated events [16].
The algorithm uses electric field profiles produced using TCAD tools to evaluate the time carriers will take to reach the collecting electrode.The electric field is deformed by the presence of defects due to radiation damage; in particular, it has a deep minimum in the sensor mid-plane [17].This has two main effects: first and most important, carriers will have reduced velocity and so more prone to be trapped; second, the Lorentz angle will become dependent on the position along the bulk.The combination of these effects determines a reduction in the CCE, a deformation of the cluster shape, and, eventually, a degradation of the space-point resolution.
In Figure 1, a comparison between the data and the simulated MC events for the CCE in IBL planar sensors is shown.Comparison of the evolution of CCE with integrated luminosity for IBL planar sensors in data (points) and MC events (bands) [18]; the corresponding radiation fluence is also indicated.The vertical bands include the estimated uncertainty affecting the input parameters of TCAD radiation damage model and of the trapping constants.
The comparison covers almost two order of magnitude of fluence and the level of agreement is excellent, down to 1 %.More results can be found in [10].
Despite the success of this algorithm, it cannot be used for the HL-LHC as it is excessively time-consuming.Therefore, a new, faster, yet equally precise algorithm is required, which will be presented in the next section.

A Radiation Damage Digitizer Based on Look-Up Tables for Silicon Detectors at HL-LHC Experiments
Starting in 2025 and extending through 2028, the LHC will be upgraded to the highluminosity LHC (HL-LHC), making the expansion of experimental particle physics research and the exploration of fundamental physics of the LHC possible.The HL-LHC will deliver an instantaneous luminosity of 7.5 × 10 34 cm −2 s −1 , a 5-7 times increase from the nominal value resulting in an average pileup of about 200 inelastic proton-proton collisions per bunch crossing.The innermost parts of the pixel detector, situated closest to the interaction point, will be exposed to unprecedented levels of radiation, reaching a non-ionizing radiation fluence of 1 − 2 × 10 16 n eq cm −2 during the expected lifetime of the ATLAS detector [19,20].
The actual ATLAS IBL, Pixel, and Strip detectors will not be able to cope with such harsh conditions.For this reason, the current detectors will be replaced by the ATLAS Inner Tracker (ITk), a full-silicon vertexing and tracking detector [19,20].The ATLAS ITk will have pixels detectors in its innermost part and strips at the outer radii.The layout of ATLAS ITk and that of its pixels part are presented in Figure 2. The ATLAS ITk pixels detector comprises five barrel layers (from L0 to L4) and a system of as many rings covering the forward region.The innermost pixel layer (L0) and ring (R0) will be equipped with 3D n + -on-p sensors.All the rest of the detector will be equipped with planar n + -on-p sensors, with thicknesses of either 100 µm (in L1 and R1) or 150 µm (in all other layers and rings).
The anticipated increase in particle density at HL-LHC calls for a faster algorithm to model the effects of radiation damage in the ATLAS MC events.A novel methodology centered on charge re-weighting from Look-Up Tables (LUTs) is under study and is outlined in the subsequent sections.

LUT Method
Similar to the actual radiation damage digitizer (see Section 2.2), the Look-Up Table (LUT) method simulates radiation damage effects in ATLAS MC events after Geant4 [12][13][14] simulates charge deposition and before signal digitization, where signal is discriminated and transformed into digital format by simulating the response of the readout electronics chip.
The task of the radiation damage digitizer is to estimate on which pixels the deposited charges will induce a signal and how large this signal will be.The LUT method is proposed to accomplish this task for the data taking of ATLAS during the high-luminosity phase of LHC (HL-LHC).ATLAS MC events will be simulated assuming an ideal, undamaged pixel sensor.The impact of radiation damage will then be incorporated as a correction, utilizing the LUTs.
The LUT method presented herein is inspired by the "template" method used by the CMS collaboration [22].The primary distinction between the CMS and LUT methods is that, while the former applies corrections to simulated clusters of pixels, the latter aims to reweigh observables of a group of drifting carriers in the silicon bulk.
The design of the LUT method to model the radiation damage effects in ATLAS MC events is based on the dynamics of charge carriers in the bulk of pixel sensors.The situation is sketched in Figure 3 for the case of planar sensors.
A schematic diagram of carrier dynamics in silicon planar pixel sensors.As a MIP crosses the sensor (at an angle φ trk ), electron-hole pairs are created and transported to their respective electrodes under the influence of electric and magnetic fields.Electrons and holes may be trapped before reaching the electrodes, but still induce a charge on the primary and neighbor electrodes.
Here, the sensor is assumed to be a planar pixel sensor of thickness w; y is the direction of the magnetic field ⃗ B and the direction of the electric field ⃗ E is along the bulk depth axis z; carriers are deflected in the x direction due to magnetic field.
Electrons and holes produced by a minimum ionizing particle (MIP) traversing the sensor drift toward collecting electrodes at an angle known as the Lorentz angle (θ LA ) with respect to the electric field.
In pristine detector-i.e., without radiation damage effects-the final position of the carriers is calculated projecting the carrier position to the respective collecting electrode.The final position is then smeared to take into account diffusion.The LUT method is intended to apply corrections to this basic method.In the following the initial position of the charge carrier will be referred to as deposited position⃗ r dep = (x dep , y dep , z dep ), and the final one the propagated position⃗ r prop = (x prop , y prop , z prop ); see also Figure 3.
In Figure 3, drifting carriers can be trapped before reaching the collecting electrode due to radiation damage in the sensor bulk.In such cases, they induce a signal that is only a fraction, k < 1, of their charge.This fraction can be calculated using the Shockley-Ramo theorem [23,24], which requires the initial and final positions of the carrier as input.
The carriers that are trapped cover a distance ∆z = z prop − z dep along the bulk depth direction, which depends on the ⃗ r dep position of the carriers.This dependence on the carrier generation position is due to the non-uniform electric field in the sensor bulk after radiation damage [10].As mentioned above, the carriers will drift at an angle θ LA relative to the direction of the electric field.This angle also depends on the ⃗ r dep position where the carrier was created, again due to the non-uniformity of the electric field.The carrier will end up at a position in the bending direction (orthogonal to the electric and magnetic field) that is displaced-on average, (the carrier movement is affected by diffusion too which adds dispersion to the final carrier position)-by ∆x = tan(θ LA ) • ∆z; of course, x prop = x dep + ∆x, if there was no diffusion.As discussed above, the carrier's final position ⃗ r prop will determine the amplitude of the signal induced on the electrodes.This final position depends on the carrier's initial position⃗ r dep and on the combination of electric and magnetic fields.The signal will be a fraction k of the carrier charge q, and this fraction k is also a function of the original carrier position⃗ r dep .
In summary, for each group of carriers, the method calculates the free path ∆z, the Lorentz angle θ LA , and the fraction of the induced signal k.These three quantities will be referred to as the observables.
To determine the values of the three observables, repeated simulations of the drift of carriers deposited at precise positions⃗ r dep within the bulk are conducted.These simulations are performed using the Allpix 2 [25] simulation framework, which will be presented in Section 3.2, using precise electric field and Ramo potential maps produced using TCAD (Technology-Computer-Aided Design) simulations.
For each simulated event, the initial position⃗ r dep , the final position⃗ r prop of the carrier will be saved, together with the signal fraction k induced on the pixel matrix.In principle, the values of three observables should be recorded as a function of the initial position⃗ r dep of the carrier, but due to the aforementioned limited computing resources, the LUTs will be a function of only z ≡ z dep .For a fixed z value, the average over all (x dep , y dep ) ≡ (x, y) positions will be carried out for each observable and assigned to that z value; to make the notation lighter, "dep" is dropped from here onward.Thus, for a single event and a single value of z, the three LUTs will be: Using the LUTs defined in 1, the propagated position⃗ r prop and the induced signal q ind of a charge, q deposited at depth z in the sensor bulk is calculated as follows: where x,y is a Gaussian-distributed random number added to simulate the effect of diffusion.
The essence of the LUT method is summarized in Equation ( 2): the precise dynamics of carrier drift are substituted with an "average" drift, and the same principle applies to the signal amplitude.
The procedure of charge deposition and drifting is repeated for each⃗ r dep several times in order to assess the dispersion of the carriers dynamics.
It is worth noting that, while this study primarily focuses on planar pixel sensors, the same methodology can also be readily applied to strip sensors.Additionally, for the sake of simplicity, the method will be presented here only for planar pixel sensors, but ongoing efforts are being made to extend it to 3D sensors as well-see, for example [10].
In the following sections, the Allpix 2 simulation framework, the TCAD simulations utilized as input to Allpix 2 simulations, and the process of calculating the Look-Up Tables (LUTs) will be presented.

Allpix-Squared for Radiation Damage Digitizer
Allpix 2 is a generic, open-source software framework for the simulation of silicon pixel detectors [25].The framework allows the user to create detailed simulations of the entire experimental chain of a testbeam, from incident radiation to digitized detector response.An extensible system of modules is responsible for executing the distinct simulation steps, such as realistic charge carrier deposition using the Geant4 toolkit and the propagation of charge carriers in silicon through a drift-diffusion model.Detailed electric field maps imported from TCAD simulations can be used to model the drift behavior of charge carriers within the silicon, introducing a higher level of realism to Monte Carlo-based simulations of particle detectors.
In this study, a planar n + -on-p pixel sensor with a thickness of 100 µm and a pitch of 50 × 50 µm 2 is simulated.This type of sensor is representative of what will be used in the second-to-innermost pixel layer (L1) of the ITk pixel detector.The sensor is simulated after irradiation corresponding to a fluence of 4 × 10 15 n eq cm −2 and operated at a voltage of 600 V.These conditions are expected for the aforementioned pixel layer at the end of the HL-LHC phase of ATLAS data taking.Figure 4 illustrates the simulation chain implemented for evaluating the LUTs, highlighting the use of electric field and Ramo potential maps from TCAD, along with the simulation of trapping effects.

Simulation Inputs
The electric field and Ramo potential maps calculated using TCAD tools can be inputted into Allpix 2 simulations to model the behavior of silicon radiation detectors, as outlined in Section 3.1.The following sections provide details of the TCAD simulations used for the calculation of LUTs in this study.All simulation results presented here have been obtained using Silvaco TCAD (https://silvaco.com/tcad/,accessed on 17 June 2024) tools.
The structure simulated in TCAD represents one-quarter of a 3 × 3 pixels matrix from a 100 µm n + -on-p planar pixel sensor with a pitch of 50 × 50 µm 2 ; refer to Figure 5a,c.The simulation of one pixel plus its neighbors is required to model charge sharing accurately via Ramo potential calculation.However, due to the symmetry of the design, it is sufficient to simulate just one-quarter of a 3 × 3 pixels matrix.To simulate the effects of radiation damage, the LHCb TCAD radiation damage model [26] was used.This model was developed for planar n + -on-p pixels intended for the LHCb Velo Upgrade [27], where irradiation fluences of several 1 × 10 15 n eq cm −2 are expected.Similarly, planar n + -on-p pixel sensors in the ATLAS ITk detector will be exposed to similar fluences during their expected operational lifetime.Therefore, the LHCb TCAD radiation damage model was chosen for this study.
In the ATLAS ITk detector, planar pixels will be utilized everywhere except in the innermost section of the pixel detector, where 3D sensors will be employed [19].
The largest fluence expected for planar pixels is Φ max ∼ 4 × 10 15 n eq cm −2 [28].Electric field maps were obtained for a sensor irradiated at the fluence Φ max once polarized at the maximal expected voltage (600 V).The extracted electric field profile along the bulk depth is shown in Figure 5a,b.In Figure 5a,c 2D projection of the weighting field is depicted; here, the result is presented for the full 3 × 3 pixels matrix as it is already in the format for Allpix 2 .
The electric field and weighting potential maps from TCAD are integrated into the Allpix 2 simulation chain, as illustrated in Figure 4, to generate the three LUTs introduced in Section 3.1 for correcting pixels response in ATLAS simulations.In the following section, the evaluation of LUTs using Allpix 2 will be presented.

Look-Up Tables Generation
Charged particles produced in LHC collisions will ionize the pixels sensors volume homogeneously.So in order to obtain a realistic representation of the pixel response, it is important to study it for all possible charge deposition locations ⃗ r dep .The LUTs are generated using the "scan" mode of the "DepositionPointCharge" module in Allpix 2 [25].This module places a specified quantity of charge carriers at a precise location within the detector's active volume.In the studies reported here, the scan mode was configured to deposit q dep = 1000 electron-hole pairs every 2 µm along the z direction and every 1 µm for the other two directions for the simulated device.
Once carriers have drifted to their final position⃗ r prop , the induced signal is evaluated for the pixel cell where they ended up and for the first eight neighbors, as illustrated in Figure 6.It is important to note that, while the propagated charge may land in the same pixel cell as the deposited charge, this is not always the case.(bottom) Illustration of the 3 × 3 matrix used to evaluate the induced signal; the central pixel is the one containing q prop ; hence,⃗ r prop .
In the following, the propagated position will always correspond to that of the electrons.This assumption is made in calculating the LUTs due to the negligible impact of holes on signal "position" compared to electrons.This simplification is based on the principle known as the "small-pixel effect" [29].In heavily segmented sensors such as pixel sensors, the weighting potential profile along the bulk has a steep gradient near the readout side.Consequently, only carriers in close proximity to the pixel electrode can induce a significant signal.Given that holes move away from the readout electrode, their influence on signal position can be safely neglected.
The pixel with the largest induced signal is identified, its position and signal q max are retained for the subsequent steps.The CCE per deposition position⃗ r dep is defined as q max divided by the deposited charge q dep : The CCE(z) is then evaluated by extracting the most probable CCE values (MPV) for each z.An example of CCE distribution for charges deposited at a z position close to the midplane is reported in Figure 7a.The ∆z(z) observable is evaluated by taking the average of the distribution of the free path traveled by the carriers as a function of the generation depth.This process is illustrated in Figure 7b.
While for the ∆z distribution, the average over repeated simulations is a good representative of the distribution itself, it is not the case for charge collection efficiency; for that, the most probable value is a better choice.
Finally, the Lorentz angle θ LA (z) value is calculated by fitting a straight line to the electron drift distribution ∆x vs. ∆z for each z dep position for various x dep and y dep positions, as shown in Figure 7c.The slope of the fit at each z dep is then extracted, and the LUT is constructed by plotting the slope (< tan(θ LA ) >) as a function of the generation depth, as illustrated in Figure 8c.
The large spread of ∆x values observed when electrons reach the electrode is an artifact of the simulation.At this position, electrons are already within the pixel implant and cease to diffuse further.To prevent biasing the results, the upper limit of the fit range for extracting θ LA (z) is constrained to a few micrometers away from the pixel implant.
Finally, Figure 8 displays the LUTs for charge collection efficiency (CCE(z)), average free path (∆z(z)), and tangent of Lorentz angle (tan θ LA (z)) for the case presented in this paper.The CCE peaks at over 85% at approximately 20 µm below the pixel implant.This phenomenon occurs because holes become trapped near the pixel, thereby partially screening a portion of the signal.However, the influence of holes becomes negligible when charge is deposited farther away from the pixel electrode.
The average distance ∆z covered by electrons becomes significantly less than the distance to the electrode rather quickly.For instance, at z = 20 µm, the average distance ∆z is approximately 15 µm, which is 25 % less.
The Lorentz angle deflection is predicted to be very small everywhere in the bulk compared to the typical value in the ATLAS sensors before irradiation-about 0.22 radians.The main reason for this is the large electric field inside the sensor, especially at the pixel side.This reduction in the Lorentz angle after irradiation has already been observed with collision data in ATLAS pixels [30].The largest deflection in x will be for charges deposited at the backside, and it will be just 2 µm.
The following section will present the first validation of the LUT algorithm.

Internal Validation of the LUTs Digitizer: Closure Test
A first validation of the new algorithm has been carried out through a closure test.As explained above in Allpix 2 , it is possible to implement different algorithms for each of the various steps of the simulation.In particular, after the energy deposition step (see Figure 4), the charge transport can be modeled in several ways, depending on the needed accuracy.
This feature was exploited for the purpose of our closure test by implementing a new charge transport algorithm based on LUTs, the "LUT propagator".This propagator is based on Equation (2).For every deposited electron the induced signal and final position is evaluated using the LUTs.
In the closure test presented, events simulated using the LUT propagator have been compared to fully simulated (FS) events.
Events were simulated using sensors of the same type and under conditions used to generate the LUTs, as discussed in Section 3. The agreement between LUT-based and full simulation (FS) events are assessed by examining cluster charge and cluster size, independently in the transverse (x) and longitudinal (y) directions (ATLAS uses a righthanded coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the z-axis coinciding with the axis of the beam pipe.The x-axis points from the IP towards the center of the LHC ring, and the y-axis points upward.Cylindrical coordinates (r,ϕ) are used in the transverse plane, ϕ being the azimuthal angle around the beam axis.The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2)).
For charge deposition, a beam of pions at three distinct energies is simulated, corresponding to transverse momentum (p T ) values of 1, 10, and 100 GeV/c.
In order to reproduce the configuration in which the simulated pixel sensors will be installed, the module is tilted around the beam axis by 0.25 radians and at different η values, from 0 (normal incidence) to 1.4 (about 1.1 radian with respect to the normal).
In the following, the results for cluster charge (Section 4.1) and for cluster size (Section 4.2) for p T = 1 GeV/c will be presented.Results will be summarized for all p T values in Section 4.3, accompanied by a brief discussion.

Results of the Closure Test on Cluster Charge
The cluster charge distributions from simulated events for p T = 1 GeV/c pions are shown in Figures 9-11 for η = 0, 1 and 1.4, respectively, for both the full simulation and the LUT-based simulation.

Results of the Closure Test on Cluster Size
The cluster size distributions in the transverse direction from simulated events for p T = 1 GeV/c pions are shown in   Despite an overall good agreement between the average of the cluster size distribution, some tension exists in the transverse distribution for η = 1.The difference is of the order of 18%, while for the other η values the agreement is within few %, as seen for the cluster charge.Similar behavior is observed for p T = 10 and 100 GeV/c.

Summary and Discussion on Closure Test
In Tables 1-3, the summary of the results of the closure tests is presented for all quantities and pseudorapidities for pions with p T = 1, 10 and 100 GeV/c, respectively.The discrepancy reported in Section 4.2 is evident for η = 1 for all p T values in transverse cluster size CS X .
In order to investigate these discrepancies between FS-and LUT-based events for η = 1 in transverse cluster size, more pseudorapidity values were explored for pions with p T = 100 GeV/c.Simulations were run at η = 0.4, 0.8 and 1.2.The results are reported in Table 4.For these values of η, the difference in the average cluster size in the transverse direction between FS-and LUT-based events is again within 1-4%.
These results indicate that the discrepancy observed at η = 1 is limited to a very narrow range around that value.In order to confirm that this the divergence of the beam for simulations at η = 1, it was increased to ±10 mrad, corresponding to a pseudorapidity range: η ∈ [0.99, 1.01].The results are reported in Table 5.
Table 5.Average cluster size in the transverse direction (CS X ) for various beam divergences at η = 1 and p T = 100 GeV/c; ϵ denotes the relative difference between FS and LUT.When the beam divergence is increased in the y direction, the level of agreement between LUT and FS events is much better than 1%.

Beam
In general, simulations with such beam divergence show an improvement in the agreement for all observables, reaching about 2 % in CS Y and less than 0.1 % in C Q .Since charged particles from LHC collisions are distributed uniformly in η from −2 to +2, the observed local discrepancy poses no problem at all.

Conclusions and Outlook
The extreme data rates and fluences expected at the HL-LHC impose stringent constraints on both detectors and computing resources.Radiation damage effects will significantly impair the performance of silicon pixel detectors in the ATLAS Inner Tracker.Thus, it is crucial to accurately replicate such degradation in MC events, as demonstrated by the current pixel detector and the existing radiation damage digitizer.
The pursuit of a faster digitizer while maintaining precision led to the development of the LUT digitizer.In a closure test, the new radiation damage digitizer demonstrated excellent agreement with fully simulated data across all ranges of pseudorapidity and transverse momentum.The discrepancy at η = 1 is confined to a range of less than 1% in pseudorapidity, as detailed in Section 4.3.
In the near future, data from ITk pixel modules equipped with the new readout chip will become available, allowing for the comparison of predictions from LUT-based simulations with testbeam data.
The LUT method will also be extended to 3D sensors [10].For these sensors, the LUT of CCE as a function of the radial distance of the carriers from the central electrode should suffice, as the Lorentz deflection is negligible.Efforts are underway to model edge areas for all types of sensors (planar, 3D, and strips).
Preliminary tests indicate that the LUT radiation damage digitizer operates as fast as the digitizer without radiation damage corrections.This was somewhat expected, as the computational load is similar for both.As shown in Equation ( 2), for unirradiated sensors, k = 1, ∆z = w − z, and tan θ LA is constant, resulting in no algorithmic difference.This also confirms that reading values from LUTs is not the dominant factor in the computing load.
In summary, the LUT method for simulating radiation damage effects is as precise as full simulation methods and as fast as algorithms that do not model radiation damage.The LUT method is expected to be used as the default for the generation of MC events for the HL-LHC phase of ATLAS.

Figure 1 .
Figure 1.Comparison of the evolution of CCE with integrated luminosity for IBL planar sensors in data (points) and MC events (bands)[18]; the corresponding radiation fluence is also indicated.The vertical bands include the estimated uncertainty affecting the input parameters of TCAD radiation damage model and of the trapping constants.

Figure 2 .
ATLAS ITk layout.(a) Full detector; pixels/strips parts are in red/blue.(b) Pixel part.The z coordinate is measured along the beam axis and r is the radial distance from the center of ATLAS, from [21].

Figure 4 .
Figure 4. Modified Allpix 2 simulation chain with a single detector for radiation damage digitizer.

Figure 5 .
(a) Depicts one-quarter of a 3 × 3 pixel matrix structure simulated in TCAD and (b) illustrates the extracted electric field distribution as a function of charge generation depth of the structure in (a) irradiated at a fluence of 4 × 10 15 n eq cm −2 at 600 V bias voltage.(c) A 2D projection of the weighting potential for the simulated structure.

Figure 6 .
Figure 6.(top) Illustration of deposited (propagated) charge q dep (q prop ) and position ⃗ r dep (⃗ r proop ).(bottom) Illustration of the 3 × 3 matrix used to evaluate the induced signal; the central pixel is the one containing q prop ; hence,⃗ r prop .

Figure 7 .
Distribution of observables for charge injected close to the sensor midplane.(a) Distribution of simulated CCE.(b) Distribution of simulated ∆z.(c) Visualization of electron drift ∆x vs. ∆z towards the pixel side.The red line represents the fit of a straight line to the electron drift, and the average tan θ LA is estimated as the slope of the fitted line.

Figure 9 .Figure 10 .
Cluster charge distribution for pions with p T = 1 GeV/c impinging at η = 0. (a) Full simulation; (b) LUT-based simulation.level of agreement between FS and LUT simulation is remarkable: the agreement between the mean values of the two distributions is well within 5%.Similar results are obtained for p T = 10 and 100 GeV/c.Cluster charge distribution for pions with p T = 1 GeV/c impinging at η = 1.(a) Full simulation; (b) LUT-based simulation.

Figure 12 .Figure 13 .Figure 14 .Figure 15 .
for η = 0, 1 and 1.4, respectively, for both full simulation and LUT-based simulation.Cluster size distribution in the transverse direction for pions with p T = 1 GeV/c impinging at η = 0. (a) Full simulation; (b) LUT-based simulation.Cluster size distribution in the transverse direction for pions with p T = 1 GeV/c impinging at η = 1.(a) Full simulation; (b) LUT-based simulation.Cluster size distribution in the transverse direction for pions with p T = 1 GeV/c impinging at η = 1.4.(a) Full simulation; (b) LUT-based simulation.The cluster size distributions in the longitudinal direction from simulated events for p T = 1 GeV/c pions are shown in Figures 15-17 for η = 0, 1 and 1.4, respectively, for both full simulation and LUT-based simulation.Cluster size distribution in the longitudinal direction for pions with p T = 1 GeV/c impinging at η = 0. (a) Full simulation; (b) LUT-based simulation.

Figure 16 .Figure 17 .
Cluster size distribution in the longitudinal direction for pions with p T = 1 GeV/c impinging at η = 1.(a) Full simulation; (b) LUT-based simulation.Cluster size distribution in the longitudinal direction for pions with p T = 1 GeV/c impinging at η = 1.4.(a) Full simulation; (b) LUT-based simulation.

Table 1 .
Average cluster size in transverse (CS X ) and bending (CS Y ) direction, and average cluster charge C Q for FS-and LUT-based events for different pseudorapidity values η and p T = 1 GeV/c; ϵ is the relative difference between FS and LUT.

Table 2 .
Average cluster size in transverse (CS X ) and bending (CS Y ) direction, and average cluster charge C Q for FS-and LUT-based events for different pseudorapidity values η and p T = 10 GeV/c; ϵ is the relative difference between FS and LUT.

Table 3 .
Average cluster size in transverse (CS X ) and bending (CS Y ) direction, and average cluster charge C Q for FS-and LUT-based events for different pseudorapidity values η and p T = 100 GeV/c; ϵ is the relative difference between FS and LUT.

Table 4 .
Average cluster size in transverse direction CS X , for different pseudorapidity values η and p T = 100 GeV/c; ϵ is the relative difference between FS and LUT.