A fast simulation method of ion scattering track clusters driven by a bypass factor

Abstract

The application discloses a kind of fast simulation method of ion scattering track cluster driven by bypass factor, including constructing basic radiation parameter table;According to the basic radiation parameter table, search discrete energy point, form lookup table;According to the variation of bypass factor of each discrete energy point in lookup table, determine the scattering angle distribution characteristic parameter corresponding to discrete energy point, and generate scattering angle sample based on scattering angle distribution characteristic parameter and store in scattering angle sample pool;Initialize three-dimensional grid, set the initial position and direction of incident ion;Carry out incident ion transport and energy deposition cycle;Parallel execution of the simulation task of multiple incident ions, respectively in independent three-dimensional grid copy incident ion transport and energy deposition cycle, finally all three-dimensional grid copy is obtained Total dose distribution.This method achieves optimal balance between calculation efficiency and accuracy, so as to generate high-fidelity three-dimensional dose distribution chart.

Description

Technical Field

[0001] This invention belongs to the field of ion radiation and matter interaction technology, and particularly relates to a rapid simulation method for ion scattering track clusters driven by the bypass factor. Background Technology

[0002] Ion therapy, as a currently advanced radiotherapy technology, demonstrates significant advantages in the field of tumor treatment due to its unique Bragg peak physical properties. Compared with traditional photon (X-ray) therapy, ion beams can deposit the maximum dose in the tumor area while significantly reducing the radiation dose to normal tissues before and after the tumor. Therefore, it has irreplaceable application value in the treatment of pediatric tumors, tumors near critical organs, and recurrent tumors.

[0003] However, the precise implementation of ion therapy faces severe challenges. The transport of ion beams in human tissues involves complex physical processes such as energy deposition, multiple Coulomb scattering, and nuclear reactions. To ensure therapeutic efficacy and maximize the protection of normal tissues, highly accurate dose calculations and distribution predictions are essential before treatment. This places extremely high demands on the accuracy and efficiency of ion transport simulation technology. Current mainstream ion transport simulation technologies suffer from a core contradiction: "precision and speed are mutually exclusive," making it difficult to meet actual clinical needs and becoming a key bottleneck restricting the development of ion therapy technology.

[0004] Existing ion transport simulation techniques are mainly divided into two categories, both of which have significant limitations. The first category is Monte Carlo simulation methods based on first principles, represented by Geant4, FLUKA, and MCNP. These methods achieve high-precision simulations of ion-matter interactions through complete physical models, but the computational cost is enormous, and each simulation is extremely time-consuming, making them unsuitable for clinical treatment planning optimization or real-time applications requiring rapid iteration. The second category is simplified pencil beam algorithms, widely used to meet clinical speed requirements. These algorithms decompose the ion beam into independent pencil beams, calculate dose contributions using empirical formulas, and then sum them to obtain the total dose distribution. While this method improves computational speed, the physical model is overly simplified and cannot accurately characterize the dose distribution at the tail of the Bragg peak caused by multiple Coulomb scattering of ions, especially at heterogeneous media or interfaces where errors are significant.

[0005] Therefore, there is an urgent need for a fast simulation method for ion scattering track clusters driven by bypass factor that can significantly improve computational efficiency, achieve three-dimensional fine dose calculation, and is easy to parallelize. Summary of the Invention

[0006] To address the problems existing in the prior art, this invention proposes a rapid simulation method for ion scattering track clusters driven by the bypass factor.

[0007] The technical solution of the present invention is as follows:

[0008] A rapid simulation method for ion scattering track clusters driven by a bypass factor, wherein the bypass factor is used to characterize the scattering degree of incident ions in a target material, comprising:

[0009] A basic radiation parameter table reflecting the transport behavior of incident ions in a target material is constructed. The basic radiation parameter table includes at least the energy, projected range, and circumduction factor of the incident ions.

[0010] Based on the aforementioned basic radiation parameter table, search for discrete energy points where the detour factor changes as energy decreases, forming several adjacent discrete energy point pairs and a lookup table recording energy loss and detour factor changes;

[0011] Based on the change in the circumference factor of each discrete energy point in the lookup table, the scattering angle distribution characteristic parameters corresponding to the discrete energy point are determined. Based on the scattering angle distribution characteristic parameters, scattering angle samples are independently sampled from the corresponding polar angle distribution and azimuth angle distribution to pre-generate scattering angle samples and store them in the scattering angle sample pool.

