Proton or heavy ion split dose calculation method, device and radiotherapy system

Abstract

The application discloses a proton or heavy ion separation type dose calculation method and device and a radiotherapy system, relates to the technical field of proton heavy ion radiotherapy, and comprises the following steps: constructing a three-dimensional voxel model of a patient and determining material physical parameters of each voxel based on CT image data of the patient; determining electromagnetic interaction parameters and strong interaction cross-section parameters based on the CT image data of the patient and particle physical parameters; decomposing electromagnetic interaction and strong interaction processes of a charged particle into independent parallel calculation processes, and respectively calculating dose contributions of the electromagnetic interaction process and dose contributions of the strong interaction process; and coupling the obtained doses by probability weighting based on probabilities of the electromagnetic interaction and the strong interaction, and obtaining a final three-dimensional dose distribution. Through decoupling of two types of physical processes at a calculation level and forced sampling by using conditional probabilities, structured parallel step-by-step calculation of a particle transport process is realized, and the calculation efficiency is effectively improved while the dose calculation precision is ensured.

Description

Technical Field

[0001] This invention relates to the field of proton and heavy ion radiotherapy technology, and in particular to a method, apparatus, and radiotherapy equipment for calculating proton or heavy ion dose separately. Background Technology

[0002] Proton and heavy ion radiotherapy, due to its ideal physical dose distribution characteristics and high biological effects, has become an important means of tumor radiotherapy. In the process of proton and heavy ion therapy planning, the accuracy of dose calculation directly affects the quality of the treatment plan; therefore, it is necessary to establish a dose calculation model that can accurately describe the transport process of charged particles in human tissues. Currently, commonly used dose calculation methods in clinical research mainly include analytical algorithms, pencil beam algorithms, and Monte Carlo simulation algorithms. Among them, the Monte Carlo method, by randomly sampling the particle-matter interaction process, can accurately describe the electromagnetic interaction, nuclear reaction, and secondary particle transport processes of particles, and is considered the "gold standard" for particle dose calculation.

[0003] However, traditional Monte Carlo dosimetry methods typically handle both electromagnetic and strong interactions of particles within the same computational framework when simulating particle transport. Specifically, during particle transport, it's necessary to determine whether a nuclear reaction occurs at each transport step, and based on random sampling, decide whether the particle should continue with electromagnetic interactions or undergo strong interactions to generate secondary particles. While this unified simulation accurately reflects the physical process, the continuous nature of electromagnetic interactions and the low-probability discrete nature of strong interactions lead to numerous branch decisions and different computational paths, resulting in reduced computational efficiency. Furthermore, when implementing Monte Carlo dosimetry on parallel computing platforms (such as GPUs), the varying probabilities of nuclear reactions in different transport steps for different particles can cause inconsistent thread execution paths, leading to thread divergence and further reducing parallel computing efficiency.

[0004] On the other hand, clinical treatment planning systems typically require large-scale particle transport simulations to be completed within a short period of time to obtain the patient's three-dimensional dose distribution. Therefore, how to improve the efficiency of particle transport calculations, especially the computational efficiency of parallel computing platforms, while ensuring the accuracy of dose calculations, has become a key technical problem that urgently needs to be solved in the field of particle therapy dose calculation. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a proton or heavy ion separation-based dose calculation method, apparatus, and radiotherapy equipment. It distinguishes between the electromagnetic and strong interactions of charged particles, decoupling these two physical processes at the computational level and employing a conditional probability forced sampling mechanism. This achieves structured parallel step-by-step calculation of the particle transport process, avoiding random branching problems caused by the coupling of different physical processes, reducing thread divergence in parallel computing, and effectively improving computational efficiency while ensuring dose calculation accuracy.

[0006] In a first aspect, the present invention provides a proton or heavy ion separation-based dose calculation method, the method comprising:

[0007] Based on the patient's CT image data, a three-dimensional voxel model of the patient was constructed, and the material physical parameters of each voxel were determined;

[0008] Electromagnetic interaction parameters and strong interaction section parameters were determined based on patient CT image data and particle physics parameters.

[0009] The electromagnetic and strong interactions of charged particles are decomposed into independent and parallel computational processes to calculate the dose contribution of electromagnetic interaction and the dose contribution of strong interaction separately.

[0010] Based on the dose contribution of the electromagnetic interaction and the dose contribution of the strong interaction, as well as the probability of the interaction occurring, the obtained doses are probability-weighted to obtain the final three-dimensional dose distribution.

[0011] Preferably, the process of decomposing the electromagnetic and strong interactions of charged particles into independent and parallel computational processes includes:

[0012] The transport process of charged particles is divided into independent and parallel electromagnetic interaction calculation processes and strong interaction calculation processes;

[0013] In the electromagnetic interaction calculation process, only the electromagnetic interaction process of particles is simulated to calculate the dose deposition size of electromagnetic interaction in each voxel;

[0014] In the strong interaction calculation, only the strong interaction process of particles is simulated, and the dose deposition size of each voxel under strong interaction is calculated.

[0015] Based on the probability of electromagnetic interaction and strong interaction between particles in each voxel, the electromagnetic dose and strong interaction dose deposition are weighted by probability to obtain the final voxel dose distribution.

[0016] Preferably, the dose deposition size of each voxel under electromagnetic interaction is calculated by the following method:

[0017] The probability that a particle will only undergo electromagnetic interaction in a voxel is calculated and multiplied by the dose contribution calculated in the electromagnetic interaction phase to obtain the electromagnetic dose deposition.

[0018] Wherein, the probability P that only electromagnetic interaction occurs EM Calculated using an exponential decay model:

[0019]

[0020] Σhad represents the macroscopic cross section of the strong interaction, and the path length of particle x in the voxel.

[0021] Preferably, the strong interaction dose calculation step includes:

[0022] The strong interaction process is divided into elastic scattering process and inelastic nuclear reaction process;

[0023] Based on the relative proportions of the elastic scattering cross section and the inelastic reaction cross section, the type of strong interaction is sampled using conditional probability to determine whether the particle undergoes elastic scattering or inelastic nuclear reaction, and the dose contribution of elastic scattering and inelastic nuclear reaction is calculated respectively.

[0024] Preferably, in the strong interaction calculation process, a forced collision method is adopted so that each particle experiences a strong interaction event once during the calculation phase; and a statistical weight correction factor is introduced to compensate for the probability deviation caused by forced collision sampling. The weight correction factor is determined based on the actual strong interaction probability of the particle in the voxel.

[0025] Preferably, the probability of strong particle interactions occurring in voxels is calculated using an exponential decay model based on macroscopic sampling cross-sections:

[0026]

[0027]

[0028] Where: P had Σhad represents the probability of strong interactions occurring between particles in a voxel; Σhad represents the macroscopic cross section of the strong interaction; Σ ela For elastic collision sections; Σ nonela is the inelastic nuclear reaction cross section; x is the path length of the particle in the voxel;

[0029] The final voxel dose is calculated using the following formula:

[0030]

[0031] Where D is the final voxel dose; D EM The dose calculated for the electromagnetic interaction phase; Dela D is the dose of elastic scattering in strong interactions; nonels The dose of inelastic nuclear reactions in strong interactions; P EM P represents the probability of electromagnetic interaction occurring. had P represents the probability of a strong interaction occurring. ela P represents the conditional probability of elastic scattering. nonela The conditional probability of an inelastic nuclear reaction; under the condition that strong interactions occur, P satisfies... ela +P nonela =1.

[0032] Preferably, the electromagnetic interaction dose is the electromagnetic dose of the particle under conditions where no strong interaction occurs; the strong interaction dose includes the electromagnetic dose generated by the electromagnetic interaction during the transport process of the particle before the strong interaction occurs, and the dose generated by the secondary particles generated by the strong interaction during subsequent transport.

[0033] The secondary particles produced by the strong interaction are determined in the following way:

[0034] Based on particle physics properties and voxel parameters, information about the reaction channels in which nuclear reactions occur is obtained;

[0035] The types of secondary particles are obtained by sampling the reaction channel using a probability distribution.

[0036] The reaction channel was obtained by sampling the probability distribution of the reaction cross section;

[0037] The energy of the secondary particles is obtained by sampling the differential cross section;

[0038] The emission angle of secondary particles is obtained by sampling double differential cross sections or angular distribution data.

[0039] Preferably, the electromagnetic interaction calculation process and the strong interaction calculation process are executed independently on a parallel computing platform, so that particles in the same calculation stage perform the same type of physical process; in order to reduce thread branch divergence and achieve parallel acceleration of the particle transport process.

[0040] Optionally, the strong interaction calculation process employs importance sampling, weight correction, and / or path splitting methods to achieve statistical enhancement of strong interaction events.

[0041] Preferably, the electromagnetic interaction calculation process is implemented using deterministic or semi-analytical methods.

[0042] Optionally, the process of decomposing the electromagnetic and strong interactions of charged particles into independent and parallel computational processes includes:

[0043] Electromagnetic interaction dose and strong interaction dose are calculated using deterministic transport methods and Monte Carlo methods, respectively, and the total dose is synthesized through a coupling model; the deterministic transport methods include continuous slowing approximation or analytical models.

[0044] In a second aspect, the present invention provides a proton or heavy ion separation dose calculation device, comprising a memory and a processor, wherein the memory stores a computer program, and the processor executes the following steps when executing the computer program:

[0045] Based on the patient's CT image data, a three-dimensional voxel model of the patient was constructed, and the material physical parameters of each voxel were determined;

[0046] Electromagnetic interaction parameters and strong interaction section parameters were determined based on patient CT image data and particle physics parameters.

[0047] The electromagnetic and strong interactions of charged particles are decomposed into independent and parallel computational processes to calculate the dose contribution of electromagnetic interaction and the dose contribution of strong interaction separately.

[0048] Based on the dose contribution of the electromagnetic interaction and the dose contribution of the strong interaction, as well as the probability of the interaction occurring, the obtained doses are probability-weighted to obtain the final three-dimensional dose distribution.

[0049] Thirdly, the present invention provides a radiotherapy system, comprising:

[0050] Particle accelerators are used to generate proton beams or heavy ion beams.

[0051] The proton or heavy ion separation dose calculation device described in this embodiment of the invention is used to calculate the three-dimensional dose distribution of a patient.

[0052] Fourthly, the present invention provides a computer-readable storage medium that stores computer instructions, wherein when a computer reads the computer instructions, the computer executes the steps of the method described in any one of the embodiments of the present invention.

[0053] Compared with existing technologies, the beneficial effects of this invention include at least the following: By distinguishing between electromagnetic interactions and strong interactions, this invention avoids the overhead of branch judgment at each transport step in the traditional Monte Carlo method. On parallel computing platforms, by centralizing the processing of similar physical processes (such as staged computation, particle swarm grouping, or path splitting), thread branch divergence is significantly reduced, and the utilization rate of parallel computing resources such as GPUs is improved.

