Correlated histogram representation of Monte Carlo derived medical accelerator photon-output phase space

[0002] This patent application claims priority to U.S. Provisional Patent Application Serial No. 60/123,242 filed Mar. 8, 1999.

[0004] 1. Field of the Invention

[0005] The present invention relates radiation therapy, and more specifically it relates to a method for condensing the photon energy and angular distributions obtained from Monte Carlo simulations of medical accelerators.

[0006] 2. Description of Related Art

[0007] Monte Carlo methods are valuable for radiation dose calculations because of their ability to predict dose accurately for the entire range of conditions encountered in radiation treatment of cancer. They can correctly characterize beams of electrons and photons emerging from a linear accelerator, as well as photons, electrons and alpha particles emerging from brachytherapy sources. They can further accurately simulate the transport of particles through beam shaping devices or, for brachytherapy, source shielding devices. The transport of the resulting particles through the patient may then similarly be simulated, taking into account surface irregularities and the internal structure of the patient, including any metallic implants and prostheses, to determine the dose deposited in the patient. This disclosure provides photon treatments with linear accelerators.

[0008] While Monte Carlo methods are the most accurate means of predicting radiation dose, the use of generalized Monte Carlo software packages for routine clinical applications is presently impractical due partly to computer constraints, in particular, the running time needed to obtain a converged calculation. The lengthy running time of generalized Monte Carlo software packages warrants the development of algorithms optimized specifically for the problem at hand. The problem may be divided into three steps: (1) Determining the characteristics of the photon radiation emanating from the accelerator head; (2) propagating this radiation down to the patient, through beam modifiers (if any) and through the air above the patient; and (3), the actual transport within the patient of the original radiation, as well as of any radiation produced or scattered between the accelerator head and the patient.

[0009] It is an object of the present invention to provide, inter alia, a method for condensing the photon energy and angular distributions obtained from Monte Carlo simulations of medical accelerators. This method represents the output as a series of correlated histograms and as such is well-suited for inclusion as the photon-source package for Monte Carlo codes used to determine the dose distributions in photon teletherapy.

[0010] The method accounts for the isocenter-plane variations of the photon energy spectral distributions with increasing distance from the beam central axis for radiation produced in the bremsstrahlung target as well as for radiation scattered by the various treatment machine components within the accelerator head. Comparison of the isocenter energy fluence computed by this algorithm with that of the underlying full-physics Monte Carlo photon phase-space indicates that energy fluence errors are less than 1% of the maximum energy fluence for a range of open-field sizes. Comparison of jaw-edge penumbrae shows that the angular distributions of the photons are accurately reproduced. The Monte Carlo sampling efficiency (the fraction of generated photons which clear the collimator jaws) of the algorithm is approximately 83 % for an open 10×10 field, rising to approximately 96% for an open 40×40 field. Data file sizes for a typical medical accelerator, at a given energy, are approximately 150 kB, compared to the 1 GB size of the underlying full-physics phase-space file.

[0011] FIG. 1 shows a high-energy electron beam strikes a thin, high-Z target.

[0012] FIG. 2 shows how backtracking of the photons to their point of origin indicates which portions of the accelerator head contribute to the output energy fluence.

[0013] FIG. 3 shows energy fluence distributions from a Varian 6MV accelerator.

[0014] FIGS. 4A-D shows the diameters of the subsources and their relative positions with respect to the jaw opening lead to different areas of illumination at the isocenter plane.

[0015] FIGS. 5A-D show the effect of jaw position on energy fluence from the primary collimator and the flattening filter.

[0016] FIGS. 6A-6F shows that the on-axis photon energy distributions depend strongly upon the treatment head component in which the photons are created and the energy distributions at the isocenter plane of these subsources at a distance of 20 cm from the central axis of the accelerator.

[0017] FIGS. 7A and 7B show closeup views of the probability distribution for photons from the target.

[0018] FIGS. 8A-8C show the effect of source position on photon enery.

[0019] FIGS. 9A-9E show closeup views of the radial and angular distributions used to generate the starting positions of the primary collimator photos.

[0020] FIGS. 10A-10F show closeup views of the radial and angular distributions used to generate the starting positions of the flattening filter photons.

[0021] FIGS. 11A-F, a modified form of FIG. 8, show the correlated histograms used to describe the subsources constituting a treatment machine.

[0022] FIGS. 12A-C show radial variations of the heights of the probability bins for several photon energies.

[0023] FIG. 13 shows the “knot points” of the piecewise-linear functions shown in FIG. 12 are the bin-centers and bin-heights obtained during histogramming of the raw data.

