Command based Scoring statistics the dose deposit in human body

[Hans]

Hello, esteemed seniors and users of Geant4. I have currently completed the modeling of a Varian Vital Beam medical accelerator using Geant4, including the modeling of a specific patient. The beam uses a 6MV phase space, and I have simulated the irradiation of over 1,000 subfields for a treatment plan. The scoring method used is Command-based Scoring to tally energy deposition (the energy deposition matrix matches the resolution of the CT). Subsequently, I divided the energy deposition matrix by the mass matrix (converted from CT values to material densities using four fitting formulas) to obtain the dose matrix. I then compared this with the dose calculation results from the AAA algorithm in the TPS (Treatment Planning System). Note that I can only compare relative values, not absolute values, due to computational limitations preventing the accumulation of 1 Gy dose in Geant4. Surprisingly, I found abnormally high dose values in the low-density lung region, which is unreasonable. This is particularly noteworthy because, compared to the AAA algorithm—which is known to overestimate dose in low-density lung regions—my Monte Carlo results are even higher. I am uncertain about the cause.

I am aware that in Command-based Scoring, the physical unit can be specified as doseDepositto tally dose directly. However, I have a question: when a scoring mesh grid contains material heterogeneity (i.e., multiple materials with different densities occupy the same scoring volume), could there be issues with the Geant4 kernel’s dose calculation? In principle, since dose = energy / mass, I would like to seek discussion and collaboration with experts, experienced users, or anyone interested in medical physics on the forum.

[Hans]

this is physic model (MLC+JAW+human body+source)

[Hans]

[Hans]

[ctwhite]

@Hans

could there be issues with the Geant4 kernel’s dose calculation?

Yes there probably will be issues. The command base scorer calculates the dose by getting the material’s density at the energy deposition step, if the dose grid contains more than one material this will result in inaccurate calculation of the dose. You need to make sure that the scoring grid is aligned with the CT voxels. If you want to score the dose at a coarser resolution than the CT scan, which is typical for TPS, its best to downsample it in a pre-processing step.

The overdose to the lung could be caused by multiple things, and there isn’t enough information from the figures to say with certainty, but it seems like your scoring grid might not be aligned correctly (in particular in the Z direction). Also double check your material definitions and ensure that it is consistent with the TPS. Possibly you might have used the wrong density, either in the simulation or in the post-processing. Actually saying that, it seems like there is a lot of dose at the skin too, including a strange wedge-like shape on the back, in the second plot you shared. This indicates an alignment issue or material definition problem. How are you loading the patient model into Geant4?

Simulating the full treatment, in terms of Gy delivered, isn’t feasibly in Geant4. You need to normalise the dose from the TPS and the Geant4 simulation. This looks like an IMRT plan, so make sure that you take the meterset weight into consideration for each control point, as they will likely be different.

[Hans]

Hello, thank you very much for your suggestions and guidance. I would like to respond with the following points:

1. I have recently rerun the simulations. In the Command-based Scoring, I explicitly specified the dose to be scored in units of Gy, and I ensured that the mesh segmentation is fully aligned with the entire CT matrix. That is, I guaranteed that each mesh voxel contains only a single CT voxel, so the density within each mesh element is homogeneous.

2. For patient modeling in Geant4, I used a combination of the G4VNestedParameterisation class and G4Replica, which allows me to efficiently build a full 1:1 patient geometry with a resolution of 512 × 512 × 50 CT voxels while keeping memory usage under control.

   Regarding material assignment, I followed the method described in Schneider, Bortfeld, Schlegel, Phys. Med. Biol. 45:459–478, 2000, which provides the conversion from HU values to mass density and elemental composition. For example, lung tissue with low density is typically in the range of 0.15–0.25 g/cm³, while normal soft tissue is close to 1 g/cm³.

   In principle, mapping all 512 × 512 × 50 voxels individually according to their HU values would require defining about 3001 materials (HU range from −1000 to 2000), which is computationally impractical and leads to severe performance issues. Therefore, following the referenced paper, we grouped the HU values into 24 intervals, resulting in only 24 materials, each defined by an average density and a corresponding elemental composition mixture (details can be found in the cited publication).

