Depth dose curves for large fields

[JAMKARAN]

Hi

For square fields (3x3, 5x5, 8x8, 10x10,12x12, 15x15, 20x20,25x25 and 30x30 cm2) I use the following code where bixel width corresponds to the side length of the square at the iso center. For 3x3, 5x5, 8x8, 10x10 cm2 PDD curves vary when side length of the square but for larger filed (12x12, 15x15, 20x20,25x25 and 30x30 cm2) , PDD curves become similar.

matRad_rc

% load patient data, i.e. ct, voi, cst
load BOXPHANTOM.mat

% meta information for treatment plan
pln.radiationMode   = 'photons';     % either photons / protons / carbon
pln.machine         = 'Generic';

pln.numOfFractions  = 30;

% beam geometry settings
pln.propStf.bixelWidth      = 100; % [mm] / also corresponds to lateral spot spacing for particles
pln.propStf.gantryAngles    = 0; % [°]
pln.propStf.couchAngles     = 0; % [°]
pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);

% dose calculation settings
pln.propDoseCalc.doseGrid.resolution.x = ct.resolution.x; % [mm]
pln.propDoseCalc.doseGrid.resolution.y = ct.resolution.y; % [mm]
pln.propDoseCalc.doseGrid.resolution.z = ct.resolution.z; % [mm]

% optimization settings
pln.propOpt.optimizer       = 'IPOPT';
pln.propOpt.bioOptimization = 'none'; % none: physical optimization;             const_RBExD; constant RBE of 1.1;
                                      % LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of RBE-weighted dose
pln.propOpt.runDAO          = false;  % 1/true: run DAO, 0/false: don't / will be ignored for particles
pln.propOpt.runSequencing   = false;  % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below

%% define steering file by hand

SSD = 800;

if pln.propStf.gantryAngles == 0 && pln.propStf.gantryAngles == 0
    stf.gantryAngle       = pln.propStf.gantryAngles;
    stf.couchAngle        = pln.propStf.couchAngles;
    stf.bixelWidth        = pln.propStf.bixelWidth;
    stf.radiationMode     = 'photons';
    stf.SAD               = 1000;
    stf.isoCenter         = [240 240+878-SSD 240];
    stf.numOfRays         = 1; 
    stf.sourcePoint_bev   = [0 -stf.SAD 0];
    stf.sourcePoint       = [0 -stf.SAD 0];
    stf.numOfBixelsPerRay = 1;
    stf.totalNumOfBixels  = 1;

    stf.ray.rayPos_bev = [0 0 0];
    stf.ray.targetPoint_bev = [0 stf.SAD 0];
    stf.ray.rayPos = [0 0 0];
    stf.ray.targetPoint = [0 stf.SAD 0];
    stf.ray.beamletCornersAtIso = nan; % not needed for pencil beam dose calc
    stf.ray.rayCorners_SCD = nan; % not needed for pencil beam dose calc
    stf.ray.energy = 6;
else
    error('not implemented');
end

%% dose calculation
resultGUI = matRad_calcDoseDirect(ct,stf,pln,cst,1);

%% start gui for visualization of result
matRadGUI

[wahln]

Hi, and what exactly is your issue / your question?

[JAMKARAN]

I think for large fields matRad understimate absorbd dose on central axis especially at higher depths. If not, why PDD curves of large fields overlap.

[markbangert]

So you'd expect a continued decline in dose with increasing field size? What physical effect could induce such behaviour.

JAMKARAN notifications@github.com schrieb am Di., 1. Sept. 2020, 07:45:

I think for large fields matRad understimate absorbd dose on central axis especiallhy at higher depths.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/e0404/matRad/issues/441#issuecomment-684362631, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNPY7L2SD6OMSSKXGKXEP3SDSC73ANCNFSM4QQNUJYQ .

[JAMKARAN]

Thanks Mark

because I have observed in commissioning for accelerators -PDD curves change as the field size increases- in fact with an increment of field size, PDD curves increase. In the following can be seen:

**Field Size and Shape**

The field size is directly related to the amount of primary radiation entering the patient and the resulting dose distribution. A dose at any point results from both the primary and scattered radiation. Larger field sizes lead to increased generation of scattered radiation, which tends to increase dose at a specified point. Conversely, small field sizes may lead to loss of electronic equilibrium and the central flat part of the beam. As field size increases:

Zmax increases (less rapidly for very large field sizes over 20 cm) PDD increases at other points along the central beam axis The ratio of the penumbra to the central portion of the field decreases

As field size decreases:

Zmax decreases (rapidly for field sizes under 6 cm) PDD decreases at points along the central beam axis The ratio of the penumbra to the central portion increases, and may lead to complete loss of the central flat section

[JAMKARAN]

Reminder @wahln @markbangert

[wahln]

Patience please.... We are not really a customer service team here.

I could not really reconstruct your images, which is also related to me not really understanding the "PDD" curves you show. PDD is dose relative to the maximum, so each curve should occupy the interval [0,1]. The y axis say "unnormalized", but in the context of a PDD I don't know what unnormalized means.

