Merge lp:~christopher-hunt08/maus/minor_tweaks into lp:maus/merge

[Christopher Hunt]

Some small changes to v2.9.1

[Adam Dobbs (ajdobbs)]

Thanks Chris, do you have a Jenkins build for this branch?

[Christopher Hunt (christopher-hunt08)]

Sure its:
MAUS_Hunt_v2.6.1_tweaks

There should be more to come, I think this just some little things I found.

Chris

On 31-08-17 10:55, Adam Dobbs wrote:
> Thanks Chris, do you have a Jenkins build for this branch?
> --
> https://code.launchpad.net/~christopher-hunt08/maus/minor_tweaks/+merge/329995
> You are the owner of lp:~christopher-hunt08/maus/minor_tweaks.

--

Room: Blackett-514 Ext: 41625
E-mail: <email address hidden>
Imperial College London

[Adam Dobbs (ajdobbs)]

OK, I juts tried to do a clean build on lx02, after merging it into a clean merge branch:

bzr branch lp:maus/merge chunt
cd chunt
bzr merge lp:~christopher-hunt08/maus/minor_tweaks
./install_build_test.bash -j 4 --use-system-gcc true

I get the following error when it reaches the start of the unit tests:

Install file: "tests/integration/test_simulation/test_tof/tof_mc_plotter" as "build/tests/tof_mc_plotter"
Install file: "tests/cpp_unit/test_cpp_unit" as "build/tests/test_cpp_unit"
scons: done building targets.

INFO: Running the unit tests

Failure: ImportError (No module named maus_cpp.globals) ... ERROR

======================================================================
ERROR: Failure: ImportError (No module named maus_cpp.globals)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/loader.py", line 390, in loadTestsFromName
    addr.filename, addr.module)
  File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/importer.py", line 39, in importFromPath
    return self.importFromDir(dir_path, fqname)
  File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/importer.py", line 86, in importFromDir
    mod = load_module(part_fqname, fh, filename, desc)
  File "/vols/build/mice/adobbs/MAUS/maus/chunt/build/tests/test_maus_cpp/test_simulation.py", line 27, in <module>
    import maus_cpp.globals
ImportError: No module named maus_cpp.globals

----------------------------------------------------------------------
Ran 1 test in 0.001s

FAILED (errors=1)

^C
adobbs@lx02:~$

Any thoughts?

review: Needs Fixing

[Christopher Hunt (christopher-hunt08)]

I have no idea.

Unit test run fine for me and Jenkins.

I've seen things argue about the python includes before but I have no idea why. Can you perform the import on the command line?

Chris

On 31-08-17 15:44, Adam Dobbs wrote:
> Review: Needs Fixing
>
> OK, I juts tried to do a clean build on lx02, after merging it into a clean merge branch:
>
> bzr branch lp:maus/merge chunt
> cd chunt
> bzr merge lp:~christopher-hunt08/maus/minor_tweaks
> ./install_build_test.bash -j 4 --use-system-gcc true
>
> I get the following error when it reaches the start of the unit tests:
>
> Install file: "tests/integration/test_simulation/test_tof/tof_mc_plotter" as "build/tests/tof_mc_plotter"
> Install file: "tests/cpp_unit/test_cpp_unit" as "build/tests/test_cpp_unit"
> scons: done building targets.
>
> INFO: Running the unit tests
>
> Failure: ImportError (No module named maus_cpp.globals) ... ERROR
>
> ======================================================================
> ERROR: Failure: ImportError (No module named maus_cpp.globals)
> ----------------------------------------------------------------------
> Traceback (most recent call last):
> File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/loader.py", line 390, in loadTestsFromName
> addr.filename, addr.module)
> File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/importer.py", line 39, in importFromPath
> return self.importFromDir(dir_path, fqname)
> File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/importer.py", line 86, in importFromDir
> mod = load_module(part_fqname, fh, filename, desc)
> File "/vols/build/mice/adobbs/MAUS/maus/chunt/build/tests/test_maus_cpp/test_simulation.py", line 27, in <module>
> import maus_cpp.globals
> ImportError: No module named maus_cpp.globals
>
> ----------------------------------------------------------------------
> Ran 1 test in 0.001s
>
> FAILED (errors=1)
>
> ^C
> adobbs@lx02:~$
>
> Any thoughts?
> --
> https://code.launchpad.net/~christopher-hunt08/maus/minor_tweaks/+merge/329995
> You are the owner of lp:~christopher-hunt08/maus/minor_tweaks.

--

Room: Blackett-514 Ext: 41625
E-mail: <email address hidden>
Imperial College London

[Christopher Hunt (christopher-hunt08)]

Interestingly I also don't have a build/tests/ directory. And I don't think it does on jenkins either.

And yours does for some reason?

Chris

On 31-08-17 15:44, Adam Dobbs wrote:
> Review: Needs Fixing
>
> OK, I juts tried to do a clean build on lx02, after merging it into a clean merge branch:
>
> bzr branch lp:maus/merge chunt
> cd chunt
> bzr merge lp:~christopher-hunt08/maus/minor_tweaks
> ./install_build_test.bash -j 4 --use-system-gcc true
>
> I get the following error when it reaches the start of the unit tests:
>
> Install file: "tests/integration/test_simulation/test_tof/tof_mc_plotter" as "build/tests/tof_mc_plotter"
> Install file: "tests/cpp_unit/test_cpp_unit" as "build/tests/test_cpp_unit"
> scons: done building targets.
>
> INFO: Running the unit tests
>
> Failure: ImportError (No module named maus_cpp.globals) ... ERROR
>
> ======================================================================
> ERROR: Failure: ImportError (No module named maus_cpp.globals)
> ----------------------------------------------------------------------
> Traceback (most recent call last):
> File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/loader.py", line 390, in loadTestsFromName
> addr.filename, addr.module)
> File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/importer.py", line 39, in importFromPath
> return self.importFromDir(dir_path, fqname)
> File "/vols/build/mice/adobbs/MAUS/maus/chunt/third_party/install/lib/python2.7/site-packages/nose/importer.py", line 86, in importFromDir
> mod = load_module(part_fqname, fh, filename, desc)
> File "/vols/build/mice/adobbs/MAUS/maus/chunt/build/tests/test_maus_cpp/test_simulation.py", line 27, in <module>
> import maus_cpp.globals
> ImportError: No module named maus_cpp.globals
>
> ----------------------------------------------------------------------
> Ran 1 test in 0.001s
>
> FAILED (errors=1)
>
> ^C
> adobbs@lx02:~$
>
> Any thoughts?
> --
> https://code.launchpad.net/~christopher-hunt08/maus/minor_tweaks/+merge/329995
> You are the owner of lp:~christopher-hunt08/maus/minor_tweaks.

--

Room: Blackett-514 Ext: 41625
E-mail: <email address hidden>
Imperial College London

[chrisrogers1234 (chris-rogers)]

FYI import maus_cpp.globals is looking for $MAUS_ROOT_DIR/build/maus_cpp/globals.so

If you do

nm build/maus_cpp/globals.so | grep init

you should see a line like

0000000000004070 T initglobals

if that doesn't exist then somewhere the compile effort failed; it is compiling

src/py_cpp/PyGlobals.cc

which is done by method

build_python_modules(env)

in src/common_py/maus_build_tools/core_builder.py

hope that helps!

[Adam Dobbs (ajdobbs)]

Dobbs's fault, on my last push I missed a : when I was setting PYTHONPATH in configure. Not sure why my push didn't break Jenkins. Tests running on your branch now Chris.

=== modified file 'bin/utilities/data_tracker_verify.py'
--- bin/utilities/data_tracker_verify.py	2017-05-11 16:12:21 +0000
+++ bin/utilities/data_tracker_verify.py	2017-08-31 10:53:10 +0000
@@ -36,6 +36,7 @@
 import event_loader
 from analysis import tools
 import ROOT
+import glob
 
 # Useful Constants and configuration
 REFERENCE_PLANE = 0
@@ -1010,7 +1011,11 @@
 
 ##### 2. Load SciFi Events ####################################################
     print "\nLoading Spills...\n"