[0054] The particle transport process in the medium is divided into electromagnetic interaction calculation steps and strong interaction calculation steps. The electromagnetic transport step independently calculates the continuous energy loss processes such as ionization, excitation, and multiple scattering; the strong interaction step independently simulates the strong interaction process of particles and calculates the energy deposited in situ during the strong interaction process and the energy deposition contribution of secondary particles generated by the reaction. By structurally decoupling different physical mechanisms, the complexity of a single transport process is reduced, and a foundation is provided for the independent optimization and parallel implementation of different interaction processes. By introducing enhanced sampling methods such as forced collision mechanisms or importance sampling, the sampling efficiency of low-probability nuclear reaction events is improved, avoiding the waste of a large amount of computation for non-interaction transport, and significantly improving computational efficiency while maintaining statistical accuracy. Modeling or treating electromagnetic interactions, elastic scattering, and inelastic nuclear reactions independently or separately isolates the optimization and extension of different physical processes, facilitating the introduction of more refined nuclear reaction models or electromagnetic transport models, and improving the maintainability and scalability of the algorithm. By obtaining the final dose distribution through probability weighting or event reconstruction, it is possible to quantitatively separate the contributions of electromagnetic interactions, elastic scattering, and inelastic nuclear reactions to the total dose, providing richer information for physical mechanism analysis and treatment plan optimization. Attached Figure Description

[0055] Figure 1 This is a schematic diagram of a proton or heavy ion separation-based dose calculation method according to an embodiment of the present invention;

[0056] Figure 2 This is a roadmap of electromagnetic-strong interaction-separated dose calculation technology according to an embodiment of the present invention. Detailed Implementation

[0057] Exemplary embodiments will now be described more fully with reference to the accompanying drawings. However, these exemplary embodiments can be implemented in many forms and should not be construed as limited to the embodiments set forth herein; rather, they are provided to make the invention more comprehensive and complete, and to fully convey the concept of the exemplary embodiments to those skilled in the art. The same reference numerals in the drawings denote the same or similar structures, and therefore repeated descriptions of them will be omitted.

[0058] The terms used to express position and direction in this invention are illustrated with the accompanying drawings, but changes can be made as needed, and all such changes are included within the scope of protection of this invention.

[0059] In traditional Monte Carlo particle transport algorithms, charged particles need to simultaneously determine whether electromagnetic and strong interactions have occurred in each transport step. Electromagnetic interactions typically exhibit continuous energy loss processes, while nuclear reactions are discrete processes with lower probabilities. These two types of physical processes are intertwined within the same computational framework, leading to a complex algorithm structure, frequent branching decisions, and thus increasing the difficulty of program implementation and reducing computational efficiency.

[0060] On parallel computing platforms such as GPUs, a large number of particle events are simulated simultaneously, but different particles have different probabilities of undergoing nuclear reactions during transport. Some particles only undergo electromagnetic transport, while others may undergo nuclear reactions at a certain location and generate secondary particles. Because different particles execute different computational paths, the execution paths of GPU threads become inconsistent, resulting in thread branch divergence and significantly reducing the efficiency of parallel computing.

[0061] When protons or heavy ions undergo nuclear reactions with a medium, various types of secondary particles are produced, such as protons, neutrons, alpha particles, and heavy ion fragments. The energy, angle, and type of these particles exhibit strong statistical fluctuations. In traditional algorithms, the generation and transport processes of these secondary particles need to be processed in real-time within the same transport framework, increasing computational complexity and potentially leading to slower statistical convergence in dose calculations. Furthermore, due to the low probability of nuclear reactions, extensive simulations are required to obtain reliable results with low uncertainty.

[0062] In proton or heavy ion therapy in the clinical energy range, most particles only undergo electromagnetic interactions during transport, while nuclear reactions occur at a relatively low probability. Traditional Monte Carlo algorithms still require continuous assessment of nuclear reaction conditions for each particle, meaning that a significant amount of computational resources are primarily used for electromagnetic transport simulation, thus reducing the overall dose calculation efficiency.

[0063] In view of the above problems, embodiments of the present invention construct a differentiated particle transport and dose calculation framework. The core idea of ​​this framework is to distinguish between the electromagnetic interaction and strong interaction of charged particles, and to decouple or enhance the two types of physical processes at the computational level through different methods, thereby reducing the complexity of a single transport process, while providing a foundation for the independent optimization and parallel implementation of different interaction processes.

[0064] Example 1, refer to Appendix Figure 1 and attached Figure 2 This embodiment provides a proton or heavy ion separation-based dose calculation method, the method comprising:

[0065] Based on the patient's CT image data, a three-dimensional voxel model of the patient is constructed, and the material physical parameters of each voxel are determined. Specifically, the patient's CT image data is discretized into a three-dimensional voxel mesh according to the spatial resolution set by the treatment planning system, with each voxel having spatial coordinates. The material physical parameters of each voxel are determined, including the material composition, density, and average excitation energy. The material composition can be obtained from a CT value-tissue equivalent material conversion table, the density can be calculated from a CT value-density calibration curve, and the average excitation energy can be obtained by looking up the material composition in a table or by calculating using an empirical formula.

[0066] Electromagnetic interaction parameters and strong interaction cross-section parameters are determined based on patient CT image data and particle physics parameters. The particle beam irradiation parameters include the energy of the primary particles produced by the accelerator and the beam spot size. Electromagnetic interaction parameters include particle ionization energy loss rate, multiple scattering angle distribution parameters, etc., which can be obtained through continuous slowing approximation or calculated based on the Bethe-Bloch formula. Strong interaction cross-section parameters include the macroscopic cross-section of elastic scattering and the macroscopic cross-section of inelastic nuclear reactions, which can be calculated based on particle energy, target nucleus type, and atomic mass number using nuclear reaction theoretical models or nuclear databases (such as ENDF, TALYS, etc.).

[0067] The electromagnetic interactions and strong interactions of charged particles are treated differently to calculate the dose contribution of electromagnetic interactions and strong interactions separately; the treatment of electromagnetic interactions and strong interactions of charged particles is divided into independent and parallel calculation processes.

[0068] As an exemplary implementation, the transport process of charged particles is divided into independent and parallel electromagnetic interaction calculation processes and strong interaction calculation processes. In the electromagnetic interaction calculation process, only the electromagnetic interaction processes such as particle ionization, excitation, and multiple scattering are simulated. The energy deposition of particles in voxels is calculated using Monte Carlo methods or analytical methods to obtain the dose contribution generated by electromagnetic interactions. Simultaneously, the probability that particles only undergo electromagnetic interactions is calculated, and this probability is combined with the dose contribution generated by electromagnetic interactions to obtain the electromagnetic dose deposition.

[0069] In the calculation of strong interactions, the probability of a particle engaging in a strong interaction is determined based on the particle's strong interaction cross-section information. Then, based on this probability, random sampling is used to determine whether a strong interaction event occurs. When a strong interaction event is determined to occur, the energy deposition of the strong interaction process is simulated, including both elastic scattering and inelastic nuclear reactions, and the dose contribution generated by the strong interaction is calculated.

[0070] It should be noted that the above-described phased calculation is only one implementation method for distinguishing between electromagnetic and strong interactions in this invention. Differential processing can also be achieved in other ways, such as: improving the sampling efficiency of strong interaction events by introducing importance sampling or bias sampling for nuclear reaction-related events during the unified transport process; recording the energy deposition of each physical process during the unified simulation process and decomposing and reconstructing the dose according to the probability of occurrence of each process in the post-processing stage; splitting particle trajectories, with one part of the particles specifically used to simulate nuclear reaction processes and the other part used for electromagnetic transport; using deterministic transport methods to calculate the dose distribution for electromagnetic interactions and Monte Carlo methods to simulate strong interactions; dividing particles into two particle groups to handle electromagnetic transport and strong interaction transport separately; or processing different interactions in segments according to the particle energy range. All of the above methods fall within the scope of protection of this invention.

[0071] Based on the electromagnetic interaction contribution and the strong interaction contribution, the three-dimensional dose distribution of the patient is obtained. The electromagnetic interaction contribution and the strong interaction contribution are combined in a certain way to form the final dose value of each voxel. As an example of the combination method, the two types of contributions can be probabilistically weighted based on the probability of strong interaction between particles; after obtaining the dose values ​​of all voxels, the three-dimensional dose distribution of the patient is obtained. It should be noted that the specific combination method of the electromagnetic interaction contribution and the strong interaction contribution can be adjusted according to actual calculation needs, as long as the final dose distribution can be obtained based on the two types of contributions.

[0072] The effects of the above technical solution are as follows:

[0073] By distinguishing between electromagnetic and strong interactions, the complex branching structure of traditional Monte Carlo methods, which simultaneously assess both interactions at each transport step, is avoided. On parallel computing platforms, by centralizing the processing of similar physical processes, thread branch divergence is significantly reduced, and the utilization rate of parallel computing resources such as GPUs is improved, thereby drastically shortening the dose calculation time.

[0074] By treating electromagnetic and strong interactions separately, the modeling and optimization of the two types of physical processes become independent. Researchers can introduce more accurate continuous slowing models for electromagnetic interactions or more refined nuclear reaction models for strong interactions without modifying the other part of the algorithm structure, thus improving the maintainability and scalability of the algorithm.

[0075] The structured design of this embodiment is naturally adapted to large-scale parallel computing platforms such as GPUs. By dividing the particle transport process into stages with similar computational characteristics, a large number of particles perform the same type of computational task in the same stage, giving full play to the advantages of parallel computing architecture and achieving efficient large-scale particle transport simulation.

[0076] See attached document Figure 2 In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0077] The transport process of charged particles is divided into independent and parallel electromagnetic interaction calculation processes and strong interaction calculation processes;

[0078] In the electromagnetic interaction calculation, only the electromagnetic interaction processes such as particle ionization, excitation, and multiple scattering are simulated to calculate the electromagnetic interaction dose. Combined with the probability that particles only undergo electromagnetic interaction, the dose deposition magnitude of electromagnetic interaction in each voxel is obtained. Specifically, the Monte Carlo method is used to simulate the electromagnetic transport of particles in voxels, recording the energy deposited by particles in voxels due to electromagnetic interaction, and obtaining the electromagnetic interaction dose contribution D. EM .

[0079] The dose deposition size of each voxel under electromagnetic interaction is calculated in the following ways:

[0080] Among them, the Monte Carlo method was used to simulate the electromagnetic transport of particles in voxels, and the energy of particles deposited in voxels due to electromagnetic interactions was recorded to obtain the electromagnetic interaction dose contribution D. EM .

[0081] Since particles may interact strongly during actual transport, it is necessary to consider the probability P that the particles only interact electromagnetically. EM For D EM Corrections are made to obtain the true electromagnetic dose deposition. .

[0082] The probability P that a particle undergoes only electromagnetic interaction within a voxel EM Calculated using an exponential decay model:

[0083]

[0084] Where Σhad is the macroscopic cross section of the strong interaction (unit: ), representing the probability that a particle will undergo any type of strong interaction within a unit path length in the medium (i.e., the total probability of elastic scattering or inelastic nuclear reaction); x is the path length of the particle in the voxel (unit: cm); the model is based on the Poisson process assumption, that is, the probability that a particle will not undergo any strong interaction within a path length x follows an exponential distribution.

[0085] The dose D calculated during the electromagnetic interaction phase EM With probability P EM Multiplication yields electromagnetic dose deposition. :

[0086] .

[0087] In the strong interaction calculation, only the strong interaction process of particles is simulated, and the dose deposition size of each voxel under strong interaction is calculated.

[0088] The probability Phad of strong interactions between particles in voxels is also calculated using an exponential decay model:

[0089] Satisfying P EM +P had =1.

[0090]

[0091] Σela: Elastic collision cross section; the probability contribution of the particle to the elastic scattering (kinetic energy conservation, direction change) between the particle and the atomic nucleus.

[0092] Σnonela: Inelastic nuclear reaction cross section, corresponding to the probability contribution of an inelastic collision between a particle and an atomic nucleus (which may produce secondary particles and violate kinetic energy conservation).

[0093] The formula states that the total macroscopic cross section of strong interaction is equal to the sum of the macroscopic cross section of elastic scattering and the macroscopic cross section of inelastic nuclear reaction. This reflects that the two types of strong interaction are probabilistically mutually exclusive and complete. Under the condition of strong interaction, a particle will either undergo elastic scattering or inelastic nuclear reaction, and the total cross section is the sum of the two.

[0094] When it is determined that a strong interaction event has occurred, it is necessary to further determine the type of strong interaction.

[0095] The strong interaction process is divided into elastic scattering process and inelastic nuclear reaction process;

[0096] Based on the relative proportions of the elastic scattering cross section and the inelastic reaction cross section, the type of strong interaction is sampled using conditional probability to determine whether the particle undergoes elastic scattering or inelastic nuclear reaction, and the dose contribution of elastic scattering and inelastic nuclear reaction is calculated respectively.

[0097] Under conditions where strong interactions occur, the conditional probability P of elastic scattering ela It equals the elastic scattering cross section divided by the total strong interaction cross section, i.e., P ela =Σela / Σhad, the conditional probability P of an inelastic nuclear reaction nonela It equals the inelastic nuclear reaction cross section divided by the total strong interaction cross section, i.e., P. nonela =Σnonela / Σhad, and satisfy P ela +Pnonela =1; the particle undergoes elastic scattering or inelastic nuclear reaction based on the conditional probability through random sampling.

[0098] Based on the sampling results, the elastic scattering dose contribution D was calculated respectively. ela Or inelastic nuclear reaction dose contribution D nonela For inelastic nuclear reactions, further simulations of the generation and transport processes of secondary particles are needed, including the energy and angular distribution of secondary particles and subsequent energy deposition.

[0099] To further improve the sampling efficiency of strong interaction events, this embodiment employs a forced collision method during the strong interaction calculation process, ensuring that each particle experiences a strong interaction event during this calculation phase. Specifically, by normalizing the strong interaction probability, the particle is forced to undergo either elastic scattering or an inelastic nuclear reaction.

[0100] To avoid statistical bias caused by forced occurrence events, a statistical weighting correction factor is introduced to weight the dose contribution. The weighting correction factor is determined based on the actual probability of strong interactions between particles in the voxel; in this embodiment, the weighting correction factor is equal to P. had That is, under forced collision sampling, the dose contribution of elastic scattering or inelastic nuclear reaction needs to be multiplied by P. had In order to restore the probabilistic and statistical properties of real physical processes.

[0101] Based on the probability of electromagnetic interaction and strong interaction between particles in each voxel, the electromagnetic dose and strong interaction dose deposition are weighted by probability to obtain the final voxel dose distribution.

[0102] The final voxel dose is calculated using the following formula:

[0103]

[0104] Where D is the final voxel dose; D EM The dose calculated for the electromagnetic interaction phase; D ela D is the dose of elastic scattering in strong interactions; nonels The dose of inelastic nuclear reactions in strong interactions; P EM P represents the probability of electromagnetic interaction occurring. had P represents the probability of a strong interaction occurring. ela P represents the conditional probability of elastic scattering. nonela The conditional probability of an inelastic nuclear reaction;

[0105] The physical meaning of this formula is that the energy deposition of a particle during transport consists of two parts: one part is the energy deposited by electromagnetic interaction alone when the particle does not have any strong interaction, and the other part is the energy deposited when the particle has a strong interaction (including the electromagnetic energy deposited before the strong interaction and the energy deposited during the strong interaction itself). The two parts are weighted and summed according to their respective probabilities, ensuring that the calculation result is statistically unbiased.

[0106] After traversing all voxels, the three-dimensional dose distribution of the patient was obtained.

[0107] The effects of the above technical solution are as follows:

[0108] In traditional Monte Carlo methods, electromagnetic and strong interactions are intertwined and judged within the same transport step, leading to a complex algorithm structure. Strong interactions are low-probability events, with many particles undergoing no nuclear reactions during transport, resulting in wasted computational resources. This embodiment separates the two physical processes, allowing the electromagnetic interaction stage to handle only continuous energy loss processes and the strong interaction stage to handle only discrete event processes, thus achieving structural decoupling of the physical processes. This decoupling ensures that the computational tasks of each stage have similar physical characteristics, facilitating algorithm optimization and parallel implementation. It avoids the complex branching structure of traditional Monte Carlo methods, which simultaneously judge both interactions at each transport step. On parallel computing platforms, the electromagnetic and strong interaction stages are executed independently, ensuring that particles within the same stage perform the same type of physical process, significantly reducing thread branch divergence problems in parallel computing architectures such as GPUs.

[0109] By introducing a forced collision sampling method, each particle is guaranteed to experience a strong interaction event during the strong interaction calculation, thus avoiding the waste of computation on a large amount of non-interaction transport. Combined with a statistical weighting correction factor, the sampling efficiency for low-probability nuclear reaction events is significantly improved while ensuring unbiased calculation results. For scenarios with a low probability of nuclear reactions in clinical proton therapy (such as high-energy proton beams), this embodiment can achieve lower statistical uncertainty within the same computation time, or significantly shorten the computation time while achieving the same statistical accuracy.

[0110] By modeling electromagnetic interactions, elastic scattering, and inelastic nuclear reactions independently, the optimization and extension of different physical processes are isolated from each other. More accurate continuous slowing models (such as the Bethe-Bloch formula considering density effects) can be introduced for electromagnetic interactions, or more refined nuclear reaction models (such as quantum molecular dynamics models) can be introduced for strong interactions, without modifying the algorithm structure of other parts. This modular design improves the maintainability and scalability of the algorithm.

[0111] Obtaining the final dose distribution through a probabilistic weighting method allows for the quantitative separation of the contributions of electromagnetic interactions, elastic scattering, and inelastic nuclear reactions to the total dose. This decompositional modeling not only contributes to a deeper understanding of the dose formation mechanism of nuclear interactions in particle therapy but also provides richer information for optimizing treatment plans. For example, clinical physicists can separately assess the impact of secondary particles generated by nuclear reactions on the target area dose and the dose to organs at risk, thereby developing more precise treatment protocols.

[0112] In one possible implementation, the electromagnetic interaction dose is the electromagnetic dose of the particle under conditions where no strong interaction occurs; the strong interaction dose includes the electromagnetic dose generated by the electromagnetic interaction during the transport process of the particle before the strong interaction occurs, and the dose generated by the secondary particles generated by the strong interaction during subsequent transport.

[0113] The secondary particles produced by the strong interaction are determined in the following way:

[0114] Based on particle physics properties and voxel parameters, information about the reaction channels in which nuclear reactions occur is obtained;

[0115] The types of secondary particles are obtained by sampling the reaction channel using a probability distribution.

[0116] The reaction channel was obtained by sampling the probability distribution of the reaction cross section;

[0117] The energy of the secondary particles is obtained by sampling the differential cross section;

[0118] The emission angle of the secondary particles is obtained by sampling the double differential cross section or angular distribution data;

[0119] The data sources include TALYS, ENDF database, experimental data, or theoretical model calculation results.

[0120] Electromagnetic interaction dose is defined as the energy deposited in a medium by a particle through electromagnetic interaction processes such as ionization, excitation, and multiple scattering, without any strong interactions occurring. This dose contribution reflects the interaction between the particle and the electrons outside the atomic nuclei in the medium and is the main pathway for particle energy loss.

[0121] In actual transport processes, particles may never engage in strong interactions along the entire path length, or they may engage in strong interactions at a specific location. This embodiment distinguishes between electromagnetic interaction dose and strong interaction dose through probability separation. Specifically, in the electromagnetic interaction calculation process, particles are assumed to only engage in electromagnetic interactions, and the calculated dose contribution is D. EM This dose contribution multiplied by the probability P that the particle only undergoes electromagnetic interaction. EMAfterwards, the actual electromagnetic dose deposition is obtained; this dose corresponds to the scenario where the particles only undergo electromagnetic interaction throughout the transport process.

[0122] Strong interaction dose consists of two parts:

[0123] (1) Electromagnetic energy deposition before strong interaction: Before a strong interaction event occurs, particles have already deposited a portion of energy through electromagnetic interaction. This portion of energy is closely related to the location of the strong interaction event and needs to be accumulated based on the transport path of the particles in voxels.

[0124] (2) Transport and energy deposition of secondary particles generated by the strong interaction: When particles undergo inelastic nuclear reactions, various types of secondary particles are produced, including protons, neutrons, alpha particles, and heavy ion fragments. These secondary particles continue to be transported in the medium and deposit energy through electromagnetic interactions or further nuclear reactions. The contribution of secondary particles is an important component of the strong interaction dose.

[0125] In practical implementation, the strong interaction dose D had This is achieved by first accumulating the energy deposited by the particle through electromagnetic interactions from its initial position to the event location (i.e., the electromagnetic energy deposited before the strong interaction). Then, the energy deposit of the event itself is calculated based on the type of strong interaction (elastic scattering or inelastic nuclear reaction). For inelastic nuclear reactions, the generation and transport of secondary particles must be further simulated, and all energy deposited by the secondary particles during subsequent transport must be accumulated.

[0126] When particles undergo inelastic nuclear reactions, the types, energies, and emission angles of the reaction products exhibit significant statistical fluctuations. This embodiment employs a random sampling method based on nuclear reaction data and theoretical models to accurately describe the physical properties of secondary particles.

[0127] The reaction channels (i.e., the types and combinations of reaction products) are obtained through sampling of the reaction cross-section probability distribution. For a given incident particle energy and target nucleus type, nuclear reactions may occur through multiple reaction channels, each with a corresponding reaction cross-section. In this embodiment, the cross-section value of each reaction channel is first calculated to construct a normalized probability distribution, and then the specific reaction channel that occurs is determined through random number sampling.

[0128] The sources of reaction cross-section data include:

[0129] Nuclear databases: such as ENDF (Evaluated Nuclear Data File), JEFF, JENDL and other evaluation nuclear databases, provide validated nuclear reaction cross-sectional data;

[0130] Nuclear reaction theory models: such as the TALYS program, which can calculate cross-sectional data for any combination of incident energy and target nuclei based on nuclear reaction theory;

[0131] Experimental data: Publicly available experimental measurements of nuclear reactions can be used to verify or replace theoretical calculations.

[0132] After determining the reaction channel, the energy of the secondary particles is obtained by sampling using the differential cross section dσ / dE. The differential cross section describes the distribution probability of the products within a specific energy range. Specifically, for a selected reaction channel, its energy differential cross section data is stored in the form of discrete points or analytical functions. In this embodiment, sampling is performed using the inverse function method or the discarding method to obtain the energy values ​​of the secondary particles.

[0133] For many-body reactions (such as reactions that produce multiple secondary particles), energy sampling must comply with the law of conservation of energy. This embodiment employs a phase space sampling method to jointly sample the energy of each secondary particle while satisfying the laws of conservation of energy and momentum.

[0134] The emission angle of the secondary particles is obtained through a double differential section. Alternatively, it can be obtained by sampling angular distribution data. The double differential section describes the probability density function of the joint distribution of particle energy and emission angle. In specific implementations, the following methods can be used:

[0135] Sampling based on angular distribution data: For a specific reaction channel, the angular distribution data is stored in the form of Legendre polynomial coefficients or discrete angle points. In this embodiment, the emission angle is obtained by sampling based on the angular distribution function using either a discarding method or an inverse function method.

[0136] Joint sampling based on dual differential cross sections: When there is a correlation between energy and angle, joint sampling is performed using dual differential cross sections, sampling energy first and then conditionally sampling angle, or vice versa.

[0137] Isotropic assumption: For some reaction channels (such as low-energy neutron emission), an isotropic angular distribution approximation can be used to simplify the sampling process.

[0138] This embodiment supports multiple data sources, which can be flexibly selected according to the requirements of calculation accuracy and efficiency:

[0139] TALYS: A nuclear reaction theory program capable of calculating nuclear reaction data for various incident particles, including neutrons, protons, deuterons, and alpha particles, covering a wide energy range (keV to hundreds of MeV). The cross-sectional data, energy distribution, and angular distribution output by TALYS can be directly used in the sampling process of this embodiment.

[0140] ENDF database: An internationally recognized evaluation nuclear database containing a large amount of data on the reaction cross sections and secondary particle distributions of nuclides. This embodiment can read ENDF format data, parse it, and construct the probability distribution required for sampling.

[0141] Experimental data: When high-precision experimental measurement data is available, experimental data should be used first to obtain a more accurate description of the reaction.

[0142] Theoretical model calculation results: For energy or target core ranges not covered by databases or experimental data, calculation results from theoretical models (such as IQMD, BUU, etc.) can be used as a supplement.

[0143] After sampling and obtaining the type, energy, and emission angle of the secondary particles, these secondary particles are used as new particle sources to continue transport simulations in the medium. The transport process of secondary particles is similar to that of primary particles, also following the phased computational framework of this embodiment:

[0144] During transport, secondary particles first undergo an electromagnetic interaction calculation process, resulting in the deposition of electromagnetic interaction doses.

[0145] When secondary particles undergo nuclear reactions, they further trigger strong interaction calculations, producing secondary particles of higher generations.

[0146] This recursive process continues until the energy of all secondary particles falls below a set threshold or escapes the computational domain. The energy deposited by particles from all generations is accumulated into the corresponding voxels according to their corresponding probability weights, constituting the complete contribution of the total dose.

[0147] The effects of the above technical solution are as follows:

[0148] The energy deposition process of a particle in a medium can be decomposed into the contributions of several independent physical events. Assume that along the transport path of the particle from its starting point to its destination, zero, one, or multiple strong interaction events may occur. According to the law of total probability, the total expected energy deposition E of the particle can be expressed as:

[0149]

[0150] Where P(n) is the probability of n strong interactions occurring. This represents the expected energy deposition when n strong interactions occur. In this embodiment, n=0 (electromagnetic interaction only) and n≥1 (at least one strong interaction occur) are treated separately, and the electromagnetic dose deposition and strong interaction dose are calculated separately. This decomposition method strictly follows probability and statistics theory, ensuring the unbiasedness of the calculation results.

[0151] Before a strong interaction event occurs, the particle has already deposited a certain amount of electromagnetic energy along its transport path. This energy is closely related to the location of the event. This embodiment accurately reflects the spatiotemporal distribution of energy loss during particle transport by accumulating the electromagnetic energy deposited before the event, avoiding the conflation of electromagnetic energy before and after the strong interaction, and improving the physical accuracy of dose calculation. Compared with traditional methods that mix electromagnetic and strong interaction data, this embodiment avoids dose statistical bias caused by ambiguous event sequence, improving the spatial accuracy and physical accuracy of dose calculation.

[0152] By clearly defining the dosage composition of electromagnetic interaction and strong interaction, and accurately simulating the complete electromagnetic energy deposition and secondary particle transport process before the generation of strong interaction, this embodiment can more accurately reflect the physical mechanism of particle energy deposition in the medium. Compared with the traditional method of mixing electromagnetic and strong interactions, this embodiment avoids the dose statistical bias caused by the ambiguity of the event sequence, and improves the spatial accuracy and physical realism of dose calculation.

[0153] In one possible implementation, the electromagnetic interaction calculation process and the strong interaction calculation process are executed independently on a parallel computing platform, so that particles in the same calculation stage perform the same type of physical process; thereby reducing thread branch divergence and achieving parallel acceleration of the particle transport process.

[0154] The effects of the above technical solution are as follows:

[0155] In traditional Monte Carlo transport methods, each thread simulates the transport process of one particle. Because different particles follow different transport paths, the physical events that occur also differ: some particles only undergo electromagnetic interactions, some undergo elastic scattering, and some undergo inelastic nuclear reactions and produce secondary particles. This difference in computational paths leads to threads within the same thread bundle executing different instruction streams, a phenomenon known as thread divergence. In GPU architectures, when thread divergence occurs within a thread bundle, instructions for different branches need to be executed sequentially, resulting in decreased thread bundle utilization and reduced computational efficiency.

[0156] This embodiment eliminates branching divergence by executing electromagnetic and strong interactions independently in stages, ensuring that all particles perform the same type of physical process within each stage. In the electromagnetic stage, all particles only process electromagnetic interactions; in the strong interaction stage, all particles process strong interaction events (whether elastic scattering or inelastic nuclear reactions, the computational flow is essentially similar). Therefore, the instruction flow of the thread bundle is completely consistent within each stage, achieving 100% thread utilization.

[0157] Optimize memory access patterns through particle grouping and batch processing:

[0158] Coalesced access: When threads within the same thread bundle access contiguous memory addresses, the GPU can merge multiple accesses into a single memory transaction, significantly improving bandwidth utilization. This embodiment employs a SoA data layout, ensuring that the access patterns of particle attributes naturally meet the requirements for coalesced access.

[0159] Shared memory reuse: For particle grouping, the state data of particles within a group is loaded into shared memory. Since the access latency of shared memory is much lower than that of global memory, and it can be reused by all threads within the same thread block, the number of accesses to global memory is significantly reduced.

[0160] GPU computation typically depends on CPU-GPU data transfer. This embodiment employs asynchronous execution and multi-streaming techniques to achieve overlap between computation and data transfer:

[0161] Asynchronous memory copy: Using CUDA streams, memory copy operations and kernel function execution are placed in different streams, enabling data transfer and computation to be performed in parallel.

[0162] Double buffering: Two buffers are allocated, one for computation of the current batch and the other for data transfer of the next batch, to achieve pipelined parallelism.

[0163] Because strong interaction events have a low probability of occurrence, the computational load of different particle batches may vary (some batches may contain more particles undergoing nuclear reactions). This embodiment employs a dynamic load balancing strategy, which dynamically adjusts the scheduling order based on the statistical characteristics of particle batches (such as the nuclear reaction rate, the number of secondary particles, etc.), allocating batches with higher computational loads to GPU cores or thread blocks with stronger computing power, thus avoiding resource waste caused by uneven load distribution.

[0164] In one possible implementation, the strong interaction calculation process employs importance sampling, weight correction, and / or path splitting methods to achieve statistical enhancement of strong interaction events.

[0165] The strong interaction calculation process employs an importance sampling method to statistically enhance strong interaction events, including: introducing biased sampling for nuclear reaction-related events, increasing the sampling probability of high-contribution events, and restoring unbiased statistics through a weight correction factor.

[0166] The strong interaction calculation process employs a weighting correction method to statistically enhance strong interaction events. This includes allowing particles to interact randomly according to their true mean free path, while introducing a statistical weighting correction factor to weight the dose contribution in order to compensate for statistical bias caused by insufficient sampling of rare events.

[0167] The strong interaction calculation process employs a path splitting method to achieve statistical enhancement of strong interaction events, including: during particle transport, when a particle meets preset conditions, the particle trajectory is split to generate multiple sub-particle copies, which respectively simulate the nuclear reaction process and the electromagnetic transport process, and statistical consistency is maintained through weight allocation.

[0168] The statistical augmentation methods can be used in combination, including the combined use of importance sampling and weight adjustment, the combined use of path splitting and weight adjustment, or the combined use of importance sampling and path splitting.

[0169] In one specific implementation, a bias factor is defined based on the particle's energy and the type of medium it inhabits. The bias factor reflects the degree of sample enhancement for nuclear reaction events; a value greater than 1 indicates an increased sampling probability. The biased probability distribution is obtained by multiplying the true probability distribution by the bias factor and then normalizing it.

[0170] The following bias strategies can be adopted:

[0171] Energy-dependent bias: For high-energy particles, the nuclear reaction cross-section is larger, so a smaller bias factor is used, usually with a value of 2 to 5; for low-energy particles, the nuclear reaction probability is extremely low, so a larger bias factor is used, usually with a value of 10 to 20, in order to increase the sampling opportunity of nuclear reactions in the low-energy region.

[0172] Region-dependent bias: For critical regions such as the target area or organs at risk, increase the sampling probability of nuclear reaction events to more accurately assess the contribution of nuclear reactions to the dose.

[0173] In the calculation of strong interactions, a biased probability distribution is used instead of the true probability distribution for sampling to determine whether a strong interaction event will occur. Specifically, the calculation is performed using a biased macroscopic cross-section, i.e., the true macroscopic cross-section is multiplied by a bias factor, and then the biased probability of strong interaction occurrence is calculated based on this biased cross-section. Random sampling is performed using this biased probability, increasing the probability of nuclear reaction events to a factor of the original bias factor.

[0174] To maintain statistical unbiasedness, a weight correction factor is introduced, initially set to 1. When biased sampling is used, the weights are updated by multiplying the original weights by the reciprocal of the bias factor. That is, after each sampling, the particle weights are reduced according to the reciprocal of the bias factor. In the final dose contribution calculation, the particle energy deposition is multiplied by the current weight to achieve statistical compensation.

[0175] To avoid statistical efficiency degradation due to excessive weight differences, when the weight is too small, the Russian roulette method is used to randomly terminate the particle; when the weight is greater than a set threshold, the particle is split to generate multiple copies and assign weights, with the weight of each copy being the original weight divided by the number of split copies.