[0012] Based on the target area geometry and the preset voxel resolution, initialize the three-dimensional mesh for energy deposition and set the initial position and orientation of the incident ions;

[0013] An incident ion transport and energy deposition cycle is performed. Based on the current energy of the incident ion, the lookup table and the scattering angle sample cell are consulted to determine the energy loss, scattering angle sample, and step size representing the difference between the density-normalized projected range corresponding to the current energy and the next energy in the discrete energy point pair to which the incident ion belongs. The energy loss is distributed to the three-dimensional mesh according to the path ratio. Based on the scattering angle sample and step size, the position and orientation of the incident ion are updated, and the energy loss is subtracted to update the energy of the incident ion, until the termination condition is met.

[0014] Multiple incident ion simulation tasks are executed in parallel, with each performing incident ion transport and energy deposition cycles in an independent 3D mesh replica. Finally, all 3D mesh replicas are merged to obtain the total dose distribution.

[0015] Output the incident ion track clusters and dose analysis results based on the total dose distribution.

[0016] Furthermore, the specific method for constructing the basic radiation parameter table reflecting the transport behavior characteristics of incident ions in the target material includes: establishing interpolation functions of the projected range and the circumduction factor with respect to energy, and constructing the basic radiation parameter table by performing density normalization processing on the projected range and combining it with the density of the target material.

[0017] Furthermore, the specific method for forming the lookup table includes: starting from the initial energy and searching with a predetermined step size, and recording the discrete energy point pair and its corresponding change in the detour factor and energy loss only when the absolute value of the change in the detour factor between adjacent discrete energy points is greater than a preset threshold.

[0018] Furthermore, the formula for calculating the characteristic parameters of the scattering angle distribution is as follows: In the formula, These are the characteristic parameters of the scattering angle distribution; This represents the change in the detour factor.

[0019] Furthermore, the specific method for independently sampling from the corresponding polar angle distribution and azimuth angle distribution, pre-generating scattering angle samples, and storing them in the scattering angle sample pool includes: for each discrete energy point E i Based on the change in its detour factor Calculate the corresponding scattering angle distribution characteristic parameters. ;

[0020] For each discrete energy point E i Pre-independently from the scale parameter is A large number of polar angles were extracted from the Rayleigh distribution. Sample and validate, from Extracting a large number of azimuth angles from a uniform distribution Samples, forming sample pairs And it is stored as a scattering angle sample in the scattering angle sample pool.

[0021] Furthermore, after the pre-generated scattering angle sample, the method further includes comparing the statistical characteristics of the generated scattering angle sample with the theoretical distribution based on the scattering angle distribution characteristic parameters. If the deviation exceeds the allowable range, the sample is regenerated.

[0022] Furthermore, the specific steps for the incident ion transport and energy deposition cycle are as follows:

[0023] S5.1) Based on the current energy of the incident ion and using a lookup table, find the corresponding discrete energy pair and energy loss, and randomly select a preset scattering angle sample from the scattering angle sample pool. As the scattering angle sample for the current step size;

[0024] S5.2) Using the basic radiation parameter table, calculate the density-normalized projected range difference between the current energy and the next energy, and use it as the step length for the current step;

[0025] S5.3) The three-dimensional digital differential analysis algorithm is used to distribute the energy loss of the incident ions deposited within the step size into the three-dimensional voxel grid according to the length ratio of the path through each voxel;

[0026] S5.4) Update the position and orientation of the incident ions based on the scattering angle sample and step size, and subtract energy loss to update the energy of the incident ions;

[0027] S5.5) Repeat S5.1) to S5.4) until the termination condition is met.

[0028] Furthermore, the allocation of energy loss to the three-dimensional mesh according to the path ratio is achieved using a three-dimensional digital difference analysis algorithm. The specific method includes: linearly allocating energy loss to the corresponding voxels according to the ratio of the length of the incident ion step segment in each voxel to the length of the bus segment.

[0029] Furthermore, the termination conditions include the energy of the incident ions being lower than the cutoff threshold or the incident ions penetrating the target area.

[0030] Furthermore, the parallel execution specifically involves: dividing the total number of ions to be simulated into multiple batches, and using a multi-process parallel computing framework, allocating the simulation task of each batch to different computing processes; each computing process has an independent copy of the three-dimensional mesh, and summing all the three-dimensional mesh copies after the simulation is completed.

