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

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.

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.

In this repo you will find an example of the simulation of ^{1}H, ^{6}He, ^{19}C, ^{27}S and ^{56}Ni ions, entering a 1x1x20 m^{3} cube filled with isobutane gas (2.67 mg/cm3 in STP conditions, composition C_{4}H_{10}), 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.

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

`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;
}
```

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:

${\overline{)\frac{\Delta H}{\Delta 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.

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; k
```SetBinContent( 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);
}

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\left(x\right)={a}_{0}+{a}_{1}x+{a}_{2}{x}^{2}+{a}_{3}{x}^{3},$

$f\text{'}\text{'}\left(x\right)=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) );
}
```

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