how to plot Bragg curve?

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

http://www.med.harvard.edu/JPNM/physics/nmltd/radprin/sect7/7.1/bragg.GIF

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!

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

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

[TheDestroyer]

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

[Ciolla]

sorry, but what do you mean?

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

[Ciolla]

i want to plot the energy-loss vs penetration energy
that is wat you see here:
http://www.med.harvard.edu/JPNM/physics/nmltd/radprin/sect7/7.1/bragg.GIF
(i suppose..)

[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.

[Bob S]

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

[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/

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

[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

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

[Bob S]

the point is: what program did you used?
can you post it!?
thanks!!
Ciolla-
Here it is. It is written in TrueBasic.
Bob S

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

[uart]

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:


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

The Bethe_formula (http://en.wikipedia.org/wiki/Bethe_formula) 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?

[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

[uart]

but this is in contraddiction with "your" definition of E=m0(gam-1)


Well actually I wrote E=m_0 c^2 (\gamma-1)[/tex] as I [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.

[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.

[Bob S]

hi bob! 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!
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

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

[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

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

[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

[Ciolla]

thank:D

[ekosulistya]

Hi,

If I want to use alfa particle and air as target, which one do I have to change in that code?

Thanks.

[Bob S]

For the Bethe Bloch equation, use Eq(1) in

http://en.wikipedia.org/wiki/Bethe_formula

For the value of the factors in ()2, look up the value of the classical electron radius in

http://pdg.lbl.gov/2009/reviews/rpp2009-rev-phys-constants.pdf

Calculate initial value of β from initial kinetic energy E

calculate dE/dx from equation

penetration step increment x + Δx (x is in grams per cm2)

new energy = E - (dE/dx)Δx (where E is kinetic energy)

calculate new β from new E

iterate until E = 0

To convert x to length (cm), divide x (in grams/cm2) by density (in grams per cm3. For air use 0.0012 grams/cm3.

Bob S

[ekosulistya]

Thanks Bob,

I rewrite your code in Python, and get same graph. Here is the code :

import math
mpc2=938.27 #proton mass in MeV
mec2=0.511 #electron mass in MeV
re=2.817e-13 #clasical electron radius
rho=1 #density of water in grams per cm3
IZ=75 #ionisation constant for water times Z-eff
Z=1 #Z proton charge
TZ=10 #Z-eff target
TA=18 #A-eff target
n=6.02e23*rho*TZ/TA #target electron density
konst=(4*math.pi*n*Z**2*mec2*re**2)
f=open('F:/proton_water.txt','w')
x=0
dx=0.01
Estart=5
E=Estart
while E>0:
x=x+dx
gam=1+E/mpc2
bgam=math.sqrt(gam**2-1)
beta=bgam/gam
faktor=(1/beta**2)*(math.log(2*mec2*1e6*beta**2/(IZ*(1-beta**2)))-beta**2)
dEdx=konst*faktor
dE=dEdx*dx
E=E-dE
xstr=str(x)
dEdxstr=str(dEdx)
data=xstr+','+dEdxstr+'\n'
f.write(data)
print x,dEdx
f.close()

Now I want to reproduce Bragg Curve for alpha particle with energy 5.49 MeV in air. Air composition is 80% N2 and 20% O2. I guess I have to change the value of Z, IZ, rho, TZ and TA. I changed Z = 2, TZ=30, TA=60, rho=0.0012, and IZ =16*TZ^0.9, and run the code, but I didn't get the correct curve.

[Bob S]

A 5 MeV alpha particle is very non relativistic, so the Bethe-Bloch equation has to be modified. A good discussion is in the paper

http://www.hep.wisc.edu/~prepost/407/alpha/alpha.pdf

Another problem is that the slow alpha particle occasionally picks up an electron from the air, so the alpha's apparent charge is e rather than 2e. This extends the range a little.

Fig. 3 in the reference is taken verbatim from page 650 in Evans The Atomic Nucleus (1955). It shows the range of a 5 MeV alpha to be about 3.5 cm of air. Evans is worth reading, from pages 637 to 650, if you can find a copy.

Bob S

[ekosulistya]

Thank you very much Bob,

I am new to this topic. I am studying the interaction of radiation with matters. I've got the paper that you suggested, and I'll try to find a copy of the book. I am sorry if I'm asking too much, but if you can give me short explanation on how to get Bragg curve for 5.49 MeV in air by programming, that will help me a lot to understand the stopping process.

Eko S

[Bob S]

A good comprehensive review of charged particle energy loss in matter is given by the LBL Particle Data Group:

http://www.google.com/url?sa=t&source=web&cd=14&ved=0CC8QFjADOAo&url=http%3A%2F%2Fwww.phy.bris.ac.uk%2Fpeople%2Fcus sans_dg%2Fphy

See especially Eqn(27.3) and text. The Bethe Bloch ionization energy loss theory is based on the fact that the major energy loss is collisions with atomic electrons and not with nuclei. Incident charged particles elestically scatter on nuclei with very little energy loss.

Here in attachment is my code for a 5.3 MeV alpha particle (Polonium alpha) in air, including a plot of the Bragg curve. It is surprisingly accurate. Note that the log on line 350 of my code is the natural log and not base 10.

I don't know of a simple discussion of the derivation of these equations.

Bob S

[ekosulistya]

Thank you very much Bob. You've been very helpful. I don't have Truebasic, so I write your code in Excel worksheet. I've tried write your code in Python, but this program also new to me, I can't figure out how to do the condition in the loop part. I attach the picture from Excel worksheet, and part of the worksheet. I can't open the page that you suggested. Anyway, thanks again.

[ekosulistya]

Hi Bob,

Is it okay if use TZ for air with (0.8*14+0.2*16) and TA with (0.8*28+0.2*32)? The graph is't change much. And if I change the alpha energy to 5.49 MeV, the graph shows that the projected range is 3,83 cm.

Eko

[ekosulistya]

Thank you. :-D

Eko S

[naiveBen]

hi .( I`m not used to using english.)

I have a question.

There isn`t any clue on distance with the Beth-Bloch Formular

I mean this formular is independent on the distance except for β
(I hope it makes sense)


so Bob sugested to insert dx

but It seems to me that there is no way we could find out

real amount of dx

it`s a differential distance, isn`t it

but how on earth we just use dx without any real value??
(if i use it with some kind of coding program. there`s no problem?)

so my queation is

what is the value of dx? and when it comes to deciding the value of dx, is there any reasonable reason

[naiveBen]

it`s needless to say that dx is supposed to be really small amount.

but i need some reasonable clue.. let me know why it has to be 0.01 or something..

T^T

[spenmurphy]

Thank you very much Bob. You've been very helpful. I don't have Truebasic, so I write your code in Excel worksheet.

What did this code look like? I'm just learning how to use Excel VBA.

Thanks!

[ekosulistya]

naiveBen,
did you mean that it has to do with the free path length of the ion? I think, dx in the equation is a small step the ion make, whatever the interaction that has happened to the ion, and in any direction. it's computational step, the smaller step we choose, the graph is more smooth.

spenmurphy,
I wasn't using VBA Excel, but worksheet formula only. You can see the snapshot of part of the worksheet in my previous message. I rewritten the code from Bob in worksheet cell, step by step, copy all the way down, and the result is a data set, from which we can choose the value of x and dE/dx, and plot it with Excel scatter chart