[0031] Compared with the prior art, the present invention has the following beneficial effects:

[0032] This invention proposes a fast simulation method for ion scattering track clusters driven by the bypass factor. This method pre-emptively and parallelizes the most time-consuming real-time sampling process of random scattering angles in traditional Monte Carlo simulations. The method pre-generates and verifies an angular distribution sample pool of all possible scattering events independently based on the change in the bypass factor. During simulation, only rapid random reads are required, completely eliminating the sequential random number generation and function call overhead of traditional Monte Carlo methods. Combined with a lock-free parallel computing framework based on independent grid replicas, the simulation time for large-scale ion clusters is significantly reduced, with efficiency improvements reaching several orders of magnitude, meeting the real-time or near-real-time computational needs of clinical applications.

[0033] This invention can directly inherit validated ion transport data (energy, projected range, and bypass factor) from databases such as NIST PSTAR. By constructing an accurate interpolation function and using the change in bypass factor to deduce the root mean square of the scattering angle distribution, it ensures that the scattering process is statistically strictly equivalent to the real physical process. Therefore, this method can accurately simulate lateral dose broadening, the true Bragg peak shape, and its tail dose while enjoying extremely high computational speed, with accuracy comparable to traditional MC methods.

[0034] The method of this invention uses a three-dimensional digital differential analysis algorithm to linearly distribute the energy loss of the incident ion at each step according to the proportion of the length of the three-dimensional voxel grid through which the path passes. This algorithm directly traverses the affected voxels, achieving computational complexity linearly related to the number of voxels. Compared to the complexity or accuracy loss caused by coarse grid averaging that may exist in some approximation algorithms, this method achieves an optimal balance between computational efficiency and accuracy, thereby generating a high-fidelity three-dimensional dose distribution map. Attached Figure Description

[0035] Figure 1 This is a flowchart of a fast simulation method for ion scattering track clusters driven by the bypass factor according to the present invention;

[0036] Figure 2(a) shows the polar angle. A schematic diagram of the Rayleigh distribution;

[0037] Figure 2(b) shows the azimuth angle. A schematic diagram of a uniform distribution;

[0038] Figure 3 This is a three-dimensional track map of a single proton;

[0039] Figure 4(a) shows the change of polar angle with the number of steps;

[0040] Figure 4(b) is a schematic diagram showing the change of azimuth angle with the number of steps;

[0041] Figure 4(c) is a schematic diagram showing the change of energy with the number of steps;

[0042] Figure 4(d) shows the Rayleigh distribution parameters. A diagram illustrating the change with the number of steps;

[0043] Figure 5 A 3D visualization of 200 proton track clusters;

[0044] Figure 6(a) is a schematic diagram of the depth-dose curve (Bragg peak);

[0045] Figure 6(b) is a schematic diagram of the lateral average dose distribution at the Bragg peak (z=17.6cm);

[0046] Figure 6(c) is a schematic diagram of the average dose distribution in the XZ plane (y=1.0cm);

[0047] Figure 6(d) shows the final energy distribution. Detailed Implementation

[0048] The present invention will be further illustrated below with reference to the accompanying drawings and specific embodiments. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. After reading this invention, any modifications of the invention in various equivalent forms by those skilled in the art will fall within the scope defined by the appended claims.

[0049] Example 1:

[0050] This invention provides a rapid simulation method for ion scattering track clusters driven by a bypass factor, whereby the bypass factor characterizes the degree of scattering of incident ions in a target material, such as... Figure 1 As shown, the specific steps include the following:

[0051] S1. Construct a basic radiation parameter table that reflects the transport behavior of incident ions in the target material. The basic radiation parameter table shall include at least the energy, projected range and circumduction factor of the incident ions.

[0052] S2. Based on the basic radiation parameter table, search for discrete energy points where the detour factor changes as energy decreases, forming several pairs of adjacent discrete energy points and a lookup table recording energy loss and detour factor changes.

[0053] S3. Based on the change in the circumference factor of each discrete energy point in the lookup table, determine the scattering angle distribution characteristic parameters corresponding to the discrete energy point, and independently sample from the corresponding polar angle distribution and azimuth angle distribution based on the scattering angle distribution characteristic parameters, pre-generate scattering angle samples and store them in the scattering angle sample pool.