Here are the curves that I find by iterating through the various field sizes (left absolute dose, right PDDs).

The PDD shows the expected effect of relatively increasing in depth with increasing field size. Also there is no equilibrium within the field size range. Note that the x-axes is not completely correct because it was just a quick hack on the full dose cube. Here's my adaptation of your script including the loop:

matRad_rc

% load patient data, i.e. ct, voi, cst
load BOXPHANTOM.mat

% meta information for treatment plan
pln.radiationMode   = 'photons';     % either photons / protons / carbon
pln.machine         = 'Generic';

pln.numOfFractions  = 30;

% beam geometry settings
pln.propStf.bixelWidth      = 100; % [mm] / also corresponds to lateral spot spacing for particles
pln.propStf.gantryAngles    = 0; % [°]
pln.propStf.couchAngles     = 0; % [°]
pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);

% dose calculation settings
pln.propDoseCalc.doseGrid.resolution.x = ct.resolution.x; % [mm]
pln.propDoseCalc.doseGrid.resolution.y = ct.resolution.y; % [mm]
pln.propDoseCalc.doseGrid.resolution.z = ct.resolution.z; % [mm]

% optimization settings
pln.propOpt.optimizer       = 'IPOPT';
pln.propOpt.bioOptimization = 'none'; % none: physical optimization;             const_RBExD; constant RBE of 1.1;
                                      % LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of RBE-weighted dose
pln.propOpt.runDAO          = false;  % 1/true: run DAO, 0/false: don't / will be ignored for particles
pln.propOpt.runSequencing   = false;  % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below

%% define steering file by hand

SSD = 800;

bixelWidths = [3 5 8 10 12 15 20 25 30];

hF = figure;
hAxAbs = subplot(1,2,1);
hAxNorm = subplot(1,2,2);

for width = bixelWidths
    pln.propStf.bixelWidth = width;
    if pln.propStf.gantryAngles == 0 && pln.propStf.gantryAngles == 0
        stf.gantryAngle       = pln.propStf.gantryAngles;
        stf.couchAngle        = pln.propStf.couchAngles;
        stf.bixelWidth        = pln.propStf.bixelWidth;
        stf.radiationMode     = 'photons';
        stf.SAD               = 1000;
        stf.isoCenter         = [240 240+878-SSD 240];
        stf.numOfRays         = 1;
        stf.sourcePoint_bev   = [0 -stf.SAD 0];
        stf.sourcePoint       = [0 -stf.SAD 0];
        stf.numOfBixelsPerRay = 1;
        stf.totalNumOfBixels  = 1;

        stf.ray.rayPos_bev = [0 0 0];
        stf.ray.targetPoint_bev = [0 stf.SAD 0];
        stf.ray.rayPos = [0 0 0];
        stf.ray.targetPoint = [0 stf.SAD 0];
        stf.ray.beamletCornersAtIso = nan; % not needed for pencil beam dose calc
        stf.ray.rayCorners_SCD = nan; % not needed for pencil beam dose calc
        stf.ray.energy = 6;
    else
        error('not implemented');
    end
    resultGUI_tmp = matRad_calcDoseDirect(ct,stf,pln,cst,1);
    plot(hAxAbs,ct.resolution.x*(1:ct.cubeDim(2)),resultGUI_tmp.physicalDose(:,80,80),'DisplayName',sprintf('%dx%d',width,width),'LineWidth',1.5); hold(hAxAbs,'on');
    plot(hAxNorm,ct.resolution.x*(1:ct.cubeDim(2)),resultGUI_tmp.physicalDose(:,80,80) ./ max(resultGUI_tmp.physicalDose(:)),'DisplayName',sprintf('%dx%d',width,width),'LineWidth',1.5); hold(hAxNorm,'on');
    resultGUI.(['physicalDose_fSize' num2str(width)]) = resultGUI_tmp.physicalDose;
end
title(hAxAbs,'absolute dose'); legend(hAxAbs);
xlim(hAxAbs,[100 400]);
title(hAxNorm,'PDD'); legend(hAxNorm);
xlim(hAxNorm,[100 400]);
ylim(hAxNorm,[0.1 1.1]);

%% start gui for visualization of result
matRadGUI

[JAMKARAN]

Very thanks @wahln and also, I'm so sorry for reminder. I have another question about your script. because isosenter of box phantom is located at center of box you set stf.isoCenter = [240 240+878-SSD 240]; if isocenter was on surface of phantom and SSD=SAD=100 cm, Is the the below correct?

SSD = 1000;
.
.
.
stf.SAD               = 1000;
stf.isoCenter         = [240 0 240];
.
.
.
stf.sourcePoint_bev   = [0 -stf.SAD 0];
stf.sourcePoint       = [0 -stf.SAD 0];

[markbangert]

Thank you @wahln for taking over :)

One addition from my side. Please note that effects regarding lateral electron equilibrium are accounted for in the kernel calculation via the output factor -> https://github.com/e0404/photonPencilBeamKernelCalc/blob/master/referenceData/of.dat.