-    file_reader = event_loader.maus_reader(namespace.maus_root_files)
+    files = []
+    for filename in namespace.maus_root_files :
+      files += glob.glob(filename)
+
+    file_reader = event_loader.maus_reader(files)
 
     try :
       while file_reader.next_event() and \

=== added file 'src/common_py/analysis/beam_sampling.py'
--- src/common_py/analysis/beam_sampling.py	1970-01-01 00:00:00 +0000
+++ src/common_py/analysis/beam_sampling.py	2017-08-31 10:53:10 +0000
@@ -0,0 +1,1025 @@
+#!/usr/bin/env python
+
+# This file is part of MAUS: http://micewww.pp.rl.ac.uk/projects/maus
+#
+# MAUS is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# MAUS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+"""
+  Used for selecting  distributions of particles from a larger
+  distribution.
+
+  The parent distribution must be known prior to invocation as both the
+  covariance matrix and number of particles is required.
+"""
+
+# pylint: disable = W0311, E1101, W0102, R0902, C0103, R0201, W0223
+# pylint: disable = R0912, R0913, R0914, R0915, W0221, C0302
+
+import numpy
+from scipy.stats import chi2
+import math
+import ROOT
+
+import analysis.hit_types
+import analysis.inspectors
+import analysis.covariances
+
+
+MAX_STATISTICAL_WEIGHT = 10.0
+
+
+################################################################################
+
+def gaussian(x, mean, var) :
+  """
+    Analytically calculate values for a Gaussian distribution with given mean
+    and variance
+  """
+  var_inv = 1.0/var
+  vector = x - mean
+  const = 1.0 / math.sqrt( 2.0*math.pi * var )
+  exp = -0.5*( vector*vector*var_inv )
+
+  return const * math.exp( exp )
+
+def multivariate_gaussian(x, mean, cov) :
+  """
+    Analytically calculate vectors for a multivariate Gaussian distribution
+    with given mean vector and covariance matrix
+  """
+  cov_inv = numpy.linalg.inv(cov)
+  vector = x - mean
+  const = 1.0 / math.sqrt( numpy.linalg.det( 2.0*math.pi * cov ) )
+  exp = -0.5*( vector.transpose().dot(cov_inv.dot(vector)) )
+
+  return const * math.exp( exp )
+
+################################################################################
+
+
+class Sampler_base(object) :
+  """
+    Base class for sampler-type object definitions.
+    Outlines the interface.
+  """
+  def __init__(self) :
+    pass
+
+
+  def accept(self, hit) :
+    """
+      Test method. Returns true if the supplied hit should be accepted into the 
+      daughter distribution.
+    """
+    raise NotImplementedError("Sample_base.accept()"+\
+        " Method not overriden by derived class")
+
+
+  def weight(self, hit) :
+    """
+      Test method. Returns the statistical weight of the hit for sampling.
+    """
+    raise NotImplementedError("Sample_base.accept()"+\
+        " Method not overriden by derived class")
+
+
+  def calculate_parent(self, x) :
+    """
+      Returns the value of the parent distribution at a position x
+    """
+    raise NotImplementedError("Sample_base.calculate_parent()"+\
+        " Method not overriden by derived class")
+
+
+  def calculate_daughter(self, x) :
+    """
+      Returns the value of the daughter distribution at a position x
+    """
+    raise NotImplementedError("Sample_base.calculate_daughter()"+\
+        " Method not overriden by derived class")
+
+  def get_plots(self) :
+    """
+      Returns a dictionary of plots for analysis.
+    """
+    return {}
+
+
+
+
+class GaussianMomentumSampler(Sampler_base) :
+  """
+    Select a gaussian momentum distribution
+  """
+  def __init__(self, parent_dist, mean, std) :
+    """
+      Init.
+      Requires parent TH1F momenutm distribution
+      Mean of daughter distribution
+      STD of daughter distribution
+    """
+    Sampler_base.__init__(self)
+    self.__mean = mean
+    self.__var = std*std
+
+    self.__parent_distribution = parent_dist
+    self.__parent_distribution.Scale(\
+                                 1.0 / self.__parent_distribution.GetEntries())
+
+    self.__max = 5.0*std
+    self.__normalisation = 1.0e+20
+    self.__weight_normalisation = 0.0
+
+
+    # ACCEPT / REJECT
+    for bin_i in range(self.__parent_distribution.GetNbinsX()) :
+      x = self.__parent_distribution.GetBinCenter(bin_i)
+      if abs(x - mean) > self.__max :
+        continue
+
+      if self.calculate_daughter(x) > 0.0 : #Divide by zero errors
+        norm = self.calculate_parent(x) / self.calculate_daughter(x)
+      else :
+        norm = self.__normalisation
+
+      if norm < self.__normalisation :
+        self.__normalisation = norm
+
+
+    # STATISTICAL WEIGHTING
+    max_daughter = 0.0
+    max_parent = 0.0
+
+    for bin_i in range(self.__parent_distribution.GetNbinsX()) :
+      x = self.__parent_distribution.GetBinCenter(bin_i)
+      if abs(x - mean) > self.__max :
+        continue
+
+      daughter = self.calculate_daughter(x)
+      parent = self.calculate_parent(x)
+      if daughter > max_daughter :
+        max_daughter = daughter
+        max_parent = parent
+        
+    self.__weight_normalisation = max_parent / max_daughter
+
+
+  def calculate_parent(self, x) :
+    """
+      Calculate parent probability for momentum
+    """
+    found_bin = self.__parent_distribution.FindBin(x)
+    return self.__parent_distribution.GetBinContent(found_bin)
+
+
+  def calculate_daughter(self, x) :
+    """
+      Calculate daughter probability for momentum
+    """
+    return gaussian( x, self.__mean, self.__var )
+
+
+  def accept(self, hit) :
+    p = hit.get_p()
+
+    if abs(p - self.__mean) > self.__max :
+      return False
+
+    u = numpy.random.sample() #Uniform random number (0,1)
+    g_y = self.calculate_daughter(p)*self.__normalisation
+    f_y = self.calculate_parent(p)
+
+    if ( u >= (g_y / f_y) ) :
+      return False
+    else :
+      return True
+
+
+  def weight(self, hit) :
+    p = hit.get_p()
+
+    if abs(p - self.__mean) > self.__max :
+      return 0.0
+
+    g_y = self.calculate_daughter(p)*self.__weight_normalisation
+    f_y = self.calculate_parent(p)
+
+    ratio = (g_y / f_y)
+    return ratio
+
+
+
+class MomentumWindowSampler(Sampler_base) :
+  """
+    Sample a square momentum selection
+  """
+  def __init__(self, parent_dist, mean, window) :
+    Sampler_base.__init__(self)
+    self.__mean = mean
+    self.__window = window/2.0
+
+    self.__parent_distribution = parent_dist
+
+
+  def calculate_parent(self, x) :
+    """
+      Nothing to do
+    """
+    return None
+
+
+  def calculate_daughter(self, x) :
+    """
+      Nothing to do
+    """
+    return None
+
+
+  def accept(self, hit) :
+    p = hit.get_p()
+
+    if abs(p - self.__mean) > self.__window :
+      return False
+    else :
+      return True
+
+
+
+
+class GaussianXYSampler(Sampler_base) :
+  """
+    Multivariate Gaussian Sampling class.
+    Assumes both the parent and daughter distributions are gaussian.
+  """
+
+  def __init__(self, parent_means, parent_cov, means, cov) :
+    """
+      Expects 
+    """
+    Sampler_base.__init__(self)
+    if parent_means.shape != (2,) :
+      raise ValueError("Parent Means must be a 2x1 vector")
+    if parent_cov.shape != (2, 2) :
+      raise ValueError("Parent Covariance must be a 2x2 matrix")
+    if means.shape != (2,) :
+      raise ValueError("Means must be a 2x1 vector")
+    if cov.shape != (2,2) :
+      raise ValueError("Covariance must be a 2x2 matrix")
+
+    self.__parent_means = parent_means
+    self.__parent_covariance = parent_cov
+    self.__means = means
+    self.__covariance = cov
+
+    self.__normalisation = math.sqrt(numpy.linalg.det(self.__covariance) / \
+                                     numpy.linalg.det(self.__parent_covariance))
+
+
+  def calculate_parent(self, x) :
+    """
+      Returns the value of the parent distribution at a position x
+    """
+    return multivariate_gaussian(x, self.__parent_means, \
+                                                      self.__parent_covariance)
+
+
+  def calculate_daughter(self, x) :
+    """
+      Returns the value of the daughter distribution at a position x
+    """
+    return multivariate_gaussian(x, self.__means, self.__covariance)
+
+
+  def accept(self, hit) :
+    """
+      Returns true if the hit-vector is to be acceted into the daughter
+      distribution, otherwise returns false.
+    """
+    vector = [ hit.get_x(), hit.get_y() ]
+
+    u = numpy.random.sample() #Uniform random number (0,1)
+    g_y = self.calculate_daughter(vector)*self.__normalisation
+    f_y = self.calculate_parent(vector)
+
+    if ( u >= (g_y / f_y) ) :
+      return False
+    else :
+      return True
+
+
+
+
+class Gaussian4DSampler(Sampler_base) :
+  """
+    Multivariate Gaussian Sampling class.
+    Assumes both the parent and daughter distributions are gaussian.
+  """
+
+  def __init__(self, parent_means, parent_cov, means, cov) :
+    """
+      Expects 
+    """
+    Sampler_base.__init__(self)
+    if parent_means.shape != (4,) :
+      raise ValueError("Parent Means must be a 4x1 vector")
+    if parent_cov.shape != (4, 4) :
+      raise ValueError("Parent Covariance must be a 4x4 matrix")
+    if means.shape != (4,) :
+      raise ValueError("Means must be a 4x1 vector")
+    if cov.shape != (4,4) :
+      raise ValueError("Covariance must be a 4x4 matrix")
+
+    self.__parent_means = parent_means
+    self.__parent_covariance = parent_cov
+    self.__means = means
+    self.__covariance = cov
+
+    self.__normalisation = math.sqrt(numpy.linalg.det(self.__covariance) / \
+                                     numpy.linalg.det(self.__parent_covariance))
+
+
+  def calculate_parent(self, x) :
+    """
+      Returns the value of the parent distribution at a position x
+    """
+    return multivariate_gaussian(x, self.__parent_means, \
+                                                      self.__parent_covariance)
+
+
+  def calculate_daughter(self, x) :
+    """
+      Returns the value of the daughter distribution at a position x
+    """
+    return multivariate_gaussian(x, self.__means, self.__covariance)
+
+
+  def accept(self, hit) :
+    """
+      Returns true if the hit-vector is to be acceted into the daughter
+      distribution, otherwise returns false.
+    """
+    vector = hit.get_as_vector()[2:6]
+
+    u = numpy.random.sample() #Uniform random number (0,1)
+    g_y = self.calculate_daughter(vector)*self.__normalisation
+    f_y = self.calculate_parent(vector)
+
+    if ( u >= (g_y / f_y) ) :
+      return False
+    else :
+      return True
+
+
+class Amplitude4DSampler(Sampler_base) :
+  """
+    1D Chi-Squared Distribution Sampler
+    Assumes both the parent and daughter distributions are Chi-Squared
+    distributed.
+  """
+
+  def __init__(self, parent_distribution, covariance, emittance, max_x=1.0e12) :
+    """
+      Setup Amplitude4D Sampler. Requires the parent amplitude distribtion,
+      parent covariance matrix, required emittance and a maximum amplitude
+      range.
+    """
+    Sampler_base.__init__(self)
+    if type(parent_distribution) != ROOT.TH1F :
+      raise ValueError("Parent distribution must be a ROOT::TH1F")
+    if covariance.shape != (4,4) :
+      raise ValueError("Covariance must be a 4x4 matrix")
+
+    self.__parent_distribution = parent_distribution
+    self.__parent_distribution.Scale(\
+                                 1.0 / self.__parent_distribution.GetEntries())
+    self.__weight_normalisation = 0.0
+    self.__normalisation = 1.0e+20
+    self.__parent_emittance =\
+                         analysis.covariances.emittance_from_matrix(covariance)
+    self.__emittance = emittance
+    self.__parent_covariance = covariance
+    self.__parent_covariance_inv = numpy.linalg.inv(covariance)
+    self.__ndf = 4
+    self.__max = max_x
+
+    # ACCEPT / REJECT
+    for bin_i in range(self.__parent_distribution.GetNbinsX()) :
+      x = self.__parent_distribution.GetBinCenter(bin_i)
+      if x > self.__max :
+        break
+
+      if self.calculate_daughter(x) > 0.0 : #Divide by zero errors
+        norm = self.calculate_parent(x) / self.calculate_daughter(x)
+      else :
+        norm = self.__normalisation
+
+      if norm < self.__normalisation :
+        self.__normalisation = norm
+
+
+    # STATISTICAL WEIGHTING
+    max_daughter = 0.0
+    max_parent = 0.0
+
+    for bin_i in range(self.__parent_distribution.GetNbinsX()) :
+      x = self.__parent_distribution.GetBinCenter(bin_i)
+      if x > self.__max :
+        continue
+
+      daughter = self.calculate_daughter(x)
+      parent = self.calculate_parent(x)
+      if daughter > max_daughter :
+        max_daughter = daughter
+        max_parent = parent
+        
+    self.__weight_normalisation = max_parent / max_daughter
+
+    print "Parent =", self.__parent_emittance
+    print "Select =", self.__emittance
+
+
+  def _get_amplitude(self, x) :
+    """
+      Return the amplitude of a particle given the state vector.
+    """
+    amplitude = self.__parent_emittance*x.transpose().dot(\
+                                         self.__parent_covariance_inv.dot(x))
+    return amplitude
+
+
+  def calculate_parent(self, x) :
+    """
+      Calculate parent probability for amplitude
+    """
+    found_bin = self.__parent_distribution.FindBin(x)
+    return self.__parent_distribution.GetBinContent(found_bin)
+
+
+  def calculate_daughter(self, x) :
+    """
+      Calculate daughter probability for amplitude
+    """
+#    c = x*self.__emittance/(self.__parent_emittance**2)
+    c = x/self.__emittance
+    return chi2.pdf(c, self.__ndf)
+
+
+  def accept(self, hit) :
+    """
+      Returns true if the hit-vector is to be acceted into the daughter
+      distribution, otherwise returns false.
+    """
+
+    x = self._get_amplitude(numpy.array(hit.get_as_vector()[2:6]))
+
+    if x > self.__max :
+      return False
+
+    u = numpy.random.sample() #Uniform random number (0,1)
+    g_y = self.calculate_daughter(x)*self.__normalisation
+    f_y = self.calculate_parent(x)
+
+    if ( u >= (g_y / f_y) ) :
+      return False
+    else :
+      return True
+
+
+  def weight(self, hit) :
+    x = self._get_amplitude(numpy.array(hit.get_as_vector()[2:6]))
+
+    if x > self.__max :
+      return 0.0
+
+    g_y = self.calculate_daughter(x)*self.__weight_normalisation
+    f_y = self.calculate_parent(x)
+
+    ratio = (g_y / f_y)
+    return ratio
+
+
+
+
+class UniformAmplitude4DSampler(Sampler_base) :
+  """
+    1D Flat Amplitude Distribution Sampler
+  """
+
+  def __init__(self, parent_distribution, covariance, amplitude) :
+    """
+
+    """
+    Sampler_base.__init__(self)
+    if type(parent_distribution) != ROOT.TH1F :
+      raise ValueError("Parent distribution must be a ROOT::TH1F")
+    if covariance.shape != (4,4) :
+      raise ValueError("Covariance must be a 4x4 matrix")
+
+    self.__parent_distribution = parent_distribution
+    self.__parent_distribution.Scale(\
+                                 1.0 / self.__parent_distribution.GetEntries())
+    self.__weight_normalisation = 0.0
+    self.__normalisation = 1.0e+20
+    self.__parent_emittance =\
+                         analysis.covariances.emittance_from_matrix(covariance)
+    self.__max_amplitude = amplitude
+    self.__square = 1.0 / amplitude
+    self.__parent_covariance = covariance
+    self.__parent_covariance_inv = numpy.linalg.inv(covariance)
+
+    # ACCEPT / REJECT
+    for bin_i in range(1, self.__parent_distribution.GetNbinsX()) :
+      x = self.__parent_distribution.GetBinCenter(bin_i)
+      if x > self.__max_amplitude :
+        break
+
+      if self.calculate_daughter(x) > 0.0 : #Divide by zero errors
+        norm = self.calculate_parent(x) / self.calculate_daughter(x)
+      else :
+        norm = self.__normalisation
+
+      if norm < self.__normalisation :
+        self.__normalisation = norm
+
+    # STATISTICAL WEIGHTING
+    x = self.__max_amplitude * 0.5
+    max_daughter = self.calculate_daughter(x)
+    max_parent = self.calculate_parent(x)
+    self.__weight_normalisation = max_parent / max_daughter
+
+
+  def _get_amplitude(self, x) :
+    """
+      Calculate the amplitude of a particle given the state vector
+    """
+    amplitude = self.__parent_emittance*x.transpose().dot(\
+                                         self.__parent_covariance_inv.dot(x))
+    return amplitude
+
+
+  def calculate_parent(self, x) :
+    """
+      Calculate parent probability for amplitude
+    """
+    found_bin = self.__parent_distribution.FindBin(x)
+    return self.__parent_distribution.GetBinContent(found_bin)
+
+
+  def calculate_daughter(self, x) :
+    """
+      Calculate daughter probability for amplitude
+    """
+    if x > self.__max_amplitude :
+      return 0.0
+    else :
+      return self.__square
+
+
+  def accept(self, hit) :
+    """
+      Returns true if the hit-vector is to be acceted into the daughter
+      distribution, otherwise returns false.
+    """
+
+    x = self._get_amplitude(numpy.array(hit.get_as_vector()[2:6]))
+
+    if x > self.__max_amplitude :
+      return False
+
+    u = numpy.random.sample() #Uniform random number (0,1)
+    g_y = self.calculate_daughter(x)*self.__normalisation
+    f_y = self.calculate_parent(x)
+
+    if ( u >= (g_y / f_y) ) :
+      return False
+    else :
+      return True
+
+
+  def weight(self, hit) :
+    x = self._get_amplitude(numpy.array(hit.get_as_vector()[2:6]))
+
+    if x > self.__max_amplitude :
+      return 0.0
+
+    g_y = self.calculate_daughter(x)*self.__weight_normalisation
+    f_y = self.calculate_parent(x)
+
+    ratio = (g_y / f_y)
+    return ratio
+
+
+
+class XYPhaseSpaceSampler(Sampler_base) :
+  """
+    1D Chi-Squared Distribution Sampler
+    Assumes both the parent and daughter distributions are Chi-Squared
+    distributed.
+  """
+
+  def __init__(self, parent_distribution_x, parent_distribution_y,\
+                                                           means, covariance) :
+    """
+      Initialise the transverse phase space sampler. Requires the x-px phase
+      space distribution, the y-py phase space distribution, the means of the 
+      distributions and the covariance matrix.
+    """
+    Sampler_base.__init__(self)
+    if type(parent_distribution_x) != ROOT.TH2F :
+      raise ValueError("Parent distributions must be a ROOT::TH2F")
+    if type(parent_distribution_y) != ROOT.TH2F :
+      raise ValueError("Parent distributions must be a ROOT::TH2F")
+    if means.shape != (2,) :
+      raise ValueError("Means must be a 2x1 vector")
+    if covariance.shape != (2, 2) :
+      raise ValueError("Covariance must be a 2x2 vector")
+
+    self.__parent_distribution_x = parent_distribution_x
+    self.__parent_distribution_y = parent_distribution_y
+    self.__parent_distribution_x.Scale(\
+                               1.0 / self.__parent_distribution_x.GetEntries())
+    self.__parent_distribution_y.Scale(\
+                               1.0 / self.__parent_distribution_y.GetEntries())
+    self.__weight_normalisation_x = 0.0
+    self.__weight_normalisation_y = 0.0
+    self.__normalisation_x = 1.0e+20
+    self.__normalisation_y = 1.0e+20
+    self.__max = 1.0e+20
+    self.__means = means
+    self.__covariance = covariance
+    self.__covariance_inv = numpy.linalg.inv(covariance)
+
+    self.__max = 9.0
+
+    # ACCEPT / REJECT
+    for bin_i in range(self.__parent_distribution_x.GetNbinsX()) :
+      for bin_j in range(self.__parent_distribution_x.GetNbinsY()) :
+        x = self.__parent_distribution_x.GetXaxis().GetBinCenter(bin_i)
+        px = self.__parent_distribution_x.GetYaxis().GetBinCenter(bin_j)
+        vec = numpy.array([x,px])
+        if not self._is_contained(vec) :
+          continue
+
+        if self.calculate_daughter(vec) > 0.0 : #Divide by zero errors
+          norm = self.calculate_parent(vec, self.__parent_distribution_x) / \
+                                      self.calculate_daughter(numpy.array(vec))
+        else :
+          norm = self.__normalisation_x
+
+        if norm < self.__normalisation_x :
+          self.__normalisation_x = norm
+
+    for bin_i in range(self.__parent_distribution_y.GetNbinsX()) :
+      for bin_j in range(self.__parent_distribution_y.GetNbinsY()) :
+        x = self.__parent_distribution_y.GetXaxis().GetBinCenter(bin_i)
+        px = self.__parent_distribution_y.GetYaxis().GetBinCenter(bin_j)
+        vec = numpy.array([x,px])
+        if not self._is_contained(vec) :
+          continue
+
+        if self.calculate_daughter(vec) > 0.0 : #Divide by zero errors
+          norm = self.calculate_parent(vec, self.__parent_distribution_y) / \
+                                      self.calculate_daughter(numpy.array(vec))
+        else :
+          norm = self.__normalisation_y
+
+        if norm < self.__normalisation_y :
+          self.__normalisation_y = norm
+
+
+    # STATISTICAL WEIGHTING
+    max_daughter = 0.0
+    max_parent = 0.0
+
+    for bin_i in range(self.__parent_distribution.GetNbinsX()) :
+      for bin_j in range(self.__parent_distribution_x.GetNbinsY()) :
+        x = self.__parent_distribution_x.GetXaxis().GetBinCenter(bin_i)
+        px = self.__parent_distribution_x.GetYaxis().GetBinCenter(bin_j)
+        vec = numpy.array([x,px])
+        if not self._is_contained(vec) :
+          continue
+
+        daughter = self.calculate_daughter(vec)
+        parent = self.calculate_parent(vec)
+        if daughter > max_daughter :
+          max_daughter = daughter
+          max_parent = parent
+          
+    self.__weight_normalisation_x = max_parent / max_daughter
+
+    max_daughter = 0.0
+    max_parent = 0.0
+
+    for bin_i in range(self.__parent_distribution_y.GetNbinsX()) :
+      for bin_j in range(self.__parent_distribution_y.GetNbinsY()) :
+        x = self.__parent_distribution_y.GetXaxis().GetBinCenter(bin_i)
+        px = self.__parent_distribution_y.GetYaxis().GetBinCenter(bin_j)
+        vec = numpy.array([x,px])
+        if not self._is_contained(vec) :
+          continue
+
+        daughter = self.calculate_daughter(vec)
+        parent = self.calculate_parent(vec)
+        if daughter > max_daughter :
+          max_daughter = daughter
+          max_parent = parent
+          
+    self.__weight_normalisation_y = max_parent / max_daughter
+
+
+
+  def _is_contained(self, x) :
+    """
+      Function to determine whether or not a particle's state vector is 
+      contained withing the amplitude maximum
+    """
+    amp = x.transpose().dot(self.__covariance_inv.dot(x))
+#    print amp
+    if amp < self.__max :
+      return True
+    else :
+      return False
+
+
+  def calculate_parent(self, x, distribution=None) :
+    """
+      Calculate parent probability for the vector x, given the parent
+      distribution
+    """
+    if distribution is None :
+      distribution = self.__parent_distribution_x
+
+    found_bin = distribution.FindBin(x[0], x[1])
+    return distribution.GetBinContent(found_bin)
+
+
+  def calculate_daughter(self, x) :
+    """
+      Calculate daughter probability for a 4D vector
+    """
+    return multivariate_gaussian(x, self.__means, self.__covariance)
+
+
+  def accept(self, hit) :
+    """
+      Returns true if the hit-vector is to be acceted into the daughter
+      distribution, otherwise returns false.
+    """
+    x = [ hit.get_x(), hit.get_px() ]
+    y = [ hit.get_y(), hit.get_py() ]
+
+    if not self._is_contained(x) or not self._is_contained(y) :
+      return False
+
+    u = numpy.random.sample() #Uniform random number (0,1)
+    g_y_x = self.calculate_daughter(x)*self.__normalisation_x
+    f_y_x = self.calculate_parent(x, self.__parent_distribution_x)
+
+    g_y_y = self.calculate_daughter(y)*self.__normalisation_y
+    f_y_y = self.calculate_parent(y, self.__parent_distribution_y)
+
+    if ( u >= (g_y_x / f_y_x) ) and ( u >= (g_y_y / f_y_y) ) :
+      return False
+    else :
+      return True
+
+
+  def weight(self, hit) :
+    x = [ hit.get_x(), hit.get_px() ]
+    y = [ hit.get_y(), hit.get_py() ]
+
+    if not self._is_contained(x) or not self._is_contained(y) :
+      return 0.0
+
+    g_y_x = self.calculate_daughter(x)*self.__weight_normalisation_x
+    f_y_x = self.calculate_parent(x, self.__parent_distribution_x)
+
+    g_y_y = self.calculate_daughter(y)*self.__weight_normalisation_y
+    f_y_y = self.calculate_parent(y, self.__parent_distribution_y)
+
+    ratio_x = (g_y_x / f_y_x)
+    ratio_y = (g_y_y / f_y_y)
+    return ratio_x*ratio_y
+
+
+
+class XY4DPhaseSpaceSampler(Sampler_base) :
+  """
+    4D Sampler for beams that require the 6 2D projections from which the sample
+    is made.
+  """
+
+  def __init__(self, x_px, x_py, y_px, y_py, x_y, px_py, means, covariance) :
+    """
+
+    """
+    Sampler_base.__init__(self)
+    if type(x_px) != ROOT.TH2F :
+      raise ValueError("Parent distributions must be a ROOT::TH2F. Not "\
+                                                              +str(type(x_px)))
+    if type(x_py) != ROOT.TH2F :
+      raise ValueError("Parent distributions must be a ROOT::TH2F. Not "\
+                                                              +str(type(x_py)))
+    if type(y_py) != ROOT.TH2F :
+      raise ValueError("Parent distributions must be a ROOT::TH2F. Not "\
+                                                              +str(type(y_py)))
+    if type(y_px) != ROOT.TH2F :
+      raise ValueError("Parent distributions must be a ROOT::TH2F. Not "\
+                                                              +str(type(y_px)))
+    if type(px_py) != ROOT.TH2F :
+      raise ValueError("Parent distributions must be a ROOT::TH2F. Not "\
+                                                             +str(type(px_py)))
+    if type(x_y) != ROOT.TH2F :
+      raise ValueError("Parent distributions must be a ROOT::TH2F. Not "\
+                                                               +str(type(x_y)))
+
+    if means.shape != (4,) :
+      raise ValueError("Means must be a 4x1 vector")
+    if covariance.shape != (4, 4) :
+      raise ValueError("Covariance must be a 4x4 vector")
+
+    self.__all_cells = [(0, 1), (0, 2), (0, 3), (2, 1), (1, 3), (2, 3)]
+
+    self.__distributions = [ [ None for _ in range(4) ] for _ in range(4) ]
+    self.__distributions[0][1] = x_px
+    self.__distributions[0][2] = x_y
+    self.__distributions[0][3] = x_py
+    self.__distributions[2][1] = y_px
+    self.__distributions[1][3] = px_py
+    self.__distributions[2][3] = y_py
+
+    self.__normalisations = [ [ None for _ in range(4) ] for _ in range(4) ]
+    for a, b in self.__all_cells:
+      self.__normalisations[a][b] = 1.0e+20
+
+    self.__weight_normalisations =\
+                               [ [ None for _ in range(4) ] for _ in range(4) ]
+    for a, b in self.__all_cells:
+      self.__weight_normalisations[a][b] = 0.0
+
+    for a, b in self.__all_cells:
+      self.__distributions[a][b].Scale(\
+                             1.0/(self.__distributions[a][b].GetEntries()**2) )
+
+    self.__means = means
+    self.__covariance = covariance
+    self.__covariance_inv = numpy.linalg.inv(covariance)
+
+    self.__max = 12.0
+    self.__parent_min = 0.5
+
+    # ACCEPT / REJECT
+    while True :
+      for a, b in self.__all_cells :
+        dist = self.__distributions[a][b]
+        self.__normalisations[a][b] = 1.0e+20
+
+        for bin_i in range(dist.GetNbinsX()) :
+          for bin_j in range(dist.GetNbinsY()) :
+            x = dist.GetXaxis().GetBinCenter(bin_i)
+            px = dist.GetYaxis().GetBinCenter(bin_j)
+
+            vec = numpy.array([0.0, 0.0, 0.0, 0.0])
+            vec[a] = x
+            vec[b] = px
+
+            if not self._is_contained(vec) :
+              continue
+
+            if self.calculate_daughter(vec) > 0.0 : #Divide by zero errors
+              norm = self.calculate_parent(numpy.array([x, px]), dist) / \
+                                                   self.calculate_daughter(vec)
+            else :
+              norm = self.__normalisations[a][b]
+
+            if norm < self.__normalisations[a][b] :
+              self.__normalisations[a][b] = norm
+
+        if self.__normalisations[a][b] < self.__parent_min :
+          self.__max = self.__max - 1.0
+          break
+
+      else :
+        break
+
+    # STATISTICAL WEIGHTING
+    for a, b in self.__all_cells :
+      dist = self.__distributions[a][b]
+      self.__normalisations[a][b] = 0.0
+      max_daughter = 0.0
+      max_parent = 0.0
+
+      for bin_i in range(dist.GetNbinsX()) :
+        for bin_j in range(dist.GetNbinsY()) :
+          x = dist.GetXaxis().GetBinCenter(bin_i)
+          px = dist.GetYaxis().GetBinCenter(bin_j)
+
+          vec = numpy.array([0.0, 0.0, 0.0, 0.0])
+          vec[a] = x
+          vec[b] = px
+
+          if not self._is_contained(vec) :
+            continue
+
+          daughter = self.calculate_daughter(vec)
+          parent = self.calculate_parent(numpy.array([x, px]), dist)
+
+          if daughter > max_daughter :
+            max_daughter = daughter
+            max_parent = parent
+
+      self.__weight_normalisation[a][b] = max_parent / max_daughter
+
+
+  def _is_contained(self, x) :
+    """
+      Function to determine whether or not a particle's state vector is 
+      contained withing the amplitude maximum
+    """
+    amp = x.transpose().dot(self.__covariance_inv.dot(x))
+    if amp < self.__max :
+#      print amp, x
+      return True
+    else :
+      return False
+
+
+  def calculate_parent(self, x, distribution=None) :
+    """
+      Calculate parent probability for a vector x, given the parent
+      distribution
+    """
+    if distribution is None :
+      distribution = self.__x_y
+
+    found_bin = distribution.FindBin(x[0], x[1])
+    return distribution.GetBinContent(found_bin)
+
+
+  def calculate_daughter(self, x) :
+    """
+      Calculate daughter probability for a 4D vector
+    """
+    return multivariate_gaussian(x, self.__means, self.__covariance)
+
+
+  def accept(self, hit) :
+    """
+      Returns true if the hit-vector is to be acceted into the daughter
+      distribution, otherwise returns false.
+    """
+
+    vec = numpy.array(hit.get_as_vector()[2:6])
+
+#    if not self._is_contained(vec) :
+#      return False
+
+    u = numpy.random.sample() #Uniform random number (0,1)
+    daughter = self.calculate_daughter(vec)
+
+    for a, b in self.__all_cells :
+      dist = self.__distributions[a][b]
+      vector = [ vec[a], vec[b] ]
+      parent_result = self.calculate_parent(vector, dist)
+
+      result = (daughter*self.__normalisations[a][b] / parent_result)
+      if (u < result ) :
+        continue
+      else :
+        break
+    else : 
+#      self.__selection_inspector.add_hit(vec)
+      return True
+
+    return False
+
+
+  def weight(self, hit) :
+    vec = numpy.array(hit.get_as_vector()[2:6])
+
+    daughter = self.calculate_daughter(vec)
+    ratio = 1.0
+
+    for a, b in self.__all_cells :
+      dist = self.__distributions[a][b]
+      vector = [ vec[a], vec[b] ]
+      parent_result = self.calculate_parent(vector, dist)
+
+      result = (daughter*self.__weight_normalisations[a][b] / parent_result)
+
+      ratio *= result
+
+    return ratio
+

=== modified file 'src/common_py/analysis/covariances.py'
--- src/common_py/analysis/covariances.py	2017-05-10 10:29:14 +0000
+++ src/common_py/analysis/covariances.py	2017-08-31 10:53:10 +0000
@@ -119,7 +119,7 @@
     self._mean_vector = numpy.zeros( self._numVars )
     self._product_matrix = numpy.zeros( ( self._numVars, self._numVars ) )
     self._mean_momentum = 0.0
-    self._num_particles = 0
+    self._num_particles = 0.0
 
     self._correction_matrix = None
 
@@ -147,7 +147,7 @@
     self._mean_vector = numpy.zeros( self._numVars )
     self._product_matrix = numpy.zeros( ( self._numVars, self._numVars ) )
     self._mean_momentum = 0.0
-    self._num_particles = 0
+    self._num_particles = 0.0
 
   def add_hit_scifi( self, scifi_hit ) :
     """