[0024] FIG. 14 shows roll-off in photon energy fluence for flattening filter photons, at field edge defined by upper jaw at the position for a 40×40 field.

[0025] FIG. 15 shows the algorithm generates photon trajectories which may intersect the edges of the jaws

[0026] FIG. 16 shows minimizing R_closestfor most of the particles of an azimuthally-anisotropic subsource leads to an off-center distribution.

[0027] FIGS. 17A and 17B show that the alias sampling technique allows any distribution to be sampled efficiently

[0028] This disclosure focuses on Step 1, as described above in the Background section, determining the characteristics of the photon radiation emanating from the accelerator head. This step may be called the “patient-independent” part of the overall calculation, since the radiation output of the accelerator head is not affected by the presence of downstream beam modifiers. Steps 2 and 3 described in the Background section will necessarily vary from treatment plan to treatment plan, and may be called “patient-dependent” portions of the overall dose calculation. Step 1 merits a detailed treatment for two reasons. The first is that it is not purposeful to re-simulate the accelerator output during each new treatment plan; the time is better spent on the patient-specific part of the problem. The second reason is practical: Monte Carlo simulations of a medical accelerator head require a prohibitively long time to perform (many tens of hours on single-CPU machines), and it is impractical to store and manipulate (on a routine basis) the resulting gigabyte-size data files.

[0029] The problem, then, becomes one of developing a method to generate photons whose energies and trajectories represent those of the underlying Monte Carlo simulation of the accelerator head. Three issues arise in this context. First, the energy distributions and angular distributions of the photons generated by any such algorithm must accurately simulate those found in the underlying gigabyte-size data files. Second, the algorithm used to generate these photons must be computationally efficient, so that the running time of the overall dose calculation is not affected. Third, if the algorithm relies on condensing the gigabyte-size files into smaller data files for use during run-time, this process of condensation should be automated so that data files may be easily produced for a number of accelerators.

[0030] Several references suggest that dose delivered to the patient should be accurate to within 5% and that the computed dose should be accurate to within 2% of the maximum dose delivered to the patient. Therefore, the phase-space information must also be at least as accurate, and any algorithm used for sampling particles from the phase-space must produce data that is equally accurate. Furthermore, the sampling algorithm must be computationally efficient so as not to affect computational efficiency and the overall run-time.

[0031] The nature of the source modeling problem is illustrated in FIG. 1. A beam of electrons strikes a thin target 10 made of a high-Z material. The resulting bremsstrahlung photons are collimated by a conical collimator 12. The photon beam passes through a flattening filter 14, which is thicker in the center and therefore attenuates the central portion of the bremsstrahlung photon distribution relative to the edges. The intended result is a flat dose distribution at a specified depth in a water phantom. Although the dose distribution is uniform, the photon energy fluence and energy spectrum vary as a function of distance from the central ray. There are two reasons for this. The first is that the energy of the bremsstrahlung photons decreases with increased scattering angle. The second is that the flattening filter removes more low-energy photons from the center of the beam than from the edge.

[0032] In addition, non-negligible amounts of radiation scatter from the collimator and the flattening filter and pass through the jaws 16 to arrive at the patient plane 18. The combined radiation field can be characterized by several distributions of bremsstrahlung and scattered photons. All of these distributions must be understood in order to develop a useful source algorithm for input into any Monte Carlo dose calculation code.

[0033] This method uses a series of correlated histograms to reproduce the phase-space information obtained from full-scale Monte Carlo simulations of accelerator treatment heads. Errors in the total energy fluence at the isocenter, when compared to that given by the underlying full-physics Monte Carlo simulation, are below 1% of the maximum energy fluence in the region of illumination. Errors for the photon energy fluence arising from any one hardware component, e.g., the target or the flattening filter, are kept well below this criterion of 1% of the maximum energy fluence in the region of illumination.

[0034] The algorithm's computational efficiency has two bases. The first is conceptual: the algorithm is designed to generate only photons which have a high probability of going through the aperture defined by the jaws. Little time is spent generating photons that will not reach the patient. The second basis is the use of alias-sampling techniques which allow the use of finely-binned distributions without a corresponding increase in the run-time of the code.

[0035] The correlated histograms used by this algorithm may be generated automatically (i.e., without user intervention) from the data files resulting from the Monte Carlo simulation of the accelerator head.

[0036] This disclosure deals with photons emanating from the accelerator head. Note that these photons, once generated, must next be transported through any patient-dependent beam modifiers as well as through the air column above the patient.

[0037] Methods and Materials