[0054] S4. Based on the target area geometry and the preset voxel resolution, initialize the three-dimensional mesh for energy deposition and set the initial position and orientation of the incident ions.

[0055] S5. Perform incident ion transport and energy deposition cycle. Based on the current energy lookup table of the incident ion and the scattering angle sample cell, determine the energy loss, scattering angle sample, and step size representing the difference between the density normalized projected range corresponding to the next discrete energy point in the discrete energy point pair to which the incident ion belongs. Distribute the energy loss to the three-dimensional mesh according to the path ratio. Update the position and direction of the incident ion according to the scattering angle sample and step size, and subtract the energy loss to update the energy of the incident ion until the termination condition is met.

[0056] S6. Simulate multiple incident ion tasks in parallel, with each performing incident ion transport and energy deposition cycles in an independent 3D mesh copy. Finally, merge all 3D mesh copies to obtain the total dose distribution.

[0057] S7. Output the incident ion track clusters and dose analysis results based on the total dose distribution.

[0058] Furthermore, the specific method for constructing a basic radiation parameter table that reflects the transport behavior characteristics of incident ions in the target material includes: establishing interpolation functions of the projected range and the circumduction factor with respect to energy, and constructing a basic radiation parameter table by performing density normalization on the projected range and combining it with the density of the target material.

[0059] Density normalization refers to converting the projected range data from surface density (g / cm²) to geometric length (cm) by dividing the projected range data by the target density.

[0060] Furthermore, the specific method for forming the lookup table includes: starting from the initial energy E0, searching at predetermined step sizes to identify all detour factors. The discrete energy point E changes i The discrete energy point pair (E) is recorded only when the absolute value of the change in the circumference factor between adjacent discrete energy points is greater than a preset threshold. i, E i+1 and the corresponding change in detour factor. and energy loss This forms an energy-bypass factor increment lookup table.

[0061] Furthermore, the formula for calculating the characteristic parameters of the scattering angle distribution is as follows: In the formula, These are the characteristic parameters of the scattering angle distribution; This represents the change in the detour factor.

[0062] Furthermore, specific methods for independently sampling from the corresponding polar angle distribution and azimuth angle distribution, pre-generating scattering angle samples, and storing them in the scattering angle sample pool include:

[0063] For each discrete energy point E i Based on the change in its detour factor Calculate the corresponding scattering angle distribution characteristic parameters. ;

[0064] Pre-generate parameters that conform to the scale. The polar angle sample set with Rayleigh distribution and the azimuth angle sample set with uniform distribution;

[0065] For each discrete energy point E i Pre-independently from the scale parameter is A large number of polar angles were extracted from the Rayleigh distribution. Sample and validate, from Extracting a large number of azimuth angles from a uniform distribution Samples, forming sample pairs And it is stored as a scattering angle sample in the scattering angle sample pool.

[0066] Furthermore, after pre-generating scattering angle samples, the statistical characteristics of the generated scattering angle samples are compared with the theoretical distribution based on the scattering angle distribution characteristic parameters. If the deviation exceeds the allowable range, the samples are regenerated.

[0067] Specifically, the theoretical mean of the Rayleigh distribution is calculated based on the scattering parameters of discrete energy points. and theoretical variance ; for the extracted polar angle Statistical analysis of the sample to determine the theoretical mean and theoretical variance ; Calculate the relative deviation of the mean between the actual statistic and the theoretical statistic. and variance relative deviation Set deviation threshold Automatic judgment will be performed: if and If the scattering angle sample quality of the current discrete energy point is deemed acceptable, the simulation process can proceed to the next step; if any deviation exceeds... If so, the scattering angle of the discrete energy point is resampled.

[0068] Furthermore, the specific steps for the incident ion transport and energy deposition cycle are as follows:

[0069] S5.1) Based on the current energy of the incident ion and using a lookup table, find the corresponding discrete energy pair and energy loss, and randomly select a preset scattering angle sample from the scattering angle sample pool. As the scattering angle sample for the current step size;

[0070] S5.2) Using the basic radiation parameter table, calculate the density-normalized projected range difference between the current energy and the next energy, and use it as the step length for the current step;

[0071] S5.3) The three-dimensional digital differential analysis algorithm is used to accurately allocate the energy loss of the incident ions deposited within the step size to the three-dimensional voxel grid according to the length ratio of the path through each voxel.

