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

[Christopher Hunt]

Some tweaks and additions to the analysis code in src/common_py, plus a bug fix in single_thread.py where using the data object wasn't being correctly deleted at the end of the pipeline.

[Paolo Franchini (p-franchini)]

Thanks Chris,

I will do it later today,

cheers,

Paolo

On 09/05/18 10:32, Christopher Hunt wrote:
> Christopher Hunt has proposed merging lp:~christopher-hunt08/maus/maus_analysis_devel into lp:maus/merge.
>
> Requested reviews:
> MAUS Maintainers (maus-maintainers)
>
> For more details, see:
> https://code.launchpad.net/~christopher-hunt08/maus/maus_analysis_devel/+merge/345284
>
> Some tweaks and additions to the analysis code in src/common_py, plus a bug fix in single_thread.py where using the data object wasn't being correctly deleted at the end of the pipeline.
>

=== modified file 'bin/scifi/tracker_resolution_plots.py'
--- bin/scifi/tracker_resolution_plots.py	2017-12-12 10:11:06 +0000
+++ bin/scifi/tracker_resolution_plots.py	2018-05-09 09:31:35 +0000
@@ -318,6 +318,11 @@
               tracker+'_p_residual_p', "p Residuals in p", \
                                      PZ_BIN, PZ_MIN, PZ_MAX, 500, -50.0, 50.0 )
 