[0176] In one specific implementation, nuclear reactions are low-probability events during strong interaction calculations. Traditional methods often fail to involve a large number of particles, leading to an underestimation of the statistical contribution of nuclear reaction doses. Weighting correction methods compensate for the statistical bias caused by insufficient sampling of rare events by assigning higher weights to particles that do participate in nuclear reactions.

[0177] Particles are allowed to interact randomly according to their true mean free path. Specifically, the free path length of particles is randomly sampled using an exponential distribution. The free path length is sampled based on uniform random numbers obtained through the inverse function of an exponential distribution. If the sampled free path length is less than the path length within the current voxel, the particle interacts strongly within that voxel; otherwise, the particle continues to be transported.

[0178] When particles interact strongly, their weights are multiplied by a weight correction factor. The weight correction factor is determined based on the event type and its importance.

[0179] Standard weight correction: The weight correction factor is the reciprocal of the true probability that a particle will interact strongly in the current voxel. This correction factor compensates for the statistical contribution of low-probability events.

[0180] Energy-dependent correction: For highly penetrating particles such as secondary neutrons produced by high-energy particles, a larger correction factor is used, that is, multiplying the standard correction factor by a coefficient greater than 1, to enhance the statistical accuracy of their dose contribution.

[0181] Region-dependent correction: For critical regions, a larger correction factor is used to improve the reliability of dose calculation in those regions.

[0182] The energy of particle deposition is multiplied by the current weight and then added to the corresponding voxel to achieve dose contribution statistics after weight correction.

[0183] To prevent excessive weighting from exacerbating statistical fluctuations, a weight cap is set. When a particle's weight exceeds the weight cap, the particle is split into multiple copies. The weight of each copy is the original weight divided by the number of copies split, and the process continues.

[0184] Path splitting is a technique that enhances the sampling efficiency of specific physical processes by replicating particle trajectories. During particle transport, when a particle enters a specific region or reaches a certain energy condition, its trajectory is split, generating multiple sub-particle copies. Some of these copies are specifically used to simulate nuclear reaction processes, while others are used for electromagnetic transport processes, thereby increasing the sampling density of nuclear reaction events.

[0185] During particle transport, it is determined in real time whether the splitting conditions are met. The splitting conditions include:

[0186] Region condition: Splitting is triggered when a particle enters a critical region;

[0187] Energy conditions: Splitting is triggered when the particle energy is below a set threshold. Because the nuclear reaction cross section of low-energy particles is small, enhanced sampling is required.

[0188] Depth condition: Splitting is triggered when the particle depth exceeds the Bragg peak position, and the nuclear reaction in this region contributes significantly to the dose tail;

[0189] Random condition: Splits are triggered randomly with a certain probability, used to enhance versatility.

[0190] When the splitting condition is met, the original particle is split into multiple sub-particle copies. The splitting factor is set according to computational resource requirements, typically ranging from 2 to 10. The task allocation method for the sub-particle copies is as follows:

[0191] Dedicated nuclear reaction replicas: Certain subparticle replicas are specifically designed to simulate nuclear reaction processes. These replicas force nuclear reaction events to be triggered during strong interaction calculations to enhance nuclear reaction sampling.

[0192] Electromagnetic transport replicas: The remaining sub-particle replicas only perform electromagnetic transport simulations and do not trigger nuclear reaction events.

[0193] Hybrid replicas: Some replicas maintain the standard transport pattern, i.e., nuclear reactions occur according to the true probability.

[0194] The weights of the original particles are proportionally distributed among the sub-particle copies. Weight distribution methods include:

[0195] Equal weight allocation: Each replica has the same weight, that is, the original weight is divided by the number of split replicas.

[0196] Importance-weighted allocation: the nuclear reaction copy is assigned a higher weight, and the electromagnetic copy is assigned a lower weight, to reflect the difference in the importance of the different copies to the dose contribution.

[0197] Each subparticle copy is independently simulated for subsequent transport. During dose accumulation, the energy deposited by each copy is multiplied by its weight and then accumulated into the corresponding voxel.

[0198] To prevent the computational load from spiraling out of control due to explosive particle growth, a maximum split depth is set, for example, allowing a maximum of 3 splits. Once this depth is exceeded, no further splits will occur.

[0199] This embodiment supports the combined use of multiple statistical augmentation methods to leverage their respective advantages:

[0200] In the importance sampling process, the particle weights already include a bias correction factor. To further enhance the contribution of rare events, weight correction can be superimposed on the importance sampling. In this case, the total weight correction factor is the product of the importance sampling weight correction factor and the weight correction factor from the weight correction method. Using them together can simultaneously obtain the dual advantages of improved sampling efficiency from importance sampling and enhanced statistics on rare events from weight correction.

[0201] The nuclear reaction replicas generated by path splitting can be further weighted using a weighting correction method to enhance their statistical contribution to the nuclear reaction event. In this case, the initial weight of the nuclear reaction replica is the original particle weight divided by the splitting factor and then multiplied by the weighting correction factor.

[0202] Before path splitting, importance sampling is used to increase the sampling probability of nuclear reaction events, and then path splitting is performed on important regions to achieve a second improvement in sampling efficiency.

[0203] By employing importance sampling, weight correction, and / or path splitting methods, this embodiment significantly improves the sampling efficiency for low-probability nuclear reaction events.

[0204] In one possible implementation, the electromagnetic interaction calculation process is carried out using deterministic or semi-analytical methods.

[0205] Deterministic methods model the average behavior of particle-matter interactions, ignoring statistical fluctuations and directly calculating the energy deposition distribution of particles in the medium through analytical formulas or numerical integration. The core assumption of this method is that the collective behavior of a large number of particles can be approximated by continuous slowing down, meaning that the energy loss of particles in the medium is continuous and deterministic, rather than a discrete stochastic process.

[0206] In proton and heavy ion therapy, electromagnetic interactions (ionization, excitation, and multiple scattering) are the primary mechanisms for particle energy loss, and their average behavior is highly predictable. Deterministic methods leverage this property to rapidly calculate particle range, energy loss distribution, and dose distribution using analytical formulas.

[0207] The specific implementation steps include:

[0208] Based on the initial energy of the particle beam and the material properties of the medium, a relationship between range and energy is established. For proton beams, the range-energy relationship can be obtained through empirical formulas or lookup tables; for heavy ion beams, charge states and relativistic effects must be considered. Range-energy data can be pre-calculated and stored as lookup tables for quick retrieval in subsequent calculations.

[0209] The energy loss of a particle in a medium is calculated using the continuous slowing-down approximation. Specifically, the medium is discretized into multiple thin layers along the particle's trajectory. Within each thin layer, the energy loss of the particle after passing through that layer is calculated based on its current energy and the stopping power of the medium. The stopping power is calculated using the Bethe-Bloch formula, taking into account density effects and shell corrections.

[0210] The energy loss is calculated using a layer-by-layer recursive approach: the initial energy is known, the energy loss of the first layer is calculated to obtain the remaining energy, and this remaining energy is then used as the input for the next layer, and so on, until the particle energy drops to the cutoff threshold or the particle escapes the calculation region. The energy deposited in each layer is the electromagnetic dose contribution of that layer's voxel.

[0211] For scanning or scattering beams, beam spot size and beam intensity distribution must be considered. The beam is spatially discretized into multiple sub-beams, each with a specific energy and weight. A continuous slowing approximation is performed on each sub-beam to obtain a one-dimensional dose distribution along the beam's central axis. Then, considering the transverse distribution of the beam spot (typically Gaussian), the one-dimensional dose distribution is convolved with the transverse distribution to obtain a three-dimensional dose distribution.

[0212] For the lateral diffusion of the particle beam, analytical models are used to approximate the multiple scattering effect. Commonly used models include the Fermi-Eyges theory, which describes the beam broadening in the medium by calculating the angular diffusion and lateral displacement distribution of the particles. Specifically, based on the particle energy and the radiation length of the medium, the mean square displacement of multiple scattering is calculated, and a Gaussian distribution is used to approximate the lateral distribution of the beam.

[0213] Because deterministic methods neglect energy loss fluctuations, the calculated results deviate somewhat from the actual physical processes. To correct this deviation, a dose correction factor is introduced. This correction factor can be obtained by comparing with Monte Carlo simulation results or given as an empirical formula. For different combinations of energy and media, the correction factor can be pre-calculated and stored as a lookup table.

[0214] Semi-analytical methods lie between deterministic and Monte Carlo methods. While retaining some analytical descriptions of the physical processes, they employ random sampling or numerical integration for key processes. This approach, while maintaining computational speed, can more accurately describe the statistical properties of electromagnetic interactions, making it suitable for scenarios requiring high precision.

[0215] The specific implementation steps include:

[0216] Similar to deterministic methods, the medium is discretized into multiple thin layers, and the average energy loss of particles is calculated using a continuous slowing-down approximation. However, in semi-analytical methods, the average energy loss is not directly used as the dose contribution of the layer; instead, it is corrected based on the statistical distribution of the energy loss.

[0217] The statistical distribution of energy loss can be described by a Landau or Vavilov distribution, the width of which depends on particle energy, medium thickness, and material properties. For each thin layer, the mean and variance of the energy loss are calculated, and then the actual energy loss of that layer is obtained by sampling according to the distribution function, or the mean is used directly with an additional statistical fluctuation correction.

[0218] The range of a particle beam in a medium exhibits statistical fluctuations, known as range deviation. Semi-analytical methods describe this range deviation using analytical models, typically employing a Gaussian distribution approximation, with its variance calculated based on particle energy and medium properties. In dose calculation, the dose distribution obtained from the continuous slowing approximation is convolved with the range deviation distribution to obtain a dose distribution that more closely approximates the true physical distribution.

[0219] For multiple scattering, the semi-analytical method employs a hybrid approach: for most particles, the analytical solution of Fermi-Eyges theory is used to describe lateral diffusion; for critical regions (such as beam edges and steep gradient regions), a small amount of random sampling is introduced to more accurately describe the distribution of scattering angles. This hybrid approach improves the dose calculation accuracy in edge regions while maintaining computational efficiency.

[0220] Although the electromagnetic interaction phase does not directly simulate nuclear reactions, their occurrence affects particle beam transport, leading to beam attenuation and energy changes. The semi-analytical method corrects for the effects of nuclear reactions using an analytical model: it calculates the attenuation probability of particles during transport based on the nuclear reaction cross-section and subtracts the energy lost due to nuclear reactions from the energy loss calculation. The contribution of secondary particles generated by nuclear reactions is calculated separately in the strong interaction phase.

[0221] The semi-analytical method can dynamically adjust the computational accuracy based on the importance of the computational region. For the target area and high-dose gradient regions, a higher-precision processing method is used (e.g., increasing the number of layers and introducing random sampling); for low-dose regions and homogeneous areas, a lower-precision analytical approximation is used. This adaptive accuracy control can further improve overall computational efficiency while ensuring computational accuracy in critical regions.

[0222] While deterministic or semi-analytical methods are used to calculate the electromagnetic interaction, the strong interaction phase is still simulated using the Monte Carlo method. The coupling between the two phases is achieved in the following way:

[0223] First, deterministic or semi-analytical methods are used to quickly calculate the dose distribution when particles only interact electromagnetically in the medium, thus obtaining the basic field of electromagnetic dose deposition.

[0224] The probability of strong interactions between particles is calculated based on the nuclear reaction cross-section, and the attenuation of the particle beam is corrected accordingly. Specifically, during electromagnetic transport, the intensity of the particle beam attenuates exponentially, with the attenuation coefficient equal to the macroscopic cross-section of the strong interaction. This correction ensures that electromagnetic dose deposition matches the probability of nuclear reactions.

[0225] In the calculation of strong interactions, the Monte Carlo method is used to simulate nuclear reaction events and the transport of secondary particles. For particles undergoing nuclear reactions, the electromagnetic energy deposited before the nuclear reaction occurs is included in the strong interaction dose, while the energy deposit of the nuclear reaction itself and the contribution of secondary particles are also simulated. The final total dose is the sum of the electromagnetic dose deposit and the strong interaction dose.

[0226] Deterministic and semi-analytical methods avoid the random sampling and extensive particle history tracing required in Monte Carlo simulations, significantly improving computational speed. For typical proton therapy plans for prostate cancer, the time required to calculate the electromagnetic dose distribution using deterministic methods can be reduced from minutes to seconds, meeting the time requirements for online dose validation and real-time plan optimization.

[0227] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0228] Instead of explicitly separating electromagnetic and strong interaction processes, the sampling efficiency of strong interaction events is improved by introducing importance sampling or bias sampling of nuclear reaction-related events in a unified transport process.

[0229] or,

[0230] The energy deposition of each physical process is recorded during the unified simulation process, and the dose is decomposed and reconstructed according to the occurrence probability of each process in the post-processing stage.

[0231] In one possible implementation, instead of explicitly separating the electromagnetic and strong interaction processes, the sampling efficiency of strong interaction events is improved by introducing importance sampling or bias sampling of nuclear reaction-related events within the unified transport process. Specifically, this includes:

[0232] In the process of particle transport, a biased probability distribution of nuclear reaction events is constructed based on the energy of the particles and the medium in which they are located.

[0233] By using a biased probability distribution instead of the true probability distribution for sampling, the sampling probability of nuclear reaction events is increased;

[0234] A weighting correction factor is introduced to weight the dose contribution in order to compensate for the probability bias caused by biased sampling.

[0235] The bias probability distribution is constructed based on the energy dependence or region dependence of the particles:

[0236] Energy-dependent approach: A larger bias factor is used for low-energy particles, and a smaller bias factor is used for high-energy particles;

[0237] Region-dependent approach: Increase the sampling probability of nuclear reaction events for critical regions.

[0238] The weight correction factor is the reciprocal of the bias factor, and the weight of a particle is multiplied by this weight correction factor after each bias sampling.

[0239] The process of recording the energy deposition of each physical process during the unified simulation and decomposing and reconstructing the dose according to the occurrence probability of each process in the post-processing stage specifically includes:

[0240] In the unified transport simulation, the energy deposited by each physical process (including electromagnetic interaction, elastic scattering, and inelastic nuclear reaction) when each particle interacts with the medium is recorded.

[0241] After the simulation is completed, the recorded energy deposition is statistically decomposed according to the actual occurrence probability of each physical process to obtain the independent dose contribution of each physical process.

[0242] The dose contributions of each process after decomposition are reconstructed to form the total dose distribution.

[0243] The energy deposition recorded for each physical process includes:

[0244] Multiple dose counters are established for each voxel, corresponding to the contributions of electromagnetic interaction, elastic scattering, inelastic nuclear reaction, and various types of secondary particles, respectively.

[0245] At each step of particle transport, the deposited energy is added to the corresponding counter, depending on the type of physical process currently occurring.

[0246] The process of decomposing and reconstructing the dose based on the probability of occurrence of each process includes:

[0247] Calculate the actual probability of each physical process occurring in each voxel, which is determined based on particle energy, medium material, and interaction cross section;

[0248] The recorded dose contribution was normalized according to the probability of occurrence to eliminate the influence of statistical fluctuations;

[0249] The dose contributions of each normalized process are summed to obtain the final voxel dose.

[0250] The dose decomposition and reconstruction also includes recursive processing of the contribution of secondary particles: the energy deposits generated by secondary particles are traced back to the initial nuclear reaction event according to their source physical process and assigned to the corresponding dose counter.

[0251] The effects of the above technical solution are as follows:

[0252] By recording the energy deposition of each physical process and decomposing and reconstructing it in the post-processing stage, the independent dose contributions of electromagnetic interactions, elastic scattering, inelastic nuclear reactions, and various types of secondary particles can be quantitatively obtained. The separated dose distribution helps to gain a deeper understanding of the contribution of different physical mechanisms to the total dose, providing a quantitative basis for treatment planning optimization and radiation risk assessment. This can be achieved without modifying the core logic of the unified transport framework, only by adding a dose counter, resulting in minimal code invasiveness and easy integration into existing Monte Carlo code.

[0253] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0254] During particle transport, when a particle enters a specific region or reaches a certain energy condition, its trajectory splits. One part of the particles are used to simulate nuclear reaction processes, while the other part is used for electromagnetic transport.

[0255] The dose contribution of each sub-particle replica is calculated separately and then superimposed according to the weight allocation to obtain the total dose distribution.

[0256] The splitting conditions include one or more of the following:

[0257] Regional conditions: Particles enter the target area, endanger organs, or are located in a pre-defined critical region;

[0258] Energy condition: Particle energy is below a set threshold;

[0259] Depth condition: The particle's depth exceeds the Bragg peak position;

[0260] Random conditions: Triggered randomly with a preset probability.

[0261] The task allocation method for the sub-particle replicas includes:

[0262] Nuclear Reaction Dedicated Copy: Specifically designed to simulate nuclear reaction processes, forcibly triggering nuclear reaction events during strong interaction calculations;

[0263] Electromagnetic transport-specific copy: Only performs electromagnetic transport simulations and does not trigger nuclear reaction events;

[0264] Hybrid replicas: Maintain standard transport patterns and allow nuclear reactions to occur with true probability.

[0265] The weight allocation method includes:

[0266] Equal weight allocation: The weights of each sub-particle copy are equal, equal to the weight of the original particle divided by the splitting factor;

[0267] Importance-weighted allocation: nuclear reaction copies are assigned higher weights, and electromagnetic copies are assigned lower weights.

[0268] The method also includes setting an upper limit for the splitting depth, and stopping splitting after the upper limit is reached to avoid an explosive increase in the number of particles.

[0269] The step of calculating the dose contribution of each sub-particle replica separately and then superimposing them according to the weight allocation includes:

[0270] Each sub-particle replica independently undergoes subsequent transport simulations. During dose accumulation, the energy deposited by each replica is multiplied by its weight.

[0271] The weighted dose contributions of all copies are summed to the corresponding voxel to obtain the final dose of that voxel.

[0272] This method can be used in conjunction with a phased computational framework. Within the phased framework, path splitting can be applied to the strong interaction computation process to further enhance the sampling of nuclear reaction events; it can also be applied to the transition between the electromagnetic and strong interaction phases to achieve synergistic optimization of the two types of processes. This synergistic effect further improves the overall computational efficiency.

[0273] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0274] Electromagnetic interaction dose and strong interaction dose are calculated using deterministic transport methods and Monte Carlo methods, respectively, and the total dose is synthesized through a coupling model; the deterministic transport methods include continuous slowing approximation or analytical models.

[0275] The deterministic transport method includes a continuously slowing approximation or analytical model, specifically including:

[0276] The continuous slowing-down approximation is used to calculate the energy loss of particles in the medium, and a continuous field of electromagnetic dose distribution is obtained.

[0277] An analytical model is used to handle the multiple scattering effect and describe the lateral diffusion of the particle beam;

[0278] The penetration depth of particles can be quickly determined using the range-energy relationship.

[0279] The Monte Carlo method was used to simulate strong interaction doses, specifically including:

[0280] Based on the information of the strong interaction cross section, determine the probability of a particle undergoing a nuclear reaction;

[0281] Whether a strong interaction event occurs is determined by random sampling;

[0282] Reactor sampling, energy sampling, and angle sampling are performed on the events that occur to simulate the generation and transport of secondary particles;

[0283] Record the energy of strong interaction processes and the deposition of secondary particles.

[0284] The synthesis of the total dose through a coupling model specifically includes:

[0285] The electromagnetic dose distribution calculated by deterministic methods is used as the basic dose field;

[0286] The dose contribution of strong interactions calculated by the Monte Carlo method is superimposed on the baseline dose field;

[0287] The electromagnetic dose is attenuated and corrected based on the probability of nuclear reaction to avoid double counting;

[0288] The statistical fluctuations in the Monte Carlo simulation were smoothed to obtain the final dose distribution.

[0289] The coupling model also includes dose weighting:

[0290] For particles that do not undergo nuclear reactions, their dose contribution is calculated entirely by deterministic methods;

[0291] For particles that undergo nuclear reactions, their dose contribution includes the electromagnetic energy deposition before the nuclear reaction and the energy deposition of the nuclear reaction itself, which are calculated by deterministic and Monte Carlo methods respectively and then superimposed.

[0292] The method further includes: cross-validating the calculation results of the deterministic method and the Monte Carlo method, and adjusting the coupling parameters according to the validation results to ensure the consistency of dose calculation.

[0293] The effects of the above technical solution are as follows:

[0294] By treating electromagnetic interactions separately from strong interactions, the electromagnetic component can be calculated independently using deterministic methods (continuous slowing approximation, analytical models), avoiding the complex branching structure of determining nuclear reactions at each step in the traditional Monte Carlo method. The electromagnetic dose calculation time is reduced from several minutes to several seconds, and the overall calculation speed is several times to an order of magnitude faster than the pure Monte Carlo method, fully demonstrating the computational efficiency advantage of this differentiated treatment.

[0295] By separating and treating the strong interaction from the electromagnetic interaction, the Monte Carlo method can be used to accurately simulate nuclear reaction events and the transport of secondary particles. This separation ensures that the nuclear reaction component is unaffected by the simplified model of the electromagnetic component, maintaining high accuracy independently and guaranteeing that the nuclear reaction dose calculation is comparable to the gold standard Monte Carlo method.

[0296] By differentiating the processing methods, computational resources can be allocated according to the characteristics of the two types of physical processes: the electromagnetic part uses a deterministic method, which has low computational resource requirements; the strong interaction part uses a Monte Carlo method, which only needs to simulate nuclear reaction events, and the number of particles required is much lower than that of the pure Monte Carlo method. This optimized resource allocation reduces the overall computational resource requirements, making it easier to deploy in a wider range of medical institutions.

[0297] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0298] The particles are divided into two groups: one group participates only in electromagnetic transport and the other group participates only in strong interaction transport. The statistical ratio of the two groups is controlled by initial particle number allocation or weight allocation, and the doses are calculated separately and then superimposed.

[0299] The particle classification is completed before the simulation begins, specifically including:

[0300] Based on a preset statistical ratio, the total number of particles is divided into electromagnetic particle groups and nuclear reaction particle groups.