[0038] Use of Monte Carlo radiation transport simulation codes to study the output of accelerator heads is well established. A schematic of the modeling process is shown in FIG. 1. Electrons with kinetic energies at or near the nominal operating energy are incident on a target material at the top of the accelerator head. The resulting bremsstrahlung photons are tracked through the accelerator head. Secondary photons and electrons produced within the accelerator head are also tracked. Particles arriving at the bottom of the head are tallied. A “tally” consists of recording the particle's phase-space information (position, direction, energy, particle type), its Monte Carlo weight (accounting for the various physics-biasing schemes used), and any additional information desired about the particle. This might include where the particle was last scattered, the type of scattering event, how many scattering events it underwent, etc.

[0039] The simulations described herein were performed using the Monte Carlo codes BEAM96 and MCNP4B based on machine drawings and materials data supplied by the accelerator manufacturers. (For the purposes of this disclosure, patches were applied to the “as-released” MCNP4B to include “BEAM-style” LATCH capability and to use BEAM-style bremsstrahlung splitting). Approximately 5×10⁶to 2×10⁷incident electrons are used in the simulations, depending on which accelerator is being modeled. (To date we have simulated eight accelerators made by Varian, Inc. These are the 600C: 4Mv and 6MV, 600CD: 6MV, 2100C: 6MV, 8MV, 10MV, 15MV, and 18MV). Given the variance reduction schemes used (e.g., forced collisions, particle splitting), the resulting phase-space files contain information for ˜3×10⁷photons (of varying weights) and occupy ˜1 GB to ˜3 GB each, depending on the code used to generate the files.

[0040] Unless otherwise specified, the plots shown in this disclosure are all based on a 6MV Varian 2100C data set calculated with BEAM96.

[0041] Analysis and Algorithm Development

[0042] Preliminary Analysis of Phase-Space Files Resulting from Monte Carlo Simulation of Treatment Heads

[0043] The initial check of the photon phase-space file involves graphing the place of creation (or of last scattering) (x_c, y_c, z_c) of the first ˜10⁵photons in the file. Such a scatter plot of x_cvs z_cfor a stylized (fictitious) accelerator head is shown in FIG. 2. This step in the phase-space analysis serves two purposes—the first being practical, the second being conceptual. First, it is a useful check on the input to the program, since the locations of the photon creations should correlate with the physical structure of the accelerator head; and second, it shows which portions of the treatment machine head contribute to the output of the machine. For the example shown in FIG. 2, photons originating from the target “target photons”) come from a well-defined spot. Photons coming from the primary collimator (“primary collimator photons”) are fewer in number, and they tend to come from the upper edge of the collimator. Thus, the inner surface of the primary collimator is not a uniform source of photons. Rather, the primary collimator appears to be more of a “ring” source. The flattening filter is also a source of photons (“flattening filter photons”). Unlike the primary collimator, however, where the photons originate in only a small fraction of the collimator, photons originate in the entire volume of the flattening filter.

[0044] The energy fluence distributions at the isocenter plane are shown in FIG. 3. Most of the energy comes directly from the target, with contributions at the several percent level from the flattening filter and the primary collimator. This gives a more quantitative description of the relative importances of the treatment machine head components. The energy fractions coming from each treatment head component and the physics mechanisms involved are given in Table I. Bremsstrahlung is the dominant production mechanism of target photons; incoherent scattering produces most of the photon energy from the primary collimator and flattening filter.

[0045] Illumination Regions at the Isocenter Plane

[0046] Next the regions on the patient plane that are illuminated by photons from the electron target, the primary collimator, and flattening filter are established.

[0047] Target Photons

[0048] FIG. 3 shows that the energy fluence for target photons decreases quite rapidly at approximately 25 cm from the central axis. The spatial extent of the energy fluence of target photons is limited by the opening angle of the primary collimator. Beyond 29 cm from the central axis, the energy fluence curve for the target photons has decreased to below 0.3% of its maximum value. Based on the accuracy requirement of <1% given in the above, 29 cm from the central axis is the maximum radius that needs to be considered “reachable” by photons from the bremsstrahlung target.

[0049] Primary Collimator and Flattening Filter Photons

[0050] FIG. 3 shows that photons from the primary collimator and flattening filter would, in the absence of the collimating jaws, illuminate very large areas of the isocenter plane. As is shown below, however, the jaws define the maximum area illuminated by these photons.