[0072] S5.4) Update the position and orientation of the incident ions based on the scattering angle sample and step size, and subtract energy loss to update the energy of the incident ions;

[0073] S5.5) Repeat S5.1) to S5.4) until the termination condition is met.

[0074] Furthermore, the energy loss is distributed to the three-dimensional mesh according to the path ratio using a three-dimensional digital difference analysis algorithm. The specific method includes: linearly distributing the energy loss to the corresponding voxel according to the ratio of the length of the incident ion step line segment in each voxel to the length of the bus segment.

[0075] Specifically, the energy loss of this step is linearly distributed to the corresponding voxels according to the proportion of the length of the line segment within each voxel to the total length of the line segment: In the formula, Energy deposited within a single voxel The length of the path within this voxel. This represents the total path length.

[0076] Furthermore, the termination conditions include the energy of the incident ions being below the cutoff threshold or the incident ions penetrating the target area.

[0077] Furthermore, parallel execution specifically involves dividing the total number of ions to be simulated into multiple batches, using a multi-process parallel computing framework to assign the simulation task of each batch to different computing processes; each computing process has an independent copy of the three-dimensional mesh, and after the simulation is completed, all three-dimensional mesh copies are summed.

[0078] Furthermore, the dose analysis results include at least one of the following: depth dose curve, lateral dose distribution, dose volume histogram, ion track cluster three-dimensional visualization, and dose distribution cross-sectional diagram.

[0079] This embodiment also applies the rapid simulation method of the present invention to a parallel proton beam containing 100,000 protons with an initial energy of 160 MeV. This parallel proton beam is incident perpendicularly on a pure water target (density 1.0 g / cm³). The dimensions of the water target in the XY plane are... The depth dimension (Z-axis) is dynamically determined based on the proton range.

[0080] The interpolation function used in this example includes ion energy, projected range, and bypass factor data that are generally available from publicly available databases. A piecewise cubic Hermite interpolation method is used on the raw data to establish an interpolation function that covers the range of ion incident energies and how these parameters change with energy.

[0081] Specifically, this embodiment loads proton transport data of water from the NIST PSTAR public database file "pstar_water.txt": raw data such as energies from 0.01 to 200 MeV, projected range (g / cm²), and bypass factor. Piecewise cubic Hermite interpolation is performed on the raw data, and the projected range data is density normalized and converted to centimeters. Three interpolation functions are constructed, and verification tests are conducted at energies of 0.01 MeV, 1.0 MeV, and 10.0 MeV. This is combined with the density of liquid water... It is encapsulated into a unified interface for subsequent simulated calls.

[0082] This example, using data from the basic radiation parameter table, starts with the incident energy E0 of the ions and searches towards lower energies in predetermined step sizes to identify all energy points E where the orbital factor changes. i Calculate the change in detour factor: Only when At that time, among them It is a preset threshold for the change in the detour factor; record this energy pair (E i E i+1 ), change in detour factor and energy loss The continuous scattering process is discretized into a series of... The scattering events characterized lay the foundation for subsequent parallel presampling.

[0083] Specifically, in this embodiment, the incident energy of each proton is... The predetermined step size is 0.01 MeV, and the preset threshold for the change in the detour factor is... The energy was set at 1e-15, starting from 160 MeV, and searching towards lower energies in 0.01 MeV steps. A total of 192 points of change in the circumduction factor were found, and the energy pairs of each scattering event (E) were recorded. i E i+1 ), change in detour factor and energy loss .

[0084] Step S3, the construction of the scattering angle sample pool, includes: for each obtained scattering energy point Ei, based on its circumference factor change... Calculate the characteristic parameters of the scattering in this step: ;

[0085] Each Ei is pre-determined independently from the scale parameter. The number of polar angle samples drawn from the Rayleigh distribution is twice the number of incident ions.

[0086] All sampled angles are examined to see if they conform to a Rayleigh distribution. Based on the characteristic parameters of the scattering energy point, the theoretical mean of the Rayleigh distribution is calculated. and theoretical variance : ; ;

[0087] For the extracted polar angle Statistical analysis of the sample: ; ; This indicates the calculation of the sample mean; This indicates the calculation of sample variance; Indicates the polar angle of extraction The sample data set;

