Bragg Curves

This is part of a simulation of an active target with the Geant4 libraries.

The code

Run this simulation in Geant4 to generate an output file that contains information about the interactions, positions, energy, energy loss per step, etc. of a particle entering a gas detector. The output file is in ROOT format.

Detector box geometry Detector tube geometry The figures above show examples of the cubic and tube detector geometries.

The macros in the bragg folder allow you to open this ROOT simulation file, read and manipulate its contents, fit data, etc. The aim of these snippets is to calculate the range of certain ions in a gaseous medium. One of the possible ways to do that is by calculating the inflection point of the Bragg curve.

Simulation

In this repo you will find an example of the simulation of ¹H, ⁶He, ¹⁹C, ²⁷S and ⁵⁶Ni ions, entering a 1x1x20 m³ cube filled with isobutane gas (2.67 mg/cm3 in STP conditions, composition C₄H₁₀), with energies from 1 to 150 MeV/u, and with the electric and magnetic fields turned off. You can also set the number of runs to be launched for each energy.

Tracks inside the detector.

The files

`range.C` (and `range.h`)

You can manipulate the simulation file with range.C, which has to be run in a ROOT prompt. This program opens the simulation file, reads all the events in it, fills in the relevant histograms and saves them into the simulation file. Then it performs a fit and writes the fit parameters to an output file.

$ root -l           // -l avoids printing the header
$root[] .L range.C  // .L Loads the macro
$root[] range t     // Creates an object
$root[] t.Loop();   // Loop on all entries
$root[] .q          // Quit

The object of interest inside the simulation file is stepSumLengthOnGasPrim1, which saves the range of the particle in the material for every run (i.e., the sum of the length of all the steps until the particle is stopped in the gas). The values of the range will slightly vary from one run to the other, but they will tend to fall around a value (so, to have decent statistics, launch around 1000 runs for each energy in the simulation). The program builds a histogram with the values of the range and fits it to a Gaussian curve:

TF1 *f1 = new TF1("f1","gaus");
Double_t par[3];
Double_t *parErr;
parErr = new Double_t[3];

hStepSumLengthOnGas1->Fit("f1");
f1->GetParameters(par);
parErr = f1->GetParErrors();

par[3] saves the mean value of the range, the sigma and the chi squared of the fit, and parErr saves the errors. In the output file (range.dat) we save the energy, the range and the sigma with their corresponding errors.

Gaussian fit of the range for an energy The figure shows the Gaussian fit for a Carbon ion entering isobutane. The fluctuations are due to the nature of the low energy interactions at the end of the track.

`readrange.C`

