C++/MATLAB coding about nuclear decay & detection

Start with an ensemble of N identical radioactive nuclei that decay at a certain point in time. To give you an idea of a possible way to go about this, you could start with an array of integers (or bool) for determining whether a particle has decayed or not. Create a time loop of Nt time steps with step size ∆t. The inverse-mean-lifetime λ of a radioactive sample is basically the probability of decay per unit time, hence for every time step i the probability of decay is pi = λ∆t. For each nucleus use the random number generator to generate a number between 0 and 1 and compare with pi to choose if a nucleus decays or not — remember you will also have to add a condition so that a particular nucleus does not decay more than once!

We will now build on this assignment to develop a virtual detector for the radiation, including angular acceptance.

1. When each nucleus decays it emits Nγ photons with different energies that travel in random directions. A spherical polar coordinate system can be defined with angles θ and φ, where θ is measured from an axis in the coordinate system as in figure 1.

Figure 1: Spherical coordinate system.

1

You should now consider a detector with a circular cross-section so that only some of these photons are incident (and can therefore be detected). To make this problem easier, we can orientate the coordinate system so that the detector is situated on the axis that θ is measured from, which means φ doesn’t matter by symmetry and we can ignore it. Any photons that are in the range 0 ≤ θ ≤ θ0, where θ0 is the collection angle of the detector, will be detected. For a circular detector with a radius rd at a distance d from the point source, the collection angle will be tan θ0 = rd/d.

First, we need to randomly assign angles θ to the emitted photons. To do this, for each photon generate a random number x between -1 and 1, and find the angle via θ = arccos x. Note that this means for a single radioactive decay releasing Nγ photons, there need to be Nγ angles generated. Write a detector function/class that allows only particles with θ < θ0 to be registered. Also include an efficiency � by using a random number to accept only � of the incident photons.

Why do we generate the angle this way and not just generate a random number between 0 and π?

For a detector of radius 10 cm with an efficiency � = 0.3, find the total number of detected photons for a radiation source that emits 2 photons per decay, starting with 100,000 nuclei with λ∆t as in the previous question, detecting for 25 years1 at the following distances from the source:

(a) 2 cm

(b) 10 cm

(c) 50 cm

Does the number of photons detected vary with distance as we might expect?

1The reason we are having to detect for such a long time is the low number of “nuclei” – which is simply due to computational limitations

2

2. Finally, we are going to make a spectrum, assuming some inherent energy detection width in the detector. For each photon that is detected by the detector, we record its energy, but add a random shift in energy that arises because of the detection process. To do this, add to the photon energy an energy given by ∆E multiplied by a random number from a normal distribution. This is not the usual random number taken from a uniform distribution (i.e. all numbers equally probable over an interval), but is a number that has a probability of occurring following the (unnormalized) normal distribution p(x) = e−x

2/2.

To do this, in matlab use the randn function (not rand), and in C/C++ you can use the following function:

// Gaussian random distribution function

double gaussian()

{

double x, y;

do

{

x = 6.0*random()/RAND_MAX-3.0; // +/-6 sigma range

y = exp(-x*x/2.0);

} while (1.0*random()/RAND_MAX > y);

return x;

}

The photon energy should therefore be recorded as Edetected = Eγ + x∆E, where x is the number sampled from the normal distribution.

For N = 105 initial nuclei emitting two photons having energies E1 = 1173 keV and E2 = 1333 keV, and for a detector resolution of ∆E = 20 keV, using � = 0.3, rd = 10 cm, d = 10 cm and detecting for 25 years, make a histogram plot of the energies with 100 bins from 1000 keV to 3000 keV.

3. Finally, we are going to account for the fact that the detector can’t distinguish pho- tons arriving within a short space of time – basically the total energy deposited in the detector gets summed. In your code, for any photons originating from the same “nucleus” and both falling within the acceptance angle θ0, instead of recording the individual energies, record the summation of the two photon energies. As before, keep- ing the same conditions as the last question, make histogram plots of the energies with 100 bins from 1000 keV to 3000 keV with this coincidence detection for the following distances d:

(a) 2 cm

(b) 10 cm

(c) 50 cm

How does the distance of the detector from the source affect the relative size of the coincidence peak?

3