[0301] Electromagnetic particle swarms only simulate electromagnetic interactions and do not judge or process any nuclear reaction events;

[0302] Nuclear reaction particle swarms are specifically designed to simulate strong interactions, using forced collisions or real free path sampling to handle nuclear reaction events.

[0303] The statistical ratio of the two types of particles is controlled by the initial particle number allocation, specifically including:

[0304] The proportion of particles in the nuclear reaction swarm is determined based on the actual probability of the nuclear reaction occurring.

[0305] The proportion of particles in a nuclear reaction swarm is equal to the probability of a nuclear reaction occurring multiplied by an amplification factor. An amplification factor greater than 1 is used to enhance nuclear reaction sampling.

[0306] The proportion of electromagnetic particles is the same as the remaining proportion, and the total number of particles remains unchanged.

[0307] The statistical ratio of the two types of particles is controlled by weight allocation, specifically including:

[0308] All particles participate in both types of transport simultaneously, but are assigned different statistical weights;

[0309] The weight of electromagnetic particles is 1, and the weight of nuclear reaction particles is the reciprocal of the probability of a nuclear reaction occurring.

[0310] By adjusting the weights, the statistical contribution of the nuclear reaction particle swarm is matched with the actual physical probability.

[0311] The process of calculating and summing the doses separately specifically includes:

[0312] Electromagnetic particle swarms are transported independently, and the electromagnetic dose deposited is recorded.

[0313] Independent transport of nuclear reaction particle swarms, recording the strong interaction doses deposited (including electromagnetic energy deposition before the nuclear reaction and energy deposition of the nuclear reaction itself).

[0314] The electromagnetic dose and the strong interaction dose are weighted and superimposed to obtain the final voxel dose.

[0315] The method also includes dynamically adjusting the particle swarm ratio:

[0316] During transport, the ratio of particle numbers in the two types of particle groups is dynamically adjusted based on changes in particle energy or regional importance.

[0317] For low-energy or critical regions, increase the proportion of nuclear reaction particle swarms and enhance nuclear reaction sampling;

[0318] For high-energy or non-critical regions, reduce the proportion of nuclear reaction particle swarms to save computational resources.

[0319] The transport of the two types of particle swarms can be performed in parallel:

[0320] The transport simulations of electromagnetic particle swarms and nuclear reaction particle swarms were executed in parallel on independent computing units;

[0321] After the simulation is completed, the dose contributions of the two parts are superimposed according to their weights to form the total dose distribution.

[0322] The effects of the above technical solution are as follows:

[0323] By dividing particles into two groups—the electromagnetic particle group, which simulates only electromagnetic interactions without requiring nuclear reaction assessments, and the nuclear reaction particle group, which specifically simulates strong interactions without handling electromagnetic transport details—this partitioning avoids the overhead of each particle simultaneously simulating two types of physical processes, as is common in traditional methods, thus improving computational efficiency.

[0324] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0325] Based on the particle energy range, the system is divided into high-energy and low-energy regions. The high-energy region is modeled using strong interactions, while the low-energy region is modeled using electromagnetic transport, with a smooth transition at the energy boundary.

[0326] The dose contributions calculated from high-energy region modeling and low-energy region modeling are superimposed according to energy segments to obtain the final dose distribution.

[0327] During the simulation, a high-energy dose counter and a low-energy dose counter were established for each voxel. The high-energy dose counter recorded the dose deposited when the particle energy was above a threshold, and the low-energy dose counter recorded the dose deposited when the particle energy was below a threshold.

[0328] For particles with energies within the transition range, their dose contributions are weighted and simultaneously recorded in both the high-energy and low-energy dose counters. The weights are determined based on the distance between the particle's energy and the threshold: the closer the energy is to the high-energy region, the greater the weight of the high-energy counter; the closer the energy is to the low-energy region, the greater the weight of the low-energy counter. The dose contributions of all particles are then accumulated to their respective counters according to their weights.

[0329] The effects of the above technical solution are as follows:

[0330] By partitioning energy levels, the high-energy region employs the Monte Carlo method to accurately simulate nuclear reactions, ensuring the accuracy of nuclear reaction dose contributions; the low-energy region uses deterministic methods to quickly calculate electromagnetic doses, avoiding the overhead of complex models. This partitioning optimizes the allocation of computing resources according to energy ranges, improving overall computing speed and achieving a balance between efficiency and accuracy.

[0331] By establishing separate high-energy and low-energy dose counters, a physical partitioning of the nuclear reaction-dominant dose and the electromagnetic-dominant dose was achieved. Outputting the high-energy dose distribution allows for the assessment of the influence of secondary particles on the dose, while outputting the low-energy dose distribution allows for the analysis of dose characteristics near the Bragg peak, providing a quantitative basis for a deeper understanding of the contribution of particles at different energy levels to the total dose.

[0332] By setting a transition interval at the energy boundary and performing weighted fusion, the modeling results of the high-energy region and the low-energy region are smoothly connected, avoiding dose discontinuity caused by zone switching. The fusion weight is dynamically adjusted according to the distance between the particle energy and the threshold, ensuring the physical rationality and visual smoothness of the dose distribution across the entire energy range.

[0333] By superimposing the dose contributions from the high-energy region and the low-energy region, a complete three-dimensional dose distribution is formed, which can be directly used for treatment plan evaluation, dose verification, and biological effect calculation. This partitioned processing not only separates the two types of physical processes but also ensures the integrity and usability of the total dose results through the superposition of the partitioned energy regions.

[0334] The separate dose calculation method described in this invention is not limited to the specific implementation described above. For those skilled in the art, without departing from the core idea of ​​this invention, other methods such as weight correction, importance sampling, path splitting, post-processing reconstruction, or deterministic and Monte Carlo coupling can be used to achieve the separate calculation of electromagnetic interaction and strong interaction and dose coupling. These equivalent alternatives should all fall within the protection scope of this invention.

[0335] Example 2: This example provides a proton or heavy ion separation dose calculation device, including a memory and a processor. The memory stores a computer program, and the processor executes the following steps when executing the computer program:

[0336] Based on the patient's CT image data, a three-dimensional voxel model of the patient was constructed, and the material physical parameters of each voxel were determined; the material physical parameters include the material composition, density, and average excitation energy of each voxel.

[0337] Electromagnetic interaction parameters and strong interaction cross-section parameters are determined based on patient CT image data and particle physics parameters; the particle beam irradiation parameters include the energy of the primary particles generated by the accelerator and the beam spot size; the electromagnetic interaction parameters include the particle ionization energy loss rate and multiple scattering angle distribution parameters.

[0338] The electromagnetic interactions and strong interactions of charged particles are treated differently to calculate the dose contribution of electromagnetic interactions and strong interactions separately; the treatment of electromagnetic interactions and strong interactions of charged particles is divided into independent and parallel calculation processes.

[0339] Based on the dose contribution of the electromagnetic interaction and the dose contribution of the strong interaction, as well as the probability of the interaction occurring, the obtained doses are probability-weighted to obtain the final three-dimensional dose distribution.

[0340] In one possible implementation, the decomposition of the electromagnetic and strong interactions of charged particles into independent and parallel computational processes includes:

[0341] The transport process of charged particles is divided into independent and parallel electromagnetic interaction calculation processes and strong interaction calculation processes;

[0342] In the electromagnetic interaction calculation process, only the electromagnetic interaction process of particles is simulated to calculate the dose deposition size of electromagnetic interaction in each voxel;

[0343] In the strong interaction calculation, only the strong interaction process of particles is simulated, and the dose deposition size of each voxel under strong interaction is calculated.

[0344] Based on the probability of electromagnetic interaction and strong interaction between particles in each voxel, the electromagnetic dose and strong interaction dose deposition are weighted by probability to obtain the final voxel dose distribution.

[0345] In one possible implementation, the dose deposition size of each voxel under the electromagnetic interaction is calculated by:

[0346] Calculate the probability P that the particle undergoes only electromagnetic interaction within the voxel. EM and the dose contribution D calculated in the electromagnetic interaction phase. EM Multiplication yields electromagnetic dose deposition;

[0347] Wherein, the probability P that only electromagnetic interaction occurs EM Calculated using an exponential decay model:

[0348]

[0349] Σ had denoted as the macroscopic cross section of the strong interaction, and x represents the path length of the particle in the voxel.

[0350] In one possible implementation, the strong interaction dose calculation step includes:

[0351] The strong interaction process is divided into elastic scattering process and inelastic nuclear reaction process;

[0352] Based on the relative proportions of the elastic scattering cross section and the inelastic reaction cross section, the type of strong interaction is sampled using conditional probability to determine whether the particle undergoes elastic scattering or inelastic nuclear reaction, and the dose contribution of elastic scattering and inelastic nuclear reaction is calculated respectively.

[0353] Under conditions where strong interactions occur, the conditional probability P of elastic scattering ela It equals the elastic scattering cross section divided by the total strong interaction cross section, i.e., P ela =Σela / Σhad, the conditional probability P of an inelastic nuclear reaction nonela It equals the inelastic nuclear reaction cross section divided by the total strong interaction cross section, i.e., P. nonela =Σnonela / Σhad, and satisfy P ela +P nonela =1; the particle undergoes elastic scattering or inelastic nuclear reaction based on the conditional probability through random sampling.

[0354] In one possible implementation, a forced collision method is used during the strong interaction calculation, so that each particle experiences a strong interaction event once during the calculation phase; and a statistical weight correction factor is introduced to compensate for the probability deviation caused by forced collision sampling. The weight correction factor is determined based on the actual strong interaction probability of the particles in the voxel.

[0355] In one possible implementation, the probability of strong particle interactions occurring in a voxel is calculated using an exponential decay model based on a macroscopic sampling cross-section:

[0356]

[0357]

[0358] Where: P had Σ represents the probability of strong interactions between particles in a voxel; had For the macroscopic cross section of strong interaction, Σ ela For elastic collision sections; Σ nonela is the inelastic nuclear reaction cross section; x is the path length of the particle in the voxel;

[0359] The final voxel dose is calculated using the following formula:

[0360]

[0361] Where D is the final voxel dose; D EM The dose calculated for the electromagnetic interaction phase; D ela D is the dose of elastic scattering in strong interactions; nonels The dose of inelastic nuclear reactions in strong interactions; P EM P represents the probability of electromagnetic interaction occurring. had P represents the probability of a strong interaction occurring. ela P represents the conditional probability of elastic scattering. nonela The conditional probability of an inelastic nuclear reaction; under the condition that strong interactions occur, P satisfies... ela +P nonela =1.

[0362] In one possible implementation, the electromagnetic interaction dose is the electromagnetic dose of the particle under conditions where no strong interaction occurs; the strong interaction dose includes the electromagnetic dose generated by the electromagnetic interaction during the transport process of the particle before the strong interaction occurs, and the dose generated by the secondary particles generated by the strong interaction during subsequent transport.

[0363] The secondary particles produced by the strong interaction are determined in the following way:

[0364] Based on particle physics properties and voxel parameters, information about the reaction channels in which nuclear reactions occur is obtained;