[JAMKARAN]

Thanks @markbangert Excuse me @markbangert , is it possible clarify The effects of lateral electron equilibrium for dose calculation?

[wahln]

Regarding the isocenter: If you set the coordinate to 0, the isocenter will be at the beginning of the CT (which is air and therefore not the same as the surface of the watercube). Note that there is also a display bug in this script, because while the isocenter is changed in the stf, it is not changed in the pln (in the script the isocenter is actually deeper than the center). I think the water surface is at 118.5 in the boxphantom.

I restructured the script a little to use (1) the correct machine parameters, (2) use a preliminary SSD calculation to compute the surface coordinates and (3) put the isocenter on the surface automatically (also for display in the GUI):

matRad_rc

% load patient data, i.e. ct, voi, cst
load BOXPHANTOM.mat

% meta information for treatment plan
pln.radiationMode   = 'photons';     % either photons / protons / carbon
pln.machine         = 'Generic';

pln.numOfFractions  = 30;

% beam geometry settings
pln.propStf.bixelWidth      = 100; % [mm] / also corresponds to lateral spot spacing for particles
pln.propStf.gantryAngles    = 0; % [°]
pln.propStf.couchAngles     = 0; % [°]
pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);

% dose calculation settings
pln.propDoseCalc.doseGrid.resolution.x = ct.resolution.x; % [mm]
pln.propDoseCalc.doseGrid.resolution.y = ct.resolution.y; % [mm]
pln.propDoseCalc.doseGrid.resolution.z = ct.resolution.z; % [mm]

% optimization settings
pln.propOpt.optimizer       = 'IPOPT';
pln.propOpt.bioOptimization = 'none'; % none: physical optimization;             const_RBExD; constant RBE of 1.1;
                                      % LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of RBE-weighted dose
pln.propOpt.runDAO          = false;  % 1/true: run DAO, 0/false: don't / will be ignored for particles
pln.propOpt.runSequencing   = false;  % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below

%% define steering file by hand

bixelWidths = [3 5 8 10 12 15 20 25 30];

load(['basedata/' pln.radiationMode '_' pln.machine '.mat']);
stf.gantryAngle       = pln.propStf.gantryAngles;
stf.couchAngle        = pln.propStf.couchAngles;
stf.bixelWidth        = pln.propStf.bixelWidth;
stf.radiationMode     = pln.radiationMode;
stf.SAD               = machine.meta.SAD;
stf.isoCenter         = pln.propStf.isoCenter;
stf.numOfRays         = 1;
stf.sourcePoint_bev   = [0 -stf.SAD 0];
stf.sourcePoint       = [0 -stf.SAD 0];
stf.numOfBixelsPerRay = 1;
stf.totalNumOfBixels  = 1;

stf_tmp = matRad_computeSSD(stf,ct);
stf.isoCenter(2) = stf.isoCenter(2) - (stf.SAD - stf_tmp.ray.SSD);
pln.propStf.isoCenter = stf.isoCenter;

hF = figure;
hAxAbs = subplot(1,2,1);
hAxNorm = subplot(1,2,2);

for width = bixelWidths
    pln.propStf.bixelWidth = width;
    stf.bixelWidth = width;   

    resultGUI_tmp = matRad_calcDoseDirect(ct,stf,pln,cst,1);
    plot(hAxAbs,ct.resolution.x*(1:ct.cubeDim(2)),resultGUI_tmp.physicalDose(:,80,80),'DisplayName',sprintf('%dx%d',width,width),'LineWidth',1.5); hold(hAxAbs,'on');
    plot(hAxNorm,ct.resolution.x*(1:ct.cubeDim(2)),resultGUI_tmp.physicalDose(:,80,80) ./ max(resultGUI_tmp.physicalDose(:)),'DisplayName',sprintf('%dx%d',width,width),'LineWidth',1.5); hold(hAxNorm,'on');
    resultGUI.(['physicalDose_fSize' num2str(width)]) = resultGUI_tmp.physicalDose;
end
title(hAxAbs,'absolute dose'); legend(hAxAbs);
xlim(hAxAbs,[100 400]);
title(hAxNorm,'PDD'); legend(hAxNorm);
xlim(hAxNorm,[100 400]);
ylim(hAxNorm,[0.1 1.1]);

%% start gui for visualization of result
matRadGUI

[JAMKARAN]

Very thanks @wahln

[markbangert]

Sorry for the delay. I was on holidays...

I have the following understanding of the term "lateral electron equilibrium": For photon irradiation, the main dose deposition mechanism is via secondary electrons that travel on the range of mm-cm for therapeutic energies. On the central axis for small fields, there is more lateral scatter away from the central axis than towards the central axis (merely because the field extend in lateral direction is small, thence less secondary electrons are generated away from the central axis). Consequently we have a lower dose for smaller fields, which is accounted for by the output factor during kernel calculation (cf https://pubmed.ncbi.nlm.nih.gov/8497215/).

I will close the issue for the time being. Please feel free to reopen if there are remaining open points on your end.