@@ -162,15 +162,16 @@
       Add a hit and extract the information into the class.
       Accepts both xboa type hits and the local emittance analysis type hits.
     """
+    weight = hit.get_weight()
     for row, rowvar in enumerate( VARIABLE_LIST ) :
       for col, colvar in enumerate( VARIABLE_LIST ) :
         self._product_matrix[row][col] += hit.get( rowvar ) *\
-                                           hit.get( colvar )
-
-      self._mean_vector[row] += hit.get( rowvar )
-
-    self._mean_momentum += hit.get_p()
-    self._num_particles += 1
+                                           hit.get( colvar ) * weight**2.0
+
+      self._mean_vector[row] += hit.get( rowvar ) * weight
+
+    self._mean_momentum += hit.get_p()*weight
+    self._num_particles += 1.0*weight
 
 
   def add_dict( self, dictionary ) :

=== modified file 'src/common_py/analysis/hit_types.py'
--- src/common_py/analysis/hit_types.py	2017-05-10 10:29:14 +0000
+++ src/common_py/analysis/hit_types.py	2017-08-31 10:53:10 +0000
@@ -39,6 +39,7 @@
   def __init__( self, x = 0.0, y = 0.0, z = 0.0, \
                       px = 0.0, py = 0.0, pz = 0.0, station = 0, \
                       time = 0.0, mass = 105.6583715, p_value = 1.0, pid = 13, \
+                      weight = 1.0, \
                       scifi_track_point = None, virtual_track_point = None ) :
     """
       Initialise the object. this can be done in three ways:
@@ -58,6 +59,7 @@
       self.__p_value = p_value
       self.__pid = pid
       self.__station = station
+      self.__weight = weight
 #      self.__reference = reference
 
     elif scifi_track_point is not None and virtual_track_point is None :
@@ -73,13 +75,12 @@
       self.__pid = pid
       self.__station = tools.calculate_plane_id(scifi_track_point.tracker(), \
           scifi_track_point.station(), scifi_track_point.plane())
+      self.__weight = weight
 #      if reference is None :
 #        self.__reference = str(scifi_track_point.tracker())+"."+\
 #            str(scifi_track_point.station())+"."+str(scifi_track_point.plane())
 #      else :
 #        self.__reference = reference
-      if math.isnan(self.__x) :
-        raise ValueError("NaN Values Received from Scifi Track Point")
 
     elif scifi_track_point is None and virtual_track_point is not None :
       self.__x = virtual_track_point.GetPosition().x()
@@ -93,6 +94,7 @@
       self.__p_value = 1.0
       self.__pid = virtual_track_point.GetParticleId()
       self.__station = virtual_track_point.GetStationId()
+      self.__weight = weight
 #      self.__reference = reference
 #      if reference is None :
 #        self.__reference = virtual_track_point.GetStationId()
@@ -104,6 +106,10 @@
       raise ValueError( 'Please supply precisely one of "virtual_track_point"'+\
                  ' or "scifi_track_point", or specify all values explicitly.' )
 