Repeat the simulation for different values of energy to have a collection of pairs of energy and range. Inside ROOT, run readrange.C (it's in the bragg/fits directory) to read this collection and fit the resulting curve:

$ root -l
$root: .x readrange.C

In this case, the macro doesn't use a ROOT fit function, but a user defined function, myfit.C:

$R = a E^{b}$

where R is the range, E is the energy, a gives information on the particle and b gives information on the medium.

Double_t myfit(Double_t *x, Double_t *par) {
  // Trying to fit to R = a * E^b
  // x      = Array of energies (E)
  // par[0] = Fist parameter    (a)
  // par[1] = Second parameter  (b)
  // fitval = Range             (R)

   Double_t fitval = par[0]*pow(x[0],par[1]);
   return fitval;
}

Fit of the range as a function of the energy The figure shows the range dependency on the energy for a proton stopped in isobutane.

These results can be easily compared with tables and calculations from SRIM, ATIMA, etc.

`bragg.C` (and `bragg.h`)

To compare the simulations with experimental data, it is easier to work with the Bragg peak, because it is easy to identify it in the data (however, the end of the track is affected by fluctuations and the detection threshold). From the tracks of the simulation file, it is possible to calculate their corresponding Bragg curves.

The range can be obtained from the Bragg curve by searching the inflection point in the region at the end of the curve. The inflection point is the point where the curvature changes from concave to convex. At such a point, the second derivative is zero and changes sign from there on.

In bragg.C, the relevant objects are loaded and the histograms built. Since a histogram is not a continuous curve but a discrete set of bins, a "discrete" derivative has to be calculated:

${\frac{Δ H}{Δ z}}_{i} = \frac{H_{i + 1} - H_{i}}{z_{i + 1} - z_{i}}$

where H_i is the content of the bin z_i of the histogram, and z_i+1 - z_i is the size of the pads in the active wall of the detector (the part of the detector that actually "detects"). The pad size is 1.5 mm.

Projection of a track over the pads When a particle enters the detector, its interaction with the gas produces tracks that are detected by the pads in the active wall.

In the code, the braggcurve histogram is built, and then the derivative histogram is manually built from it, as explained above. The second derivative is calculated in the same way and stored in the derivative2 histogram:


// MANUAL CALCULATION OF THE DERIVATIVES OF THE BRAGG CURVE
derivative->SetBinContent( 1,(braggCurve->GetBinContent(2) - braggCurve->GetBinContent(1))/pad);
derivative->SetBinContent( 2,(braggCurve->GetBinContent(3) - braggCurve->GetBinContent(2))/pad);
derivative2->SetBinContent(1,(derivative->GetBinContent(2) - derivative->GetBinContent(1))/pad);

for(Int_t k=1; kSetBinContent( k+2,(braggCurve->GetBinContent(k+3) - braggCurve->GetBinContent(k+2))/pad);
derivative2->SetBinContent(k+1,(derivative->GetBinContent(k+2) - derivative->GetBinContent(k+1))/pad);
derivative3->SetBinContent( k, (derivative2->GetBinContent(k+1)- derivative2->GetBinContent(k)) /pad);
}

Bragg curves and derivative First panel: Bragg curve. Second panel: First derivative. Third panel: Second derivative. Fourth panel: All together.

Then the program fits the last part of the Bragg curve to a third order function. Imposing the condition that the second derivative has to be zero, the inflection point x can be obtained from the parameters of the fit.:

$f (x) = a_{0} + a_{1} x + a_{2} x^{2} + a_{3} x^{3},$
$f'' (x) = 0 \Rightarrow x = \frac{a_{2}}{3 a_{3}}$

The program first performs a smoothing of the histogram, and then takes the last 5 pads (bins) to make the fit. Then, it writes the fit parameters to fit.dat.

TF1 *myfit = 0;
myfit = new TF1("myfit", "pol3", minPos, maxPos); // Fit to a function in the interval [minPos, maxPos]
h->Smooth(50, -1, -1); // Smooth the entire histo 50 times before fitting
h->Fit("myfit", "RQ");

myfit->GetParameters(par);
parErr = myfit->GetParErrors();

if (par[3] != 0 && par[2] != 0) {
	inflexionPoints    = - ( par[2] / (3*par[3]) );
	inflexionPointsErr = (inflexionPoints / 3)*sqrt( pow(parErr[2] / par[2], 2) + pow(parErr[3] / par[3], 2) );
}

Third order fit Fit to a third order function to calculate the inflection point.

This is done for each run, and the statistics of the inflection point for that ion / gas are saved in another histogram. Then, it performs a Gaussian fit to this histogram, and the results of the fit with their errors are saved in inflectionP.dat.

TF1 *gauss = new TF1("gauss", "gaus");
Double_t param[3];
Double_t *paramErr;
paramErr = new Double_t[3];

inflectionPoint->Fit("gauss");
gauss->GetParameters(param);
paramErr = gauss->GetParErrors();

The inflection point of the Bragg curve represents the range of the particle.

The inflection point is what can be easily obtained with an emperiment. The range is what can be easily obtained with a simulation. It would be ideal to find a relationship between the two.

Published: 2007-02-05 21:07:34