e0404 / matRad

An open source multi-modality radiation treatment planning sytem developed by e0404 @ DKFZ
http://www.matRad.org
Other
223 stars 176 forks source link

Question on the Proton_generic machine data #336

Closed KW-Wong closed 3 years ago

KW-Wong commented 5 years ago

I am a new user to MatRad and I have a question which I cant find the info on the internet. I would like to know how matRad determine the range and peak position of the proton beam for a particular energy? In selecting the proton energies for the SOBP, which definition of treatment depth is used? (it seems to be it is d_100)

markbangert commented 5 years ago

Hi KW-Wong, you are right. The energies are selected according to the peak position (or d_100/d_max) in the function https://github.com/e0404/matRad/blob/master/matRad_generateStf.m around line 250. The peak position is stored along with the base data for every energy... Let me know if that clears up things for you?! Cheers / Mark

KW-Wong commented 5 years ago

Hello Mark,

Thanks for your reply on the treatment depth. For the range and peak position in the machine data, how do you determine them? via a specific equation or from experiment? I am a bit puzzled about the definition of range, is that a CSDA range?

I tried with some depth-dose curve equation derived from the Bethe law, the peak position that I found from the equation matches the range value in the machine data. I am not sure whether I did it incorrectly, so I want to know what definition of range is used in the machine data table.

Thank you very much for your time.

Regards, Kin Wing

markbangert commented 5 years ago

Hi Kin Wing,

you just need to make sure that the data you use for planning is consistent. Therefore make sure that

[~,ix] = max(machine.data(i).Z)
machine.data(i).peakPos = machine.data(i).depths(ix)

for all i.

That's it!

If you have a non-zero offset in range for the measured IDD data in machine.data(i).offset (e.g. if the standard measurement protocol does not include some ionization chambers that are in the beam for treatment) this automatically taken care of in matRad_generateStf.

Still unclear?

KW-Wong commented 5 years ago

Hi Mark,

I understand the implementation in matRad. My question is to ask about the way of creating the machine data for peak position. And I found that you guys created the machine data using the expressions from some literature which mentioned in matRad wiki.

Thank you for your help.

markbangert commented 5 years ago

Argh now I finally get it. Yes, for the generation of the generic base data we have an implementation for the proton Bragg curve according to https://www.ncbi.nlm.nih.gov/pubmed/9434986 and for the lateral scattering we use https://www.sciencedirect.com/science/article/pii/0168583X9395944Z. Currently this is not open source. Do you need this code? What for? For us it would probably mean some work to put it in a form that we can offer it for download...

For the clinical machines we use measured and Monte Carlo simulated data!

markbangert commented 5 years ago

As I have not heard from you in a while I am going to close this issue...

Please note that I put the code to compute a multiple Coulomb scattering sigma at https://de.mathworks.com/matlabcentral/fileexchange/70416-multiple-coulomb-scattering for the time being. If I find the time to clear up the code I will also put the Bragg curve computation at Matlab Central.

justingm commented 3 years ago

Hello!

I also am wondering about the same thing as Kin Wing. What is the definition of the range in the base data? I find that it's a bit redundant if you already have the Bragg peak position which you use to select the energy. Do you use the range information in some other part in matRad?

Also, is the code to compute the Bragg curve also available? Thanks!

Best, Justin

wahln commented 3 years ago

For the generic base data, range should be equal to the CSDA range. But as you said, the field is not necessary. For the positioning of the Bragg-Peaks only the peakPos is necessary. There are, I think, a few computation where the range field can be used, but is not needed (we also have some internal base data that does not have a range field).

Regarding the code for generation of pencil-beams: It is not within matRad, but we have it internally. @markbangert since that was your work, would you have any issues with us sharing the code, e.g., by mail? Or I could also clean it up a little and add it to matRad.

markbangert commented 3 years ago

We can def put this function online. Originally this was not written by myself, though, but by Thomas Bortfeld - who told me that he would be fine with that code being published :)

The range in the base data corresponds to the range parameter within this very function that we used to generate the base data. Consequently, we also put it in as a reference for ourselves. And if I am not mistaken his corresponds to the 80% dose at the distal falloff - which is independent from range straggeling effects.

justingm commented 3 years ago

I see. Thanks for the clarification. We are planning to create our own base data file for our beam line and I am trying to understand all the information that I'd need to gather.

Anyway, it would be a tremendous help if you can share the code for the calculation of the Bragg curves. In the meantime, I will read the reference paper. :)

I hope both of you will have a great Easter weekend!

Best, Justin

wahln commented 3 years ago

OK I will take a look after easter if I can prepare the code snippet.

I also want to refer to my answer on #500 :

What data do you want to use to create your own base-data in matRad? do you have phase space parameterizations? do you have IDDs + beam widths? Maybe then I could tell you more exactly what parts of the code could be usefeul.

justingm commented 3 years ago

Perfect! Thank you so much!

justingm commented 3 years ago

Hello Niklas,

I am wondering if you've taken a look at Bortfeld's Bragg curve code. But I understand if you are busy rt now. :)

Best, Justin

justingm commented 3 years ago

Also, I tried out the the MCS code from the link above and compared them to my simulations in Geant4. Are the values of sigma in cm or mm? In the plot, it says cm but when I compared it to my results, it agrees if it is in mm.

stale[bot] commented 3 years ago

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

wahln commented 3 years ago

sorry for the late response @justingm , but I found the function and pushed it to a special branch specials/analyticalBaseDataGeneration. Direct link: https://github.com/e0404/matRad/blob/specials/analyticalBaseDataGeneration/basedata/generation/matRad_getDDDfromAnalyCalc.m Feel free to play around with it but don't take its output for granted!