+    tracker_dict['dpp_residual_p'] = ROOT.TH2F( \
+              tracker+'_dpp_residual_p', "dp/p Residuals in p", \
+                                     PZ_BIN, PZ_MIN, PZ_MAX, 500, -1.0, 1.0 )
+
+
 
     tracker_dict['x_residual_pt'] = ROOT.TH2F( \
                tracker+'_x_residual_pt', "X Residuals in p_{t}", \
@@ -1016,6 +1021,7 @@
     tracker_plots['pt_residual_p'].Fill( P_mc, Pt_res )
     tracker_plots['pz_residual_p'].Fill( P_mc, res_mom[2] )
     tracker_plots['p_residual_p'].Fill( P_mc, P_res )
+    tracker_plots['dpp_residual_p'].Fill( P_mc, P_res/P_mc )
 
 
     tracker_plots['x_residual_pz'].Fill( Pz_mc, res_pos[0] )
@@ -1230,6 +1236,45 @@
         plot_dict[tracker][component+plot_axis+'_bias'] = bias_graph
 
 
+  for tracker in [ "upstream", "downstream" ] :
+    for a_plot in [ "dpp_residual_p" ] :
+      plot = plot_dict[tracker][a_plot]
+
+      rms_error = array.array( 'd' )
+      bin_size = array.array( 'd' )
+      bins = array.array( 'd' )
+      rms = array.array( 'd' )
+      mean = array.array( 'd' )
+      mean_error = array.array( 'd' )
+
+      width = plot.GetXaxis().GetBinWidth(1)
+      for i in range( 0, plot.GetXaxis().GetNbins() ) :
+        projection = plot.ProjectionY( a_plot+'_pro_'+str(i), i, (i+1) )
+
+        plot_mean = plot.GetXaxis().GetBinCenter( i ) + width
+        pro_mean, pro_mean_err, pro_std, pro_std_err = \
+                                      analysis.tools.fit_gaussian(projection)
+
+        bin_size.append( width*0.5 )
+        bins.append( plot_mean )
+        rms.append( pro_std )
+        rms_error.append( pro_std_err )
+        mean.append( pro_mean )
+        mean_error.append( pro_mean_err )
+
+      if len(bins) != 0 :
+        resolution_graph = ROOT.TGraphErrors( len(bins), \
+                                             bins, rms, bin_size, rms_error )
+        bias_graph = ROOT.TGraphErrors( len(bins), \
+                                           bins, mean, bin_size, mean_error )
+      else :
+        resolution_graph = None
+        bias_graph = None
+
+      plot_dict[tracker][a_plot+'_resolution'] = resolution_graph
+      plot_dict[tracker][a_plot+'_bias'] = bias_graph
+
+
 
   for tracker in [ "upstream", "downstream" ] :
 #    for component in [ "pt_", "pz_", ] :
 
=== modified file 'src/common_py/analysis/__init__.py'
--- src/common_py/analysis/__init__.py	2017-02-17 12:38:12 +0000
+++ src/common_py/analysis/__init__.py	2018-05-09 09:31:35 +0000
@@ -31,4 +31,5 @@
 import analysis.hit_types
 import analysis.inspectors
 import analysis.scifitools # pylint: disable = E0401
+import analysis.systematic_error_estimation
 
 
=== added file 'src/common_py/analysis/analysis_track.py'
--- src/common_py/analysis/analysis_track.py	1970-01-01 00:00:00 +0000
+++ src/common_py/analysis/analysis_track.py	2018-05-09 09:31:35 +0000
@@ -0,0 +1,141 @@
+#!/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/>.
+#
+
+"""
+  Defines the AnalysisTrack class.
+
+  A generic track class that can be used to simplify the analysis of track data.
+"""
+
+# pylint: disable = W0311, R0902, R0904, R0913, C0103, W0102
+
+class AnalysisTrack() :
+  """
+    A simple, unified class that holds the essential information, when copied
+    from MAUS data types, during various analyses.
+
+    Some syntax was borrowed from XBoa for simplicity and because I like the
+    interface!
+  """
+  def __init__( self, trackpoints={}, chisq=0.0, ndf=0, p_value=1.0, status=0 ) :
+    """
+      Initialise object. 
+    """
+    self.__trackpoints = trackpoints
+    self.__chisq = chisq
+    self.__ndf = ndf
+    self.__p_value = p_value
+    self.__status = status
+
+
+  def __getitem__(self, index) :
+    """
+      Return a trackpoint from the track at index ``index''
+    """
+    return self.__trackpoints[index]
+    
+
+  def __setitem__(self, index, item) :
+    """
+      Set the value of the trackpoint in the track at index ``index''
+    """
+    self.__trackpoints[index] = item
+
+
+  def __len__(self) :
+    """
+      Return the number of trackpoints in the track>
+    """
+    return len(self.__trackpoints)
+    
+
+  def get_trackpoint( self, index ) : 
+    """
+      Return a trackpoint from the track at index ``index''
+    """
+    return self.__trackpoints[index]
+
+
+  def set_trackpoint( self, index, item ) : 
+    """
+      Set the value of the trackpoint in the track at index ``index''
+    """
+    self.__trackpoints[index] = item
+
+
+  def get_p_value( self ) :
+    """
+      Return the track p-value
+    """
+    return self.__p_value
+
+
+  def set_p_value( self, pvalue ) :
+    """
+      Set the track p-value
+    """
+    self.__p_value = pvalue
+
+
+  def get_chisq( self ) :
+    """
+      Return the track chisq squared
+    """
+    return self.__chisq
+
+
+  def set_chisq( self, chisq ) :
+    """
+      Set the track chisq squared
+    """
+    self.__chisq = chisq
+
+
+  def get_ndf( self ) :
+    """
+      Get the track No. Degrees of Freedom
+    """
+    return self.__ndf
+
+
+  def set_ndf( self, ndf ) :
+    """
+      Set the track No. Degrees of Freedom
+    """
+    self.__ndf = ndf
+
+
+  def get_chisq_ndf(self) :
+    """
+      Get the track Chisquared / No. Degrees of Freedom
+    """
+    return self.__chisq / self.__ndf
+
+
+  def set_status(self, status) :
+    """
+      Set a undefined status flag
+    """
+    self.__status = status
+
+
+  def get_status(self) :
+    """
+      Get the status flag
+    """
+    return self.__status
+
 
=== modified file 'src/common_py/analysis/hit_types.py'
--- src/common_py/analysis/hit_types.py	2017-08-15 16:19:42 +0000
+++ src/common_py/analysis/hit_types.py	2018-05-09 09:31:35 +0000
@@ -23,7 +23,7 @@
   Note that this class is used by the add_hit meth
 """
 
-# pylint: disable = W0311, R0902, R0904, R0913, C0103
+# pylint: disable = W0311, R0902, R0904, R0913, C0103, R0914, R0915
 
 import math
 import analysis.tools as tools
@@ -40,14 +40,15 @@
                       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 ) :
+                      scifi_track_point = None, virtual_track_point = None, \
+                      global_track_point = None ) :
     """
       Initialise the object. this can be done in three ways:
       1. Specify all components by hand
       2. Build from a sci-fr trackpoint
       3. Build from a virtual trackpoint
     """
-    if scifi_track_point is None and virtual_track_point is None :
+    if scifi_track_point is None and virtual_track_point is None and global_track_point is None :
       self.__x = x
       self.__y = y
       self.__z = z
@@ -62,7 +63,7 @@
       self.__weight = weight
 #      self.__reference = reference
 
-    elif scifi_track_point is not None and virtual_track_point is None :
+    elif scifi_track_point is not None and virtual_track_point is None and global_track_point is None :
       self.__x = scifi_track_point.pos().x()
       self.__y = scifi_track_point.pos().y()
       self.__z = scifi_track_point.pos().z()
@@ -82,7 +83,7 @@
 #      else :
 #        self.__reference = reference
 
-    elif scifi_track_point is None and virtual_track_point is not None :
+    elif scifi_track_point is None and virtual_track_point is not None and global_track_point is None :
       self.__x = virtual_track_point.GetPosition().x()
       self.__y = virtual_track_point.GetPosition().y()
       self.__z = virtual_track_point.GetPosition().z()
@@ -101,6 +102,20 @@
 #      else :
 #        self.__reference = reference
 
+    elif scifi_track_point is None and virtual_track_point is None and global_track_point is not None :
+      self.__x = global_track_point.get_position().X()
+      self.__y = global_track_point.get_position().Y()
+      self.__z = global_track_point.get_position().Z()
+      self.__px = global_track_point.get_momentum().X()
+      self.__py = global_track_point.get_momentum().Y()
+      self.__pz = global_track_point.get_momentum().Z()
+      self.__time = global_track_point.get_position().T()
+      self.__mass = mass
+      self.__p_value = p_value
+      self.__pid = pid
+      self.__station = station
+      self.__weight = weight
+
     else :
       print "WTF!"
       raise ValueError( 'Please supply precisely one of "virtual_track_point"'+\
@@ -388,3 +403,200 @@
 
   for key, val in __hit_get_variables.iteritems() :
     __hit_get_keys.append( key )
+
+
+################################################################################
+
+
+class AnalysisSpacePoint() :
+  """
+    A simple, unified class that holds the essential information, when copied
+    from MAUS data types, during various analyses.
+
+    Some syntax was borrowed from XBoa for simplicity and because I like the
+    interface!
+  """
+  def __init__( self, x = 0.0, y = 0.0, z = 0.0, station = 0, \
+                      time = 0.0, weight = 1.0, tof = None, scifi = None ) :
+    """
+      Initialise the object. this can be done in three ways:
+      1. Specify all components by hand
+      2. Build from a sci-fr trackpoint
+      3. Build from a virtual trackpoint
+    """
+    if tof is None and scifi is None :
+      self.__x = x
+      self.__y = y
+      self.__z = z
+      self.__time = time
+      self.__station = station
+      self.__weight = weight
+
+    elif scifi is not None and tof is None :
+      self.__x = scifi.pos().x()
+      self.__y = scifi.pos().y()
+      self.__z = scifi.pos().z()
+      self.__time = time
+      self.__station = tools.calculate_plane_id(scifi.tracker(), \
+          scifi.station(), 1)
+      self.__weight = weight
+
+    elif scifi is None and tof is not None :
+      self.__x = tof.GetGlobalPosX()
+      self.__y = tof.GetGlobalPosY()
+      self.__z = tof.GetGlobalPosZ()
+      self.__time = tof.GetTime()
+      self.__station = tof.GetStation()
+      self.__weight = weight
+
+    else :
+      print "WTF!"
+      raise ValueError( 'Please supply precisely one of "virtual_track_point"'+\
+                 ' or "scifi_track_point", or specify all values explicitly.' )
+
+
+  def __str__( self ) :
+    """
+      Return a string of the parameters.
+      Mostly useful for debuging.
+    """
+    return '(' + str( self.__x ) + ',' +\
+                 str( self.__y ) + ',' +\
+                 str( self.__z ) + '):[' +\
+                 str( self.__time ) +\
+                 str( self.__station ) + ']'
+
+
+  def __repr__( self ) :
+    """
+      Return a string of the parameters.
+      Mostly useful for debuging.
+    """
+    return '(' + str( self.__x ) + ',' +\
+                 str( self.__y ) + ',' +\
+                 str( self.__z ) + '):[' +\
+                 str( self.__time ) +\
+                 str( self.__station ) + ']'
+
+
+
+  def set_x( self, val ) :
+    """
+      Set the x position
+    """
+    self.__x = val
+
+  def set_y( self, val ) :
+    """
+      Set the y position
+    """
+    self.__y = val
+
+  def set_z( self, val ) :
+    """
+      Set the z position
+    """
+    self.__z = val
+
+  def set_time( self, val ) :
+    """
+      Set the time
+    """
+    self.__time = val
+
+  def set_station( self, station ) :
+    """
+      Set the station ID
+    """
+    self.__station = station
+
+  def set_weight( self, w ) :
+    """
+      Set the statistical weight of the hit
+    """
+    self.__weight = w
+
+
+  def get_x( self ) :
+    """
+      Get the x position
+    """
+    return self.__x
+
+  def get_y( self ) :
+    """
+      Get the y position
+    """
+    return self.__y
+
+  def get_z( self ) :
+    """
+      Get the z position
+    """
+    return self.__z
+
+  def get_time( self ) :
+    """
+      Get the time
+    """
+    return self.__time
+
+  def get_r( self ) :
+    """
+      Get the position radius
+    """
+    return math.sqrt( self.__x**2 + self.__y**2 )
+
+  def get_phi( self ) :
+    """
+      Get the Phi angle (Position)
+    """
+    return math.atan2( self.__y, self.__x )
+
+  def get_station( self ) :
+    """
+      Get the particle ID
+    """
+    return self.__station
+
+
+  def get_weight( self ) :
+    """
+      Return the statistical weight of the hit
+    """
+    return self.__weight
+
+
+  def get( self, string ) :
+    """
+      Return the value of the variable described in the string
+    """
+    if not string in AnalysisSpacePoint.__hit_get_keys :
+      return None
+    else :
+      return AnalysisSpacePoint.__hit_get_variables[ string ]( self )
+
+
+  def get_as_vector( self ) :
+    """
+      Returns the 6D phase-space vector corresponding to the particle.
+      [t, x, y, z]
+    """
+    return [self.__time, self.__x, self.__y, self.__z]
+
+
+  def get_variable_list() :
+    """
+        Return the list of variables that the user can request
+    """
+    return AnalysisSpacePoint.__hit_get_keys
+  get_variable_list = staticmethod( get_variable_list )
+
+
+
+  __hit_get_keys = []
+  __hit_get_variables = { 'x':get_x, 'y':get_y, 'z':get_z,
+                              't':get_time, 'r':get_r, 'station':get_station }
+
+  for key, val in __hit_get_variables.iteritems() :
+    __hit_get_keys.append( key )
 
=== modified file 'src/common_py/analysis/inspectors.py'
--- src/common_py/analysis/inspectors.py	2017-08-15 16:19:42 +0000
+++ src/common_py/analysis/inspectors.py	2018-05-09 09:31:35 +0000
@@ -82,13 +82,13 @@
     
     self.emittance_plot = ROOT.TH1F(\
               'inspected_emittance_{0}'.format(self.plane), 'Emittance', \
-                                                               1000, 0.0, 20.0)
+                                                               5000, 0.0, 20.0)
     self.emittance_x_plot = ROOT.TH1F(\
            'inspected_emittance_x_{0}'.format(self.plane), 'X Emittance', \
-                                                            1000, 0.0, 20.0)
+                                                            5000, 0.0, 20.0)
     self.emittance_y_plot = ROOT.TH1F(\
          'inspected_emittance_y_{0}'.format(self.plane), 'Y Emittance', \
-                                                            1000, 0.0, 20.0)
+                                                            5000, 0.0, 20.0)
     self.alpha_plot = ROOT.TH1F(\
                'inspected_alpha_{0}'.format(self.plane), 'Alpha Function', \
                                                              2000, -5.0, 5.0)
@@ -109,11 +109,11 @@
                                                              3000, 0.0, 3000.0)
     self.total_momentum_plot = ROOT.TH1F(\
               'inspected_total_momentum_{0}'.format(self.plane), 'Momentum', \
-                                                              1000, 0.0, 500.0)
+                                                              2000, 0.0, 500.0)
 
     self.amplitude_plot = ROOT.TH1F(\
             'single_particle_amplitudes_{0}'.format(self.plane), 'Amplitude', \
-                                                              200, 0.0, 100.0)
+                                                              500, 0.0, 100.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 )
@@ -170,15 +170,15 @@
     if self.ensemble_size > 1 :
       self.ensemble_covariance.add_hit(hit)
 
-      if self.covariance.length() > self.ensemble_size :
+      if self.ensemble_covariance.length() > self.ensemble_size :
         self.fill_plots()
 
 
   def fill_plots(self) :
-    if ( self.covariance.length() < (self.ensemble_size / 2) ) or self.ensemble_size < 2 :
+    if ( self.ensemble_covariance.length() < (self.ensemble_size / 2) ) or self.ensemble_size < 2 :
       return
 
-    self.emittance_plot.Fill(self.covariance.get_emittance(\
+    self.emittance_plot.Fill(self.ensemble_covariance.get_emittance(\
                                                           ['x','px','y','py']))
     self.emittance_x_plot.Fill(self.ensemble_covariance.get_emittance(['x','px']))
     self.emittance_y_plot.Fill(self.ensemble_covariance.get_emittance(['y','py']))
@@ -446,3 +446,86 @@
 
     return ins_data
 
+
+####################################################################################################
+
+
+class EmittanceSystematicErrorInspector() :
+
+  def __init__(self, plane_id, ensemble_size) :
+    if ensemble_size < 2 :
+      raise ValueError("Ensemble size is too small.")
+
+    self.plane = plane_id
+    self.ensemble_size = ensemble_size
+    self.covariance = analysis.covariances.CovarianceMatrix()
+    self.true_covariance = analysis.covariances.CovarianceMatrix()
+    self.systematic_error_calculator = analysis.systematic_error_estimation.chi_squared_estimator()
+    self.true_emittances = []
+
+
+  def add_hit(self, hit, true_hit) :
+    self.covariance.add_hit(hit)
+    self.true_covariance.add_hit(true_hit)
+
+    if self.covariance.length() > self.ensemble_size :
+      self.fill_calculator()
+
+
+  def fill_calculator(self) :
+    if ( self.covariance.length() < (self.ensemble_size / 2) ) :
+      return
+
+    emittance = self.covariance.get_emittance(['x','px','y','py'])
+    true_emittance = self.true_covariance.get_emittance(['x','px','y','py'])
+
+    self.true_emittances.append( ( true_emittance, true_emittance / math.sqrt( 2.0*(self.true_covariance.length() - 1) ) ) )
+    self.systematic_error_calculator.add_value( emittance, emittance / math.sqrt( 2.0*(self.covariance.length() - 1) ) )
+
+    self.covariance.clear()
+    self.true_covariance.clear()
+
+
+  def get_current_ensemble_size(self) :
+    return self.covariance.length()
+
+
+  def get_systematic_error(self) :
+    num_measurements = len(self.systematic_error_calculator)
+
+    if num_measurements > 0 :
+      numer = 0.0
+      denom = 0.0
+      for meas, weight in self.true_emittances :
+        numer += meas / (weight**2)
+        denom += 1.0 / (weight**2)
+      true_emittance = numer / denom
+
+      bias = self.systematic_error_calculator.calculate_mean() - true_emittance
+      error = math.sqrt(self.systematic_error_calculator.calculate_variance()/num_measurements)
+      bias_lower = bias - error
+      bias_upper = bias + error
+
+      result = self.systematic_error_calculator.calculate_systematic_error()
+      lower, upper = self.systematic_error_calculator.calculate_bounds( result )
+    else :
+      bias = 0.0
+      error = 0.0
+      bias_lower = 0.0
+      bias_upper = 0.0
+
+      result = 0.0
+      lower = 0.0
+      upper = 0.0
+
+    return bias, bias_lower, bias_upper, result, lower, upper
+
+
+  def get_data_dictionary(self) :
+    return {}
+
+
+  def get_plot_dictionary(self) :
+    return self.systematic_error_calculator.get_plots()
+
+
 
=== added file 'src/common_py/analysis/systematic_error_estimation.py'
--- src/common_py/analysis/systematic_error_estimation.py	1970-01-01 00:00:00 +0000
+++ src/common_py/analysis/systematic_error_estimation.py	2018-05-09 09:31:35 +0000
@@ -0,0 +1,247 @@
+#!/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/>.
+#
+
+"""
+  Some classes for the calculation of systematic errors.
+
+  Inspiration is from:
+    L Lyons 1992 J. Phys. A: Math. Gen. 25 1967
+"""
+
+# pylint: disable = W0311, E1101, W0102, R0902, C0103, R0904
+
+import math
+import ROOT
+
+import scipy
+import scipy.optimize
+
+
+class chi_squared_estimator(object) :
+  """
+    Performs a chi-squared estiamtion of the systematic error based on pairs of values + errors.
+  """
+
+  def __init__(self, unique_name="chisq_error_estimation") :
+    """
+      Initialise empty object.
+      Option to set the convergence criteria for the chisq-minimisation routine
+    """
+    self.__unique_name = unique_name
+    self.__measurements = []
+    self.__errors = []
+    self.__systematic_error = 0.0
+    self.__mean = 0.0
+
+
+  def __len__(self) :
+    """
+      Return the number of measurements
+    """
+    return len(self.__measurements)
+
+
+  def add_value(self, measurement, standard_error) :
+    """
+      Add a measurement and its standard error to the object
+    """
+    self.__measurements.append(measurement)
+    self.__errors.append(standard_error)
+
+
+  def _calculate_chisquared(self) :
+    """
+      Private function: calculate the chisquared of the measurements
+    """
+    total = 0.0
+    s_sq = self.__systematic_error**2
+
+    for meas, error in zip(self.__measurements, self.__errors) :
+      total += (meas - self.__mean)**2 / (error*error + s_sq)
+
+    return total
+
+
+  def _calculate_mean(self) :
+    """
+      Private function: calculate the weighted mean of the measurements
+    """
+    numerator = 0.0
+    denominator = 0.0
+#    s_sq = self.__systematic_error**2
+
+    for meas, error in zip(self.__measurements, self.__errors) :
+      numerator += meas / (error*error)
+      denominator += 1.0 / (error*error)
+
+    return numerator / denominator
+
+
+  def calculate_mean(self) :
+    """
+      Calculate the weighted mean of the measurements assuming no systematic error.
+    """
+    self.__systematic_error = 0.0
+    return self._calculate_mean()
+
+
+  def calculate_variance(self) :
+    """
+      Calculate the variance of the measurements assuming no systematic error.
+    """
+    mean = self.calculate_mean()
+    summation = 0.0
+    for meas in self.__measurements :
+      summation += (meas - mean)**2
+
+    return summation / len(self.__measurements)
+
+
+  def _constant_error_systematic(self, get_mean, NDOF, stat_error) :
+    """
+      Private function: Calculate the systematic error assuming every measurement has the
+      same statistical error, and the mean is already assumed.
+
+      This permits a simpler equation to solve.
+    """
+    self.__mean = get_mean()
+    total = 0.0
+
+    for meas in self.__measurements :
+      total += (meas - self.__mean)**2 
+
+    diff = ( total / NDOF ) - stat_error**2
+
+    if diff <= 0.0 :
+      return 0.0
+    else :
+      return math.sqrt(diff)
+
+
+  def calculate_bounds(self, systematic_estimate, mean=None) :
+    """
+      Estimate an asymmetric error bar based on the chisquared distribution.
+
+      Changes the estimate of the systematic error until the delta chi-squared reaches the
+      equivalent of a 1-sigma deviation (central 68%).
+    """
+    lower_bar = -1.0e300 # Ridonkulous starting values
+    upper_bar =  1.0e300
+
+    get_mean = None
+
+    if mean is None : # Estimate mean from the measurements
+      get_mean = self._calculate_mean
+      NDOF = len(self.__measurements) - 1
+
+    else : # If the mean is specified we gain a degree of freedom.
+      get_mean = lambda : mean
+      NDOF = len(self.__measurements)
+
+
+    bounds = [(0.0, None)]
+
+    limit = scipy.stats.chi2.ppf(0.16, NDOF)
+    upper_bar = scipy.optimize.minimize(self._minimize_me, [systematic_estimate], args=(limit, get_mean), method="TNC", bounds=bounds)
+
+    limit = scipy.stats.chi2.ppf(0.84, NDOF)
+    lower_bar = scipy.optimize.minimize(self._minimize_me, [systematic_estimate], args=(limit, get_mean), method="TNC", bounds=bounds)
+
+    return lower_bar.x[0], upper_bar.x[0]
+
+
+  def calculate_systematic_error(self, mean=None, error_assumption=None) :
+    """
+      The main function - returns the systematic error given the data and measurements stored in
+      the class.
+
+      Specifying the mean stops the mean from being calculated on the fly from the data, this
+      slightly changes th behaviour - the mean is no longer a function of the systematic error and
+      the number of degrees of freedom is different.
+
+      Specifying the error_assumption, assumes that every measurement has the same (specified) 
+      error. This allows for a faster calculation of the systematic error.
+
+      Specifying both allows the for systematic error to be calculated trivially.
+    """
+    get_mean = None
+    NDOF = 0
+    self.__systematic_error = 0.0
+
+    if mean is None : # Estimate mean from the measurements
+      get_mean = self._calculate_mean
+      NDOF = len(self.__measurements) - 1
+
+    else : # If the mean is specified we gain a degree of freedom.
+      get_mean = lambda : mean
+      NDOF = len(self.__measurements)
+
+
+      if error_assumption is not None : # Assume all statistical errors are the 
+                                        # same. Shortcut!
+        return self._constant_error_systematic(get_mean, NDOF, error_assumption)
+
+
+    bounds = [(0.0, None)]
+    estimate = 0.0
+    result = scipy.optimize.minimize(self._minimize_me, [estimate], args=(NDOF, get_mean), method="TNC", bounds=bounds)
+
+    self.__systematic_error = result.x[0]
+
+    return self.__systematic_error
+
+
+  def _minimize_me(self, parameters, ndof, get_mean) :
+    """
+      Minimize function for SciPy.Optimize.
+  
+      Calculates the squared residual of the Chi-squared value and number of degrees of freedom.
+    """
+    self.__systematic_error = parameters[0]
+    self.__mean = get_mean()
+    chisq = self._calculate_chisquared()
+    return (chisq - ndof)**2
+
+
+  def get_plots(self) :
+    """
+      Return a dictionary containing histograms of the data that was stored.
+    """
+    if len(self.__measurements) == 0 :
+      return {}
+
+    min_val = min(self.__measurements)
+    max_val = max(self.__measurements)
+    total_weight = 0.0
+
+    for error in self.__errors :
+      total_weight += 1.0 / (error**2)
+
+    hist = ROOT.TH1F(self.__unique_name, "", 100, min_val, max_val)
+    weighted_hist = ROOT.TH1F(self.__unique_name+'-weighted', "", 100, min_val, max_val)
+
+    for value, error in zip( self.__measurements, self.__errors ) :
+      hist.Fill(value)
+      weighted_hist.Fill(value, (1.0/(error**2))/total_weight)
+
+    plots = {}
+    plots['measurements'] = hist
+    plots['weighted_measurements'] = weighted_hist
+
+    return plots
+
 
=== modified file 'src/common_py/framework/single_thread.py'
--- src/common_py/framework/single_thread.py	2018-02-14 14:17:17 +0000
+++ src/common_py/framework/single_thread.py	2018-05-09 09:31:35 +0000
@@ -201,15 +201,16 @@
         # done with tranform-merge, now write it out
         self.outputer.save(redevent)
         try:
-            del redevent
-        except:
-            pass
-        # if we converted to a different representation, delete the old one
-        try:
             maus_cpp.converter.del_data_repr(event)
+        except TypeError: # pylint: disable = W0702
+            print "Failed to delete MAUS::Data object."
+        try:
             maus_cpp.converter.del_data_repr(redevent)
-        except: # pylint: disable = W0702
-            pass
+        except TypeError: # pylint: disable = W0702
+            try: # If we returned a different datatype, try to delete it anyway.
+                del redevent
+            except:
+                print "Failed to delete Reduced Data object."
 
     def start_of_run(self, new_run_number):
         """