3. Regarding the skin dose you observed, this involves several steps, including calibration, resampling, spatial alignment, and masking:

   • Calibration: I need to determine how many phase-space photons correspond to 1 MU. This was done by measuring the maximum dose point in water and back-calculating proportionally. I obtained approximately 1 MU = 1 cGy ≈ 8.7 × 10¹¹ photons. Since it is not feasible to simulate such a large number of particles directly, I currently assume 1 MU corresponds to 2 × 10⁷ photons for practical simulations.

   • Resampling: The original 3D energy or dose matrix from Monte Carlo scoring has the same resolution as the CT voxels, while the TPS dose matrix has a resolution of 2.5 × 2.5 × 2.5 mm³ and covers only part of the spatial extent. Therefore, I resample the high-resolution Monte Carlo dose matrix to the TPS dose grid using the Python SimpleITK library. This is why, before resampling, the Monte Carlo dose distribution (on the right) appears with higher resolution and clearer details.

   • Masking: Since I am only interested in the dose inside the patient body, I read the contoured Body structure and apply its index set to each dose slice: voxels inside the mask are kept, while those outside are set to zero. You made a very careful observation—the initially high “skin dose” was actually caused by an error in the left–right orientation when applying the Body mask for the first time. This issue has now been corrected.

Thank you again for your detailed feedback. I look forward to further discussion and collaboration with experts who are interested in this topic. I have a background in nuclear physics and software development.

[Hans]

Regarding control point simulation, I have completed accurate mechanical modeling and motion design of the entire treatment head system (JAW + MLC). This allows me to generate arbitrary combinations of gantry angle, collimator angle, X1/X2/Y1/Y2 jaw projections, as well as opposing motion of the 120-leaf MLC, including rotation of the phase-space source plane (you may refer to the model diagram). All mechanical parameters are taken from the official Varian Monte Carlo data support package.

I read the patient’s treatment plan file, which specifies the mechanical configuration and delivered monitor units for each control point (for example, 1.561 MU). I then determine the required number of input photons for the corresponding subfield simulation as 2 × 10⁷ × 1.561. A single IMRT plan typically contains more than 1000 control points, and I simulate all of them. With the patient geometry kept fixed, the scorer also remains stationary, and the three-dimensional dose data are accumulated over all control points.

The total computational cost for this process was approximately 12 days using 20 threads.

[ctwhite]

Double check your material definitions. If one of those 24 materials is incorrectly a lung volume but with a lower than expected density (i.e 0.02 g/cm3 instead of 0.20 g/cm3) then you can have higher than expected dose in those materials.

[Hans]

ok i will do that check in the day after tomorrow

[Hans]

HU_BIN1 (-1000,-950) density = 0.026979

HU_BIN2 (-950,-900) density = 0.078517

HU_BIN3 (-900,-850) density = 0.130055

HU_BIN4 (-850,-800) density = 0.181593

HU_BIN5 (-800,-750) density = 0.233131

HU_BIN6 (-750,-700) density = 0.284669

HU_BIN7 (-700,-650) density = 0.336207

HU_BIN8 (-650,-600) density = 0.387745

HU_BIN9 (-600,-550) density = 0.439283

HU_BIN10 (-550,-500) density = 0.490821

HU_BIN11 (-500,-450) density = 0.542359

HU_BIN12 (-450,-400) density = 0.593897

HU_BIN13 (-400,-350) density = 0.645435

HU_BIN14 (-350,-300) density = 0.696973

HU_BIN15 (-300,-250) density = 0.748511

HU_BIN16 (-250,-200) density = 0.800049

HU_BIN17 (-200,-150) density = 0.851587

HU_BIN18 (-150,-120) density = 0.887663

HU_BIN19 (-120,-83) density = 0.924119

HU_BIN20 (-82,-53) density = 0.953717

HU_BIN21 (-52,-23) density = 0.980302

HU_BIN22 (-22,7) density = 1.009900

HU_BIN23 (8,18) density = 1.029165

HU_BIN24 (19,80) density = 1.079065

HU_BIN25 (80,120) density = 1.119900

HU_BIN26 (120,200) density = 1.111720

HU_BIN27 (200,300) density = 1.165000

HU_BIN28 (300,400) density = 1.224200

HU_BIN29 (400,500) density = 1.283400

HU_BIN30 (500,600) density = 1.342600

HU_BIN31 (600,700) density = 1.401800

HU_BIN32 (700,800) density = 1.461000

HU_BIN33 (800,900) density = 1.520200