[0365] The types of secondary particles are obtained by sampling the reaction channel using a probability distribution.

[0366] The reaction channel was obtained by sampling the probability distribution of the reaction cross section;

[0367] The energy of the secondary particles is obtained by sampling the differential cross section;

[0368] The emission angle of secondary particles is obtained by sampling double differential cross sections or angular distribution data.

[0369] In one possible implementation, the electromagnetic interaction calculation process and the strong interaction calculation process are executed independently on a parallel computing platform, so that particles in the same calculation stage perform the same type of physical process; thereby reducing thread branch divergence and achieving parallel acceleration of the particle transport process.

[0370] In one possible implementation, the strong interaction calculation process employs importance sampling, weight correction, and / or path splitting methods to achieve statistical enhancement of strong interaction events.

[0371] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0372] Instead of explicitly separating electromagnetic and strong interaction processes, the sampling efficiency of strong interaction events is improved by introducing importance sampling or bias sampling of nuclear reaction-related events in a unified transport process.

[0373] or,

[0374] The energy deposition of each physical process is recorded during the unified simulation process, and the dose is decomposed and reconstructed according to the occurrence probability of each process in the post-processing stage.

[0375] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0376] During particle transport, when a particle enters a specific region or reaches a certain energy condition, its trajectory splits. One part of the particles are used to simulate nuclear reaction processes, while the other part is used for electromagnetic transport.

[0377] In one possible implementation, the decomposition of the electromagnetic and strong interactions of charged particles into independent and parallel computational processes includes:

[0378] Electromagnetic interaction dose and strong interaction dose are calculated using deterministic transport methods and Monte Carlo methods, respectively, and the total dose is synthesized through a coupling model; the deterministic transport methods include continuous slowing approximation or analytical models.

[0379] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0380] The particles are divided into two groups: one group participates only in electromagnetic transport and the other group participates only in strong interaction transport. The statistical ratio of the two groups is controlled by initial particle number allocation or weight allocation, and the doses are calculated separately and then superimposed.

[0381] In one possible implementation, the distinction between the electromagnetic interactions and strong interactions of charged particles includes:

[0382] Based on the particle energy range, the system is divided into high-energy and low-energy regions. The high-energy region is modeled using strong interaction, while the low-energy region is modeled using electromagnetic transport, with a smooth transition at the energy boundary.

[0383] The specific implementation methods and effects of the above technical solutions are the same as those described in Example 1, and will not be repeated here.

[0384] This invention also provides a radiotherapy system, comprising:

[0385] Particle accelerators are used to generate proton beams or heavy ion beams.

[0386] The proton or heavy ion separation dose calculation device described in Example 2 is used to calculate the three-dimensional dose distribution of a patient.

[0387] This invention also provides an electronic device, which includes a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the steps of the method described in Embodiment 1 of this invention or the function of the device described in Embodiment 2 of this invention.

[0388] This invention also provides a computer-readable storage medium for storing a computer program. When the computer program is executed, it implements the steps of the method in Embodiment 1 of this invention. The specific implementation method is consistent with the implementation method and the technical effects achieved in the above method embodiments, and some contents will not be repeated.

[0389] In this invention, a readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device. The program product can take the form of any combination of one or more readable media. A readable medium can be a readable signal medium or a readable storage medium. A readable storage medium can be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: an electrical connection having one or more wires, a portable disk, a hard disk, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage device, magnetic storage device, or any suitable combination thereof.

[0390] Computer-readable storage media may include data signals propagated in baseband or as part of a carrier wave, carrying readable program code. Such propagated data signals may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. The readable storage medium may also be any readable medium capable of sending, propagating, or transmitting a program for use by or in conjunction with an instruction execution system, apparatus, or device. The program code contained on the readable storage medium may be transmitted using any suitable medium, including but not limited to wireless, wired, optical fiber, RF, or any suitable combination thereof. Program code for performing operations of the present invention may be written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Java and C++, as well as conventional procedural programming languages ​​such as C or similar programming languages. The program code may be executed entirely on a user computing device, partially on an associated device, as a standalone software package, partially on a user computing device and partially on a remote computing device, or entirely on a remote computing device or server. In cases involving remote computing devices, the remote computing devices can be connected to user computing devices via any type of network, including local area networks (LANs) or wide area networks (WANs), or they can be connected to external computing devices (e.g., via the Internet using an Internet service provider).

[0391] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the invention without departing from the principles and spirit of the invention, and all such changes should fall within the protection scope of the claims of the present invention.

Claims

1. A proton or heavy ion separation-based dose calculation method, characterized in that, The method includes: Based on the patient's CT image data, a three-dimensional voxel model of the patient was constructed, and the material physical parameters of each voxel were determined; Electromagnetic interaction parameters and strong interaction section parameters were determined based on patient CT image data and particle physics parameters. The electromagnetic and strong interactions of charged particles are decomposed into independent and parallel computational processes to calculate the dose contribution of electromagnetic interaction and the dose contribution of strong interaction separately. Based on the dose contribution of the electromagnetic interaction and the dose contribution of the strong interaction, as well as the probability of the interaction occurring, the obtained doses are probability-weighted to obtain the final three-dimensional dose distribution.

2. The method according to claim 1, characterized in that, The process of decomposing the electromagnetic and strong interactions of charged particles into independent and parallel computational processes includes: The transport process of charged particles is divided into independent and parallel electromagnetic interaction calculation processes and strong interaction calculation processes; In the electromagnetic interaction calculation process, only the electromagnetic interaction process of particles is simulated to calculate the dose deposition size of electromagnetic interaction in each voxel; In the strong interaction calculation, only the strong interaction process of particles is simulated, and the dose deposition size of each voxel under strong interaction is calculated. Based on the probability of electromagnetic interaction and strong interaction between particles in each voxel, the electromagnetic dose and strong interaction dose deposition are probability-weighted to obtain the final voxel dose distribution.

3. The method according to claim 2, characterized in that, The dose deposition size of each voxel under electromagnetic interaction is calculated in the following ways: The probability that a particle will only undergo electromagnetic interaction in a voxel is calculated and multiplied by the dose contribution calculated in the electromagnetic interaction phase to obtain the electromagnetic dose deposition. Wherein, the probability P that only electromagnetic interaction occurs EM Calculated using an exponential decay model: Σ had denoted as the macroscopic cross section of the strong interaction, and x represents the path length of the particle in the voxel.

4. The method according to claim 2, characterized in that, The strong interaction dose calculation steps include: The strong interaction process is divided into elastic scattering process and inelastic nuclear reaction process; Based on the relative proportions of the elastic scattering cross section and the inelastic reaction cross section, the type of strong interaction is sampled using conditional probability to determine whether the particle undergoes elastic scattering or inelastic nuclear reaction, and the dose contribution of elastic scattering and inelastic nuclear reaction is calculated respectively.

5. The method according to claim 2, characterized in that, In the strong interaction calculation process, a forced collision method is adopted to ensure that each particle experiences a strong interaction event once during the calculation phase; and a statistical weight correction factor is introduced to compensate for the probability deviation caused by the forced collision. The weight correction factor is determined based on the actual strong interaction probability of the particle in the voxel.

6. The method according to claim 5, characterized in that, The probability of strong interactions between particles in voxels is calculated using an exponential decay model based on macroscopic sampling cross-sections: Where: P had Σ represents the probability of strong interactions between particles in a voxel; had For the macroscopic cross section of strong interaction, Σ ela For elastic collision sections; Σ nonela is the inelastic nuclear reaction cross section; x is the path length of the particle in the voxel; The final voxel dose is calculated using the following formula: Where D is the final voxel dose; D EM The dose during the electromagnetic interaction phase; D ela D is the dose of elastic scattering in strong interactions; nonela The dose of inelastic nuclear reactions in strong interactions; P EM P represents the probability of electromagnetic interaction occurring. had P represents the probability of a strong interaction occurring. ela P represents the conditional probability of elastic scattering. nonela The conditional probability of an inelastic nuclear reaction; under the condition that strong interactions occur, P satisfies... ela +P nonela =1.

7. The method according to claim 2, characterized in that, The electromagnetic interaction dose is the electromagnetic dose of the particle under conditions where no strong interaction occurs; the strong interaction dose includes the electromagnetic dose generated by the electromagnetic interaction during the transport process of the particle before the strong interaction occurs, and the dose generated by the secondary particles generated by the strong interaction during subsequent transport. The secondary particles produced by the strong interaction are determined in the following way: Based on particle physics properties and voxel parameters, information about the reaction channels in which nuclear reactions occur is obtained; The types of secondary particles are obtained by sampling the reaction channel using a probability distribution. The energy of the secondary particles is obtained through differential cross-section sampling; The emission angle of secondary particles is obtained by sampling double differential cross sections or angular distribution data.

8. The method according to claim 2, characterized in that, The electromagnetic interaction calculation process and the strong interaction calculation process are executed independently on a parallel computing platform, so that particles in the same calculation stage perform the same type of physical process; This reduces thread branch divergence and accelerates the parallel processing of particle transport.

9. The method according to claim 2, characterized in that, The strong interaction calculation process employs importance sampling, weight correction, and / or path splitting methods to achieve statistical enhancement of strong interaction events.

10. The method according to claim 2, characterized in that, The electromagnetic interaction calculation process is implemented using deterministic or semi-analytical methods.

11. The method according to claim 1, characterized in that, The process of decomposing the electromagnetic and strong interactions of charged particles into independent and parallel computational processes includes: Electromagnetic interaction dose and strong interaction dose are calculated using deterministic transport methods and Monte Carlo methods, respectively, and the total dose is synthesized through a coupling model; the deterministic transport methods include continuous slowing approximation or analytical models.

12. The method according to claim 1, characterized in that, The process of decomposing the electromagnetic and strong interactions of charged particles into independent and parallel computational processes includes: The particles are divided into two groups: one group participates only in electromagnetic transport and the other group participates only in strong interaction transport. The statistical ratio of the two groups is controlled by initial particle number allocation or weight allocation, and the doses are calculated separately and then superimposed.

13. A proton or heavy ion separation-type dose calculation device, characterized in that, Includes a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to perform the following steps: Based on the patient's CT image data, a three-dimensional voxel model of the patient was constructed, and the material physical parameters of each voxel were determined; Electromagnetic interaction parameters and strong interaction section parameters were determined based on patient CT image data and particle physics parameters. The electromagnetic and strong interactions of charged particles are decomposed into independent and parallel computational processes to calculate the dose contribution of electromagnetic interaction and the dose contribution of strong interaction separately. Based on the dose contribution of the electromagnetic interaction and the dose contribution of the strong interaction, as well as the probability of the interaction occurring, the obtained doses are probability-weighted to obtain the final three-dimensional dose distribution.

14. A radiotherapy system, characterized in that, include: Particle accelerators are used to generate proton beams or heavy ion beams. The proton or heavy ion separation dose calculation device as described in claim 13 is used to calculate the three-dimensional dose distribution of a patient.

15. A computer-readable storage medium, characterized in that, The storage medium stores computer instructions, and when the computer reads the computer instructions, the computer implements the steps of the method according to any one of claims 1-12.

Images