# -*- coding: utf-8 -*-
"""
Created on Sat Oct  3 11:44:28 2020

@author: jbobowsk
"""


# This tutorial will use Monte Carlo methods to simulate the output of a
# Photomultiplier Tube (PMT) once triggered by the detection of one or more
# photons.

# This example follows notes posted online which you can find at the following
# url: http://superk.physics.sunysb.edu/~mcgrew/phy310/lectures/phy310-lecture-06-2007.pdf
import numpy as np
import matplotlib.pyplot as plt
import statistics

# In all Monte Carlo simulations it is necessary to generate random or
# pseudo-random numbers.  The Python command 'np.random.uniform()'
# from the NumPy module generates random numbers uniformly distributed 
# between zero and one. 
print(np.random.uniform())

# First, suppose N photons hit the photocathode.  Determine the number of
# photoelectrons (pe) that are generated.  Assume that the incoming light is
# 400 nm and that the quantum efficiency of the photocathode is 0.23.
#(Following  http://superk.physics.sunysb.edu/~mcgrew/phy310/lectures/phy310-lecture-06-2007.pdf)
N = 2363
QE = 0.23
pe = 0
for i in range(N):
    if np.random.uniform() < QE:
        pe += 1
print(pe)

# When an electron with energy E hits a Dynode the average number of
# secondary electrons liberated is alpha*E (i.e. number of secondary electrons
# is proportional to the energy of the incoming electron) where alpha is a
# contant.  Since we are counting electrons, the distribution of liberated
# electrons follows the Poisson distribution.  In Python random numbers
# drawn from a Poisson distribution with mean mu are generated using
# poissrnd(mu). As an example, below we generate random integers drawn from
# a Poisson parent distribtuion with mean 2000.
mu = 2000
Ydist = []
for i in range(1000):
    Ydist = Ydist + [np.random.poisson(mu)]
print('mean:', statistics.mean(Ydist))
print('variance:', statistics.stdev(Ydist)**2)
nbins = 25
plt.hist(Ydist, nbins, color = 'c', edgecolor = 'k')
plt.xlabel('Y')
plt.ylabel('counts')

# A photomultiplier tube (PMT) consists of a photocathode followed by a
# series of dynodes maintained at different electric potentials and then
# finally an anode.  When a single photon is incident on the photocathode
# it either produces a photoelectron or it doesn't.  The probability that
# it produces a photoelectron is determined by the quantum efficiency
# QE of the photocathode.  For this exercise we will assume that QE = 0.23.
# If a photoelectron is produced, it is accelerated towards the first dynode by
# means of a potential difference.  We assume that all electrons, whether
# produced at the photocathode or one of the dynodes, start with zero kinetic
# energy.  Therefore, the energy an electron gains is simply its charge times
# the potential difference between its starting and final positions.  When
# an electron collides with a dynode, secondary electrons are produced.  The
# average number of secondary electrons produced is proportional to the energy
# of the incoming electron.  In this problem, we assume that an electron
# accelerated through 20 V will, on average, produce one secondary electron
# upon colliding with the dynode.  Because we are "counting" electrons, the
# distribution of secondary electrons generated will follow a Poisson
# distribution.  

QE = 0.23; # Set the quantum efficeincy of the photocathode
maxi = int(1e4) # Set the number of Monte Carlo iterations
dynode = [0, 150, 300, 450, 600, 750, 850] # Set the number of dynodes and
# their voltages
dist = []
from datetime import datetime
print(datetime.now()) # Print the time at the simulation was started
for i in range(maxi): # This loop runs the Monte Carlo simulation maxi times
    # (10,000 in this case)
    if np.random.uniform() < QE: # Generate a random number between [0, 1] and check to see
        # if it's greater than the quantum efficiency
        electrons = 1 # If true, then a photoelectron is produced
        for j in range(len(dynode)-1): # This loop will examine the secondary
            # electrons produced at each of the dynodes
            mu = electrons*(dynode[j+1]-dynode[j])/20 # Mean number of
            # electrons generated at a dynode determined from the number of
            # incoming electons and the energy of the incoming electrons.
            # One electron accelerated through 20 V will on anverage produce
            # one electron
            electrons = electrons + np.random.poisson(mu)
        dist = dist + [electrons] # Add the total number of electrons detected at
        # the anode for this trial to a 'dist' list.  
print(datetime.now()) 

print('The fraction of trials that produced photoelectrons is', len(dist)/maxi)
 # Determines what fraction of maxi trials actually
# produced electrons at the PMT anode.  It should come out to be close the
# QE = 0.23.
import statistics
print('The mean of the distribution was', statistics.mean(dist))
# Calculate the mean of the distribution
plt.figure() # Generate a new figure
xbins = np.linspace(0, 2e6, 101) # Set bin ranges for the histogram.
# 100 equally-sized bins fron 0 to 2e6.
plt.hist(dist, xbins, color = 'c', edgecolor = 'k') # Plot the histgram of the 
# simulated number of anode electrons
plt.xlim((0, 2e6))
plt.ylim((0, 500))
plt.xlabel('Number of Anode Electons')
plt.ylabel('Counts')

