\section{Survey of physics in secondary schools}

A survey was placed online\footnote{http://brendan.sdf-eu.org/misc/detector\_survey.php} and advertised on the Institue Of Physics PTNC (Teachers) mailing list, Times Educational Supplement forums and a mailing list of local physics teachers.  In all we got $35$ responses (in fact we got $36$ but one was from a Canadian schoolteacher and so was not included in the final results).

	\subsection{Survey outline}
\footnotesize
\begin{verbatimtab}[4]
Questionnaire for schools teaching particle physics 

I am an undergraduate at Bristol University developing an affordable
particle detector for schools. This questionnaire will hopefully give a
better insight into what schools want from a project such as ours.

Specific information collected will not be passed on to commercial
parties, it is purely for our own research.

------------------------------------------------------------------------

Q1.Would you consider purchasing apparatus for the detection of cosmic
rays at a low enough price? 

Q2. What would be the maximum you would be willing to spend on such a
piece of apparatus? 

Q3. Would you use the equipment primarily for investigation or
demonstration purposes? 
 
Q4.Would your establishment consider collborating with others in the
area to form a 'super detector' similar to the NALTA n the USA? 
 
Q5. It would be a great help if we knew what equipment is already
avaiable to your establishment. Out of the following, which does your
establishment have access to?  

	dry ice 
	liquid nitrogen 
	local source of radiation (i.e. radioactive samples) 
	high voltage supply. (Please specify max voltage)
	oscilloscope.  (Please specify: digital, analogue, both)
	a gas trained technician
	laboratory PC
	video camera. (Please specify: digital, analogue, both)
	photographic darkroom
	fume cupboard 
 
Q6. Any further comments or clarifications? 

Finally could you enter your name, the name of your establishment and
the level of education which you teach there.  

Name:
Establishment:
Level of education (i.e. GCSE, A-Level etc.):  
\end{verbatimtab}

\normalsize

\subsection{Survey results}

% PIE CHARTS, HISTOGRAMS ETC

\newpage

\section{Pion and muon decays}

The following decays with associated probabilities are decays for pions and muons, both are typical components of an extended air shower. These decays were taken from http://www.cosmicrays.org/muon-rays2.php

\bigskip

\begin{center}
\begin{tabular}[h]{rclrcl}
$\mu^-$ 	& -- 	& muon 		& $\mu^+$ 	& -- 	& antimuon \\
$e^-$ 		& -- 	& electron 	& $e^+$ 	& -- 	& positron (antielectron) \\
$\pi^-$ 	& -- 	& pion 		& $\pi^+$ 	& -- 	& antipion \\
$\pi^0$ 	& -- 	& neutral pion & $\gamma$ & -- 	& photon \\
$\nu_e$		& -- 	& electron neutrino & $\bar{\mu_e}$ & -- & antielectron neutrino \\
$\nu_{\mu}$	& -- 	& muon neutrino & $\bar{\nu_{\mu}}$ & -- & antimuon netrino 
\end{tabular}
\end{center}

\bigskip


\begin{center}
\begin{tabular}[h]{rcll}
\multicolumn{3}{l}{Particle decay} & Probability (percent)\\
$\pi^+$	& $\rightarrow$	& $\mu^+ + \nu_{\mu}$ 	& $99.98770\pm0.00004$ \\ 
$\pi^+$	& $\rightarrow$	& $\mu^+ + \nu_{\mu} + \gamma$	& $2.00\pm0.25 \times 10^{-4}$ \\ 
$\pi^+$	& $\rightarrow$	& $e^+ + \nu_e$	& $1.23\pm0.004\times10^{-4}$\\ 
$\pi^+$	& $\rightarrow$	& $e^+ + \nu_e + \gamma$	& $1.61\pm0.23\times 10^{-7}$ \\ 
$\pi^+$	& $\rightarrow$	& $e^+ + \nu_e + \pi^0$	& $1.025\pm0.034\times 10^{-8}$ \\ 
$\pi^+$	& $\rightarrow$	& $e^+ + \nu_e + e^+ + e^-$ & $3.2\pm0.5\times 10^{-9}$ \\ 
$\pi^+$	& $\rightarrow$	& $e^+ + \nu_e + \nu + \bar{\nu}$	& $< 5\times 10^{-6}$ \\ 
&&& \\
$\pi^-$	& $\rightarrow$	& $\mu^- + \bar{\nu_{\mu}}$	& $99.98770\pm 0.00004$ \\ 
$\pi^-$	& $\rightarrow$	& $\mu^- + \bar{\nu_{\mu}} + \gamma$	& $2.00\pm 0.25\times 10^{-4}$ \\
$\pi^-$	& $\rightarrow$	& $e^- + \bar{\nu_e}$	& $1.230\pm 0.004\times 10^{-4}$ \\ 
$\pi^-$	& $\rightarrow$	& $e^- + \bar{\nu_e} + \gamma$	& $1.61\pm 0.23\times 10^{-7}$\\ 
$\pi^-$	& $\rightarrow$	& $e^- + \bar{\nu_e} + \pi^0$	& $1.025\pm 0.034\times 10^{-8}$\\ 
$\pi^-$	& $\rightarrow$	& $e^- + \bar{\nu_e} + e^- + e^+$	& $3.2\pm 0.5\times10^{-9}$\\ 
$\pi^-$	& $\rightarrow$	& $e^- + \bar{\nu_e} + \bar{\nu} + \nu$	& $<5\times 10^{-6}$ \\ 
&&& \\
$\pi^0$	& $\rightarrow$	& $\gamma + \gamma$ 	& $98.8$\\ 
&&& \\
$\mu^+$	& $\rightarrow$ & $e^+ + \bar{\nu_{\mu}} + \nu_e$ 	& $\sim 100$ \\
$\mu^+$	& $\rightarrow$ & $e^+ + \bar{\nu_{\mu}} + \nu_e + \gamma$	& $1.4\pm0.4$ \\
$\mu^+$	& $\rightarrow$ & $e^+ + \bar{\nu_{\mu}} + \nu_e + e^+ + e^-$	& $3.4\pm 0.4\times 10^{-5}$ \\
$\mu^-$	& $\rightarrow$ & $e^- + \nu_{\mu} + \bar{\nu_e}$ 	& $\sim100$\\
$\mu^-$	& $\rightarrow$ & $e^- + \nu_{\mu} + \bar{\nu_e} + \gamma$	& $1.4\pm 0.4$\\
$\mu^-$	& $\rightarrow$ & $e^- + \nu_{\mu} + \bar{\nu_e} + e^- + e^+$	& $3.4\pm 0.4\times 10^{-5}$
\end{tabular}
\end{center}
\newpage

\section{Example readout electronics}
\label{readout-electronics-appendix}

The following was found at http://www.cosmicrays.org/ and although is designed for Geiger Muller tubes, could easily be adapted for our purposes. It also includes a concidence trigger so two chambers could be used to reduce noise or a double gap chamber could be used.
\begin{figure*}[h]
\begin{center}
\includegraphics[scale=0.45]{../pics/readout-electronics}
%\caption{}
\label{readout-electronics}
\end{center}
\end{figure*}

\newpage

\section{Computer modelling using python}
\label{computer-modelling}

I chose to use Python\footnote{See http://www.python.org/} to model the detector since it is easy to learn and highly reusable. The module \verb"ParticleModelClasses.py" can be used again by future students to quickly and easily write further models for the detector. There is a brief guide to the module in Appendix \ref{python:reusing}

All the computers in the main computer laboratory are currently (2005) equipped with ActiveState's Python interpreter.

	\subsection{Monte Carlo detector model for determining coincidence rates}
\footnotesize
\begin{verbatimtab}[4]
# coincidence.py
# ----------------------------------------------------------------------
# A Program for determining expected coincidence rates for different
# arrangements of RPC and scintillator paddles

# IMPORT CLASSES

# See seperate file ParticleModelClasses.py for classes imported in the
# following statement
from ParticleModelClasses import *

# MODEL VARIABLES

# Takes 'num_sets' number of sets of data. each set has 'events_per_set'
# number of particles run through the model. Take the mean of each set,
# then takes the mean of the set means as the output
num_sets = 30 # Gives generally consistent results
events_per_set = 30

# Define the limits of the model. Should enclose all detectors but also
# be as small as possible to improve speed.
model_limits = Box([600,300,397.5],[0,0,163.75])

# Create the detectors in a list for easy access and expandability
dets = {}
dets['coincidence1'] = Detector([110,100,5],[0,0,360]) # Small paddle
dets['rpc']   = Detector([100,100,2],[0,0,100])        # RPC
dets['coincidence2'] = Detector([600,300,70],[0,0,0])  # Base scintillator

# MAIN ROUTINE

# Creates a particle inside first coincidence detector and passes it
# through the array of detectors logging hits every mm if it is inside a
# detector

# Initiate some vars to store data from each set
ratios = []

# Iterate over each set
for set in range(num_sets):
	# Iterate over each particle
	for id in range(events_per_set):
		# Create a particle inside first paddle
		muon = CosmicRay(dets['coincidence1'])
		# Create a place to store all hits from this particle in each
		# detector
		for det in dets.itervalues():
			det.hits[id] = 0
		# Trace particle path through detector group until it leaves
		while model_limits.inside(muon.pos):
			# Check to see if particle is currently inside a detector
			for det in  dets.itervalues():
				# Since detectors tend to be narrow in z direction,
				# check
				# this first to avoid unecessary function calls
				# (improves
				# speed)
				if muon.pos[2] < det.nr[2] and muon.pos[2] > det.far[2]:
					# Check other co-ords
					if det.inside(muon.pos):
						# If particle is inside of detector then 
						# register a hit
						det.hits[id] = det.hits[id] + 1 
			# Update particle position for next iteration
			muon.pos = muon.pos + muon.vel
			
	# Extract data stored in our detector instances
	rpc_hits = 0
	scrn_captures = 0
	for c1, rpc, c2 in zip(dets['coincidence1'].hits.itervalues(), \
						dets['rpc'].hits.itervalues(), \
						dets['coincidence2'].hits.itervalues()):
		if c1 and c2:
			scrn_captures = scrn_captures + 1
			if rpc:
				rpc_hits = rpc_hits + 1

	# Place data for this set in array
	if scrn_captures != 0:	
		ratios.append(1.0 * rpc_hits/scrn_captures)
	

# Print out the data
print 'Calculated from ',num_sets,' sets of ',events_per_set,' ionising cosmic rays:'
print '\t Can expect ', round(100.0* sum(ratios)/len(ratios)) , \
			'% of triggers to have a signal from RPC'
\end{verbatimtab}
\normalsize

	\subsection{Monte Carlo detector model for determining average path lengths}
\footnotesize
\begin{verbatimtab}[4]
# path_length.py
# ----------------------------------------------------------------------
# A program which determines typical path lengths for particles in our
# RPC

# IMPORT CLASSES

# See seperate file ParticleModelClasses.py for classes imported in the
# following statement
from ParticleModelClasses import *
from dump import *

# MODEL VARIABLES

# Set number of event want to model for the histogram
num_events = 5000

# Set bin increments for histogram (mm)
bin_inc = 0.02

# RPC dimensions x, y, z (mm)
rpc_dims = [180,200,2]

# Reolution of model in mm
resolution = 0.02

# MAIN ROUTINE

# Define model limits so it encloses RPC by least 1 mm on all sides.
# This is because CosmicRay traceback mathod places particle up to 1mm
# inside of model
model_limits = Box([l+2 for l in rpc_dims], [0,0,0])

# Create RPC
rpc = Detector(rpc_dims, [0,0,0])

for id in range(num_events):
	# Create particle
	muon = CosmicRay(model_limits)
	# Create a place to store path length from this particle in RPC
	rpc.hits[id] = 0
	# Trace particle path through model until it leaves
	while model_limits.inside(muon.pos):
		# Check to see if particle is inside RPC
		if rpc.inside(muon.pos):
			# If is inside add 'resolution' mm to path length in RPC
			rpc.hits[id] = rpc.hits[id] + resolution
		# Update particle posistion for next iteration
		muon.pos = muon.pos + resolution * muon.vel

# PROCESS DATA

# Compiles frequency data for a histogram for use in Excel

# Get lower and upper limits for the histogram
min_val = min(rpc.hits.values())
max_val = max(rpc.hits.values())

# Initialise array to store frequency data
bins = {}
# Not necessary, but makes program easier to read
curr_bin_boundary = min_val
# Iterates over all bins in histogram
while curr_bin_boundary < max_val:
	# Gets a list of values which lie between 'curr_bin_boundary' and
	# 'curr_bin_boundary + bin_inc', then stores the length of this in
	# list
	bins[curr_bin_boundary] = len([x for x in rpc.hits.values() if x >= \
		curr_bin_boundary and x < curr_bin_boundary + bin_inc])
	# Move the bin boundary to the next histogram bin
	curr_bin_boundary = curr_bin_boundary + bin_inc

# OUTPUT

# Outputs a tab seperated list of pair values for use in Excel. Use UNIX
# redirect (i.e. 'path_length.py > file.dat') to output to a file.

# Have to use this cludgey code for sorting a dictionary
bin_vals = bins.keys()
bin_vals.sort()
for bin_val in bin_vals:
	print bin_val, '\t', bins[bin_val]
\end{verbatimtab}
\normalsize

	\subsection{Common detector classes}
\footnotesize
\begin{verbatimtab}[4]
# ParticleModelClasses.py
# ---------------------------------------------------------------------
# Defines the following classes, useful for creating models of particle
# detectors. See each classes docstrign for more information.
# 
# Vector    - A list of numbers that behave mathematically as a vector
# Box       - A cuboid in the model, Detector inherits from this
# Detector  - A cuboid detector in the model
# Particle  - A particle with velocity and position variables
# CosmicRay - A particle created in a random position in a space with
#             velocity according to the azimuthal Cos^2 distribution 

# IMPORT SOME NECESSARY FUNCTIONS

from math import *
from random import random

# CLASSES

class Vector(list):
	"""Mathematical vector of any size

	Changes behaviour of lists so as to behave more like a vector w.r.t
	addition etc. Methods are described seperately. Builds on code by A.
	Pletzer in the ActiveState Python Cookbook"""
	def __add__(self, other):
		"Returns an element by element sum in a Vector"
		return Vector([i+j for i,j in zip(self, other)])
	def __neg__(self):
		"Returns Vector with all elements negatived"
		return Vector([-i for i in self])
	def __sub__(self, other):
		"Returns Vector of the difference of two vectors"
		return Vector([i-j for i, j in zip(self,other)])
	def __mul__(self, other):
		"Returns element by element multiplication in a Vector"
		try: return Vector([i*j for i,j in zip(self,other)])
		# Must be a scalar
		except: return Vector([i*other for i in self])
	def __rmul__(self, other):
		"Since multiplcation commutes, __rmul__() does as __mul__()"
		return self*other
	def __div__(self, other):
		"Divides, element by element"
		try: return Vector([i/j for i,j in zip(self, other)])
		# Must be scalar
		except: return Vector([i/other for i in self])
	def magnitude(self):
		"Returns magnitude of the vector"
		return sqrt(sum([x**2 for x in self]))
	def unit(self):
		"Returns unit vector in direction" 
		return Vector(self / self.magnitude)
	def dot(self, other):
		"Returns dot product of this and another vector"
		return sum(self * other)
	def cross(self, other):
		"""Returns the cross product of this and another vector.
		 Currently only works for Vectors with 3 values"""
		if len(other) != 3 or len(self) != 3:
			raise ValueError, 'cross method currently only works for 3 valued Vectors'
		else:
			return Vector(self[1]*other[2]-self[2]*other[1], \
					self[0]*other[2]-self[2]*other[0], \
					self[0]*other[1]-self[1]*other[0])
		

class Box:
	"""Defines a cuboid in 3D

	First argument when a Box is created is a 3 valued vector giving the
	dimensions of the cuboid, the second argument gives the position of the
	centre. If no arguments given a unit box centred on origin is created"""
	def __init__(self, dims=[1,1,1], pos=[0,0,0]):
		# Note: Do not change the dimensions and position directly,
		# since the the corner positions will not be updated and vise
		# versa. Use the 'set' methods below 
		self.set_dims_pos(dims,pos)
	def set_corners(self, nr, far):
		"Use to redefine corners of box"
		self.nr = Vector(nr)
		self.far = Vector(far)
		self.dims = Vector([n-f for n,f in zip(nr,far)])
		self.pos = Vector([0.5*d+f for d,f in zip(self.dims,far)])
	def set_dims_pos(self, dims, pos):
		"Use to redefine position and dimensions of box"
		self.dims = Vector(dims)
		self.pos = Vector(pos)
		self.nr = Vector([0.5*d+p for d,p in zip(dims,pos)])
		self.far = Vector([-0.5*d+p for d,p in zip(dims,pos)])
	def inside(self, pnt):
		"Returns True if a point is inside box, False otherwise"
		if pnt[0] > self.nr[0] or pnt[0] < self.far[0]: return False
		if pnt[1] > self.nr[1] or pnt[1] < self.far[1]: return False
		if pnt[2] > self.nr[2] or pnt[2] < self.far[2]: return False
		return True

class Detector(Box):
	"""Defines a cuboid shaped detector in 3D

	Inheriting from the Box class, the dimensions and position are
	defined as for a Box."""
	def __init__(self, dims, pos):
		# Run Box constructor
		Box.__init__(self, dims, pos)
		# Define an array to store all particle hits
		self.hits = {}

class Particle:
	"""Defines a particle

	Initialises a particle with a position and velocity vector. If no
	arguments given a stationary particle at the origin is created"""
	def __init__(self, pos=[0,0,0], vel=[0,0,0]):
		self.pos = Vector(pos)
		self.vel = Vector(vel)

class CosmicRay(Particle):
	"""Defines a cosmic ray particle

	 Creates a particle in a cuboid with random direction and
	displacement and traces back to point of entry. Parameter vol is
	instance of Box or Detector class. If no argument is passed, CosmicRay
	is created at origin"""
	def __init__(self, vol=Box([0,0,0],[0,0,0])):
		# Azimuthal angle according to Cos^2 probability distribution.
		# This uses a Von Neumann acceptance, rejection algorithm, see
		# Section 33.3 Phys. Rev. 2004
		while 1:
			self.phi = random() * pi / 2
			y = random()
			if y <= cos(self.phi)**2: break
		self.theta = random() * 2.0 * pi # Rotational angle of particle
		# Calculate unit velocity vector
		self.vel = Vector([ sin(self.phi) * sin(self.theta), \
				sin(self.phi) * cos(self.theta), \
				-cos(self.phi) ])
		# Place particle in random place in vol
		self.pos = Vector([ random() * vol.dims[0] + vol.far[0], \
				random() * vol.dims[1] + vol.far[1], \
				random() * vol.dims[2] + vol.far[2] ])
		# Traceback the particle position until it is just inside of the
		# vol
		self.traceback(vol)
	# Define traceback function which places particle at entry point to
	# vol
	def traceback(self, vol):
		# Start by skipping 10 mm at a time
		while 1:
			self.pos = self.pos - 10 * self.vel
			if not vol.inside(self.pos): break
		# Place particle back inside detectors vol
		while not vol.inside(self.pos):
			self.pos = self.pos + self.vel
\end{verbatimtab}
\normalsize

	\subsection{Probability program to predict noise levels}
\footnotesize
\begin{verbatimtab}[4]
# noise_coincidence.py
# ----------------------------------------------------------------------
# A program  for determining how many noise blips are needed to get a
# specified coincidence ratio with a signal of a specified amount.

from random import *
from math import *

# SET VARIABLES

# Signal blip length (S)
sig_dur = 60E-9

# Noise blip duration (S)
noise_dur = 10E-9

# Approx number of signals per second
sig_freq = 1.5

# Coincidence ratio (ratio between signals:coincidences)
coinc_ratio = input("Enter coincidence ratio: ")

# START MAIN ROUTINE

# We calculate the cumulative probabilities, so we can say 'if we are
# getting 50% coincidences, this could be due to a noise level of N
# blips per second' The maths is as follows,

# P1(cumul) = P1
# P2(cumul) = !P1*P2 + P1
# P3(cumul) = !P1*!P2*P3 + !p1*P2 + P1
# P4(cumul) = !P1*!P2*!P3*P4 + !P1*!P2*P3 + !P1*P2 + P1
# etc.

# Where Pn is the probability of the nth blip coinciding, !Pn is the
# probability of the nth blip NOT coinciding.

# We can use an iterative method (if you look Pn is P(n-1) with an extra
# term) and stop when Pn(cumul) is more than the coincidence ratio.

# ALGORITHM

# Initialise with probability of one blip coinciding 
n = 1
Pcumul = (noise_dur+sig_dur)*sig_freq

# We see that the new term in each iteration is comprised of all the !Pn
# from the previous iteration with the extra !P(n-1)*Pn factor. We
# therefore store the series up to !P(n-2) in an array to avoid
# recalculating. This is initialised as 1
Pnot_series = 1

# Keep iterating until we break past the ratio. n will then be set to
# the number we want
while Pcumul < coinc_ratio:
	# Bump up the amount of noise blips 1 at a time
	n = n + 1
	# Calculate !P(n-1)
	Pprev_not = 1 - ((n-1)*noise_dur+sig_dur)*sig_freq
	# Calculate prob of nth particle coinciding
	Pn = (n*noise_dur+sig_dur)*sig_freq
	# Extend array
	Pnot_series = Pnot_series * Pprev_not
	# Calculate current cumulative probability
	Pcumul = Pnot_series * Pn + Pcumul

# OUTPUT

print "Signal blip size: \t", sig_dur
print "Signal frequency: \t", sig_freq
print "Noise blip size:  \t", noise_dur
print "Coincidence ratio:\t", coinc_ratio
print "Number of noise blips:\t", n
\end{verbatimtab}
\normalsize

	\subsection{Re-using the detector classes}
	\label{python:reusing}

In order to use the classes in the module \verb"ParticleModelCLasses.py", you will need to import them into your script. Provided a copy of \verb"ParticleModelClasses.py" is in the same directory as your script the following line should work,
\footnotesize
\begin{verbatimtab}[4]
from ParticleModelClasses import *
\end{verbatimtab}
\normalsize
After this has been done you now have access to all the objects defined in the file. The first object you will probably want to define is the box which contains your model.  To do this you assign it to a variable to create a new instance, along with two lists of numbers, the $x$,$y$,$z$ dimensions of the box and the $x$,$y$,$z$ co-ordinates of the centre of the model. So for example,
\footnotesize
\begin{verbatimtab}[4]
model_limits = Box([100,100,100],[0,0,0])
\end{verbatimtab}
\normalsize
 After this you will probably want to define a \verb"Detector" in a similar way,
\footnotesize
\begin{verbatimtab}[4]
rpc = Detector([100,100,5],[0,0,0])
\end{verbatimtab}
\normalsize
Now you can create a \verb"Particle" by specifying $x$,$y$,$z$ position co-ordinates and $v_x$,$v_y$,$v_z$ velocity components in two lists,
\footnotesize
\begin{verbatimtab}[4]
electron = Particle([0,0,3],[0,0,-1])
\end{verbatimtab}
\normalsize
This creates a particle just above the `rpc' with a downwards speed of `$1$'. To move it we do the following,
\footnotesize
\begin{verbatimtab}[4]
for t in range(10):
    electron.pos = electron.pos + electron.vel
\end{verbatimtab}
\normalsize
which tracks the electron for $10$ seconds along its velocity vector. You will want to include some more interesting code in this loop before it becomes useful.

Python allows you to easily add your own methods and build upon the existing classes using inheritance to suit your needs, for example, say you want to model a particle with spin, you would simply inherit the \verb"Particle" class and add some extra attributes,
\footnotesize
\begin{verbatimtab}[4]
class ParticleWithSpin(Particle):
	def __init__(self, pos, vel):
		# Initiate as for Particle
		Particle.__init__(pos, vel)
		# Randomly chose up or down (assume random() has been imported)
		if random() > 0.5:
			self.spin = 1
		else:
			self.spin = -1
\end{verbatimtab}
\normalsize
Current attributes cover the needs of the models that I wrote and are documented in the code.