+#    if math.isnan(self.__x) or math.isnan(self.__px) or math.isnan(self.__y) or math.isnan(self.__py) :
+    if math.isnan(self.__x) or math.isnan(self.__px) or math.isinf(self.__x) or math.isinf(self.__px) :
+      raise ValueError("NaN Values Received from Scifi Track Point")
+
 
   def __str__( self ) :
     """
@@ -202,9 +208,16 @@
     """
     self.__station = station
 
+  def set_weight( self, w ) :
+    """
+      Set the statistical weight of the hit
+    """
+    self.__weight = w
+
 #  def set_reference( self, string ) :
 #    self.__reference = string
 
+
   def get_x( self ) :
     """
       Get the x position
@@ -253,6 +266,7 @@
     """
     if math.isnan(math.sqrt( self.__px**2 + self.__py**2 + self.__pz**2 )) :
       print "ISNAN!!!", self
+      raise ValueError("NaN Values Received from Scifi Track Point")
     return math.sqrt( self.__px**2 + self.__py**2 + self.__pz**2 )
 
   def get_r( self ) :
@@ -322,6 +336,13 @@
     return self.__station
 
 
+  def get_weight( self ) :
+    """
+      Return the statistical weight of the hit
+    """
+    return self.__weight
+
+
 #  def get_reference( self ) :
 #    return self.__reference
 