counts, edges = np.histogram(dist, xbins) # Here we use 'np.histogram' to have 
# Python tell us how many counts there are in each of our 100 bins.

# Calculate the centre positions of the bins from the edge locations.
spacing = (edges[1] - edges[0])/2
centres = edges[1:] - spacing

plt.figure() # Start a new figure.
plt.semilogy(centres, counts, 'bo', fillstyle = 'none') # Make a semi-log plot 
# of the counts versus the number of detected anode electrons
plt.xlim((0, 2e6))
plt.ylim((.8, 500))
plt.xlabel('Number of Anode Electrons')
plt.ylabel('Counts')

# This next block of code is very similar to the previous.  This time,
# however, we imagine that there are 8 photons incident on the photocathode
# instead of one. Now there can be 0, 1, 2, ..., or 8 photoelectrons
# generated and then accelerated towards the first dynode.  How does the
# distribution of electrons arriving at the anode change?  We only require
# a relatively minor change to the code to study this problem.  Plot the
# histograms using the same scale and the same binwidths to make comparisons
# with the previous results easy.
QE = 0.23 # Set the quantum efficeincy of the photocathode
maxi = int(1e4) # Set the number of Monte Carlo iterations
dynode = [0, 150, 300, 450, 600, 750, 850] # Set the number of dynodes and
# their voltages
numPhotons = 8 # Set the number of photons incident on the photocathode.
dist = []
print(datetime.now()) # Print the time at the simulation was started
for i in range(maxi): # This loop runs the Monte Carlo simulation maxi times
    # (10,000 in this case)
    pe = 0 # Set the initial number of photoelectrons to zero.
    for k in range(numPhotons): # This loop individually steps through each of the
        # photons incident on the photocathode
        if np.random.uniform() < QE: # Generate a random number between [0, 1] and check to see
        # if it's greater than the quantum efficiency
            pe += 1 # Increment pe by one every time a photoelectron is
            # produced
    if pe > 0: # If there are any photoelectrons, then do something.  If there
        # are not photoelectrons then do nothing
        electrons = pe # Set the initial number of electrons equal to the
        # number of photoelectrons
        for j in range(len(dynode)-1): # This loop will examine the secondary
            # electrons produced at each of the dynodes
            mu = electrons*(dynode[j+1]-dynode[j])/20; # Mean number of
            # electrons generated at a dynode determined from the number of
            # incoming electons and the energy of the incoming electrons.
            # One electron accelerated through 20 V will on anverage produce
            # one electron
            electrons = electrons + np.random.poisson(mu)
        dist = dist + [electrons] # Add the total number of electrons detected at
        # the anode for this trial to a 'dist' list.
print(datetime.now()) # Print the time at the simulation completed.  
print('The fraction of trials that produced photoelectrons is', len(dist)/maxi)
 # Determines what fraction of maxi trials actually
# produced electrons at the PMT anode.  
print('The mean of the distribution was', statistics.mean(dist))
# Calculate the mean of the distribution
plt.figure() # Generate a new figure
xbins = np.linspace(0, 2e6, 101) # Set bin ranges for the histogram.
# 100 equally-sized bins fron 0 to 2e6.
plt.hist(dist, xbins, color = 'c', edgecolor = 'k') # Plot the histgram of the 
# simulated number of anode electrons
plt.xlim((0, 2e6))
plt.ylim((0, 500))
plt.xlabel('Number of Anode Electons')
plt.ylabel('Counts')
counts, edges = np.histogram(dist, xbins) # Here we use 'np.histogram' to have 
# Python tell us how many counts there are in each of our 100 bins.
# Calculate the centre positions of the bins from the edge locations.
spacing = (edges[1] - edges[0])/2
centres = edges[1:] - spacing
plt.figure() # Start a new figure.
plt.semilogy(centres, counts, 'bo', fillstyle = 'none') # Make a semi-log plot 
# of the counts versus the number of detected anode electrons
plt.xlim((0, 2e6))
plt.ylim((.8, 500))
plt.xlabel('Number of Anode Electrons')
plt.ylabel('Counts')

