Merge pull request #86 from senkin/unfolding-tests

Unfolding tests

.travis.yml

@@ -66,4 +66,5 @@ script:
   - time python test/TMinuitFit.py
   - time python test/tools_Calculation.py
   - time python test/tools_hist_utilities.py
+  - time python test/tools_Unfolding.py
   - time python test/Integral_GetBinContent_consistency.py

src/unfolding_tests/unfolding.py → legacy/unfolding_test.py

@@ -1,5 +1,7 @@
 from rootpy.io import File
 from rootpy.plotting import Hist, Canvas
+import matplotlib
+matplotlib.use('agg')
 import rootpy.plotting.root2matplotlib as rplt
 import matplotlib.pyplot as plt
 from rootpy import asrootpy
@@ -67,7 +69,7 @@ def checkOnMC(unfolding, method):
     global bins, nbins
     RooUnfold.SVD_n_toy = 1000
     pulls = []
-    for sub in range(2,9):
+    for sub in range(1,9):
         inputFile2 = File('../data/unfolding_merged_sub%d.root' % sub, 'read')
         h_data = asrootpy(inputFile2.unfoldingAnalyserElectronChannel.measured.Rebin(nbins, 'measured', bins))
         nEvents = inputFile2.EventFilter.EventCounter.GetBinContent(1)
@@ -87,7 +89,7 @@ def checkOnMC(unfolding, method):
     filling = h_allpulls.Fill
     for entry in allpulls:
         filling(entry)
-    fit = h_allpulls.Fit('gaus', 'WW')
+    fit = h_allpulls.Fit('gaus', 'WWS')
     h_fit = asrootpy(h_allpulls.GetFunction("gaus").GetHistogram())
     canvas = Canvas(width=1600, height=1000)
     canvas.SetLeftMargin(0.15)
