1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

How to plot Bragg curve?

  1. Aug 27, 2009 #1
    Goodmorning everyone, I'm new to this forum!
    I'm graduing in physics in Italy, and I have this problem.
    I want use root to plot the Bragg curve for different kinds of particles, with different initial conditions. the problem is, wich function have I to plot?


    the Beth Bloch formula? how can I make it "x dependent", where x is the penetration depth?

    sorry for my english, and my silly question..
    I hope you'll help me!
    thanks a lot!
  2. jcsd
  3. Aug 27, 2009 #2
    I've never seen this function before, although I've been studying physics for 6 years. Please try to explain in what subject this is studied.

    I can give a way to determine such a function, but you'll have to work a little.

    Get some program for curve fitting like "curve expert 1.3", it's a free share ware you can find it on the internet for free, and then pick up some points in the function, and try to interpolate it.

    Then, the program would interpolate so many functions, try to integrate the best one that fits or the one that makes sense to you (indefinite integral), and then determine the constant by some condition you have for T.

    Good luck :)
  4. Aug 27, 2009 #3
    uhu.. so complicated?
    i try to explain a little more..
    i'm studying Hadrontherapy, i.e. the use of light ions or protons against tumors. this kind of particles are characterized by a particular form of loosing-energy curve. as you se here
    http://www.gsi.de/forschung/bio/e-deposition_small.jpg [Broken]
    the advantage is that they loose energy manily at a particular depth. leaving the healty tissue untouched.
    so, the problem is how to express the dipendence from x in the bethe bloch formula, that gives the energy loss for ionization of a particle.

    i think that there is a solution! and a simple one! simpler of the one you suggested!

    btw, thank you!
    Last edited by a moderator: May 4, 2017
  5. Aug 27, 2009 #4
    Welcome :), hehehe... good that you figured it out by yourself :)
  6. Aug 27, 2009 #5
    sorry, but what do you mean?
  7. Aug 27, 2009 #6
    Do you want to plot the penetration depth (range) vs. particle energy for various particles (protons, alphas, etc.), by integrating the Bethe-Bloch dE/dx formula? If not, what are the ordinate and abscissa of what you want to plot?
  8. Aug 27, 2009 #7
  9. Aug 27, 2009 #8
    Here is the bethe-Bloch equation.
    It is correct except for very low energies (e.g., slow protons), where positive moving charged particles attract atomic electrons, and negative charged repel them. Barkas et al. compared emulsion tracks of stopping positive and negative muons to help resolve this.
  10. Aug 27, 2009 #9
    Here is a thumbnail plot of dE/dx vs. Energy, and the program for it.

    Attached Files:

  11. Aug 27, 2009 #10
    Here in thumbnail I think is a proton dE/dx plot vs. range in water for a 200 MeV proton. There is no range or energy straggling included. The max range was 25.9 grams per cm^2 (25.9 cm). Please check against your tables.

    [Note: the plot has been updated]

    Note 2: Look up Geant at http://wwwasd.web.cern.ch/wwwasd/geant/ [Broken]

    Attached Files:

    Last edited by a moderator: May 4, 2017
  12. Aug 28, 2009 #11
    hi bob, first, thank you, really so much!
    even if i don't know the program you used to plot the first graph, that is the "classical" BB formula, i'd already achieved that result and that graph.
    the problem i have now is to obtain the second graph you posted.
    have you obtained it using geant (i'm trying it!)
    or hav you written a program similar to the previous one?
    with which forumula?

    thank you so much!
  13. Aug 28, 2009 #12
    The second plot was done using my own program similar to the first one, by numerical integration of the Bethe Bloch equation. The only tricky part is getting the correct ionization constant, because the approximation that I =10 for every material is only approximate. Furthermore, my curve is for only a single particle, so the range spread from incident beam energy spread and range straggling (including multiple scattering) is not included. Here is another site (besides GEANT) for proton and alpha particle range data for literally 100's of materials (National Institue of Science and Technology):
    Bob S
  14. Aug 28, 2009 #13
    so bob, the fact that the plot represent only one particle is not a problem. I've just to show the difference in behaviour of different particles.

    the point is: what program did you used?
    can you post it!?
  15. Aug 28, 2009 #14
    Here it is. It is written in TrueBasic.
    Bob S

    Attached Files:

  16. Aug 29, 2009 #15
    hi bob! step by step, you are saving my week! :D
    thanks for the program. however, I don't know true basic. btw, i'm trying to understand the process to translate it in C. let's try to sum up:
    - you define all the stuff..
    - you define the starting point for x and E.
    - an iteration starts, from 0 step by step in x
    - you define gam, bgam, beta as functions of total energy E
    - you define the function to plot dE/dx
    - you plot a single point (right!?)
    - the cicle ends when E reaches 0

    it's ok?
    my remaining doubts:
    -- the command: "PLOT x,0" -- what is it for?
    -- the command: "xx=round(x,2) ??
    -- isn't gam defined as E/mc2? why do you use 1+E/mpc2?

    thanks so much bob!
    Last edited: Aug 29, 2009
  17. Aug 29, 2009 #16


    User Avatar
    Science Advisor

    Perhaps it would be helpful if I just summarized the algorithm, without the fine details.

    The http://en.wikipedia.org/wiki/Bethe_formula" [Broken] can be thought of as being in the form of :

    [tex]-\frac{d E}{d x} = f(E)[/tex]

    So if you choose some suitably small iteration interval [itex]\Delta x[/itex] you can numerically iterate this formula as per :

    [tex]|\Delta E| = f(E) \Delta x[/tex]

    [tex] E_{\mbox{next}} = E - |\Delta E|[/tex]

    [tex] x_{\mbox{next}} = x + \Delta x[/tex]

    You get the idea, you just iterate the above steps while ever E is still positive.

    Before you can get started however you'll have to get the equation into the form shown above. Currently you have the equation in the form :

    [tex]-\frac{d E}{d x} = f(\beta)[/tex]

    where [itex]\beta = v/c[/itex]. So you'll need to use the relativistic energy equation to express [tex]\beta[/tex] in terms of E. Start with,

    [tex] E = m_0 c^2 \left(\frac{1}{\sqrt{1-\beta^2}} - 1 \right)[/tex]

    Re-arrange it to get [itex]\beta[/itex] in terms of E and then you're good to go. Do you follow that Ciolla?
    Last edited by a moderator: May 4, 2017
  18. Aug 29, 2009 #17
    thanks uart, i think I got your arguments. But now I'm worried of getting lost in confusion: i've been tought that gam=E/mc2, where E is the total eneregy of the particile, E^2=p^2+m0^2, but this is in contraddiction with "your" definition of E=m0(gam-1)
    I apologize for my insistence, but I'think i'm loosing some points!
    for now, thanks you so much
  19. Aug 29, 2009 #18


    User Avatar
    Science Advisor

    Well actually I wrote [itex]E=m_0 c^2 (\gamma-1)[/tex] as I assumed that the relevant "E" in the equation was the relativistic KE, so I didn't include the rest mass energy (that's the only difference between this and the other equation that you're quoting btw).

    Anyway that seems like a fair assumption to me but yeah you'll have to check the physics as that's not really my area, I was just trying to give you an overview of the maths/algorithm.
    Last edited: Aug 29, 2009
  20. Aug 29, 2009 #19


    User Avatar
    Science Advisor

    Ok thinking about it some more I see that it doesn't really matter whether you use the total energy or just the KE in that equation as they only differ by a constant (m0 c^2). You just have to make sure that you are consistent throughout with whichever one you choose (with things such as initial and final energy and the translation between energy and velocity especially). In particular you'll have to end the simulation when the energy is reduced to m0 c^2 if you express everything in terms of total energy but you end when E=0 if you just use KE. I think just using KE is going to be more straight forward here.
    Last edited: Aug 29, 2009
  21. Aug 29, 2009 #20
    Hi Ciolla-
    1) Looks good.
    2) "Plot x.0" terminates the "plot x,dEdx" commands by drawing a vertical line down to the horizontal axis.
    3) "xx=round(x,2)" rounds off to two places (but don't really need it here).
    4) I am using E as kinetic energy in the program. Sorry for the confusion.
    Bob S
  22. Aug 29, 2009 #21
    thank you two guys, now i'think i'm on!
    now is my turn to turn it into an executable program!
    i will keep update you!
  23. Aug 29, 2009 #22
    Here is a plot (see thumbnail) of the energy loss ionization constant for the Bethe-Bloch equation plotted vs. Z of the target material. I told you earlier that IZ (in my program) was about 10*Z, but it does vary a little bit, depending on Z.
    Bob S

    Attached Files:

  24. Sep 2, 2009 #23
    just one question, Bob.
    in your plot you've written that the x assis is grams per cm2. the x assis' unit is centimeter, isn't it!?
  25. Sep 2, 2009 #24
    Hi Ciolla-
    The x axis in my program in previous post is in units of grams/cm2. To plot the result in cm, you have to divide x by the density, as folows. See changes in bold:

    310 DO
    320 x=x+dx
    325 y=x/rho
    330 gam=1+E/mpc2
    340 bgam=sqr(gam^2-1)
    350 beta=bgam/gam
    360 WHEN error in
    370 factor=(1/beta^2)*(log(2*mec2*1E6*beta^2/(IZ*(1-beta^2)))-beta^2)
    380 dEdx=Const*factor
    390 PLOT y,dEdx;
    400 USE
    410 EXIT DO
    420 END WHEN
    430 dE=dEdx*dx
    440 E=E-dE
    450 LOOP until E<=0
    460 PLOT y,0
    470 yy=round(y,2)
    480 SET COLOR "black"
    490 PLOT TEXT, AT y,2.5: "range = " & str$(yy) & " cm"
    500 END

    Bob S
  26. Sep 2, 2009 #25
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook