CN117873963A - Method for converting phase space file into input file usable by other programs - Google Patents

Abstract

The invention discloses a method for converting a phase space file into an input file usable by other programs, which comprises the following steps: reading phase space files generated by different programs, and generating an accurate geometric model; defining physical properties for each tissue and substance in the geometric model; coupling calculation is carried out through a definite theory method and a probability theory method, and a Boltzmann neutron transport equation is solved to obtain flux distribution; in the simulation process, recording the interaction condition of each particle and the substance, and calculating the distribution condition of the dose in different tissues and positions; verifying consistency of simulation results and experimental data, and calibrating; and outputting input files available for other programs. The invention can read the phase space files generated by different programs, couple the transportation programs according to the needed determinism and probability theory, output the needed particle information, and also output the files into the phase space file format defined by IAEA, and use the files as the input files of the programs.

Description

A method to convert phase space files into input files that can be used by other programs

Technical Field

The invention relates to the field of radiation therapy dose calculation, and in particular to a method for converting a phase space file into an input file available for other programs.

Background technique

Boron Neutron Capture Therapy (BNCT) is a tumor treatment method that achieves the purpose of tumor treatment through the high capture probability of neutrons by the boron-10 isotope. The basic principle of this therapy is that the patient is first injected with a compound containing the boron-10 isotope. Then, the patient is irradiated with a neutron beam, which will react with the boron-10 to release an alpha particle and a lithium ion, which has a high local killing effect on tumor tissue. Due to the high capture cross-section (capture probability) of the boron-10 isotope for neutrons and the high energy release of the reaction products, BNCT has potential advantages in tumor treatment. Because the neutron beam only reacts with tumor cells carrying boron-10, more precise tumor treatment can be achieved compared to surrounding normal tissue, thereby reducing damage to healthy tissue.

Phase space files are files used to store state information of particles or systems, which usually contain state information such as the position, momentum, energy, weight, and other possible related information of each particle in the system. These files are used to record the motion trajectory and state changes of particles during the simulation so that they can be analyzed and visualized after the simulation. The content of the phase space file can vary according to the needs of the system and simulation. According to the phase space data of the International Atomic Energy Agency (IAEA), the phase space file at least includes information such as the type of particle, the energy of the corresponding particle, direction, position, statistical weight, and history number. By using phase space files, simulation calculations can be performed under the same conditions of the particle source to obtain the dose received by the patient, which is very helpful for formulating treatment plans, optimizing dose distribution, and evaluating the effects of different treatment plans.

Although the functionality to create phase space files is widely implemented in most probabilistic computational programs, each code uses its own specific format, which limits compatibility and reproducibility between different codes. To address this issue, the IAEA has defined a standard phase space file format that can be read and written using a database it provides. Using this format, the IAEA has created a common database (IAEA NAPCNuclear Data Section 2020) for phase space files for linear accelerators used in external radiation therapy by compiling existing data that has been appropriately validated.

Different equipment manufacturers have different input formats for phase space files. When hospitals or researchers need to use equipment from different manufacturers for different purposes such as calibration or dose estimation in radiotherapy planning, and considering that the manufacturer's information is usually subject to strict confidentiality agreements, it is necessary to convert the phase space files into input files that can be used by other programs.

Summary of the invention

In order to overcome the deficiencies of the prior art, the present invention provides a method for converting a phase space file into an input file that can be used by other programs.

The technical solution of the present invention is implemented as follows: a method for converting a phase space file into an input file that can be used by other programs, comprising the steps of:

S1, reading phase space files generated by different programs to generate an accurate geometric model, wherein the geometric model includes a neutron source, capture agent distribution, tumor tissue, and surrounding normal tissue;

S2, for each tissue and substance in the geometric model, define the physical properties of nuclear reaction cross section, interaction probability and energy transfer;

S3, by coupling the deterministic method with the probabilistic method, solve the Boltzmann neutron transport equation and obtain the flux distribution, including the following steps:

S31, uses deterministic methods to quickly calculate the distribution of neutron flux;

S32, using the distribution of the neutron flux as a value function to guide the Monte Carlo program simulation, and using a Monte Carlo method based on probability theory to simulate the neutron motion state;

S4, during the simulation, records the interaction between each particle and the material and calculates the distribution of dose in different tissues and locations;

S5, verify the consistency of simulation results with experimental data or other calculation methods, and perform calibration;

S6, outputs input files that can be used by other programs.

Furthermore, the phase space files generated by the different programs read in S1 include but are not limited to MCNP files, PHITS files, IAEAheader and IAEAphsp files.

Furthermore, the step S31 includes the steps of:

S311, approximation of the Boltzmann neutron transport equation;