[0088] Calculate the relative deviation of the mean between the actual statistic and the theoretical statistic. and variance relative deviation : ; ;

[0089] Set deviation threshold ,like and If the scattering angle sample quality at the current energy point is deemed acceptable, the simulation process can proceed to the next step; if any deviation exceeds... If so, the scattering angle of that energy point is resampled.

[0090] A large number of azimuth angle samples are extracted from the uniform distribution and paired with the verified polar angle samples to form sample pairs, which are then stored in the scattering angle sample pool.

[0091] Specifically, in this embodiment, the characteristic parameters calculated at an energy of 160 MeV are... If the number of incident protons is N=100,000, then the total number of samples is There are 192 points of change in the orbital factor. For each energy point, 200,000 polar angle samples are extracted. The theoretical and actual mean and variance are calculated for each energy point's polar angle samples. The relative deviation of the mean and variance between the actual and theoretical values ​​is also calculated, and a deviation threshold is set. The value is 10%. The theoretical mean and variance calculated for 160MeV energy are 0.002097 and 0.000001, respectively. The actual sample mean and variance are 0.002099 and 0.000001, respectively. The relative deviation of the mean is 0.09% and the relative deviation of the variance is 0.32%, both less than the deviation threshold. If the result is positive, the verification is successful. A large number of azimuth angle samples are extracted from the uniform distribution and paired with the verified polar angle samples to form sample pairs, which are then stored in the scattering angle sample pool.

[0092] Please refer to Figures 2(a) and 2(b), which are verification diagrams showing that the polar angle sample with a proton energy of 160 MeV conforms to a Rayleigh distribution and the azimuth angle sample conforms to a uniform distribution in this embodiment. In Figure 2(a), the vertical axis represents the probability density and the horizontal axis represents the polar angle. (radians), red curve ( The blue bar represents the theoretical distribution curve, and the blue bar represents the sampling results. In Figure 2(b), the vertical axis represents the probability density, and the horizontal axis represents the azimuth. (radians), the red dashed line represents the theoretically uniform distribution (1 / (2π)). 0.1592), the green bars represent the sampling results.

[0093] In this example, a three-dimensional zero-value array is initialized as an energy deposition grid based on the target area geometry and the preset voxel resolution, and the incident point and initial direction of the ions are defined.

[0094] Specifically, in this embodiment, the water target is set to X = 2.0 cm, Y = 2.0 cm, and the thickness Z is 1.5 times the density-normalized projected range data corresponding to the proton incident energy: 26.45 cm; the ion incident point is (1.00, 1.00, 0.00) cm, and the incident direction is along the positive Z-axis; the preset voxel resolution is 0.02, and the spatial range is X [0.00, 2.00] cm, Y [0.00, 2.00] cm, Z [0.00, 26.45] cm, then the grid latitude is... Individual factors.

[0095] This example of single-ion transport and energy deposition simulation includes:

[0096] S5.1) Based on the current energy of the incident ion and using a lookup table, find the corresponding discrete energy pair and energy loss, and randomly select a preset scattering angle sample from the scattering angle sample pool. As the scattering angle sample for the current step size;

[0097] S5.2) Using the basic radiation parameter table, calculate the density-normalized projected range difference between the current energy and the next energy, and use it as the step length for the current step;

[0098] S5.3) The three-dimensional digital differential analysis algorithm is used to accurately allocate the energy loss of the incident ions deposited within the step size to the three-dimensional voxel grid according to the length ratio of the path through each voxel.

[0099] S5.4) Update the position and orientation of the incident ions based on the scattering angle sample and step size, and subtract the energy loss to update the energy of the incident ions;

[0100] S5.5) Repeat S5.1) to S5.4) until the termination condition is met.

[0101] This example of a 3D digital differential analysis algorithm includes: traversing all voxels traversed by the proton stepping segment, and linearly distributing the energy loss of this step to the corresponding voxel according to the proportion of the segment's length within each voxel to the total length of the entire segment. Specifically, the 3D motion path of the proton is parameterized to determine the parameters when the path enters each voxel. and parameters at departure ,according to Calculate the relative length of the path within the voxel; compare the parameters of the path reaching the voxel boundaries in each direction. , , The next voxel boundary to be traversed and the corresponding voxel index update direction are determined; energy allocation is performed using a linear scaling model. In the formula, Energy deposited within a single voxel The length of the path within this voxel. This represents the total path length.