=== modified file 'src/common_py/analysis/inspectors.py'
--- src/common_py/analysis/inspectors.py	2017-05-10 10:29:14 +0000
+++ src/common_py/analysis/inspectors.py	2017-08-31 10:53:10 +0000
@@ -17,7 +17,7 @@
 #
 
 # pylint: disable = W0311, E1101, W0613, C0111, R0911, W0621, C0103, R0902
-# pylint: disable = W0102, R0915
+# pylint: disable = W0102, R0915, R0201
 
 
 import analysis
@@ -71,13 +71,9 @@
       'inspected_phi_phasespace_{0}'.format(self.plane), '#phi-Phasespace', \
                                                                100, -4.0, 4.0 )
     self.theta_phasespace_plot = ROOT.TH1F(\
-      'inspected_phi_phasespace_{0}'.format(self.plane), '#theta-Phasespace', \
+      'inspected_theta_phasespace_{0}'.format(self.plane), '#theta-Phasespace', \
                                                                100, -4.0, 4.0 )
 
-    self.amplitude_momentum_plot = ROOT.TH2F(\
-        'inspected_A_p_phasespace_{0}'.format(self.plane), 'A-p-Phasespace' , \
-                                           200, 0.0, 100.0, 260, 130.0, 260.0 )
-
     self.pz_plot = ROOT.TH1F(\
               'inspected_pz_{0}'.format(self.plane), 'p_{z}', 400, 0.0, 400.0 )
     self.p_plot = ROOT.TH1F(\
@@ -112,7 +108,7 @@
                'inspected_beta_y_{0}'.format(self.plane), 'Y Beta Function', \
                                                              3000, 0.0, 3000.0)
     self.total_momentum_plot = ROOT.TH1F(\
-                'inspected_momentum_{0}'.format(self.plane), 'Momentum', \
+              'inspected_total_momentum_{0}'.format(self.plane), 'Momentum', \
                                                               1000, 0.0, 500.0)
 
     self.amplitude_plot = ROOT.TH1F(\
@@ -153,22 +149,22 @@
     if self.parent_covariance_inv is not None :
       vector = numpy.array(hit.get_as_vector()[2:6])
       amplitude = self.parent_emittance*vector.transpose().dot(self.parent_covariance_inv.dot(vector))
-      self.amplitude_plot.Fill(amplitude)
-      self.amplitude_momentum_plot.Fill(amplitude, hit.get_p())
+      self.amplitude_plot.Fill(amplitude, hit.get_weight())
+      self.amplitude_momentum_plot.Fill(amplitude, hit.get_p(), hit.get_weight())
 
     self.covariance.add_hit(hit)
 
-    self.position_plot.Fill(hit.get_x(), hit.get_y())
-    self.momentum_plot.Fill(hit.get_px(), hit.get_py())
-    self.x_phasespace_plot.Fill(hit.get_x(), hit.get_px())
-    self.y_phasespace_plot.Fill(hit.get_y(), hit.get_py())
-    self.xy_phasespace_plot.Fill(hit.get_x(), hit.get_py())
-    self.yx_phasespace_plot.Fill(hit.get_y(), hit.get_px())
-    self.rpt_phasespace_plot.Fill(hit.get_r(), hit.get_pt())
-    self.phi_phasespace_plot.Fill(hit.get_phi())
-    self.theta_phasespace_plot.Fill(hit.get_theta())
-    self.pz_plot.Fill(hit.get_pz())
-    self.p_plot.Fill(hit.get_p())
+    self.position_plot.Fill(hit.get_x(), hit.get_y(), hit.get_weight())
+    self.momentum_plot.Fill(hit.get_px(), hit.get_py(), hit.get_weight())
+    self.x_phasespace_plot.Fill(hit.get_x(), hit.get_px(), hit.get_weight())
+    self.y_phasespace_plot.Fill(hit.get_y(), hit.get_py(), hit.get_weight())
+    self.xy_phasespace_plot.Fill(hit.get_x(), hit.get_py(), hit.get_weight())
+    self.yx_phasespace_plot.Fill(hit.get_y(), hit.get_px(), hit.get_weight())
+    self.rpt_phasespace_plot.Fill(hit.get_r(), hit.get_pt(), hit.get_weight())
+    self.phi_phasespace_plot.Fill(hit.get_phi(), hit.get_weight())
+    self.theta_phasespace_plot.Fill(hit.get_theta(), hit.get_weight())
+    self.pz_plot.Fill(hit.get_pz(), hit.get_weight())
+    self.p_plot.Fill(hit.get_p(), hit.get_weight())
     self.position = hit.get_z()
 
     if self.ensemble_size > 1 :
@@ -312,3 +308,141 @@
     return ins_data
 
 
+####################################################################################################
+
+
+class CoolingInspector() :
+
+  def __init__(self) :
+    self.position = 0.0
+    self.number_particles = 0
+
+#    self.momentum_bias = None
+#    self.momentum_correlation = None
+
+    self.upstream_parent_covariance = None
+    self.upstream_parent_covariance_inv = None
+    self.upstream_parent_emittance = None
+
+    self.downstream_parent_covariance = None
+    self.downstream_parent_covariance_inv = None
+    self.downstream_parent_emittance = None
+
+    self.upstream_amplitude_plot = ROOT.TH1F(\
+                 'upstream_reference_amplitude', 'Amplitude', 200, 0.0, 100.0)
+    self.downstream_amplitude_plot = ROOT.TH1F(\
+                'downstream_reference_amplitude', 'Amplitude', 200, 0.0, 100.0)
+
+    self.upstream_amplitude_transmitted_plot = ROOT.TH1F(\
+                 'upstream_reference_amplitude_transmitted', 'Amplitude', 200, 0.0, 100.0)
+    self.downstream_amplitude_transmitted_plot = ROOT.TH1F(\
+                'downstream_reference_amplitude_transmitted', 'Amplitude', 200, 0.0, 100.0)
+
+    self.upstream_amplitude_scraped_plot = ROOT.TH1F(\
+                 'upstream_reference_amplitude_scraped', 'Amplitude', 200, 0.0, 100.0)
+    self.downstream_amplitude_scraped_plot = ROOT.TH1F(\
+                'downstream_reference_amplitude_scarped', 'Amplitude', 200, 0.0, 100.0)
+
+    self.amplitude_change_plot = ROOT.TH1F(\
+          'reference_amplitude_change', 'Amplitude Change', 200, -100.0, 100.0)
+    self.amplitude_change_amplitude_plot = ROOT.TH2F(\
+                  'reference_amplitude_change_amplitude', 'Amplitude Change', \
+                                          200, 0.0, 100.0, 200, -100.0, 100.0 )
+    self.amplitude_mobility_plot = ROOT.TH2F(\
+                        'reference_amplitude_mobility', 'Amplitude Mobility', \
+                                             100, 0.0, 100.0, 100, 0.0, 100.0 )
+
+
+#  def set_covariance_correction(self, correction, axes=['x', 'px', 'y', 'py']) :
+#    self.covariance.set_covariance_correction( correction, axes )
+#    self.ensemble_covariance.set_covariance_correction( correction, axes )
+
+
+  def set_upstream_covariance(self, matrix) :
+    self.upstream_parent_covariance = numpy.array(matrix)
+    self.upstream_parent_covariance_inv = numpy.linalg.inv(self.upstream_parent_covariance)
+    self.upstream_parent_emittance = analysis.covariances.emittance_from_matrix(self.upstream_parent_covariance)
+
+
+  def set_downstream_covariance(self, matrix) :
+    self.downstream_parent_covariance = numpy.array(matrix)
+    self.downstream_parent_covariance_inv = numpy.linalg.inv(self.downstream_parent_covariance)
+    self.downstream_parent_emittance = analysis.covariances.emittance_from_matrix(self.downstream_parent_covariance)
+
+
+#  def set_momentum_correction(self, constant, gradient) :
+#    self.momentum_bias = constant
+#    self.momentum_correlation = gradient
+
+
+  def add_hits(self, up_hit, down_hit) :
+#    if self.momentum_bias is not None :
+#      p = hit.get_p()
+##      correction = 1.0 + (self.momentum_bias + self.momentum_correlation*p) / p
+#      correction = ((p - self.momentum_bias) / \
+#                                         (1.0 - self.momentum_correlation)) / p
+#      hit.set_px(hit.get_px()*correction)
+#      hit.set_py(hit.get_py()*correction)
+#      hit.set_pz(hit.get_pz()*correction)
+
+    if (self.upstream_parent_covariance_inv is not None) and (self.downstream_parent_covariance_inv is not None) :
+      self.number_particles += 1
+
+      if up_hit is not None :
+        vector = numpy.array(up_hit.get_as_vector()[2:6])
+        up_amplitude = self.upstream_parent_emittance*vector.transpose().dot(self.upstream_parent_covariance_inv.dot(vector))
+        self.upstream_amplitude_plot.Fill(up_amplitude)
+        
+        if down_hit is None :
+          self.upstream_amplitude_scraped_plot.Fill(up_amplitude)
+        else :
+          self.upstream_amplitude_transmitted_plot.Fill(up_amplitude)
+
+      if down_hit is not None :
+        vector = numpy.array(down_hit.get_as_vector()[2:6])
+        down_amplitude = self.downstream_parent_emittance*vector.transpose().dot(self.downstream_parent_covariance_inv.dot(vector))
+        self.downstream_amplitude_plot.Fill(down_amplitude)
+
+        if up_hit is None :
+          self.downstream_amplitude_scraped_plot.Fill(down_amplitude)
+        else :
+          self.downstream_amplitude_transmitted_plot.Fill(down_amplitude)
+
+
+      if (up_hit is not None) and (down_hit is not None) :
+        DA = down_amplitude - up_amplitude
+        self.amplitude_change_plot.Fill( DA )
+        self.amplitude_change_amplitude_plot.Fill( up_amplitude, DA )
+        self.amplitude_mobility_plot.Fill( up_amplitude, down_amplitude )
+
+
+  def fill_plots(self) :
+    pass
+
+
+  def get_length(self) :
+    return self.number_particles
+
+
+  def get_plot_dictionary(self) :
+    ins_plots = {}
+
+    ins_plots['upstream_amplitude'] = self.upstream_amplitude_plot
+    ins_plots['downstream_amplitude'] = self.downstream_amplitude_plot
+    ins_plots['upstream_amplitude_transmitted'] = self.upstream_amplitude_transmitted_plot
+    ins_plots['downstream_amplitude_transmitted'] = self.downstream_amplitude_transmitted_plot
+    ins_plots['upstream_amplitude_scraped'] = self.upstream_amplitude_scraped_plot
+    ins_plots['downstream_amplitude_scraped'] = self.downstream_amplitude_scraped_plot
+
+    ins_plots['amplitude_change'] = self.amplitude_change_plot
+    ins_plots['amplitude_change_amplitude'] = self.amplitude_change_amplitude_plot
+    ins_plots['amplitude_mobility'] = self.amplitude_mobility_plot
+
+    return ins_plots
+
+
+  def get_data_dictionary(self) :
+    ins_data = {}
+
+    return ins_data
+

=== modified file 'src/common_py/analysis/tools.py'
--- src/common_py/analysis/tools.py	2017-02-17 12:38:12 +0000
+++ src/common_py/analysis/tools.py	2017-08-31 10:53:10 +0000
@@ -216,6 +216,10 @@
       new_location = os.path.join( location, key )
       print_plots( plot_dict[key], new_location )
       
+    elif type( plot_dict[key] ) is ROOT.TH2F :
+      canvas = ROOT.TCanvas( key+'_canvas' )
+      plot_dict[key].Draw("COL")
+      canvas.SaveAs( os.path.join(location, key ) + ".pdf", "pdf" )
     else :
       canvas = ROOT.TCanvas( key+'_canvas' )
       plot_dict[key].Draw()

=== modified file 'src/common_py/event_loader.py'
--- src/common_py/event_loader.py	2017-05-04 14:09:47 +0000
+++ src/common_py/event_loader.py	2017-08-31 10:53:10 +0000
@@ -49,6 +49,7 @@
       self.__filenames = [ filename ]
 
     self.__current_filenumber = -1
+    self.__current_statistical_weight = 1.0
 
     self.__tree = None
     self.__data = None
@@ -120,6 +121,13 @@
                                                   str(print_progress) )
 
 
