Fractional Counts - The Simulation of Low Probability events

R.L. Coldwell

Department of Physics, University of Florida, Gainesville, FL 32611

G. P. Lasche, and A. Jadczyk

Constellation Technology Corporation, 7887 Bryan Dairy Road, Suite 100, Largo, FL 33777


Abstract. The code RobSim has been added to RobWin1. It simulates spectra resulting from gamma rays striking an array of detectors made up of different components. These are frequently used set coincidence and anti-coincidence windows that decide if individual events are part of the signal. The first problem is the construction of the detector. Then owing to statistical nature of the responses of these elements there is a random nature in the response that can be taken into account by including fractional counts in the output spectrum somewhat complicating the error analysis as Poisson statistics are no longer applicable. The “dos” version of the code presented here along with a very brief set of directions in the readme.txt is in



The code began life as BSIMUL. It was built to aid in the construction of the BGO shielded GRAD2 detector that was eventually used in an Antarctic balloon borne observation of supernova 1987A3.

FIGURE 1. Sample detector similar to that used by the NEAR-X and Gamma Ray group at NASA Goddard.


The simulation works by assigning materials and detector properties to a series of zones. There are five zones in figure 1, which is the detector as drawn on the screen. Output spectra are assigned by almost a computer language. After an input photon has been completely absorbed or has left the detector, the energy absorbed in zone n is assigned the name En. For the spectrum consisting of the energy detected in NaI, zone 3, the direction lines are the name of the spectrum on one line followed by E3 always on the next.

The anti-coincidence spectrum direction lines are a bit more complex

E3 if not (10 < E1+E2 < 10000) or (50 < E4 < 10500)

The first line is the name of the file in which the spectrum will be placed. The second line tells the code to add the probability of energy absorbed in E3 to the spectrum only if the energies in zones 1, 2 and 4 are "essentially" zero. The code spreads the energies according to the input full width at half-maximum and input photo-peak response function. The above set of directions with 50000 input photons produce the results shown in figure 2 in 2-5 minutes on a laptop computer with a 366 MHz Celeron processor.



FIGURE 2. Input photon is 2614 KeV. Top spectrum is NaI alone. Bottom spectrum is NaI in anti-coincidence with any event in plastic or BGO.


The code follows photons through matter. These photons have 3 possible interactions. They can be completely absorbed in a single encounter, they can Compton scatter and they can engage in pair production. The experimental cross sections for these three reactions are stored for most materials. In addition it is possible to use information about the density and charge dependence to deduce the relationship of non-stored cross sections to stored ones4 -- this last has been implemented only for HgI2 at this time. Electrons are not followed they are absorbed with the photons.



The detector in figure 1 is defined by a set of lines in rho and z. A subroutine returns the material present at any value of (rho,z). The basic method is to calculate the distance to the nearest defining line in rho and in z. If the material on both sides of the line is the same, the distance is extended to the next line. Stored experimental cross sections are then used to find the probability of interaction. A random number generator5 determines if, where, and how the photon interacts. If the interaction is photo absorption, the photon tracking stops. If the interaction is a Compton scattering, its direction is changed randomly in accordance with the Klein-Nishina formula6. Pair production emits both a pair of 511 KeV photons and a series of bremsstrahlung7 photons. The raw and anti-coincidence spectra are given in figure 2 for 2614 KeV input photons. The same spectra for 6 MeV photons are shown in figure 3.

FIGURE 3. Input photon is 6000 KeV. Top spectrum is NaI alone. Bottom spectrum is NaI in anti-coincidence with any event in plastic or BGO.

The largest peak in this spectrum is the single escape peak in which all the incoming energy except that in one of the 511 KeV gamma rays produced by a pair production event is absorbed in the NaI in zone 3. The anti-coincidence spectrum is still a very clean representation of the case in which all of the energy is absorbed in the NaI, but in the presence of background from other sources, it is not as easily found as is the first escape peak.


The full width at half-maximum response of each detector is used to calculate the probability of the energy absorbed in the combination of zones describing a detector being reported as occurring in any one of the channels ranging over the full photo-response of the detector. The sums of these probabilities over many input photons are the spectra plotted in figures 2 and 3. This could have been accomplished, however, by reporting the energy in the single channel corresponding to it, followed by a convolution with the appropriate photo-response function. The earlier versions of the code did exactly this, although the convolution step was frequently not made.

Fractional counts begin to be worthwhile in simulating the attempts made to suppress the photo-peak in figure 3. The following method for isolating peaks from high-energy lines is due to the NASA Goddard X and Gamma- Ray group8. The lines defining the spectrum are

NaiCoBgo.sp 2 //NaI in coincidence with 511 or 1022 MeV in BGO //in anti-coincidence with a plastic front shield

//Coin window for BGO is 491 KeV to 531 KeV
//and from 1002 KeV to 1042

E3 if (491 < E1 + E2 < 531) or (1002 < E1 + E2 < 1042) not (10 < E4 < 10000)