Note: It requires the parabolic cylinder function from the Matlab file Exchange. check the links in the code!

justingm commented 3 years ago

@wahln Thanks a lot!

justingm commented 3 years ago

@wahln I have been looking at the code while reading Bortfeld's paper. How did you calculate relEnergySigma?

wahln commented 3 years ago

That's just your assumption on the spread of the beam's initial energy spectrum parametrized by its relative standard deviation. If I remember correctly, the energy spectrum is characterized by two parameters in Bortfeld's paper, one for the width and one for the tail (which in the code is close to zero, so we assume an almost symmetric, gaussian-like energy distribution). Thus, relEnergySigma it's just a parameter you can play around with. I think the original code it used an absolute value in MeV instead.

justingm commented 3 years ago

Indeed, there are two parameters for the energy spectrum. First is the peak which is modelled as a Gaussian with sigmaE = relEnergySigma*E0 and second is the tail. I saw in the paper that relEnergySigma = 0.01 and in the code it is 0.014. I was just wondering how that was derived. Thinking about it, this parameter would depend on our beam and is something that I need to find out.

By the way, I will be joining the masterclass next week. I don't know if you are involved in that. We have plans on using matRad extensively so I think it would be nice to participate on that.

Thanks again for you help!!

wahln commented 3 years ago

Yes, exactly, this would refer to the energy spectrum of your machine. For the generic dataset that's public in matRad, there's no specific machine considered, but it's quite unfortunate that these parameters are not given. I should look up / find out which ones were used for the proton_Generic basedata in the first place.

justingm commented 3 years ago

Hello @wahln! I have been playing around the code. I am also wondering which parameters you used on the proton_Generic basedata. When I compare it with the results from the code you sent, the plateau is higher and the Bragg peak position is earlier for the proton_Generic basedata (for E=90.5576). I also did some simulations in Geant4 but with a pencil beam (zero width, zero divergence) and the Bragg peak position that I get is at 61.8 mm which is closer to the peak position in proton_Generic basedata.

wahln commented 3 years ago

@markbangert or @becker89 can you comment on the exact parameters and setup for the proton_Generic calculation? However I think that differences can also come from the numerical side, i.e., the implementation of the parabolic cylinder function, I think we used a different one for computation of the data. As said, I am not completely sure that the uploaded code can exactly reproduce the proton_Generic data.

markbangert commented 3 years ago

To be completely honest with you: I think we did not document this process very well - unfortunately...

Am Di., 18. Mai 2021 um 18:02 Uhr schrieb Niklas Wahl < @.***>:

@markbangert https://github.com/markbangert or @becker89 https://github.com/becker89 can you comment on the exact parameters and setup for the proton_Generic calculation? However I think that differences can also come from the numerical side, i.e., the implementation of the parabolic cylinder function, I think we used a different one for computation of the data. As said, I am not completely sure that the uploaded code can exactly reproduce the proton_Generic data.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/e0404/matRad/issues/336#issuecomment-843299512, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNPY7PEZ2GKUKZLEQMDGD3TOKFRPANCNFSM4GY4RQBQ .

wahln commented 3 years ago

@markbangert do you know which implementation of the parabolic cylinder function was used? That might probably already help.

markbangert commented 3 years ago

No clue anymore. Sorry.

Niklas Wahl @.***> schrieb am Mi., 19. Mai 2021, 19:10:

@markbangert https://github.com/markbangert do you know which implementation of the parabolic cylinder function was used? That might probably already help.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/e0404/matRad/issues/336#issuecomment-844301267, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNPY7O4MYYO57WOGL74QLLTOPWH5ANCNFSM4GY4RQBQ .

justingm commented 3 years ago

That's fine. I saw that the original code in python was also linked so I am going to try that out. Besides, we also need to do some rigorous experiments to confirm the values. (currently waiting for beam time 🙏🏻)

wahln commented 3 years ago

I think it might then be time for a new & reproducible generic data set at some point while keeping the old for repredocuibility of older results ;-)

stale[bot] commented 3 years ago

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

stale[bot] commented 3 years ago

This issue has been automatically closed because it has not seen any activity in four weeks. This happens usually when the issue has already been solved or it is no longer relevant. If that's not the case, feel free to reopen the issue.

WeiWangHUST commented 3 years ago

That's fine. I saw that the original code in python was also linked so I am going to try that out. Besides, we also need to do some rigorous experiments to confirm the values. (currently waiting for beam time 🙏🏻)

Dear, I also want to know how to get the machine data, but I can't find the original code you said. Can you shall the link with me? Thanks!

wahln commented 3 years ago

I am quite sure you'll get closest to the original code using the snippets from here: https://gray.mgh.harvard.edu/software/293-analytical-bragg-curve

Erhushenshou commented 2 years ago

Got it. I'll have a try and then discuss. Many thanks.

Erhushenshou commented 2 years ago

Argh now I finally get it. Yes, for the generation of the generic base data we have an implementation for the proton Bragg curve according to https://www.ncbi.nlm.nih.gov/pubmed/9434986 and for the lateral scattering we use https://www.sciencedirect.com/science/article/pii/0168583X9395944Z. Currently this is not open source. Do you need this code? What for? For us it would probably mean some work to put it in a form that we can offer it for download...

For the clinical machines we use measured and Monte Carlo simulated data!

Hi Mark, Do you have code for modeling the lateral using double gaussian model? Thanks.