Hi,
Stopping powers
I am looking at the material database in openTPS/MCsquare and I can see that there are two stopping power tables (PSTAR and Geant4). Which one is the default table that MCsquare is using? If it is not Geant4, how can I force MCsquare to use the Geant4 tables instead? I am having a bit of a problem in comparing dose distributions in a heterogeneous phantom (i.e., lung and bone slab inserted in water) similar to Fig 2. of the MCsquare paper. While I do get the correct range for protons passing through the lung slab, I see a difference of 1.3 mm in the range of protons going through the bone as compared to TOPAS. My initial guess is this is a stopping power issue.
Isocentre
Do I understand it correctly that OpenTPS automatically sets the isocentre at the centre of the CT or target (if it exists)? What if I want to start the calculation at the surface of my CT? I have a 20x20x20 mm3 water phantom. What I did was I put the the Nozzle exit to Isocentre distance to 10 mm (half of the size of my CT). Does this mean that the calculation then starts at the surface of my CT with no air gap in between?
Also, is there an easy way to specify in the plan that I want to calculate just one spot at a certain energy? What I do is I just modify the input file to contain only one spot and load that. It seems like OpenTPS creates a spot scan plan for the target automatically.
Thanks!
Best,
Justin
Dear Justin,
For the computation of a single spot, you can use the following example code:
from opentps.core.data.plan._rtPlan import RTPlan
from opentps.core.io.scannerReader import readScanner
from opentps.core.io import mcsquareIO
from opentps.core.processing.doseCalculation.doseCalculationConfig import DoseCalculationConfig
from opentps.core.processing.doseCalculation.mcsquareDoseCalculator import MCsquareDoseCalculator
from opentps.core.data.plan._planIonBeam import PlanIonBeam
from opentps.core.data.plan._planIonLayer import PlanIonLayer
from opentps.core.data.images._ctImage import CTImage
import numpy as np
import matplotlib.pyplot as plt
ctCalibration = readScanner(DoseCalculationConfig().scannerFolder)
bdl = mcsquareIO.readBDL(DoseCalculationConfig().bdlFile)
plan = RTPlan()
plan.appendBeam(PlanIonBeam())
plan.beams[0].gantryAngle = 0
plan.beams[0].appendLayer(PlanIonLayer(100)) # 100 MeV
plan[0].layers[0].appendSpot([75], [75], [1]) # 75 75 are the x and y coordinates and 1 is the number of MU's
ctSize = 150
ct = CTImage()
ct.name = 'CT'
huAir = -1024.
huWater = ctCalibration.convertRSP2HU(1.)
data = huWater * np.ones((ctSize, ctSize, ctSize))
data[:, :50, :] = huAir
ct.imageArray = data
mc2 = MCsquareDoseCalculator()
mc2.beamModel = bdl
mc2.ctCalibration = ctCalibration
mc2.nbPrimaries = 1e5
dose = mc2.computeDose(ct, plan)
plt.imshow(dose.imageArray[75, :, :])
plt.show()
You can see that instead of building a plan (through planDesign) we can directly create a layer and a spot for the RTplan.
Regarding the stopping powers, the default ones used in MCsquare are from the Geant4 tables.
A colleague of mine will come back to you shortly with a more detailed answer on the calculation between the nozzle and the surface of the CT and on the stopping powers tables for MCsquare.
I hope this partially answers your question, do not hesitate if you have more questions.
Best regards,
Eliot
Dear Eliot,
I actually figured out my problem with the different ranges in MCsquare and TOPAS. It seems like there is an artefact in the dose calculation in TOPAS when I use the CTs. When I tried with just creating the geometry in TOPAS and putting a parallel scoring grid on top of it, I get more or less identical results with MCsquare.
I am still a bit confused about the isocentre issue so I guess I’ll wait for that one. Thanks!
Best regards,
Justin
1 Like
Hi Justin,
Your understanding of the Isocenter placement is actually pretty spot on. OpenTPS places the isocenter in the center of the target. You can check what your isocenter position is if you look in the file PlanPencil.txt in the opentps workspace after the files for the MCsquare simulation have been generated. The position of the nozzle is defined in the bdl file. If you want to start the calculation at the surface of your CT you can set the parameter bdl.nozzle_isocenter to your isocenter distance, which seems to be what you’re already doing.
Just in case it is relevant for you: The actual Monte Carlo simulation begins at the surface of the CT either way. If the nozzle is placed some distance from it, the particles are transported across the resulting air gap in a single step while taking the stopping power of air into account, but without simulating individual interactions.
I hope that answers your remaining questions. Let us know if things are still unclear or you have further issues!
Best,
Danah
Hi Danah,
Thanks for the explanation. I think everything’s clearer now.
About that nozzle being placed some distance away, do you also consider air scattering? If there’s some distance, then the spot size becomes larger. I guess that is less of a problem for humans but I am doing calculations on mice and since they’re very small, the dose calculations are quite sensitive with slight blow up of the beam.
Best regards,
Justin
Hi Justin,
A widening of the beam due to beam properties and initial velocity is taken into account, but air scattering is not. If you want to simulate it, you could add the air gap to your CT data, by defining rows of voxels around it corresponding to your gap and then set the nozzle on the surface of the modified CT. Then you’ll have a full MC simulation all the way.
Best,
Danah
Hi Danah,
Thanks! I figured that air scattering is not taken into account which is why I scored the beam parameters at the position of the animal container (i.e. after the air gap from the collimator) and did the calculations from the surface of the CT. I think that’s a better solution than adding air in the CT data.
Best regards,
Justin