# How to plot Bragg curve?

1. Aug 27, 2009

### Ciolla

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. Aug 27, 2009

### TheDestroyer

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 :)

3. Aug 27, 2009

### Ciolla

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
4. Aug 27, 2009

### TheDestroyer

Welcome :), hehehe... good that you figured it out by yourself :)

5. Aug 27, 2009

### Ciolla

sorry, but what do you mean?

6. Aug 27, 2009

### Bob S

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?

7. Aug 27, 2009

### Ciolla

8. Aug 27, 2009

### Bob S

Here is the bethe-Bloch equation.
http://en.wikipedia.org/wiki/Bethe_formula
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.

9. Aug 27, 2009

### Bob S

Here is a thumbnail plot of dE/dx vs. Energy, and the program for it.

#### Attached Files:

File size:
38.2 KB
Views:
824
• ###### DEDX_Proton_prog.jpg
File size:
20.6 KB
Views:
926
10. Aug 27, 2009

### Bob S

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:

• ###### DEDX_Proton_plot2.jpg
File size:
42.7 KB
Views:
711
Last edited by a moderator: May 4, 2017
11. Aug 28, 2009

### Ciolla

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!

12. Aug 28, 2009

### Bob S

Ciolla-
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):
http://physics.nist.gov/PhysRefData/Star/Text/contents.html
Bob S

13. Aug 28, 2009

### Ciolla

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!?
thanks!!

14. Aug 28, 2009

### Bob S

Ciolla-
Here it is. It is written in TrueBasic.
Bob S

#### Attached Files:

File size:
23.6 KB
Views:
601
• ###### DEDX_Proton_prog2_sub.jpg
File size:
12.6 KB
Views:
472
15. Aug 29, 2009

### Ciolla

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
16. Aug 29, 2009

### uart

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 :

$$-\frac{d E}{d x} = f(E)$$

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

$$|\Delta E| = f(E) \Delta x$$

$$E_{\mbox{next}} = E - |\Delta E|$$

$$x_{\mbox{next}} = x + \Delta x$$

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 :

$$-\frac{d E}{d x} = f(\beta)$$

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

$$E = m_0 c^2 \left(\frac{1}{\sqrt{1-\beta^2}} - 1 \right)$$

Re-arrange it to get $\beta$ in terms of E and then you're good to go. Do you follow that Ciolla?

Last edited by a moderator: May 4, 2017
17. Aug 29, 2009

### Ciolla

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

18. Aug 29, 2009

### uart

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
19. Aug 29, 2009

### uart

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
20. Aug 29, 2009

### Bob S

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

21. Aug 29, 2009

### Ciolla

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!
thanks!!!!

22. Aug 29, 2009

### Bob S

Ciolla-
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:

• ###### DEDX_Ionization_const.jpg
File size:
51.7 KB
Views:
359
23. Sep 2, 2009

### Ciolla

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!?

24. Sep 2, 2009

### Bob S

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

25. Sep 2, 2009

thank:D