HU_BIN34 (900,1000) density = 1.579400

HU_BIN35 (1000,1100) density = 1.638600

HU_BIN36 (1100,1200) density = 1.697800

HU_BIN37 (1200,1300) density = 1.757000

HU_BIN38 (1300,1400) density = 1.816200

HU_BIN39 (1400,1500) density = 1.875400

HU_BIN40 (1500,1600) density = 1.934600 //i finish my check and This is my final result for refining the density values in the Schneider −950 to −120 HU range. For each density bin, the midpoint of the corresponding HU interval is used to calculate the density.

[Hans]

[Hans]

[Hans]

G4Material* couch_surface_mat=new G4Material(“couch_surface”,1.4 * g/ cm3 ,2);
couch_surface_mat->AddElement(elC,0.92);
couch_surface_mat->AddElement(elO,0.08);

G4Material* couch_interior_mat=new G4Material(“couch_interior”,0.06 * g/ cm3 ,2);
couch_interior_mat->AddElement(elC,0.9);
couch_interior_mat->AddElement(elH,0.1);

G4Material* body_board=new G4Material(“body_board”,0.052 * g/ cm3 ,4);
body_board->AddElement(elC,0.50);
body_board->AddElement(elH,0.07);
body_board->AddElement(elN,0.11);
body_board->AddElement(elO,0.32); //In addition, I have completed the modeling, replacement, and material definition for the treatment couch surface, its internal structure, and the patient support board. The material properties were determined based on references from the literature and should be reliable.

[Hans]

In the first simulation, I only defined 24 materials, because the range from −950 to −120 HU was grouped into a single interval with an average density of about 0.5 g/cm³. This density was too coarse and biased toward a higher value.

Now, I have subdivided the −950 to −120 HU range into 50-HU bins, resulting in a total of 40 materials. I have rerun the simulation and expect the results to improve. At the same time, I configured the scoring to record both energy deposition and dose quantities.

[ctwhite]

One problem is that you have multiple bins with different densities in the -950 to -120 range, but likely your elemental composition is the same for all of them (I have to assume because the rest of the bins are assigned to that Schneider table). A -850 HU lung voxel does not contain the same amount of carbon that a -300 HU lung voxel has. From the Schneider paper:

”For the remaining range between air and soft tissues, we assign the composition of lung tissue, for which we have calculated −741 H.”

There might also be some artifacting in the CT scan that could effect the HU, which your simulation would be very sensitive to because you have 50 HU bins. Try plotting the assigned voxel densities

Most of the work simulating patients, either in G4 or EGS, that I’m familiar with seem to have a single “lung” material. I know of one TPS that defines materials as a mixture of the ICRP defined materials for its MC code, and so lung ends up being a mixture of air and soft tissue, but that importantly fixes the problem of elemental compositions. I recommend simplifying the lung model here and making it more homogeneous with fewer material definitions, or making sure that you properly calculate the elemental compositions for each density.

EDIT:
Ah, I apologise, I missed that you said that you used a single lung definition with a density of 0.5 g/cm3 initially. That value is likely too large, especially if you used the chemical composition at -741 HU. The ICRP110 example phantom which has a density of 0.382 g/cm3, for example.

[Hans]

This represents the initial approach in Geant4, where patient materials are defined rather than roughly approximated as only 10 materials, and density variations can be controlled using densitydiff. In any case, the current set of 40 materials is a reasonably achievable solution for me at this stage, and it feels more realistic than approaches based purely on ICRU or NIST materials.

My goal is first to obtain a Monte Carlo dose distribution for an IMRT plan. With sufficient reference data available, I can then further refine the model, since my research is not intended for a commercial product. The current simulation will likely take about 10 more days to complete.G4Element* elC = new G4Element(name = “Carbon”, symbol = “C”, z = 6.0, a = 12.011 * g / mole);
G4Element* elH = new G4Element(name = “Hydrogen”, symbol = “H”, z = 1.0, a = 1.008 * g / mole);
G4Element* elN = new G4Element(name = “Nitrogen”, symbol = “N”, z = 7.0, a = 14.007 * g / mole);
G4Element* elO = new G4Element(name = “Oxygen”, symbol = “O”, z = 8.0, a = 16.00 * g / mole);
G4Element* elNa =
new G4Element(name = “Sodium”, symbol = “Na”, z = 11.0, a = 22.98977 * g / mole);
G4Element* elMg =
new G4Element(name = “Magnesium”, symbol = “Mg”, z = 12.0, a = 24.3050 * g / mole);
G4Element* elP =
new G4Element(name = “Phosphorus”, symbol = “P”, z = 15.0, a = 30.973976 * g / mole);
G4Element* elS = new G4Element(name = “Sulfur”, symbol = “S”, z = 16.0, a = 32.065 * g / mole);
G4Element* elCl =
new G4Element(name = “Chlorine”, symbol = “Cl”, z = 17.0, a = 35.453 * g / mole);
G4Element* elK =
new G4Element(name = “Potassium”, symbol = “K”, z = 19.0, a = 30.0983 * g / mole);

G4Element* elFe = new G4Element(name = “Iron”, symbol = “Fe”, z = 26, a = 56.845 * g / mole);

G4Element* elCa = new G4Element(name = “Calcium”, symbol = “Ca”, z = 20.0, a = 40.078 * g / mole);

G4Element* elZn = new G4Element(name = “Zinc”, symbol = “Zn”, z = 30.0, a = 65.382 * g / mole);
// 需要先定义 12 个元素，包括 Ar（例子代码里没建 Ar，要补上）
G4Element* elAr = new G4Element(name=“Argon”,symbol=“Ar”, z=18.0,a=39.948 * g/mole);

// Creating Materials :
G4int numberofElements;

// Air
fAir = new G4Material(“Air”, 1.290 * mg / cm3, numberofElements = 2);
fAir->AddElement(elN, 0.7);
fAir->AddElement(elO, 0.3);

// Soft tissue (ICRP - NIST)
G4Material* softTissue = new G4Material(“SoftTissue”, 1.00 * g / cm3, numberofElements = 13);
softTissue->AddElement(elH, 10.4472 * perCent);
softTissue->AddElement(elC, 23.219 * perCent);
softTissue->AddElement(elN, 2.488 * perCent);
softTissue->AddElement(elO, 63.0238 * perCent);
softTissue->AddElement(elNa, 0.113 * perCent);
softTissue->AddElement(elMg, 0.0113 * perCent);
softTissue->AddElement(elP, 0.113 * perCent);
softTissue->AddElement(elS, 0.199 * perCent);
softTissue->AddElement(elCl, 0.134 * perCent);
softTissue->AddElement(elK, 0.199 * perCent);
softTissue->AddElement(elCa, 0.023 * perCent);
softTissue->AddElement(elFe, 0.005 * perCent);
softTissue->AddElement(elZn, 0.003 * perCent);

// Lung Inhale
G4Material* lunginhale =
new G4Material(“LungInhale”, density = 0.217 * g / cm3, numberofElements = 9);
lunginhale->AddElement(elH, 0.103);
lunginhale->AddElement(elC, 0.105);
lunginhale->AddElement(elN, 0.031);
lunginhale->AddElement(elO, 0.749);
lunginhale->AddElement(elNa, 0.002);
lunginhale->AddElement(elP, 0.002);
lunginhale->AddElement(elS, 0.003);
lunginhale->AddElement(elCl, 0.002);
lunginhale->AddElement(elK, 0.003);

// Lung exhale
G4Material* lungexhale =
new G4Material(“LungExhale”, density = 0.508 * g / cm3, numberofElements = 9);
lungexhale->AddElement(elH, 0.103);
lungexhale->AddElement(elC, 0.105);
lungexhale->AddElement(elN, 0.031);
lungexhale->AddElement(elO, 0.749);
lungexhale->AddElement(elNa, 0.002);
lungexhale->AddElement(elP, 0.002);
lungexhale->AddElement(elS, 0.003);
lungexhale->AddElement(elCl, 0.002);
lungexhale->AddElement(elK, 0.003);

// Adipose tissue
G4Material* adiposeTissue =
new G4Material(“AdiposeTissue”, density = 0.967 * g / cm3, numberofElements = 7);
adiposeTissue->AddElement(elH, 0.114);
adiposeTissue->AddElement(elC, 0.598);
adiposeTissue->AddElement(elN, 0.007);
adiposeTissue->AddElement(elO, 0.278);
adiposeTissue->AddElement(elNa, 0.001);
adiposeTissue->AddElement(elS, 0.001);
adiposeTissue->AddElement(elCl, 0.001);