+  def get_current_statistical_weight( self ) :
+    """
+      Returns the statistical weight of the last loaded event.
+    """
+    return self.__current_statistical_weight
+
+
   def select_events( self, selected_events_dict ) :
     """
       Only loads events that are present in the supplied dictionary.
@@ -139,7 +147,7 @@
       return False
 
 
-  def save_event( self ) :
+  def save_event( self, weight=1.0 ) :
     """
       Add the event location to an internal list of good events.
     """
@@ -149,9 +157,9 @@
     if filename not in self.__saved_events :
       self.__saved_events[filename] = {}
     if spill not in self.__saved_events[filename] :
-      self.__saved_events[filename][spill] = []
+      self.__saved_events[filename][spill] = {}
 
-    self.__saved_events[filename][spill].append(event)
+    self.__saved_events[filename][spill][event] = weight
 
 
   def get_saved_events( self ) :
@@ -252,10 +260,13 @@
     else :
       while self.next_event():
         spill = str(self.__current_spill_num)
+        event = str(self.__current_event_num)
         if self.__current_filename in self.__selected_events :
           if spill in self.__selected_events[self.__current_filename] :
-            if self.__current_event_num in \
+            if event in \
                        self.__selected_events[self.__current_filename][spill] :
+              weight = self.__selected_events[self.__current_filename][spill][event]
+              self.__current_statistical_weight = weight
               return True
       else:
         return False
@@ -290,7 +301,8 @@
 
       self.__tree.GetEntry( self.__current_spill_num )
       self.__spill = self.__data.GetSpill()
-      self.__current_num_events = self.__spill.GetReconEvents().size()
+      self.__current_num_events = max( self.__spill.GetReconEvents().size(), \
+                                            self.__spill.GetMCEvents().size() )
       self.__current_event_num = 0
 
       if self.__spill.GetDaqEventType() == "physics_event" :

=== modified file 'src/legacy/Simulation/FillMaterials.cc'
--- src/legacy/Simulation/FillMaterials.cc	2017-08-22 21:21:20 +0000
+++ src/legacy/Simulation/FillMaterials.cc	2017-08-31 10:53:10 +0000
@@ -385,6 +385,25 @@
   TufsetPU->AddElement(elO, 20.60*perCent);
   materials_list->addMaterial( TufsetPU, name );
 
+  // POLYVINYL_TOLUENE_
+  density = 1.032 * g/cm3;
+  name = "POLYVINYL_TOLUENE";
+  G4Material* PVT = new G4Material(name, density, nComp=2, kStateSolid);
+  PVT->AddElement(elC, 91.5023*perCent);
+  PVT->AddElement(elH,  8.4977*perCent);
+  materials_list->addMaterial( PVT, name );
+
+  G4Element* elCl = man->FindOrBuildElement("Cl");
+
+  // POLYVINYL_CHLORIDE_LOWD
+  density = 0.12 * g/cm3;
+  name = "POLYVINYL_CHLORIDE_LOWD";
+  G4Material* PVC_LD = new G4Material(name, density, nComp=3, kStateSolid);
+  PVC_LD->AddElement(elC,  38.436*perCent);
+  PVC_LD->AddElement(elH,   4.838*perCent);
+  PVC_LD->AddElement(elCl, 56.726*perCent);
+  materials_list->addMaterial( PVC_LD, name );
+
   // TUFNOL 
 
   // now fill in the other maps of the material properties

=== modified file 'third_party/bash/40python_extras.bash'
--- third_party/bash/40python_extras.bash	2017-02-16 23:33:58 +0000
+++ third_party/bash/40python_extras.bash	2017-08-31 10:53:10 +0000
@@ -25,7 +25,7 @@
 
 # The packages to download for the MAUS third party tarball
 download_package_list="\
- numpy anyjson python-dateutil>=1.5,<2.0 amqplib>=1.0 six>=1.4.0 \
+ numpy scipy anyjson python-dateutil>=1.5,<2.0 amqplib>=1.0 six>=1.4.0 \
  logilab-common logilab-astng suds validictory nose==1.1 nose-exclude \
  coverage doxypy astroid isort pylint==1.6.3 bitarray \
  pymongo==2.3 readline matplotlib==2.0.0 \
@@ -35,7 +35,7 @@
 # The packages to install
 # Note: pil is no longer on PyPI so we only use our existing tarball
 package_list="\
- numpy anyjson python-dateutil amqplib six \
+ numpy scipy anyjson python-dateutil amqplib six \
  logilab-common logilab-astng  suds validictory nose==1.1 nose-exclude \
  coverage doxypy pylint==1.6.3 bitarray \
  pymongo readline matplotlib \
@@ -43,7 +43,7 @@
 "
 
 # The modules to test
-module_test_list="numpy suds validictory nose coverage \
+module_test_list="numpy scipy suds validictory nose coverage \
  pylint bitarray matplotlib pymongo \
  Image django magickwand celery" #Image is pil bottom level