S312, find the approximate or analytical solution to the Boltzmann neutron transport equation.

Furthermore, the method for approximating the Boltzmann neutron transport equation in S311 includes but is not limited to multi-group approximation of the interface and Legendre expansion.

Furthermore, the method for obtaining an approximate solution or an analytical solution to the Boltzmann neutron transport equation in S312 includes but is not limited to a discrete coordinate method, a spherical harmonic function method and a finite element method.

Furthermore, the S32 includes the step of: adding a weight window parameter in the Monte Carlo simulation, specifically including the steps of:

S321, calculating an importance map based on the grid divided in the Monte Carlo simulation, wherein the importance map is defined as Wherein, (i,j,k) is the index of the grid divided in the Monte Carlo simulation, g is the particle energy group, I(i,j,k,g) is the importance of the grid (i,j,k,g), R _score (i,j,k,g) is the response of all detectors after the particle enters the grid (i,j,k,g), W _total (i,j,k,g) is the weighted sum of all particles entering the grid (i,j,k,g), and R _score (i,j,k,g) is the flux obtained by the deterministic calculation in step S31/> It is obtained that W _total (i,j,k,g) is obtained from the particle information recorded by the deterministic calculation in step S31;

S322, define the expected weight of each grid as Where/> is the average response of each particle. A window of desired weights is generated for each grid;

S323, compare the actual particle weight with the weight window boundary of splitting and roulette. If the actual weight W _actual is greater than the threshold f _upper W _expect , the particle performs splitting simulation; if the actual weight W _actual is less than the threshold f _lower W _expect , the particle needs to perform roulette to reduce simulation; if the actual weight W _actual is between the thresholds f _upper W _expect and f _lower W _expect , the particle continues to simulate, where f _upper and f _lower are constants used for threshold setting, selected from (1,+∞) and (0,1) respectively.

The beneficial effect of the present invention is that, compared with the prior art, the present method can read phase space files generated by different programs, and can also read files provided by the IAEA database, read the particle information stored in the files, output the required particle information according to the required deterministic and probabilistic coupled transport programs, and can also output the files into the phase space file format defined by the IAEA for use as input files of the program.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart of a method for converting a phase space file into an input file that can be used by other programs according to the present invention.

Detailed ways

The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

Please refer to FIG1 . The method for converting a phase space file into an input file available for other programs of the present invention comprises the following steps: A method for converting a phase space file into an input file available for other programs comprises the following steps:

S1, read the phase space files generated by different programs to generate an accurate geometric model, which includes neutron source, capture agent distribution, tumor tissue and surrounding normal tissue, etc. This model needs to reflect the actual situation as accurately as possible;

S2, for each tissue and substance in the geometric model, defining its physical properties such as nuclear reaction cross section, interaction probability and energy transfer, and this information is used to simulate the interaction between neutrons and matter;

S3, by coupling the deterministic method with the probabilistic method, solve the Boltzmann neutron transport equation and obtain the flux distribution, including the following steps:

S31, uses deterministic methods to quickly calculate the distribution of neutron flux;

S32, using the distribution of the neutron flux as a value function to guide the Monte Carlo program simulation, and using a Monte Carlo method based on probability theory to simulate the neutron motion state;

S4, during the simulation, records the interaction between each particle and matter, including energy deposition and transfer, and based on this information, calculates the distribution of dose in different tissues and locations;

S5, verify the consistency of simulation results with experimental data or other calculation methods, and perform calibration to ensure the accuracy of simulation results;

S6, outputs input files that can be used by other programs.

The core of dose calculation is to solve the Boltzmann neutron transport equation to obtain the flux distribution and provide a basis for treatment planning. The steady-state form of the Boltzmann neutron transport equation is:

in, is the angular flux, /> is the total cross section, /> is the scattering cross section, In step S3, the present invention obtains the Boltzmann neutron transport equation by coupling calculation using a deterministic method and a probabilistic method, including the steps of:

S31, uses deterministic methods to quickly calculate the distribution of neutron flux;

S32, using the distribution of the neutron flux as a value function to guide the Monte Carlo program simulation, and using the Monte Carlo method of probability theory to simulate the neutron motion state.

The step S31 comprises the steps of:

S311, approximation of the Boltzmann neutron transport equation;

S312, find the approximate or analytical solution to the Boltzmann neutron transport equation.

In step S311, because the neutron transport equation has no clear analytical solution, it is necessary to approximate and convert it to obtain an analytical solution. In step S31, the approximate methods that can be used include multi-group approximation of the cross section and Legendre expansion. Among them, Legendre expansion is: If the scattering is only related to the angle, the scattering cross section can be expanded by the L-order Legendre polynomial:

The multigroup approximation of the cross section is: the continuous energy cross section is divided into some discrete energy group structures, where each group has a different energy width.

The group flux expressed with the index g becomes

The total cross section of the group represented by the index g is

The total cross section of the group represented by the index g is

in

It can be seen that if we want to find the flux In practical applications, in order to correctly estimate the average group constant, it is necessary to use a certain approximate weight spectrum.

The approximated Boltzmann neutron transport equation becomes

By treating the scattering of other groups as valid source terms, the multigroup equations can be solved sequentially as a set of valid single-group problems, namely:

The Boltzmann neutron transport equation in the energy level group represented by g is

In step S312, the methods for approximate solution or analytical solution of Boltzmann neutron transport equation include discrete coordinate method, spherical harmonic function method and finite element method. In some embodiments of the present invention, the discrete coordinate method, also known as _SN method, is used to discretize the angle into a small number of specific directions, and then solve the particle transport equation in each direction to obtain the spatial distribution of flux. Specifically:

If there is no previous solution for the angular flux, we can assume that the direction cosines are symmetrically distributed, so that we only need to consider the distribution of the unit sphere in one quadrant. Discretization into n directions requires selecting on the coordinate axis points, and the direction cosine sets are the same, that is, {μ _n } = {η _n } = {ξ _n }, and the sum of the squares of the direction cosines is 1, so it can be deduced that in the set {α _n }, there is

Subtracting the two equations, we can get

Right now

There is a Make/>

Lian Li When C is obtained from formula (1.13), we can substitute it into formula (1.13) to get

This shows that the set {α _n } is mainly determined by α _1. By setting an appropriate α ₁ , the accuracy of the flux distribution can be improved.

The application of the Monte Carlo method of probability theory in radiotherapy can take into account complex physical processes and geometric structures, thereby providing more accurate predictions of dose distribution. The motion of particles can be directly simulated through probability theory, and the distribution of particle flux can be obtained after statistics. When performing the simulation, an initial state must be set first, that is, the spatial position, energy and direction distribution of the neutron source must be obtained, that is, And sample it as the initial state. Knowing the initial state, to obtain the change of the particle's motion state over time, it is necessary to calculate the interaction between particles. In the process represented by n-1, the transport length l of the neutron obeys the following distribution

The location where the particles interact can be obtained. Next, determine the type of particles that interact. If a medium is composed of two types of atoms, A and B, its total macroscopic cross section is

in, and/> is the total macroscopic cross section of two atoms A and B. The probability of a neutron colliding with A and B is

Then, the type of interaction between particles and the corresponding cross-sectional size are determined, such as elastic scattering, inelastic scattering, fission reaction and capture. The probability of a certain interaction is similar to the probability of the atomic type, which is the ratio of its cross-sectional size to the total cross-sectional size. Finally, the motion state of the particle after the interaction is determined. If the neutron is absorbed, there is no subsequent motion state; if elastic scattering occurs, the energy of the neutron can be calculated.

in A is the mass ratio of the interacting atom and neutron, and θ _C is the deflection angle. After these parameters are calculated, the motion state of the particle per unit time can be obtained. Repeating this step can obtain the desired particle motion state.

The approximate solution calculated by the deterministic method has a certain deviation from the exact solution. The probability method has higher calculation accuracy, but it takes a lot of time to simulate. If the deterministic method is used to quickly calculate the distribution of neutron flux and use it as a value function to guide the Monte Carlo program simulation, the sampling scheme can be changed to sample more from phase space coordinates that contribute significantly to the results, and the efficiency of the simulation can be improved while ensuring reliable accuracy, that is, increasing the weight window parameter in the Monte Carlo simulation.

In step S32, a weight window parameter is added in the Monte Carlo simulation, specifically comprising the steps of:

S321, calculating an importance map based on the grid divided in the Monte Carlo simulation, wherein the importance map is defined as Wherein, (i,j,k) is the index of the grid divided in the Monte Carlo simulation, g is the particle energy group, I(i,j,k,g) is the importance of the grid (i,j,k,g), R _score (i,j,k,g) is the response of all detectors after the particle enters the grid (i,j,k,g), W _total (i,j,k,g) is the weighted sum of all particles entering the grid (i,j,k,g), and R _score (i,j,k,g) is the flux obtained by the deterministic calculation in step S31/> It is obtained that W _total (i,j,k,g) is obtained from the particle information recorded by the deterministic calculation in step S31;

S322, define the expected weight of each grid as (0.25), where/> is the average response of each particle. A window of desired weights is generated for each grid;