# Finally, this last block of code runs the same PMT simulation above four
# times, each time using a different number of incoming photons [1, 2, 4, 8].
# This chunk of code doesn't rely on anything that came before it.  These 44 lines
# of code will simulate the expected distributions of anode electrons for 4
# difference sets of incoming photons.  There real virtue of the Monte Carlo
# simulation is that we can now vary properties of the PMT with trivial
# modifications to the code below and systematically study the effects.
# For example, what if there were mode dynodes?  All we have to do is modify
# the line dynode:=[0,150,300,450,600,750,800].  Alternatively, we could
# keep the number of dynodes fixed and modify the potential applied to the
# dynodes.  Of course, after making this relatively simple simulation work,
# we could make modifications to make it more sophisticated.  What if 8
# photons are directed towards the photocathode, but they arrive at slightly
# different times.  What does the current pulse at the anode look like?
# That's not a problem that we'll tackle here, but it does demonstrate the
# versitility of the Monte Carlo method.
QE = 0.23 # Set the quantum efficeincy of the photocathode
maxi = int(1e4) # Set the number of Monte Carlo iterations
dynode = [0, 150, 300, 450, 600, 750, 850] # Set the number of dynodes and
# their voltages
numPhotons = [1, 2, 4, 8] # Set the number of photons incident on the
# photocathode.  We will repeat the simulation for 4 different numners of
# incoming photons
print(datetime.now()) # Print the time at the simulation was started
eventsList = []
distList = []
for m in range (len(numPhotons)):
    dist = []
    plt.figure() # Generate a new figure
    for i in range(maxi): # This loop runs the Monte Carlo simulation maxi times
    # (10,000 in this case)
        pe = 0 # Set the initial number of photoelectrons to zero.
        for k in range (numPhotons[m]): # This loop individually steps through each of the
        # photons incident on the photocathode
            if np.random.uniform() < QE: # Generate a random number between [0, 1]
            # and check to see if it's greater than the quantum efficiency
                pe += 1 # Increment pe by one every time a photoelectron is
                # produced
        if pe > 0: # If there are any photoelectrons, then do something.  If there
        # are not photoelectrons then do nothing
            electrons = pe # Set the initial number of electrons equal to the
            # number of photoelectrons
            for j in range(len(dynode)-1): # This loop will examine the secondary
            # electrons produced at each of the dynodes
                mu = electrons*(dynode[j+1]-dynode[j])/20 # Mean number of
                # electrons generated at a dynode determined from the number of
                # incoming electons and the energy of the incoming electrons.
                # One electron accelerated through 20 V will on anverage produce
                # one electron
                electrons = electrons + np.random.poisson(mu)
            dist = dist + [electrons] # Add the total number of electrons detected at
            # the anode for this trial to a 'dist' list.
    # This code can take some time to execute, so provide some intermediate
    # outputs to track progress
    print('Iteration', m)
    print('The fraction of trials that produced photoelectrons is', len(dist)/maxi)
    # Determines what fraction of maxi trials actually
    # produced electrons at the PMT anode.
    print('The mean of the distribution was', statistics.mean(dist))
    # Calculate the mean of the distribution
    print(datetime.now()) # Print the time that this iteration of  the simulation completed
    eventsList = eventsList + [[numPhotons[m], len(dist)/maxi]] # Keep track of
    # how many sets of incident photons actually resulted in a signal at
    # the anode.
    distList = distList + [dist] # Store the distributions recorded into a set
    xbins = xbins = np.linspace(0, 2e6, 101) # Set bin ranges for the histogram.
    # 100 equally-sized bins fron 0 to 2e6.
    plt.hist(dist, xbins, color = 'c', edgecolor = 'k') 
    #Plot the histgram of the simulated number of anode electrons
    plt.xlim((0, 2e6))
    plt.ylim((0, 500))
    plt.xlabel('Number of Anode Electons')
    plt.ylabel('Counts')
    counts, edges = np.histogram(dist, xbins) # Here we use 'np.histogram' to have 
    # Python tell us how many counts there are in each of our 100 bins.
    # Calculate the centre positions of the bins from the edge locations.
    spacing = (edges[1] - edges[0])/2
    centres = edges[1:] - spacing
    plt.figure() # Start a new figure.
    plt.semilogy(centres, counts, 'bo', fillstyle = 'none') # Make a semi-log plot 
    # of the counts versus the number of detected anode electrons    
    plt.xlim((0, 2e6))
    plt.ylim((0.8, 500))
    plt.xlabel('Number of Anode Electrons')
    plt.ylabel('Counts')

# Using the data stored in 'eventsSet', generate a plot of the average number
# of electrons detected at the anode versus the number of photons incident
# on the photocathode.
y = []
for i in range(len(eventsList)):
    y = y + [eventsList[i][1]]
plt.figure()
plt.plot(numPhotons, y, 'ro', linewidth = 2, fillstyle = 'none')
plt.xlabel('Number of Incident Photons')
plt.ylabel('Average Number of Electrons at Anode')

# Finally, using the distributions stored in the set 'distSet', plot the
# simulated distributions for the different numbers of incident photons all on
# the same semi-log plot.
plt.figure(figsize=(12,10))
colours = ['ko', 'bo', 'ro', 'go']
for i in range(len(numPhotons)):
     counts, edges = np.histogram(distList[i], xbins) 
     spacing = (edges[1] - edges[0])/2
     centres = edges[1:] - spacing
     plt.semilogy(centres, counts, colours[i], fillstyle = 'none', linewidth = 2)
plt.xlim((0, 2e6))
plt.ylim((0.8, 500))
plt.xlabel('Number of Anode Electrons')
plt.ylabel('Counts')
plt.legend(('1 photon', '2 photons', '4 photons', '8 photons'))