// Brain (ICRP - NIST)
G4Material* brainTissue = new G4Material(“BrainTissue”, 1.03 * g / cm3, numberofElements = 13);
brainTissue->AddElement(elH, 11.0667 * perCent);
brainTissue->AddElement(elC, 12.542 * perCent);
brainTissue->AddElement(elN, 1.328 * perCent);
brainTissue->AddElement(elO, 73.7723 * perCent);
brainTissue->AddElement(elNa, 0.1840 * perCent);
brainTissue->AddElement(elMg, 0.015 * perCent);
brainTissue->AddElement(elP, 0.356 * perCent);
brainTissue->AddElement(elS, 0.177 * perCent);
brainTissue->AddElement(elCl, 0.236 * perCent);
brainTissue->AddElement(elK, 0.31 * perCent);
brainTissue->AddElement(elCa, 0.009 * perCent);
brainTissue->AddElement(elFe, 0.005 * perCent);
brainTissue->AddElement(elZn, 0.001 * perCent);

// Breast
G4Material* breast = new G4Material(“Breast”, density = 0.990 * g / cm3, numberofElements = 8);
breast->AddElement(elH, 0.109);
breast->AddElement(elC, 0.506);
breast->AddElement(elN, 0.023);
breast->AddElement(elO, 0.358);
breast->AddElement(elNa, 0.001);
breast->AddElement(elP, 0.001);
breast->AddElement(elS, 0.001);
breast->AddElement(elCl, 0.001);

// Spinal Disc
G4Material* spinalDisc = new G4Material(“SpinalDisc”, 1.10 * g / cm3, numberofElements = 8);
spinalDisc->AddElement(elH, 9.60 * perCent);
spinalDisc->AddElement(elC, 9.90 * perCent);
spinalDisc->AddElement(elN, 2.20 * perCent);
spinalDisc->AddElement(elO, 74.40 * perCent);
spinalDisc->AddElement(elNa, 0.50 * perCent);
spinalDisc->AddElement(elP, 2.20 * perCent);
spinalDisc->AddElement(elS, 0.90 * perCent);
spinalDisc->AddElement(elCl, 0.30 * perCent);

// Water
G4Material* water = new G4Material(“Water”, density = 1.0 * g / cm3, numberofElements = 2);
water->AddElement(elH, 0.112);
water->AddElement(elO, 0.888);

// Muscle
G4Material* muscle = new G4Material(“Muscle”, density = 1.061 * g / cm3, numberofElements = 9);
muscle->AddElement(elH, 0.102);
muscle->AddElement(elC, 0.143);
muscle->AddElement(elN, 0.034);
muscle->AddElement(elO, 0.710);
muscle->AddElement(elNa, 0.001);
muscle->AddElement(elP, 0.002);
muscle->AddElement(elS, 0.003);
muscle->AddElement(elCl, 0.001);
muscle->AddElement(elK, 0.004);

// Liver
G4Material* liver = new G4Material(“Liver”, density = 1.071 * g / cm3, numberofElements = 9);
liver->AddElement(elH, 0.102);
liver->AddElement(elC, 0.139);
liver->AddElement(elN, 0.030);
liver->AddElement(elO, 0.716);
liver->AddElement(elNa, 0.002);
liver->AddElement(elP, 0.003);
liver->AddElement(elS, 0.003);
liver->AddElement(elCl, 0.002);
liver->AddElement(elK, 0.003);

// Tooth Dentin
G4Material* toothDentin = new G4Material(“ToothDentin”, 2.14 * g / cm3, numberofElements = 10);
toothDentin->AddElement(elH, 2.67 * perCent);
toothDentin->AddElement(elC, 12.77 * perCent);
toothDentin->AddElement(elN, 4.27 * perCent);
toothDentin->AddElement(elO, 40.40 * perCent);
toothDentin->AddElement(elNa, 0.65 * perCent);
toothDentin->AddElement(elMg, 0.59 * perCent);
toothDentin->AddElement(elP, 11.86 * perCent);
toothDentin->AddElement(elCl, 0.04 * perCent);
toothDentin->AddElement(elCa, 26.74 * perCent);
toothDentin->AddElement(elZn, 0.01 * perCent);