S323, compare the actual particle weight with the weight window boundary of splitting and roulette. If the actual weight W _actual is greater than the threshold f _upper W _expect , the particle performs splitting simulation; if the actual weight W _actual is less than the threshold f _lower W _expect , the particle needs to perform roulette to reduce simulation; if the actual weight W _actual is between the thresholds f _upper W _expect and f _lower W _expect , the particle continues to simulate, where f _upper and f _lower are constants used for threshold setting, selected from (1,+∞) and (0,1) respectively.

The present invention can read the phase space files generated by different programs, such as MCNP, PHITS, etc., and can also read the files (IAEAheader and IAEAphsp files) provided by the IAEA database, and read the particle information stored in the files. If it is necessary to rotate, translate, etc. the surface source, it is necessary to input the rotation angle and the translation displacement, otherwise the original data is retained. Finally, according to the required deterministic and probabilistic coupled transport programs, the information of the required particles, such as type, energy, direction, position and statistical weight, etc., can be output, and the file can also be exported into the phase space file format defined by the IAEA, as the input file of the program. The above is a preferred embodiment of the present invention. It should be pointed out that for those of ordinary skill in the art, without departing from the principle of the present invention, some improvements and embellishments can also be made, and these improvements and embellishments are also considered as the protection scope of the present invention.

Claims (6)

1. A method for converting a phase space file into an input file usable by other programs, comprising the steps of:

s1, reading phase space files generated by different programs, and generating an accurate geometric model, wherein the geometric model comprises neutron sources, capturing agent distribution, tumor tissues and surrounding normal tissues;

s2, defining the nuclear reaction section, interaction probability and physical properties of energy transmission of each tissue and substance in the geometric model;

s3, solving a Boltzmann neutron transport equation by coupling calculation of a definite theory method and a probability theory method to obtain flux distribution, wherein the method comprises the following steps:

s31, rapidly calculating neutron flux distribution by using a definite theory method;

s32, taking the distribution of neutron flux as a cost function for guiding Monte Carlo program simulation, and simulating neutron motion state by using a Monte Carlo method of probability theory;

s4, in the simulation process, recording the interaction condition of each particle and the substance, and calculating the distribution condition of the dose in different tissues and positions;

s5, verifying consistency of simulation results and experimental data or other calculation methods, and calibrating;

s6, outputting input files available for other programs.

2. The method of claim 1, wherein the phase space files generated by the different programs read in S1 include, but are not limited to, MCNP files, bits files, iaeahead files, and IAEAphsp files.

3. The method for converting a phase space file into an input file usable by other programs according to claim 1, wherein said step S31 comprises the steps of:

s311, approximating a Boltzmann neutron transport equation;

s312, solving an approximate solution or an analytical solution of the Boltzmann neutron transport equation.

4. The method of converting a phase space file into an input file usable by other programs according to claim 3, wherein the method of approximating the Boltzmann neutron transport equation in S311 includes, but is not limited to, a multi-cluster approximation of an interface and Legendre expansion.

5. A method of converting a phase space file into an input file usable by other programs as claimed in claim 3, wherein the method of approximating or resolving the Boltzmann neutron transport equation in S312 includes, but is not limited to, a discrete coordinate method, a spherical harmonic method, and a finite element method.

6. The method for converting a phase space file into an input file usable by other programs according to claim 1, wherein said S32 comprises the steps of: the method for adding the weight window parameters in the Monte Carlo simulation specifically comprises the following steps:

s321, calculating an importance map based on grids divided in Monte Carlo simulation, wherein the importance map is defined asWhere (I, j, k) is the index of the grid divided in the Monte Carlo simulation, g is the particle energy group, I (I, j, k, g) is the importance of the grid (I, j, k, g), R _score (i, j, k, g) is the response of all detectors after the particle enters the grid (i, j, k, g), W _total (i, j, k, g) is the weight addition of all particles into the grid (i, j, k, g), R _score (i, j, k, g) flux +.>Obtained, W _total (i, j, k, g) from the particle information recorded by the definitive calculation in step S31;

s322, defining the expected weight of each grid asWherein->Is the average response of each particle. Generating a window of desired weight for each grid;

s323, comparing the actual particle weight with the weight window boundary of the split and roulette, if the actual weight W _actual Greater than threshold f _upper W _expect Then the particles undergo fragmentation simulation; if the actual weight W _actual Less than threshold f _lower W _expect The particle would need to make roulette to reduce the simulation; if the actual weight W _actual At threshold f _upper W _expect And f _lower W _expect In the middle, the particles continue to simulate, where f _upper And f _lower Is a constant for the threshold setting, in the respective steps (1), + -infinity) and (0, 1 is selected from the group consisting of.