[0051] FIGS. 4A-4C show, for a 10 cm×10 cm open field, photons originating in different components of the accelerator head (“subsources”) illuminate different but distinct regions of the isocenter plane. It should be understood that many other subsources may be used, depending upon the type of accelerator head used. As shown in FIG. 2, target photons (FIG. 4A) originate in a small (˜2 mm diameter) area, well removed from the jaws. Primary collimator photons (FIG. 4B) originate in a slightly larger source, (˜1.5 cm diameter) area which is slightly closer to the jaws. Finally, flattening filter photons (FIG. 4C) originate in an even larger (˜3 cm diameter) area which is fairly close to the jaws. FIGS. 4D includes a schematic that shows that this difference in vertical location leads to different, but still distinct areas of illumination at the isocenter plane. The positions of the edges of these areas are due to the relative locations of the treatment head components and the movable jaws. The positions of these edges move as each jaw moves. The upper jaws give slightly different edge positions than do the lower jaws because of the difference in vertical position variation of edge position as the jaw position is varied as shown in FIGS. 5A-D for both the primary collimator and the flattening filter. (These figures were generated by reading the photon trajectories from the phase-space file for this accelerator and seeing if the trajectory would intersect a jaw. The jaws were assumed to be infinitely attenuating.) The blocking effect of the jaw is clearly visible for each jaw position.

[0052] Radial Tracking Cutoffs

[0053] In keeping with the stated goal of keeping errors in energy fluence to below ˜1% of the maximum energy fluence, it is arbitrarily declared that a subsource only illuminates those areas where it contributes more than 0.2% of the maximum in-field energy fluence. This level thus defines the edges of the field beyond which Monte Carlo tracking of photons from the given subsource will not be performed. This 0.2% cut-off level is indicated by the dashed lines in each panel of FIGS. 5A-D. The resulting field-size limits are tabulated in Table II.