// Trabecular Bone
G4Material* trabecularBone =
new G4Material(“TrabecularBone”, density = 1.159 * g / cm3, numberofElements = 12);
trabecularBone->AddElement(elH, 0.085);
trabecularBone->AddElement(elC, 0.404);
trabecularBone->AddElement(elN, 0.058);
trabecularBone->AddElement(elO, 0.367);
trabecularBone->AddElement(elNa, 0.001);
trabecularBone->AddElement(elMg, 0.001);
trabecularBone->AddElement(elP, 0.034);
trabecularBone->AddElement(elS, 0.002);
trabecularBone->AddElement(elCl, 0.002);
trabecularBone->AddElement(elK, 0.001);
trabecularBone->AddElement(elCa, 0.044);
trabecularBone->AddElement(elFe, 0.001);

// Trabecular bone used in the DICOM Head

G4Material* trabecularBone_head =
new G4Material(“TrabecularBone_HEAD”, 1.18 * g / cm3, numberofElements = 12);
trabecularBone_head->AddElement(elH, 8.50 * perCent);
trabecularBone_head->AddElement(elC, 40.40 * perCent);
trabecularBone_head->AddElement(elN, 2.80 * perCent);
trabecularBone_head->AddElement(elO, 36.70 * perCent);
trabecularBone_head->AddElement(elNa, 0.10 * perCent);
trabecularBone_head->AddElement(elMg, 0.10 * perCent);
trabecularBone_head->AddElement(elP, 3.40 * perCent);
trabecularBone_head->AddElement(elS, 0.20 * perCent);
trabecularBone_head->AddElement(elCl, 0.20 * perCent);
trabecularBone_head->AddElement(elK, 0.10 * perCent);
trabecularBone_head->AddElement(elCa, 7.40 * perCent);
trabecularBone_head->AddElement(elFe, 0.10 * perCent);

// Dense Bone
G4Material* denseBone =
new G4Material(“DenseBone”, density = 1.575 * g / cm3, numberofElements = 11);
denseBone->AddElement(elH, 0.056);
denseBone->AddElement(elC, 0.235);
denseBone->AddElement(elN, 0.050);
denseBone->AddElement(elO, 0.434);
denseBone->AddElement(elNa, 0.001);
denseBone->AddElement(elMg, 0.001);
denseBone->AddElement(elP, 0.072);
denseBone->AddElement(elS, 0.003);
denseBone->AddElement(elCl, 0.001);
denseBone->AddElement(elK, 0.001);
denseBone->AddElement(elCa, 0.146);

// Cortical Bone (ICRP - NIST)
G4Material* corticalBone = new G4Material(“CorticalBone”, 1.85 * g / cm3, numberofElements = 9);
corticalBone->AddElement(elH, 4.7234 * perCent);
corticalBone->AddElement(elC, 14.4330 * perCent);
corticalBone->AddElement(elN, 4.199 * perCent);
corticalBone->AddElement(elO, 44.6096 * perCent);
corticalBone->AddElement(elMg, 0.22 * perCent);
corticalBone->AddElement(elP, 10.497 * perCent);
corticalBone->AddElement(elS, 0.315 * perCent);
corticalBone->AddElement(elCa, 20.993 * perCent);
corticalBone->AddElement(elZn, 0.01 * perCent);

// Tooth enamel
G4Material* toothEnamel = new G4Material(“ToothEnamel”, 2.89 * g / cm3, numberofElements = 10);
toothEnamel->AddElement(elH, 0.95 * perCent);
toothEnamel->AddElement(elC, 1.11 * perCent);
toothEnamel->AddElement(elN, 0.23 * perCent);
toothEnamel->AddElement(elO, 41.66 * perCent);
toothEnamel->AddElement(elNa, 0.79 * perCent);
toothEnamel->AddElement(elMg, 0.23 * perCent);
toothEnamel->AddElement(elP, 18.71 * perCent);
toothEnamel->AddElement(elCl, 0.34 * perCent);
toothEnamel->AddElement(elCa, 35.97 * perCent);
toothEnamel->AddElement(elZn, 0.02 * perCent);

[Hans]

I would also like to ask you how I should consider and calculate the uncertainty of my dose results, since 1 MU requires 10^11 phase-space photons, but I only simulated 10^7 due to hardware limitations. Currently, my resulting dose matrix can only be compared to the TPS after being normalized simultaneously, in terms of relative dose.