66# rootpy
77from rootpy .io import File
88from rootpy .plotting import Hist2D
9+ from rootpy .matrix import Matrix
10+ from rootpy import asrootpy
11+ import root_numpy as rnp
12+ import numpy as np
13+ from ROOT import TH2F
914# DailyPythonScripts
1015import config .RooUnfold as unfoldCfg
1116from config .variable_binning import bin_widths , bin_edges
1217from config import XSectionConfig
1318from tools .Calculation import calculate_xsection , calculate_normalised_xsection , \
14- combine_complex_results
19+ combine_complex_results , calculate_covariance_for_normalised_xsection
1520from tools .hist_utilities import hist_to_value_error_tuplelist , \
1621value_error_tuplelist_to_hist
1722from tools .Unfolding import Unfolding , get_unfold_histogram_tuple
@@ -32,7 +37,20 @@ def unfold_results( results, category, channel, k_value, h_truth, h_measured, h_
3237 unfoldCfg .Hreco = options .Hreco
3338
3439 h_unfolded_data = unfolding .unfold ( h_data )
35-
40+
41+ # unfolding.unfoldObject.GetCov()
42+
43+ covariance_matrix = None
44+ if category == 'central' :
45+ covariance_matrix = asrootpy ( unfolding .unfoldObject .Ereco (options .Hreco ) ).to_numpy ()
46+ # print list( h_unfolded_data.y() )
47+ # unfolding.unfoldObject.ErecoV(options.Hreco).Draw()
48+ # raw_input()
49+ # cov_matrix =
50+ # cov_matrix.Draw('COLZ text')
51+ # print cov_matrix
52+ # raw_input('...')
53+
3654 if options .write_unfolding_objects :
3755 # export the D and SV distributions
3856 SVD_path = path_to_JSON + '/unfolding_objects/' + channel + '/kv_' + str ( k_value ) + '/'
@@ -73,7 +91,7 @@ def unfold_results( results, category, channel, k_value, h_truth, h_measured, h_
7391 unfoldingObjectFile .Close ()
7492
7593 del unfolding
76- return hist_to_value_error_tuplelist ( h_unfolded_data )
94+ return hist_to_value_error_tuplelist ( h_unfolded_data ), covariance_matrix
7795
7896def data_covariance_matrix ( data ):
7997 values = list ( data )
@@ -265,16 +283,16 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ):
265283 scaledown_results = hist_to_value_error_tuplelist ( h_truth_scaledown )
266284 scaleup_results = hist_to_value_error_tuplelist ( h_truth_scaleup )
267285
268- TTJet_fit_results_unfolded = unfold_results ( TTJet_fit_results ,
269- category ,
270- channel ,
271- k_value ,
272- h_truth ,
273- h_measured ,
274- h_response ,
275- h_fakes ,
276- method
277- )
286+ TTJet_fit_results_unfolded , covariance_matrix = unfold_results ( TTJet_fit_results ,
287+ category ,
288+ channel ,
289+ k_value ,
290+ h_truth ,
291+ h_measured ,
292+ h_response ,
293+ h_fakes ,
294+ method
295+ )
278296
279297 normalisation_unfolded = {
280298 'TTJet_measured' : TTJet_fit_results ,
@@ -295,7 +313,9 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ):
295313 normalisation_unfolded ['powheg_v1_herwig' ] = powheg_v1_herwig_results
296314 normalisation_unfolded ['powheg_v1_pythia' ] = powheg_v1_pythia_results
297315
298- return normalisation_unfolded
316+ return normalisation_unfolded , covariance_matrix
317+
318+
299319
300320def calculate_xsections ( normalisation , category , channel , k_value = None ):
301321 global variable , met_type , path_to_JSON
@@ -396,7 +416,7 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value =
396416 write_data_to_JSON ( normalised_xsection , filename )
397417
398418if __name__ == '__main__' :
399- set_root_defaults ( msg_ignore_level = 3001 )
419+ set_root_defaults ( set_batch = False , msg_ignore_level = 3001 )
400420 # setup
401421 parser = OptionParser ()
402422 parser .add_option ( "-p" , "--path" , dest = "path" , default = 'data/' ,
@@ -579,11 +599,15 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value =
579599 filename = ''
580600
581601 # get unfolded normalisation
582- unfolded_normalisation_electron = get_unfolded_normalisation ( TTJet_fit_results_electron , category , 'electron' , k_value_electron )
583- unfolded_normalisation_muon = get_unfolded_normalisation ( TTJet_fit_results_muon , category , 'muon' , k_value_muon )
602+ unfolded_normalisation_electron , covariance_electron = get_unfolded_normalisation ( TTJet_fit_results_electron , category , 'electron' , k_value_electron )
603+ unfolded_normalisation_muon , covariance_muon = get_unfolded_normalisation ( TTJet_fit_results_muon , category , 'muon' , k_value_muon )
604+ covariance_combined = None
584605 if combine_before_unfolding :
585- unfolded_normalisation_combined = get_unfolded_normalisation ( TTJet_fit_results_combined , category , 'combined' , k_value_combined )
606+ unfolded_normalisation_combined , covariance_combined = get_unfolded_normalisation ( TTJet_fit_results_combined , category , 'combined' , k_value_combined )
586607 else :
608+ covariance_combined = covariance_electron
609+ covariance_combined += covariance_muon
610+ # covariance_combined = asrootpy( covariance_combined )
587611 unfolded_normalisation_combined = combine_complex_results ( unfolded_normalisation_electron , unfolded_normalisation_muon )
588612
589613 filename = path_to_JSON + '/xsection_measurement_results/electron/kv%d/%s/normalisation_%s.txt' % ( k_value_electron_central , category , met_type )
@@ -592,6 +616,9 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value =
592616 write_data_to_JSON ( unfolded_normalisation_muon , filename )
593617 filename = path_to_JSON + '/xsection_measurement_results/combined/%s/normalisation_%s.txt' % ( category , met_type )
594618 write_data_to_JSON ( unfolded_normalisation_combined , filename )
619+
620+
621+
595622
596623# if measurement_config.include_higgs:
597624# # now the same for the Higgs
@@ -613,15 +640,44 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value =
613640# write_data_to_JSON( unfolded_normalisation_combined_higgs, filename )
614641
615642 # measure xsection
616- calculate_xsections ( unfolded_normalisation_electron , category , 'electron' , k_value_electron_central )
617- calculate_xsections ( unfolded_normalisation_muon , category , 'muon' , k_value_muon_central )
618- calculate_xsections ( unfolded_normalisation_combined , category , 'combined' )
643+ # calculate_xsections( unfolded_normalisation_electron, category, 'electron', k_value_electron_central )
644+ # calculate_xsections( unfolded_normalisation_muon, category, 'muon', k_value_muon_central )
645+ # calculate_xsections( unfolded_normalisation_combined, category, 'combined' )
619646
620- calculate_normalised_xsections ( unfolded_normalisation_electron , category , 'electron' , k_value_electron_central )
621- calculate_normalised_xsections ( unfolded_normalisation_muon , category , 'muon' , k_value_muon_central )
622- calculate_normalised_xsections ( unfolded_normalisation_combined , category , 'combined' )
647+ # calculate_normalised_xsections( unfolded_normalisation_electron, category, 'electron', k_value_electron_central )
648+ # calculate_normalised_xsections( unfolded_normalisation_muon, category, 'muon', k_value_muon_central )
649+ # calculate_normalised_xsections( unfolded_normalisation_combined, category, 'combined' )
623650
624- normalise_to_one = True
625- calculate_normalised_xsections ( unfolded_normalisation_electron , category , 'electron' , k_value_electron_central , normalise_to_one )
626- calculate_normalised_xsections ( unfolded_normalisation_muon , category , 'muon' , k_value_muon_central , normalise_to_one )
627- calculate_normalised_xsections ( unfolded_normalisation_combined , category , 'combined' , normalise_to_one )
651+ # normalise_to_one = True
652+ # calculate_normalised_xsections( unfolded_normalisation_electron, category, 'electron', k_value_electron_central, normalise_to_one )
653+ # calculate_normalised_xsections( unfolded_normalisation_muon, category, 'muon', k_value_muon_central, normalise_to_one )
654+ # calculate_normalised_xsections( unfolded_normalisation_combined, category, 'combined', normalise_to_one )
655+
656+ if category == 'central' :
657+ normalised_xsection_covariance_electron = calculate_covariance_for_normalised_xsection ( covariance_electron , unfolded_normalisation_electron ['TTJet_unfolded' ], bin_widths [variable ] )
658+ normalised_xsection_covariance_muon = calculate_covariance_for_normalised_xsection ( covariance_muon , unfolded_normalisation_muon ['TTJet_unfolded' ], bin_widths [variable ] )
659+ normalised_xsection_covariance_combined = calculate_covariance_for_normalised_xsection ( covariance_combined , unfolded_normalisation_combined ['TTJet_unfolded' ], bin_widths [variable ] )
660+
661+ # Write covariance matrices to file
662+ filename = path_to_JSON + '/xsection_measurement_results/electron/kv%d/%s/covariance.txt' % ( k_value_electron_central , category )
663+ np .savetxt ( filename , normalised_xsection_covariance_electron , delimiter = ',' )
664+
665+ # cov_file = File( filename, 'recreate' )
666+ # covariance_electron.Write()
667+ # normalised_xsection_covariance_electron.(filename,',')
668+ # cov_file.Close()
669+
670+ filename = path_to_JSON + '/xsection_measurement_results/muon/kv%d/%s/covariance.txt' % ( k_value_muon_central , category )
671+ np .savetxt ( filename , normalised_xsection_covariance_muon , delimiter = ',' )
672+
673+ # cov_file = File( filename, 'recreate' )
674+ # covariance_muon.Write()
675+ # normalised_xsection_covariance_muon.Write()
676+ # cov_file.Close()
677+
678+ filename = path_to_JSON + '/xsection_measurement_results/combined/%s/covariance.txt' % ( category )
679+ np .savetxt ( filename , normalised_xsection_covariance_combined , delimiter = ',' )
680+ # cov_file = File( filename, 'recreate' )
681+ # covariance_combined.Write()
682+ # normalised_xsection_covariance_combined.Write()
683+ # cov_file.Close()
0 commit comments