[0102] Specifically, in this embodiment, the initial state of a single proton is: incident energy of 160 MeV, incident position of (1.0, 1.0, 0.0), and a random group of samples is selected from the scattering angle sample pool based on the incident energy. Using the sample as the scattering angle for the first step, and employing the basic radiation parameter table from step S1, the density-normalized range difference between 160 MeV and 155.1 MeV, calculated as 17.63 cm - 16.70 cm = 0.93 cm, is taken as the step length for the first step. The energy loss of protons deposited within the first step length is then considered. The ion positions are precisely allocated to the three-dimensional voxel grid according to the length ratio of each voxel through which the path passes; the ion positions and energies are updated according to the scattering angle and step length; S5.1) to S5.4) are repeated until the ion energy is below the cutoff threshold or passes out of the target area.

[0103] This example divides the total number of ions to be simulated into multiple batches, including: a preset batch size BATCH_SIZE, and the total number of ions N divided into multiple batches based on the preset batch size. Batch size... for: Each batch contains no more than BATCH_SIZE of ions, and the last batch may contain fewer than BATCH_SIZE of ions.

[0104] Create a process pool. The number of worker processes in the process pool is dynamically adjusted based on the number of CPU cores of the computing device. It is usually set to the number of CPU cores minus 1 to make full use of computing resources while avoiding system overload.

[0105] Each batch of simulation tasks is assigned to worker processes in the process pool. Each worker process independently performs the following operations: initializes a local energy grid with the same size as the global energy grid, with an initial value of zero; iteratively simulates each ion within the batch, performing a single-ion 3D Monte Carlo simulation for each ion, including processes such as scattering and energy deposition, and records the energy deposition in the local energy grid, while also recording the ion's final energy, final depth, total path length, etc.; after completing the simulation of all ions in the batch, it returns the local energy grid and statistical information for the ions in the batch.

[0106] The main process collects the results returned by all worker processes, adds up the local energy grids returned by each worker process to obtain the global energy grid, and merges the final energy, final depth, total path length, etc. of all ions into global statistics.

[0107] Specifically, in this embodiment, the total number of protons N = 100,000, and the preset batch size BATCH_SIZE = 5,000, then... =100,000 / 5,000=20, a total of 20 batches. The computing device simulating this embodiment has 24 CPU cores, therefore the number of worker processes is 23.

[0108] The output and visualizations in this example include a single-ion 3D track map, a scattering angle evolution map, a 3D visualization of ion track clusters, a 3D dose distribution map, and a scattering angle distribution verification map; the scattering angle evolution map includes the scale parameters of polar angle, azimuth angle, energy, and Rayleigh distribution. The dose distribution varies with the number of steps; the three-dimensional dose distribution map includes depth dose curves, lateral dose distribution, and dose volume histograms. Please refer to [link to relevant documentation]. Figure 3 This is a three-dimensional track map of a single proton in this embodiment (actual scale), where the water target size is... Please refer to Figures 4(a) to 4(d), which show the evolution of the scattering angle in this embodiment. Figure 5 Figure 6(a) to Figure 6(d) is a three-dimensional visualization of 200 proton track clusters in this embodiment, with 200 protons. Please refer to Figures 6(a) to 6(d), which are three-dimensional dose distribution maps of this embodiment. In Figure 6(a), the vertical axis represents energy deposition (Mev), and the horizontal axis represents depth (cm).