[0054] Table II and FIG. 5 show that when the jaws are moved to their maximum extent—corresponding to a 40×40 cm²field—no photons from the primary collimator will be tracked more than {square root over (22.631.3 cm from the beam axis on the isocenter plane. Similarly, at this same field size, no photons from the flattening filter will be tracked more than {square root over (30.540.5 cm from the beam axis on the isocenter plane. This defines the maximum radii at which the photon output from these components will be characterized.

[0055] Application to a Source Algorithm

[0056] The tabulated extents of illumination listed in Table II may be used in a source algorithm to sample x and y on the isocenter plane. For the case of a symmetric 10 cm×10 cm field (i.e., all jaws moved to their respective +5 cm setting), primary collimator photons illuminate a ∓7.0 x∓6.4≈180 cm²area. Similarly, for this same 10×10 cm field, photons from the flattening filter illuminate a ∓11.4 x∓9.2 ≈420 cm²area. This same analysis may be applied to fields that are asymmetric in x or y (or both). Specifying the field size to be used in a given treatment determines the jaw positions, and these in turn—via FIG. 5 and Table II—determine the minimum and maximum x and y positions for sampling photon positions at the isocenter plane. In addition, using the information in Table II to select a photon's (x,y) position on the isocenter plane ensures, a priori, a high probability that the photon will reach the patient without intersecting a jaw.

[0057] Energy Distributions

[0058] Photon energy distributions of the various components, as tallied at the isocenter plane and on the central axis of the beam, are shown in FIGS. 6A-C. The photon energy distributions vary strongly with the treatment head component in which they are created. Photons from the target have energies ranging from the energy of the initial incident electrons down to low, but not quite zero, energy. This is consistent with the flattening filter's removal of the lower energy photons. The energy distribution from the primary collimator reflects both this filtering process (on the low energy side) and the fact that the photons are Compton scattered through a large angle (thus setting an upper bound on the high energy side). The photon distribution from the flattening filter reflects both the lack of low-energy filtering (since this is the last treatment machine head component transited by the photons) as well as the possibility for small angle scattering (and consequently little energy loss) of the high energy photons coming from the target.

[0059] The photon energy distributions vary with the distance from the central axis of the beam. This is shown in FIGS. 6D-6F which shows the distributions as tallied at the isocenter plane and at 20 cm from the central axis of the beam. The distributions for all three subsources show an increase in the proportion of low energy (here, ˜1 MeV) photons with increasing distance from the central axis, which correlates with the decrease in the thickness of the flattening filter traversed by these photons.

[0060] Photon Origin Distributions

[0061] The subsources, in addition to illuminating different size areas of the patient, also have markedly different source distributions. The target photon source is a flat disk, the primary collimator photons come from a ring-like source, and the flattening filter photons come from a broad and centrally-peaked source.

[0062] While these distributions are quite different, each may be accurately described by a radial distribution and an angular distribution, as discussed below.

[0063] Target Photons

[0064] Target photons come from a well-defined disk whose radius matches that of the incident electron beam. A cross section through the spot is shown in FIGS. 7A and 7B. The source distribution is relatively flat. There is also a low-intensity tail at larger R, due to those photons which scatter within the target before heading towards the patient. (The photons in the target which are scattered outside of the electron beam have low energy, and strictly speaking should not be combined with the photons created within the radius defined by the electron beam. Given the small radial extent of the tail, in practice we treat target photons as coming from one source.) The energy fluence plot shown in FIG. 7A is equivalently expressed as the radial-probability function shown in FIG. 7B. The azimuthal distribution of target photons is uniform. For this reason, a uniformly distributed (cos Φ, sin Φ) pair is generateed and then obtain:

[0065] (x_start, y_start)=R_start(cos Φ, sin Φ).  (1)

[0066] Scattered Photons

[0067] Scattered photons are handled in a fashion similar to the photons from the target, in that the source position is obtained from an (R_start, Φ_start) pair. The difference here is that Φ_startis also chosen from a distribution, instead of being chosen uniformly on the interval (0,2π).

[0068] The motivation for using an angular distribution is given in Table I and FIGS. 8A and 8B. At any point on the isocenter plane 20 not on the symmetry axis of the machine, there is a distinct difference in the scattering angle for events occurring at the “near” edge and the “far” edge of the scatterer, whether it be the flattening filter as shown in the figures or the primary collimator. Calculations indicate that incoherent scattering is the dominant production mechanism for photons in both the primary collimator and the flattening filter (Table I). Since the energy of the scattered photon varies inversely with the scattering angle, more photon energy will come from the “near” portion of the treatment head component.

[0069] This leads us to the coordinate system used to describe the starting position of a photon trajectory. As shown in FIG. 8C, a polar coordinate system may be defined at the position of the source plane of the particles. The angular origin of this system may be defined by the plane containing the beam axis and the photon's position of intersection at the isocenter plane. This yields a simple decomposition of the starting position in terms of relative azimuthal angle φ_rel, where φ_relranges from the value of 0° for the “near” edge to the value, given the reflection symmetry involved, of 180° for the “far” edge of the scatterer.

[0070] Primary Collimator

[0071] The source is fairly isotropic for small radii. This isotropy disappears off-axis, with the distribution in the last panel showing the effect of moving beyond the collimator angle: the “near” edge is blocked and only photons from the back surface of the collimator can reach the isocenter plane (FIGS. 9A-9F).

[0072] Flattening Filter

[0073] This same effect appears in the case of the flattening filter. Here, too, the source distributions favor the “near” edge of the flattening filter. The azimuthal anisotropy increases with increasing R_Iso(FIGS. 10A-10F).

[0074] Sampling Algorithm

[0075] Overview

[0076] Turning now to the development of a source algorithm for a Monte Carlo photon teletherapy code. The “source problem,” as described above, was to be able to generate photons in a manner that was both accurate and efficient. Two general patterns are seen to emerge during the preceding discussions. The first trend is that the energy spectra, radial energy fluence variations, and starting distributions differ strongly for different subsources. Second, for a given subsource, the energy distributions, the radial energy fluence, and the starting distributions all vary smoothly and relatively slowly as a function of the photon's radial position at the isocenter plane.

[0077] Separate the overall photon phase-space into individual phase-spaces, each designated by the treatment-head component in which the photons originate. Next subdivide the phase-space of each subsource into “tiles”: equal-width annular regions defined at the isocenter plane. It was shown above that given the maximum ranges of travel of the collimating jaws—the maximum distance from the central axis at which a photon has a significant probability of landing is ˜29 cm for target photons, ˜31 cm for primary collimator photons, and ˜41 cm for flattening filter photons. Currently ten tiles are used for each subsource. Each tile has its own photon energy spectrum; the radially slowly-varying energy spectrum of the subsource appears as a tile-dependent energy spectrum. In addition to its energy spectrum, each tile has its own R_startand Φ_reldistributions to describe photon starting positions.

[0078] Implementation

[0079] The correlated histograms used to describe the photon phase space are shown in FIGS. 11A-11F. Consider the problem for a treatment consisting of one field size for one machine; the generalization to multiple field sizes is straightforward. The use of the correlated histograms in calculating a treatment plan requires several preprocessing steps, in addition to the steps performed at runtime:

[0080] (1) Preprocessing steps:

[0081] (a) read the data file(s) for the machine(s) used in the calculation;

[0082] (b) convert each of the energy, angle, and radial probability histograms into the form required for alias sampling;

[0083] (c) use the jaw positions specified in the treatment plan to determine the limits of the (x_Iso, y_Iso) sampling areas for each subsource of each beam;

[0084] (d) for each subsource, use the (x_Iso, y_Iso) sampling areas to calculate the energy delivered to the isocenter plane; and

[0085] (e) use the relative energy contributions as sampling probabilities for each subsource.

[0086] (2) Generating photons during run-time:

[0087] (a) use random sampling to determine which type of photon (e.g., target, primary collimator, or flattening filter) will be generated. As indicated in Table I, these photon subsources will be selected with relative probabilities of approximately 93:3:4. The actual relative probabilities will depend on the jaw settings of the particular treatment plan.

[0088] (b) for this photon type, and given the jaw settings, look up the x and y limits of illumination (FIG. 4) on the isocenter plane. Generate a random, uniformly distributed (x_Iso, y_Iso) coordinate pair within this area.

[0089] (c) given this (x_Iso, y_Iso) pair, calculate R_Iso. Adjust the weight of the particle to account for the slowly-varying energy fluence (FIG. 3) of this photon type. Given R_Iso, determine the tile index. Calculate Φ_Iso.

[0090] (d) based on the tile index, sample the particle's energy from the energy distribution appropriate for this R_Iso(FIG. 6).

[0091] (e) sample a starting radius R_startfrom the appropriate distribution;

[0092] (f) sample a relative starting angle φ_rel. For target photons, use Eq. (1).

[0093] For other photons, sample the appropriate distribution.

[0094] The particle's energy and weight have thus been determined (steps c and d), as have two points defining the particle's trajectory (steps b, e, and f). The trajectory-defining points define the particle's direction cosines, and the required phase-space information needed to start tracking this particle is now complete. The photon is backtracked to the bottom of the monitor chamber, and then transported through beam modifiers (if any) to the patient. A PEREGRINE-specified plane is a plane where PEREGRINE starts to track a particle through air towards a patient, wherein the step of calculating the phase-space information for a specific radiation accelerator head from the correlated histograms comprises using the values of X_start, Y_start, Z_start, X_Iso, y_Isoand the known z_Isofor a specific accelerator to determine a particle trajectory's intersection with the PEREGRINE-specified plane.

[0095] Accuracy

[0096] The next question is that of the accuracy of the algorithm described above. Accuracy as defined here refers to the degree to which our algorithm reproduces the underlying photon phase space. This is a two-part question. The first question deals with the energy distribution(s) of the photons at the isocenter plane. The second question deals with the trajectories of the photons.

[0097] Energy Distribution at the Isocenter Plane

[0098] These distributions are accurately reproduced because (a) the binning of the energy distributions is actually performed at this plane, and (b) because the distributions vary slowly enough with distance from the axis (FIG. 12) that a piecewise-linear distribution can track the change (FIG. 13).

[0099] This variation with radial distance is slow. FIGS. 12A-12C show curves of bin probability for several energies are plotted as functions of radius at the isocenter plane. These curves are all slowly varying and are all accurately represented by piecewise-linear functions. In addition, only a small number of pieces—on the order of 10—is sufficient.

[0100] Because the knots of the piecewise-linear functions shown in FIG. 12 are actually the bin-centers of histograms obtained by binning the photon energies over this range (FIG. 13), the radial variation in the energy distribution over the area of interest may be reproduced by creating energy histograms in a small number of annular regions (the “tiles”) at the isocenter plane. A photon energy distribution of the form shown in FIG. 6 is then created within each tile. The maximum tile radius is determined by the maximum radius of interest for the particular photon distribution. For the photons coming from the bremsstrahlung target, the maximum radius is ˜29 cm. Similarly, as noted above, for the primary collimator and flattening filter, maximum tile radii of 33.6 and 40 cm are required.

[0101] Trajectory Information

[0102] The second question deals with the trajectories of the photons. Our source algorithm determines the trajectory of the photons by “connecting the dots,” one of which is at the isocenter plane and the other of which is either at or within a few centimeters of the position of the bremsstrahlung target. The radial extent of the source distributions is on the order of 1 mm for the target photons, and on the order of 1 cm to 3 cm for photons scattered from the primary collimator or the flattening filter. The maximum expected trajectory error, therefore is expected to be no larger than δθ≈arctan 3/100≈2 deg. The actual error will be smaller because of the use of the starting R and Φ distributions.

[0103] The accuracy of the trajectories may be checked by comparing the location and width of the jaw-edge penumbra for each source as a function of jaw position with the location and width obtained with the underlying phase-space file (a “knife-edge” test). While the position of the photon on the isocenter plane is selected at random, an “error” in the photon's direction cosines will cause an error in its assumed trajectory. This trajectory error would affect, for example, whether or not a photon near the edge of a treatment field would be blocked by a jaw edge. This type of test is shown in FIG. 14 for flattening filter photons, with the upper jaw moved to its maximum distance from the centerline of the beam. The roll-off in energy fluence is well-matched by the source algorithm presented in this paper, indicating that the use of correlated R and φ distributions accurately reproduces the photons' trajectories. FIG. 14 also shows the fluence error that results when φ is assumed to be uniformly distributed.

[0104] Computational Efficiency

[0105] This algorithm's computational efficiency has two bases. The first of these is the conceptual structure, in which we first pick only those photon trajectories that enter that portion of the isocenter plane which is illuminated for the specified jaw settings. Only then do we generate a starting position to determine “where they came from.” The second basis is the efficiency of the “alias method” (see below) to sample the various probability distributions used to describe the photon phase space. This method requires two random numbers to select a bin from a distribution of arbitrary length and complexity. A third random number is required if uniform sampling within the bin is desired.

[0106] The computational burden of the algorithm is detailed in Table III. Note that the source algorithm determines the six needed pieces of the particle phase space information (x, y, u, v, E, q, weight) using approximately 14 random numbers and 3 or 4 table lookups. (Recall that, given the context of the problem, z is the isocenter plane location and w={square root over (1−u²−v²)}. Step 1, which uses alias sampling to determine which treatment head component from which machine will be selected, makes the running time of this algorithm independent of the number of beams used in the treatment plan.

[0107] Steps 2 and 6 have the most impact on the algorithm's efficiency. A number of the photon trajectories calculated via steps 2, 4, and 5 will intersect the jaw corners and subsequently not reach the patient. This behavior is intended, since this is how the penumbrae are determined. This effect is most pronounced for small field sizes, and becomes more visible as the radial extent of the treatment head component increases. This is shown in FIG. 15 for a number of square field sizes. For most of the fields, trajectories for “target” photons have a better than 90% chance of clearing the jaws. Photons from the primary collimator and flattening filter have a lower survival probability. Because most of the photons originate in the target, however, the overall efficiency is better than 80% for most fields in the range studied.

[0108] Reduction of Phase-Space Files

[0109] It was stated above that an automated procedure is useful for converting the gigabyte-sized data files obtained from Monte Carlo simulations of an accelerator head into the correlated histograms described herein. The conversion requires five passes through the 1-GB phase space file and requires approximately 23 minutes to complete. (The computer used is a Digital AlphaServer 8400 Model 5/440 running OSF1 V4.0.)

[0110] The five passes perform the following operations:

[0111] (1) This step calculates the curves shown in FIG. 5. This step also calculates the data given in Table II. In so doing, the maximum tile radius needed for each subsource, and hence the radial limits of each tile of each subsource are also calculated.

[0112] (a) For each particle in the phase-space file:

[0113] (i) read each particle and identify the hardware component in which it last scattered.

[0114] (ii) based on the particle's position at the tally plane (x_tally, y_tally, z_tally) and on the particle's direction cosines (u, v, w), determine where the particle will intersect the isocenter plane. At the same time, determine the jaw settings of Table II for which the particle would miss the jaws. Add this particle to the curves shown in FIG. 5.

[0115] (b) At this point, the curves shown in FIG. 5 have been calculated. Each curve is then searched to determine the position at which the height of the curve drops below (the nominal value of) 0.2% of the maximum energy fluence on the isocenter plane.

[0116] (c) The location of this crossing point is entered into Table II.

[0117] (d) R_maxis calculated for each subsource, based on the maximum allowed jaw positions, thus determining the radial extent for each tile of each subsource.

[0118] (2) This step makes an initial estimate of the position of the “Source Plane” (shown in FIG. 8) for each tile of each subsource. Recall that the source particle algorithm calculates the trajectory of a particle by first choosing a position at the isocenter plane and then asking “where did the particle come from?” Trajectory errors arise because of errors in determining the starting position of a particle. Minimizing the extent of the radial distributions should in turn minimize errors in the trajectories calculated during run-time. To this end, determine the value of z at which each particle comes closest to the central axis of the beam.

[0119] (a) For each particle:

[0120] (i) determine its subsource and its tile index;

[0121] (ii) using the particle's position and direction as it intersects the tally plane, trace its path back to where it was closest to the central axis of the beam: $Z_{\mathrm{closest}}=Z_{\mathrm{tally}}-\frac{w\left({\mathrm{ux}}_{\mathrm{tally}}+{\mathrm{uy}}_{\mathrm{tally}}\right)}{1-w^2}$

[0122] (i) based on this value of z_closest, add the particle to the appropriate bin of the N(z) histogram for this tile of this subsource.

[0123] (b) For each tile of each subsource, set Z_guessequal to the value of z at which the N(z) histograms has its maximum value.

[0124] (3) This step determines the high-end cutoff value for R_startfor each tile of each subsource.

[0125] (a) For each particle:

[0126] (i) determine its subsource and its tile index;

[0127] (ii) determine its distance R from the central axis of the beam, at the position z=z_guessdetermined above;

[0128] (iii) add the particle to the appropriate bin of the R_starthistogram for the appropriate tile and subsource.

[0129] (b) For each tile of each subsource, find the minimum value of R that exceeds Rstart for 99% of the particles.

[0130] (4) Fine-tune the z_startfor each tile of each subsource. The z_guessdetermined by Step 2 above yields a good, but not optimal, value of z at which to create R_startand φ_reldistributions. An azimuthally anisotropic starting distribution, shown e.g., for the primary collimator in FIG. 16, appears “off-center” when the photons belonging to this distribution are traced back to Z_guess. Such a pattern may, however, be searched for that value of ΔR which maximizes the accuracy of a (R, Φ) representation.

[0131] (a) For each particle:

[0132] (i) determine its subsource, tile index, and z_guess;

[0133] (ii) calculate its (x,y) position at z_guess; and

[0134] (iii) add the particle to the appropriate bin of a two-dimensional histogram for the appropriate tile and subsource.

[0135] (b) For each two-dimensional histogram so calculated, step along the x axis and at each point:

[0136] (i) assuming this point to be the origin, create an (R, Φ) representation of the (x,y) distribution calculated above;

[0137] (ii) find that value of x, which when used as the origin, minimizes the difference between the (x,y) distribution above and the one calculated from the (R, Φ) representation based on that origin.

[0138] (c) Determine the z_startat which no offset-would have been required.

[0139] (5) This step creates the correlated histograms needed to describe the accelerator output. For each particle:

[0140] (a) Determine its subsource, tile index, and Φ_Iso;

[0141] (b) Determine R_startand φ_relat the z_bestfor this tile and update the appropriate histograms;

[0142] (c) Update the energy distribution of this tile.

[0143] The alias sampling method is an efficient method for sampling probability distributions comprising a large number of bins, which in turn may have widely-varying heights. Also, while alias sampling is a very good way to sample distributions, it is not the only method available (other methods are more-easily coded, but tend to be slower). The distributions shown in this typically have both these properties. The alias method has been found to be faster than other methods of sampling distributions, e.g., searching cumulative distributions, while avoiding loss-of-accuracy problems associated with using equal probability distributions.

[0144] The alias-sampling method is a two-stage process. The first stage is performed as a preprocessing step, and is the remapping of the original distribution, shown in FIG. 17A, into the tableau form shown in FIG. 17B. The second stage occurs during the runtime of the Monte Carlo calculation. Two random numbers ξ1 and ξ2 are generated. ξ1 is used to select one of the columns of the tableau. ξ2 is compared to the cutoff-level of this column, and depending on the result, either the lower or upper bin index is returned. Thus, for the example shown, the pair of numbers (ξ1, ξ2)=(0.5, 0.9) indicates that “Bin 1” of FIG. 17ais selected. Similarly, the pair of numbers (ξ1, ξ2)=(0.9, 0.07) indicates that “Bin 5” of FIG. 17A should be selected.

[0145] As seen in FIG. 17B, the tableau-form of the probability distribution requires three arrays: one array for the bin edges, one array to hold the cutoff-levels, and one array to hold the “upper-bin” values. The latter two arrays are obtained during the preprocessing stage of the alias method. The preprocessing stage is an iterative process, in which “height” is transferred from “high” bins to “low” ones. The steps are described below and are summarized in Table IV.

[0146] 1. find the highest remaining bin, and the first bin whose height is less than the average. For the distribution shown in FIG. 17A, these are bins 1 and 3, respectively, for the first iteration.

[0147] 2. add enough “height” from the “high” bin to bring the “low” bin up to the average. The anti-alias index for what was the “low” bin has now been defined.

[0148] 3. reduce the “height of the “high” bin by the amount used in step 2.

[0149] 4. go to Step 1, and repeat until all bins have either been “raised” or “lowered” to the average “height”, as the case may be.

[0150] While the example used in this section was a histogram distribution with uniform bin widths, the alias sampling method has been extended to piecewise-linear functions with non-uniformly spaced abscissae.

[0151] The foregoing description of the invention has been presented for purposes of illustration and description and is not intended to be exhaustive or to limit the invention to the precise form disclosed. Many modifications and variations are possible in light of the above teaching. The embodiments were chosen and described to best explain the principles of the invention and its practical application to thereby enable others skilled in the art to best use the invention in various embodiments and with various modifications suited to the particular use contemplated. The scope of the invention is to be defined by the following claims.