@@ -160,7 +162,7 @@ def doUnfoldingSequence(unfolding, h_data, method, outputfile_suffix = '', doClo
     h_measured = asrootpy(inputFile.unfoldingAnalyserElectronChannel.measured.Rebin(nbins, 'measured', bins))
     h_fakes = asrootpy(inputFile.unfoldingAnalyserElectronChannel.fake.Rebin(nbins, 'fake', bins))
     h_response = inputFile.unfoldingAnalyserElectronChannel.response_withoutFakes_AsymBins #response_AsymBins
-    h_measured_new = h_measured - h_fakes
+    # h_measured_new = h_measured - h_fakes
 
 #    h_response = inputFile.unfoldingAnalyserElectronChannel.response_AsymBins #response_AsymBins
     nEvents = inputFile.EventFilter.EventCounter.GetBinContent(1)
@@ -169,7 +171,7 @@ def doUnfoldingSequence(unfolding, h_data, method, outputfile_suffix = '', doClo
     h_measured.Scale(lumiweight)
     h_fakes.Scale(lumiweight)
     h_response.Scale(lumiweight)
-    unfolding = Unfolding(h_truth, h_measured_new, h_response, method = method)
+    unfolding = Unfolding(h_truth, h_measured, h_response, method = method)
     #should be identical to
 #    unfolding = Unfolding(h_truth, h_measured, h_response, h_fakes, method = method)
 

setup_standalone.sh

@@ -32,8 +32,6 @@ echo "Installing argparse"
 pip install -U argparse
 echo "Installing pudb <-- interactive debugging"
 pip install -U pudb
-echo "Installing pypm"
-pip install -U pypm
 echo "Installing numpy"
 pip install -U numpy
 echo "Installing dateutil"

src/unfolding_tests/compare_unfolding_parameters.py

@@ -23,10 +23,10 @@
 from config import latex_labels
 from config.cross_section_measurement_common import translate_options
 from config.variable_binning import bin_widths, bin_edges
+from tools.plotting import Histogram_properties, compare_measurements
 from tau_value_determination import get_tau_from_global_correlation
 from tau_value_determination import get_tau_from_L_shape, get_data_histogram
 from tools.Unfolding import get_unfold_histogram_tuple, Unfolding
-from tools.plotting import Histogram_properties, compare_measurements
 from tools.Calculation import calculate_normalised_xsection
 from tools.hist_utilities import hist_to_value_error_tuplelist
 from tools.hist_utilities import value_error_tuplelist_to_hist, spread_x
@@ -198,8 +198,6 @@ def compare( central_mc, expected_result = None, measured_result = None, results
 
     input_filename_central = measurement_config.unfolding_madgraph
     input_filename_bias = measurement_config.unfolding_mcatnlo
-    input_filename_central = '/Users/phzss/work/workspace_juno/git/DailyPythonScripts/unfolding_merged.root'
-    input_filename_bias = '/Users/phzss/work/workspace_juno/git/DailyPythonScripts/unfolding_TTJets_8TeV_mcatnlo.root'
 
     variables = ['MET', 'WPT', 'MT', 'ST', 'HT']
 

src/unfolding_tests/k_value_determination.py

@@ -15,6 +15,7 @@
 from optparse import OptionParser
 from rootpy.io import File
 import matplotlib
+matplotlib.use('agg')
 
 from uncertainties import ufloat
 import numpy
@@ -36,7 +37,6 @@
 from config.cross_section_measurement_common import translate_options
 
 
-matplotlib.use('agg')
 matplotlib.rc('font',**CMS.font)
 matplotlib.rc('text', usetex = True)
 
@@ -136,7 +136,7 @@ def draw_d_i( d_i ):
     set_root_defaults()
 
     parser = OptionParser()
-    parser.add_option("-p", "--path", dest="path", default='../cross_section_measurement/data/',
+    parser.add_option("-p", "--path", dest="path", default='../../data/',
                       help="set path to JSON files")
     parser.add_option("-o", "--output_folder", dest = "output_folder", default = 'plots_k_values/',
                       help = "set path to save plots" )

test/tools_Calculation.py

@@ -139,7 +139,7 @@ def test_decombine_result_multiple_backgrounds(self):
 
         N_total_prime = N_total * 2
         N_signal_plus_bkg_2_prime, N_background_1_prime = decombine_result((N_total_prime, 0), ratio_signal_bkg_1)
-        N_signal_prime, N_background_2_prime = decombine_result((N_signal_plus_bkg_2_prime, 0), ratio_signal_bkg_2)
+        N_signal_prime, N_background_2_prime = decombine_result(N_signal_plus_bkg_2_prime, ratio_signal_bkg_2)
 
         self.assertEqual(N_signal_prime[0], N_signal * 2)
         self.assertEqual(N_background_1_prime[0], N_background_1 * 2)

test/tools_Unfolding.py

@@ -0,0 +1,71 @@
+'''
+Created on 15 May 2014
+
+@author: senkin
+'''
+from __future__ import division
+import unittest
+import importlib
+from rootpy.io import File
+from tools.Unfolding import Unfolding, get_unfold_histogram_tuple
+from tools.hist_utilities import hist_to_value_error_tuplelist, value_error_tuplelist_to_hist
+from tools.ROOT_utililities import set_root_defaults
+from tools.file_utilities import read_data_from_JSON
+from config.variable_binning import bin_edges
+
+class Test( unittest.TestCase ):
+
+    def setUp( self ):
+        # load histograms
+        self.input_file = File('data/unfolding_merged_asymmetric.root')
+        self.k_value = 3
+        self.unfold_method = 'RooUnfoldSvd'
+        self.met_type = 'patType1CorrectedPFMet'
+        self.variables = ['MET', 'WPT', 'MT' , 'ST', 'HT']
+        self.channels = ['electron', 'muon', 'combined']
+        self.dict = {}
+        for channel in self.channels:
+            self.dict[channel] = {}
+            for variable in self.variables:
+                self.dict[variable] = {}
+                h_truth, h_measured, h_response, _ = get_unfold_histogram_tuple(
+                                                            inputfile = self.input_file,
+                                                            variable = variable,
+                                                            channel = channel,
+                                                            met_type = self.met_type)
+
+                unfolding_object = Unfolding( h_truth,
+                                       h_measured,
+                                       h_response,
+                                       k_value = self.k_value,
+                                       method = self.unfold_method
+                                       )
+
+                self.dict[channel][variable] = {'h_truth' : h_truth,
+                                                'h_measured' : h_measured,
+                                                'h_response' : h_response,
+                                                'unfolding_object' : unfolding_object
+                                                }
+
+    def tearDown( self ):
+        pass
+
+    def test_closure( self ):
+        for channel in self.channels:
+            for variable in self.variables:
+                # closure test
+                unfolded_result = hist_to_value_error_tuplelist( self.dict[channel][variable]['unfolding_object'].closureTest() )
+                truth = hist_to_value_error_tuplelist( self.dict[channel][variable]['h_truth'] )
+                # the difference between the truth and unfolded result should be within the unfolding error
+                for (value, error), (true_value, _) in zip(unfolded_result, truth):
+                    self.assertAlmostEquals(value, true_value, delta = error)
+
+    def test_invalid_zero_data( self ):
+        variable = 'MET'
+        channel = 'electron'
+        pseudo_data = value_error_tuplelist_to_hist( [(0,0)]*( len( bin_edges[variable] ) - 1 ), bin_edges[variable] )
+        self.assertRaises(ValueError,  self.dict[channel][variable]['unfolding_object'].unfold, (pseudo_data))
+
+if __name__ == "__main__":
+    set_root_defaults()
+    unittest.main()

tools/Unfolding.py

@@ -146,6 +146,9 @@ def test_regularisation ( self, data, k_max ):
             raise ValueError( 'Unfolding method "%s" is not supported for regularisation parameter tests. Please use RooUnfoldSvd.' % ( self.method ) )
 
     def unfold( self, data ):
+        have_zeros = [value == 0 for value,_ in hist_to_value_error_tuplelist( data )]
+        if not False in have_zeros:
+            raise ValueError('Data histograms contains only zeros')
         self.setup_unfolding( data )
         if self.method == 'TSVDUnfold' or self.method == 'TopSVDUnfold':
             self.unfolded_data = asrootpy( self.unfoldObject.Unfold( self.k_value ) )
@@ -235,7 +238,7 @@ def pull_inputErrorOnly( self ):
         if self.unfolded_data and self.truth:
             # set unfolded_data errors to stat errors from data
             temp = self.unfolded_data.Clone()
-            temp_list = list( temp )
+            temp_list = list( temp.y() )
 #            data_list = list(self.data)
             unfolded_errors = self.get_unfolded_data_errors()
             for bin_i in range( len( temp_list ) ):
@@ -247,7 +250,7 @@ def pull_inputErrorOnly( self ):
 
             diff = temp - temp_truth
             errors = []
-            values = list( diff )
+            values = list( diff.y() )
             for bin_i in range( len( values ) ):
                 errors.append( diff.GetBinError( bin_i + 1 ) )
             result = [value / error for value, error in zip( values, errors )]
@@ -269,7 +272,7 @@ def get_unfold_histogram_tuple(
                 inputfile,
                 variable,
                 channel,
-                met_type,
+                met_type = 'patType1CorrectedPFMet',
                 centre_of_mass = 8,
                 ttbar_xsection = 245.8,
                 luminosity = 19712,