[0109] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A rapid simulation method for ion scattering track clusters driven by a bypass factor, wherein the bypass factor is used to characterize the scattering degree of incident ions in a target material, characterized in that, include: A basic radiation parameter table reflecting the transport behavior of incident ions in a target material is constructed. The basic radiation parameter table includes at least the energy, projected range, and circumduction factor of the incident ions. Based on the aforementioned basic radiation parameter table, search for discrete energy points where the detour factor changes as energy decreases, forming several adjacent discrete energy point pairs and a lookup table recording energy loss and detour factor changes. Based on the change in the circumference factor of each discrete energy point in the lookup table, the scattering angle distribution characteristic parameters corresponding to the discrete energy point are determined. Based on the scattering angle distribution characteristic parameters, scattering angle samples are independently sampled from the corresponding polar angle distribution and azimuth angle distribution to pre-generate scattering angle samples and store them in the scattering angle sample pool. Based on the target area geometry and the preset voxel resolution, initialize the three-dimensional mesh for energy deposition and set the initial position and orientation of the incident ions; An incident ion transport and energy deposition cycle is performed. Based on the current energy of the incident ion, the lookup table and the scattering angle sample cell are consulted to determine the energy loss, scattering angle sample, and step size representing the difference between the current energy and the density-normalized projected range corresponding to the next discrete energy point in the discrete energy point pair to which the incident ion belongs. The energy loss is distributed to the three-dimensional mesh according to the path ratio. Based on the scattering angle sample and step size, the position and orientation of the incident ion are updated, and the energy loss is subtracted to update the energy of the incident ion, until the termination condition is met. Multiple incident ion simulation tasks are executed in parallel, with each performing incident ion transport and energy deposition cycles in an independent 3D mesh replica. Finally, all 3D mesh replicas are merged to obtain the overall dose distribution. Output the incident ion track clusters and dose analysis results based on the total dose distribution.

2. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 1, characterized in that, The specific method for constructing a basic radiation parameter table reflecting the transport behavior of incident ions in the target material includes: establishing interpolation functions of the projected range and the circumduction factor with respect to energy, and constructing the basic radiation parameter table by performing density normalization processing on the projected range and combining it with the density of the target material.

3. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 2, characterized in that, The specific method for forming the lookup table includes: starting from the initial energy and searching with a predetermined step size; only when the absolute value of the change in the bypass factor between adjacent discrete energy points is greater than a preset threshold, recording the discrete energy point pair and its corresponding change in bypass factor and energy loss.

4. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 3, characterized in that, The formula for calculating the characteristic parameters of the scattering angle distribution is as follows: In the formula, These are the characteristic parameters of the scattering angle distribution; This represents the change in the detour factor.

5. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 4, characterized in that, The specific method for independently sampling from the corresponding polar angle distribution and azimuth angle distribution, pre-generating scattering angle samples, and storing them in the scattering angle sample pool includes: for each discrete energy point E i Based on the change in its detour factor Calculate the corresponding scattering angle distribution characteristic parameters. ; For each discrete energy point E i Pre-independently from the scale parameter is A large number of polar angles were extracted from the Rayleigh distribution. Sample and validate, from Extracting a large number of azimuth angles from a uniform distribution Samples, forming sample pairs And it is stored as a scattering angle sample in the scattering angle sample pool.

6. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 5, characterized in that, After the pre-generated scattering angle sample is generated, the statistical characteristics of the generated scattering angle sample are compared with the theoretical distribution based on the scattering angle distribution characteristic parameters. If the deviation exceeds the allowable range, the sample is regenerated.

7. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 6, characterized in that, The specific steps for the incident ion transport and energy deposition cycle are as follows: S5.1) Based on the current energy of the incident ion and using a lookup table, find the corresponding discrete energy pair and energy loss, and randomly select a preset scattering angle sample from the scattering angle sample pool. As the scattering angle sample for the current step size; S5.2) Using the basic radiation parameter table, calculate the density-normalized projected range difference between the current energy and the next energy, and use it as the step length for the current step; S5.3) The three-dimensional digital differential analysis algorithm is used to distribute the energy loss of the incident ions deposited within the step size into the three-dimensional voxel grid according to the length ratio of the path through each voxel; S5.4) Update the position and orientation of the incident ions based on the scattering angle sample and step size, and subtract energy loss to update the energy of the incident ions; S5.5) Repeat S5.1) to S5.4) until the termination condition is met.

8. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 7, characterized in that, The method of allocating energy loss to the three-dimensional mesh according to the path ratio is implemented by a three-dimensional digital difference analysis algorithm. Specifically, the method includes: linearly allocating energy loss to the corresponding voxel according to the ratio of the length of the incident ion step segment in each voxel to the length of the bus segment.

9. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 8, characterized in that, The termination conditions include the energy of the incident ions being lower than the cutoff threshold or the incident ions penetrating the target area.

10. The rapid simulation method for ion scattering track clusters driven by the bypass factor according to claim 9, characterized in that, The parallel execution specifically involves dividing the total number of ions to be simulated into multiple batches, using a multi-process parallel computing framework to assign the simulation task of each batch to different computing processes; each computing process has an independent copy of the three-dimensional mesh, and after the simulation is completed, all three-dimensional mesh copies are summed.

Images