FIGURE 4. Input photon is 6000 KeV. Spectrum is NaI in coincidence with a 511 or 1022 KeV photon in the BGO and nothing detected in the plastic.

The spectrum reported in figure 4 is


The k refers to the type of spectrum or equivalently to the set of zones in which energy E has been absorbed. The subscript i refers to energy of channel i. The probability P includes the spread due to the photo-response of the detector and also that in the other zones. That is the probability is calculated that the event would be vetoed or accepted for every event. An energy higher than the window around the 511 could be accepted owing to the fact that the detector in charge of vetoing responds as though it were in the proper range or it could be vetoed even though it is in the correct range. This is a straightforward calculation from the photo-response functions in each zone on an event by event basis, but it is not one that would be easy to make after the events are summed. If there were an energy for which a series of unlikely events need to conspire such as two or more detectors responding at the edge of their usual response functions, this event will be seen as a small value of Pk for many photons in equation 1. It is of course possible to randomly select a single channel in which to report a 1. In this case it will take 10000 photons to see a single channel with an average probability of 10-4.

The photo-response function of the BGO in figure 4 is that of the NaI, which is approximately that of either peak in figure 4. Actually BGO is usually not that energy sensitive. Quadrupling the width of the peak for the BGO leaves "essentially" the same spectrum but divides the strength of the lines by about 4. This is plotted on a logarithmic scale in figure 5.

FIGURE 5. Input photon is 6000 KeV. Spectrum is NaI in coincidence with a 511 or 1022 KeV photon in the BGO -- wide response function -- and nothing detected in the plastic.

The vertical lines in figure 5 are the error estimates in the determination of the spectrum. They are the error of the entire curve and not that of the individual points.




Include zero as a perfectly acceptable probability. Define Pj to be the probability reported for channel i of the k th spectrum after the j the photon has been stopped or exited the system. Then denote the average value of the i'th channel in the k'th spectrum as <P>. The standard deviation estimate begins with the definition


Where j is a sum over the N photons tracked by RobSim.

The variance in this value is defined by

Owing to the fact that the P's are totally uncorrelated, it is assumed that the off diagonal terms average to zero so that


This means that

The standard deviation in <P> is the square root of this quantity divided by the number of observations N


The spectrum of course is N times <P>. Thus the standard deviation in the spectrum is

To see how the usual Poisson statistics come from this, imagine that there is a channel that receives a single count M<<N times. In this case


Since by assumption M/N << 1, the second term is much less than the first, this becomes


The vertical lines at each point figure 5 extend from one standard deviation below the value of the spectrum to one above. A second simulation of the spectrum made with different random numbers has a 65% probability of finding the j th channel within this vertical line. The adjacent points are highly correlated with this point. If the j th spectral value is above the line in the second calculation, all values with approximately the same j will also be above the line and vice versa. If the standard deviations were multiplied by the square root of the full-width at half maximum, then the visual estimate based on adjacent points would seem better, but there would be other problems. FSPDISD


The code RobSim estimates a detector response function . Typically this is calculated by repeatedly allowing a photon of energy E' to enter the system of zones. The code allows a user to imagine detectors connected to each of a number of various zones and to use these to find the "best" possible response function. An attempt has been made to make this process as fast and intuitive as possible.


1 R.L. Coldwell, "Robust Fitting of Spectra to Splines with Variable Knots", CP475, Application of Accelerators in Research and Industry, edited by J.L. Duggan and I.L Morgan AIP Press, New York (1999)

2 A.C. Rester, R.B. Piercey, G. Eichhorn, R.L. Coldwell, J. McKisson, D.W. Ely, H.M. Mann and D.A. Jenkins, "The GRAD Gamma-Ray Spectrometer", IEEE Trans. on Nucl. Sci. 33, 732 (1986).

3 A.C. Rester, R.L. Coldwell, F.E. Dunnam, G. Eichhorn, J.I. Trombka, R. Starr and G.P. Lasche, "Gamma-Ray Observations on Supernova 1987A from Antarctica", Ap. J. 342, L71 (1989).

4 Robley D. Evans, The Atomic Nucleus, McGraw Hill, 1955 pp.686-715

5 R.L. Coldwell, "Correlational Defects in the Standard IBM 360 Random Number Generator and the Classical Ideal Gas Correlation Function",, J. Computational Phys. 14, 223(1974)

6 E.U. Condon, "X Rays",in Handbook of Physics, ed E.U. Condon and Hugh Odishaw, McGraw-Hill (1958)pp 7-128-7-129

7 R.L. Coldwell, M.W. Katoot, and P.S. Haskins, "Interactions of Multi-Mev Gamma Rays with Matter", pp. 125-131,High Energy Radiation Background in Space, AIP Conference Proceedings #186, Edited by A.C. Rester and J.I. Trombka (AIP, New York, 1989)

8 L.G. Evans, R. Starr, J.I. Trombka, T. McClanahan, S.H. Bailey, I. Mikheeva, J. Bhangoo, J. Bruckner, and J.O. Goldsten, "Calibration of the NEAR Gamma-Ray Spectrometer", Icarus (in press)