Fractional Counts - The Simulation of Low Probability events
R.L.
Coldwell
Department of Physics, University of Florida, Gainesville, FL 32611
ufboc@ufl.edu
G.
P. Lasche, and A.
Jadczyk
Constellation Technology Corporation, 7887 Bryan Dairy Road, Suite 100,
Largo, FL 33777 GLasche@aol.com lark1@ozline.net
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 simul.zip.
Introduction
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
NAIAC.SP
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.
COMPUTATIONAL
METHOD
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.
FRACTIONAL COUNTS
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.
STANDARD DEVIATION
OF THE SPECTRUM
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
http://www.phys.ufl.edu/~coldwell
CONCLUSION
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.