diff --git a/bin/plot b/bin/plot index 60240453..34807cdd 100755 --- a/bin/plot +++ b/bin/plot @@ -84,6 +84,8 @@ from tools.ROOT_utils import set_root_defaults from tools.file_utilities import write_data_to_JSON, read_data_from_JSON from tools.HistSet import HistSet from copy import deepcopy +from tools.ROOT_utils import get_histogram_from_file +from math import sqrt supported_commands = ['compare-files', 'compare-hists'] @@ -95,7 +97,7 @@ plot_options = {} global_options = ['files', 'file-aliases', 'histograms', 'labels', 'plot_type', 'output_folder', 'output_format', 'command', 'data_index', 'normalise', 'show_ratio', 'show_stat_errors_on_mc', 'colours', 'name_prefix', - 'fill_area', 'alpha'] + 'fill_area', 'alpha', 'subtractFiles'] def main(): options, input_values_sets, json_input_files = parse_options() @@ -173,7 +175,63 @@ def _exit_( msg ): def prepare_inputs( input_values ): input_values = check_and_fix_inputs( input_values ) global files, histograms, global_options + + subtractFiles = [] + if input_values.has_key('subtractFiles'): + subtractFiles = input_values['subtractFiles'] + sets = group_by_command( input_values['command'] ) + + for s in sets: + print s.histograms + # subtractFiles_sets = getBackgroundHists( subtractFiles, ) + + # if input_values['command'] == 'compare-hists' and len( input_values['subtractFiles'] ) > 0: + + # subtract_background( input_values['subtractFiles'] ) + if len( subtractFiles ) > 0: + hist_sets = [] + i = 0 + + for h in input_values['histograms']: + print h + # name = h.replace( '/', '_' ) + # name = name.replace( ' ', '' ) + + uncDict = {'VJets' : 0.3, + 'TTJet' : 0.06, + 'SingleTop' : 0.3 } + + for f in subtractFiles: + hist = get_histogram_from_file( histogram_path = h, + input_file = f ) + print f + print hist + unc = 0 + if 'VJets' in f or 'SingleTop' in f: + unc = 0.3 + elif 'TTJet' in f: + unc = 0.3 + + for bin in range(1,hist.GetNbinsX()): + binContent = hist.GetBinContent(bin) + binError = hist.GetBinError(bin) + newError = sqrt( binError ** 2 + ( binContent * unc ) ** 2 ) + # print 'Errors :',binError,newError + hist.SetBinError( bin, newError ) + + sets[0].histograms[i] -= hist + i += 1 + # hist_sets.append( HistSet( h, hist_inputs = hist_set, output_hist_labels = labels ) ) + + + # print hist_sets + # print 'Looping sets' + # for s in sets: + # for h in s.histograms: + # print s.inputs + # print s.histograms + if input_values.has_key('name_prefix'): for s in sets: s.name = input_values['name_prefix'] + s.name @@ -234,6 +292,7 @@ def group_by_file(): [{name:{f:h}, ...] where f = file and h = histogram path. The name is taken from the file name ''' global files, file_aliases, histograms, labels + hist_sets = [] for f,name in zip(files, file_aliases): hist_set = [] diff --git a/bin/plots2tdr b/bin/plots2tdr index 8515758a..6aba32ad 100755 --- a/bin/plots2tdr +++ b/bin/plots2tdr @@ -15,24 +15,22 @@ from glob import glob plot_map_common = {} plot_map_paper = { - # : , - # results 8 TeV - '7TeV/MET/central/normalised_xsection_combined_*.pdf': '7TeV/MET/', - '7TeV/HT/central/normalised_xsection_combined_*.pdf': '7TeV/HT/', - '7TeV/ST/central/normalised_xsection_combined_*.pdf': '7TeV/ST/', - '7TeV/WPT/central/normalised_xsection_combined_*.pdf': '7TeV/WPT/', - # results 8 TeV - '8TeV/MET/central/normalised_xsection_combined_*.pdf': '8TeV/MET/', - '8TeV/HT/central/normalised_xsection_combined_*.pdf': '8TeV/HT/', - '8TeV/ST/central/normalised_xsection_combined_*.pdf': '8TeV/ST/', - '8TeV/WPT/central/normalised_xsection_combined_*.pdf': '8TeV/WPT/', - # control plots 8 TeV - 'control_plots/after_fit/8TeV/*_HT_2orMoreBtags_with_ratio.pdf': 'control/', - 'control_plots/after_fit/8TeV/*_patType1CorrectedPFMet_2orMoreBtags_with_ratio.pdf': 'control/', - # fit variables 8 TeV - 'control_plots/after_fit/8TeV/*_angle_bl_2orMoreBtags_with_ratio.pdf': 'control/fit_variables/8TeV/', - 'control_plots/after_fit/8TeV/*_AbsEta_2orMoreBtags_with_ratio.pdf': 'control/fit_variables/8TeV/', - 'control_plots/after_fit/8TeV/*_M3_2orMoreBtags_with_ratio.pdf': 'control/fit_variables/8TeV/', + # # : , + # # results 8 TeV + # 'background_subtraction/7TeV/MET/central/normalised_xsection_combined_*.pdf': '7TeV/MET/', + # 'background_subtraction/7TeV/HT/central/normalised_xsection_combined_*.pdf': '7TeV/HT/', + # 'background_subtraction/7TeV/ST/central/normalised_xsection_combined_*.pdf': '7TeV/ST/', + # 'background_subtraction/7TeV/WPT/central/normalised_xsection_combined_*.pdf': '7TeV/WPT/', + # # # results 8 TeV + # 'background_subtraction/8TeV/MET/central/normalised_xsection_combined_*.pdf': '8TeV/MET/', + # 'background_subtraction/8TeV/HT/central/normalised_xsection_combined_*.pdf': '8TeV/HT/', + # 'background_subtraction/8TeV/ST/central/normalised_xsection_combined_*.pdf': '8TeV/ST/', + # 'background_subtraction/8TeV/WPT/central/normalised_xsection_combined_*.pdf': '8TeV/WPT/', + # # control plots 8 TeV + # 'control_plots/before_fit/8TeV/*_patType1CorrectedPFMet_2orMoreBtags_with_ratio.pdf': 'control/', + # 'control_plots/before_fit/8TeV/*_*T_2orMoreBtags_with_ratio.pdf': 'control/', + + '../tables/background_subtraction/*/*/*_normalised_xsection_*combined.tex': '../tables/' } plot_map_AN = {} diff --git a/bin/x_01_all_vars b/bin/x_01_all_vars index 39a3d201..fba8de03 100755 --- a/bin/x_01_all_vars +++ b/bin/x_01_all_vars @@ -1,14 +1,17 @@ #!/bin/bash echo "This will take a while ... grab a coffee/tea/water/cocktail" mkdir -p logs -fit_var="absolute_eta,M3,angle_bl" +# fit_var="absolute_eta,M3,angle_bl" +# fit_var="absolute_eta" +# fit_var="M3" +fit_var="angle_bl" nice_fit_var=`echo $fit_var | sed 's/,/_/g'` N_JOBS=6 echo "Using the fit variable(s): $fit_var" i=0 -for var in MET HT ST WPT MT; do +for var in MET HT ST WPT; do echo "Fitting distribution: $var" nohup time python src/cross_section_measurement/01_get_fit_results.py -V -v $var --no_combined_signal -c 7 --fit-variables $fit_var --disable-constraints &> logs/01_${var}_fit_7TeV_${nice_fit_var}.log & let i+=1 diff --git a/bin/x_01b_all_vars b/bin/x_01b_all_vars index 403deaab..65817d51 100755 --- a/bin/x_01b_all_vars +++ b/bin/x_01b_all_vars @@ -1,11 +1,11 @@ #!/bin/bash echo "This will take a while ... grab a coffee/tea/water/cocktail" mkdir -p logs -N_JOBS=5 +N_JOBS=4 i=0 -for var in MET HT ST WPT MT; do +for var in MET HT ST WPT; do echo "Fitting distribution: $var" nohup time python src/cross_section_measurement/01_get_ttjet_normalisation.py -v $var -c 7 -i config/measurements/background_subtraction &> logs/01b_${var}_bgs_7TeV.log & let i+=1 diff --git a/bin/x_02_all_vars b/bin/x_02_all_vars index f23010cb..65ad197e 100755 --- a/bin/x_02_all_vars +++ b/bin/x_02_all_vars @@ -2,6 +2,9 @@ echo "This will take a while ... grab a coffee/tea/water" mkdir -p logs fit_var="absolute_eta,M3,angle_bl" +# fit_var="absolute_eta" +# fit_var="M3" +# fit_var="angle_bl" nice_fit_var=`echo $fit_var | sed 's/,/_/g'` N_JOBS=6 diff --git a/bin/x_02b_all_vars b/bin/x_02b_all_vars index d44b6f16..64f517e7 100755 --- a/bin/x_02b_all_vars +++ b/bin/x_02b_all_vars @@ -5,9 +5,9 @@ N_JOBS=5 i=0 -for var in MET HT ST WPT MT; do +for var in MET HT ST WPT; do echo "Unfolding distribution: $var" - nohup time python src/cross_section_measurement/02_unfold_and_measure.py -v $var -c 7 -p data/normalisation/background_subtraction/ &> logs/02b_${var}_unfold_7TeV.log & + nohup time python src/cross_section_measurement/02_unfold_and_measure.py -v $var -c 7 -p data/normalisation/background_subtraction/ --test &> logs/02b_${var}_unfold_7TeV.log & let i+=1 if (( $i % N_JOBS == 0 )) then @@ -16,7 +16,7 @@ for var in MET HT ST WPT MT; do fi - nohup time python src/cross_section_measurement/02_unfold_and_measure.py -v $var -c 8 -p data/normalisation/background_subtraction/ &> logs/02b_${var}_unfold_8TeV.log & + nohup time python src/cross_section_measurement/02_unfold_and_measure.py -v $var -c 8 -p data/normalisation/background_subtraction/ --test &> logs/02b_${var}_unfold_8TeV.log & let i+=1 if (( $i % N_JOBS == 0 )) then diff --git a/bin/x_03_all_vars b/bin/x_03_all_vars index 4eba0c74..b436ef13 100755 --- a/bin/x_03_all_vars +++ b/bin/x_03_all_vars @@ -2,6 +2,9 @@ echo "This will take a while ... grab a coffee/tea/water" mkdir -p logs fit_var="absolute_eta,M3,angle_bl" +# fit_var="absolute_eta" +# fit_var="M3" +# fit_var="angle_bl" nice_fit_var=`echo $fit_var | sed 's/,/_/g'` N_JOBS=6 diff --git a/bin/x_03b_all_vars b/bin/x_03b_all_vars index 8f29ad8d..98cac863 100755 --- a/bin/x_03b_all_vars +++ b/bin/x_03b_all_vars @@ -4,7 +4,7 @@ mkdir -p logs N_JOBS=5 i=0 -for var in MET HT ST WPT MT; do +for var in MET HT ST WPT; do echo "Calculating diff. x-section for distribution: $var" nohup time python src/cross_section_measurement/03_calculate_systematics.py -s -v $var -c 7 -p data/normalisation/background_subtraction/ &> logs/03b_${var}_calculate_7TeV.log & let i+=1 diff --git a/bin/x_04_all_vars b/bin/x_04_all_vars index 4b705c10..facab725 100755 --- a/bin/x_04_all_vars +++ b/bin/x_04_all_vars @@ -3,6 +3,9 @@ echo "This will take a while ... grab a coffee/tea/water" mkdir -p logs mkdir -p plots fit_var="absolute_eta,M3,angle_bl" +# fit_var="absolute_eta" +# fit_var="M3" +# fit_var="angle_bl" nice_fit_var=`echo $fit_var | sed 's/,/_/g'` N_JOBS=5 @@ -11,7 +14,7 @@ echo "Using the fit variable(s): $fit_var" i=0 for var in MET HT ST MT WPT; do echo "Plotting diff. x-section for distribution: $var" - nohup time python src/cross_section_measurement/04_make_plots_matplotlib.py -v $var -c 7 -p data/$nice_fit_var -o plots/fit &> logs/04_${var}_plot_7TeV_${nice_fit_var}.log & # -a -s <--add -a option for additional plots, -s option to include the final results from background subtraction as well as from fitting + nohup time python src/cross_section_measurement/04_make_plots_matplotlib.py -v $var -c 7 -p data/$nice_fit_var -o plots/fit/ &> logs/04_${var}_plot_7TeV_${nice_fit_var}.log & # -a -s <--add -a option for additional plots, -s option to include the final results from background subtraction as well as from fitting let i+=1 if (( $i % N_JOBS == 0 )) then @@ -19,7 +22,7 @@ for var in MET HT ST MT WPT; do wait; fi - nohup time python src/cross_section_measurement/04_make_plots_matplotlib.py -v $var -c 8 -p data/$nice_fit_var -o plots/fit &> logs/04_${var}_plot_8TeV_${nice_fit_var}.log & # -a -s <--add -a option for additional plots, -s option to include the final results from background subtraction as well as from fitting + nohup time python src/cross_section_measurement/04_make_plots_matplotlib.py -v $var -c 8 -p data/$nice_fit_var -o plots/fit/ &> logs/04_${var}_plot_8TeV_${nice_fit_var}.log & # -a -s <--add -a option for additional plots, -s option to include the final results from background subtraction as well as from fitting let i+=1 if (( $i % N_JOBS == 0 )) then diff --git a/bin/x_04b_all_vars b/bin/x_04b_all_vars index 5e5013ca..b1fa0367 100755 --- a/bin/x_04b_all_vars +++ b/bin/x_04b_all_vars @@ -6,7 +6,7 @@ N_JOBS=5 i=0 echo "Visible phase space" -for var in MET HT ST WPT MT; do +for var in MET HT ST WPT; do echo "Plotting diff. x-section for distribution: $var" nohup time python src/cross_section_measurement/04_make_plots_matplotlib.py -v $var -c 7 -p data/normalisation/background_subtraction/ -o plots/background_subtraction &> logs/04b_${var}_plot_7TeV.log & # -a <--add this option for additional plots let i+=1 diff --git a/bin/x_05_all_vars b/bin/x_05_all_vars index 4ebe498a..6560c279 100755 --- a/bin/x_05_all_vars +++ b/bin/x_05_all_vars @@ -19,13 +19,13 @@ for var in MET HT ST MT WPT; do wait; fi - nohup time python src/cross_section_measurement/05_make_tables.py -v $var -c 8 -p data/$nice_fit_var -a -o tables/fit &> logs/05_${var}_table_8TeV_${nice_fit_var}.log & - let i+=1 - if (( $i % N_JOBS == 0 )) - then - echo "Waiting on the above to finish." - wait; - fi + # nohup time python src/cross_section_measurement/05_make_tables.py -v $var -c 8 -p data/$nice_fit_var -a -o tables/fit &> logs/05_${var}_table_8TeV_${nice_fit_var}.log & + # let i+=1 + # if (( $i % N_JOBS == 0 )) + # then + # echo "Waiting on the above to finish." + # wait; + # fi done wait; diff --git a/bin/x_05b_all_vars b/bin/x_05b_all_vars index 4c741393..30f5637c 100755 --- a/bin/x_05b_all_vars +++ b/bin/x_05b_all_vars @@ -2,26 +2,24 @@ echo "This will take a while ... grab a coffee/tea/water" mkdir -p logs mkdir -p plots -N_JOBS=5 +N_JOBS=1 i=0 -for var in MET HT ST WPT MT; do +for var in MET HT ST WPT; do echo "Tabulating diff. x-section for distribution: $var" - nohup time python src/cross_section_measurement/05_make_tables.py -v $var -c 7 -p data/normalisation/background_subtraction/ -a -o tables/background_subtraction &> logs/05b_${var}_table_7TeV.log & - let i+=1 - if (( $i % N_JOBS == 0 )) - then - echo "Waiting on the above to finish." - wait; - fi + # nohup time + echo $var + echo "7" + python src/cross_section_measurement/05_make_tables.py -v $var -c 7 -p data/normalisation/background_subtraction/ -a -o tables/background_subtraction + # &> logs/05b_${var}_table_7TeV.log & + # wait; - nohup time python src/cross_section_measurement/05_make_tables.py -v $var -c 8 -p data/normalisation/background_subtraction/ -a -o tables/background_subtraction &> logs/05b_${var}_table_8TeV.log & - let i+=1 - if (( $i % N_JOBS == 0 )) - then - echo "Waiting on the above to finish." - wait; - fi + # nohup time + echo $var + echo "8" + python src/cross_section_measurement/05_make_tables.py -v $var -c 8 -p data/normalisation/background_subtraction/ -a -o tables/background_subtraction + # &> logs/05b_${var}_table_8TeV.log & + wait; done wait; diff --git a/config/CMS.py b/config/CMS.py index 5b97410e..de83fd3e 100644 --- a/config/CMS.py +++ b/config/CMS.py @@ -15,7 +15,7 @@ } y_axis_title = { - 'fontsize':50, + 'fontsize':40, 'position' : (0., 1.), 'verticalalignment': 'bottom', 'horizontalalignment': 'right' @@ -43,7 +43,7 @@ 'pad': 12 } -legend_properties = {'size':35} +legend_properties = {'size':33} figsize = (16,16) dpi = 400 diff --git a/config/__init__.py b/config/__init__.py index 47f97dd4..72f25ece 100644 --- a/config/__init__.py +++ b/config/__init__.py @@ -42,8 +42,9 @@ class XSectionConfig(): 'unfolding_matching_up', 'unfolding_matching_up_raw', 'unfolding_mass_down', 'unfolding_mass_up', 'unfolding_mcatnlo', 'unfolding_mcatnlo_raw', - 'unfolding_powheg_pythia', 'unfolding_powheg_pythia_raw', - 'unfolding_powheg_herwig', 'unfolding_powheg_herwig_raw', + 'unfolding_powheg_v2_pythia', 'unfolding_powheg_v2_pythia_raw', + 'unfolding_powheg_v1_herwig', 'unfolding_powheg_v1_herwig_raw', + 'unfolding_powheg_v2_herwig', 'unfolding_powheg_v2_herwig_raw', 'unfolding_scale_down', 'unfolding_scale_down_raw', 'unfolding_scale_up', 'unfolding_scale_up_raw', 'unfolding_ptreweight', 'unfolding_ptreweight_raw', @@ -193,7 +194,12 @@ def __fill_defaults__( self ): } self.generator_systematics = [ 'matchingup', 'matchingdown', 'scaleup', 'scaledown' ] self.topMass_systematics = [ 'TTJets_massup', 'TTJets_massdown'] - self.topMasses = [169.5, 172.5, 173.5] + + if self.centre_of_mass_energy == 7: + self.topMasses = [166.5, 172.5, 178.5] + elif self.centre_of_mass_energy == 8: + self.topMasses = [169.5, 172.5, 173.5] + self.topMassUncertainty = 1.0 # GeV from https://twiki.cern.ch/twiki/bin/view/LHCPhysics/TtbarNNLO self.central_general_template = path_to_files + 'central/%s' + middle + '.root' self.generator_systematic_vjets_templates = { systematic: path_to_files + 'central/VJets-%s_%dpb_PFElectron_PFMuon_PF2PATJets_PFMET.root' % ( systematic, self.luminosity ) for systematic in self.generator_systematics} @@ -217,8 +223,10 @@ def __fill_defaults__( self ): } self.unfolding_madgraph_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV.root' % self.centre_of_mass_energy - self.unfolding_powheg_pythia_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_powheg.root' % self.centre_of_mass_energy - self.unfolding_powheg_herwig_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_powhegherwig.root' % self.centre_of_mass_energy + self.unfolding_powheg_v1_pythia_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_powhegV1pythia.root' % self.centre_of_mass_energy + self.unfolding_powheg_v2_pythia_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_powhegV2pythia.root' % self.centre_of_mass_energy + self.unfolding_powheg_v1_herwig_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_powhegV1herwig.root' % self.centre_of_mass_energy + self.unfolding_powheg_v2_herwig_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_powhegV2herwig.root' % self.centre_of_mass_energy self.unfolding_mcatnlo_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_mcatnlo.root' % self.centre_of_mass_energy self.unfolding_ptreweight_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_withTopPtReweighting.root' % self.centre_of_mass_energy @@ -228,8 +236,10 @@ def __fill_defaults__( self ): self.unfolding_matching_up_raw = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_matchingup.root' % self.centre_of_mass_energy self.unfolding_madgraph = self.unfolding_madgraph_raw.replace( '.root', '_asymmetric.root' ) - self.unfolding_powheg_pythia = self.unfolding_powheg_pythia_raw.replace( '.root', '_asymmetric.root' ) - self.unfolding_powheg_herwig = self.unfolding_powheg_herwig_raw.replace( '.root', '_asymmetric.root' ) + self.unfolding_powheg_v1_pythia = self.unfolding_powheg_v1_pythia_raw.replace( '.root', '_asymmetric.root' ) + self.unfolding_powheg_v2_pythia = self.unfolding_powheg_v2_pythia_raw.replace( '.root', '_asymmetric.root' ) + self.unfolding_powheg_v1_herwig = self.unfolding_powheg_v1_herwig_raw.replace( '.root', '_asymmetric.root' ) + self.unfolding_powheg_v2_herwig = self.unfolding_powheg_v2_herwig_raw.replace( '.root', '_asymmetric.root' ) self.unfolding_mcatnlo = self.unfolding_mcatnlo_raw.replace( '.root', '_asymmetric.root' ) self.unfolding_ptreweight = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_asymmetric_withTopPtReweighting.root' % self.centre_of_mass_energy @@ -240,6 +250,11 @@ def __fill_defaults__( self ): self.unfolding_mass_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_massdown_asymmetric.root' % self.centre_of_mass_energy self.unfolding_mass_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_massup_asymmetric.root' % self.centre_of_mass_energy + if self.centre_of_mass_energy == 7: + self.unfolding_mass_down = path_to_unfolding_histograms + 'unfolding_TTJets_8TeV_massdown_8To7_asymmetric.root' + self.unfolding_mass_up = path_to_unfolding_histograms + 'unfolding_TTJets_8TeV_massup_8To7_asymmetric.root' + + self.unfolding_pdfweights = {index : path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_asymmetric_pdfWeight_%d.root' % (self.centre_of_mass_energy, index) for index in range( 1, 46 )} self.histogram_path_templates = { @@ -269,6 +284,23 @@ def __fill_defaults__( self ): self.luminosity_scale = self.new_luminosity / self.luminosity + self.typical_systematics_order = [ + 'typical_systematics_theoretical', + 'typical_systematics_hadronisation', + 'typical_systematics_PDF', + 'typical_systematics_top_mass', + 'typical_systematics_pt_reweight', + 'typical_systematics_electron', + 'typical_systematics_muon', + 'typical_systematics_btagging', + 'typical_systematics_JES', + 'typical_systematics_JER', + 'typical_systematics_MET', + 'typical_systematics_PU', + 'typical_systematics_background_other', + 'typical_systematics_QCD_shape', + ] + self.typical_systematics = { "typical_systematics_electron": ['Electron_down', 'Electron_up'], @@ -419,11 +451,15 @@ def __fill_defaults_8TeV__( self ): 'JES_down': path_to_files + 'JES_down/SingleElectron' + middle + self.categories_and_prefixes['JES_down'] + '.root' } -fit_var_inputs = ['absolute_eta', 'M3', 'M_bl', 'angle_bl', - 'absolute_eta_angle_bl', - 'absolute_eta_M3', - 'absolute_eta_M_bl', - 'absolute_eta_M_bl_angle_bl', +fit_var_inputs = [ + 'absolute_eta', 'M3', +# 'M_bl', + 'angle_bl', + # 'absolute_eta_angle_bl', + # 'absolute_eta_M3', + # 'absolute_eta_M_bl', + # 'absolute_eta_M_bl_angle_bl', 'absolute_eta_M3_angle_bl', - 'absolute_eta_M_bl_M3', - 'absolute_eta_M_bl_M3_angle_bl' ] + # 'absolute_eta_M_bl_M3', + # 'absolute_eta_M_bl_M3_angle_bl' + ] diff --git a/config/latex_labels.py b/config/latex_labels.py index 088bc68c..9ceae5ba 100644 --- a/config/latex_labels.py +++ b/config/latex_labels.py @@ -3,48 +3,69 @@ @author: kreczko ''' -b_tag_bins_latex = {'0btag':'0 b-tags', '0orMoreBtag':'$\geq$ 0 b-tags', '1btag':'1 b-tag', - '1orMoreBtag':'$\geq$ 1 b-tags', - '2btags':'2 b-tags', '2orMoreBtags':'$\geq$ 2 b-tags', - '3btags':'3 b-tags', '3orMoreBtags':'$\geq$ 3 b-tags', - '4orMoreBtags':'$\geq$ 4 b-tags'} +b_tag_bins_latex = {'0btag':'0 b tags', '0orMoreBtag':'$\geq$ 0 b tags', '1btag':'1 b tag', + '1orMoreBtag':'$\geq$ 1 b tags', + '2btags':'2 b tags', '2orMoreBtags':'$\geq$ 2 b tags', + '3btags':'3 b tags', '3orMoreBtags':'$\geq$ 3 b tags', + '4orMoreBtags':'$\geq$ 4 b tags'} variables_latex = { - 'MET': 'E_{\mathrm{T}}^{\mathrm{miss}}', - 'HT': 'H_{\mathrm{T}}', - 'ST': 'S_{\mathrm{T}}', - 'MT': 'M^{\mathrm{W}}_{\mathrm{T}}', - 'WPT': 'p^\mathrm{W}_{\mathrm{T}}'} + 'MET': '\ensuremath{E_{\mathrm{T}}^{\mathrm{miss}}}', + 'HT': '\ensuremath{H_{\mathrm{T}}}', + 'ST': '\ensuremath{S_{\mathrm{T}}}', + 'MT': '\ensuremath{M^{\mathrm{W}}_{\mathrm{T}}}', + 'WPT': '\ensuremath{p^\mathrm{W}_{\mathrm{T}}}' + } + +variables_latex_macros = { + 'MET': '\met', + 'HT': '\HT', + 'ST': '\st', + 'WPT': '\wpt' + } + +variables_hepdata = { + 'MET': 'E_{T}^{miss}', + 'HT': 'H_{T}', + 'ST': 'S_{T}', + 'WPT': 'p^{W}_{T}' + } + +ttBarLatex = '$\mathrm{t}\\bar{\mathrm{t}}$' measurements_latex = {'unfolded': 'unfolded', 'measured': 'measured', - 'MADGRAPH': '$t\\bar{t}$ (MADGRAPH+Pythia)', - 'MADGRAPH_ptreweight': '$t\\bar{t}$ (MADGRAPH+$p_\mathrm{T}^\mathrm{reweight}$)', - 'MCATNLO': '$t\\bar{t}$ (MC@NLO+Herwig)', - 'POWHEG_PYTHIA': '$t\\bar{t}$ (POWHEG+Pythia)', - 'POWHEG_HERWIG': '$t\\bar{t}$ (POWHEG+Herwig)', - 'matchingdown': '$t\\bar{t}$ (matching down)', - 'matchingup': '$t\\bar{t}$ (matching up)', - 'scaledown': '$t\\bar{t}$ ($Q^{2}$ down)', - 'scaleup': '$t\\bar{t}$ ($Q^{2}$ up)', - 'TTJets_matchingdown': '$t\\bar{t}$ (matching down)', - 'TTJets_matchingup': '$t\\bar{t}$ (matching up)', - 'TTJets_scaledown': '$t\\bar{t}$ ($Q^{2}$ down)', - 'TTJets_scaleup': '$t\\bar{t}$ ($Q^{2}$ up)', - 'TTJets_massdown': '$t\\bar{t}$ (top mass down)', - 'TTJets_massup': '$t\\bar{t}$ (top mass up)', + 'MADGRAPH': r'\textsc{MadGraph \raisebox{.2ex}{+} pythia}', + 'MADGRAPH_ptreweight': r'\textsc{MadGraph} \raisebox{.2ex}{+} $p_\mathrm{T}$ reweighting', + 'MCATNLO': r'\textsc{mc@nlo \raisebox{.2ex}{+} herwig}', + 'powheg_v1_pythia': r'\textsc{powheg v1 \raisebox{.2ex}{+} pythia}', + 'powheg_v2_pythia': r'\textsc{powheg v2 \raisebox{.2ex}{+} pythia}', + 'powheg_v1_herwig': r'\textsc{powheg v1 \raisebox{.2ex}{+} herwig}', + 'powheg_v2_herwig': r'\textsc{powheg v2 \raisebox{.2ex}{+} herwig}', + 'matchingdown': 'Matching down', + 'matchingup': 'Matching up', + # 'scaledown': '$Q^{2}$ down', + # 'scaleup': '$Q^{2}$ up', + 'scaledown': '$\mu_{R}, \mu_{F}$ down', + 'scaleup': '$\mu_{R}, \mu_{F}$ up', + 'TTJets_matchingdown': 'Matching down', + 'TTJets_matchingup': 'Matching up', + 'TTJets_scaledown': '$Q^{2}$ down', + 'TTJets_scaleup': '$Q^{2}$ up', + 'TTJets_massdown': 'Top mass down', + 'TTJets_massup': 'Top mass up', 'VJets_matchingdown': 'V+jets (matching down)', 'VJets_matchingup': 'V+jets (matching up)', 'VJets_scaledown': 'V+jets ($Q^{2}$ down)', 'VJets_scaleup': 'V+jets ($Q^{2}$ up)', - 'BJet_down':'b-tagging efficiency $-1\sigma$', - 'BJet_up':'b-tagging efficiency $+1\sigma$', + 'BJet_down':'b tagging efficiency $-1\sigma$', + 'BJet_up':'b tagging efficiency $+1\sigma$', 'JES_down':'Jet energy scale $-1\sigma$', 'JES_up':'Jet energy scale $+1\sigma$', 'JER_down':'Jet energy resolution $-1\sigma$', 'JER_up':'Jet energy resolution $+1\sigma$', - 'LightJet_down':'b-tagging mis-tag rate $-1\sigma$', - 'LightJet_up':'b-tagging mis-tag rate $+1\sigma$', + 'LightJet_down':'b tagging mis-tag rate $-1\sigma$', + 'LightJet_up':'b tagging mis-tag rate $+1\sigma$', 'PU_down':'Pile-up $-1\sigma$', 'PU_up':'Pile-up $+1\sigma$', 'central':'central', @@ -55,8 +76,8 @@ 'QCD_shape' : 'QCD shape uncertainty', 'luminosity+' : 'Luminosity $+1\sigma$', 'luminosity-' : 'Luminosity $-1\sigma$', - 'TTJet_cross_section+' : '$t\\bar{t}$ cross section $+1\sigma$', - 'TTJet_cross_section-' : '$t\\bar{t}$ cross section $-1\sigma$', + 'TTJet_cross_section+' : ttBarLatex + ' cross section $+1\sigma$', + 'TTJet_cross_section-' : ttBarLatex + ' cross section $-1\sigma$', 'SingleTop_cross_section+' : 'Single top cross section $+1\sigma$', 'SingleTop_cross_section-' : 'Single top cross section $-1\sigma$', 'V+Jets_cross_section+': 'V+jets cross section \ensuremath{+1\sigma}', @@ -87,40 +108,41 @@ } samples_latex = { - 'data':'data', + 'Data':'Data', 'QCD':'QCD', 'WJets':'W $\\rightarrow \ell\\nu$', 'ZJets':'Z/$\gamma^*$ + jets', - 'TTJet':'$\mathrm{t}\\bar{\mathrm{t}}$', - 'SingleTop':'Single-Top' , - 'V+Jets' : 'W/Z + jets' + 'TTJet':ttBarLatex, + 'SingleTop':'Single Top' , + 'V+Jets' : r'W/Z \raisebox{.2ex}{+} jets' } fit_variables_latex = { 'absolute_eta' : r'lepton $|\eta|$', - 'M3' : r'$M3$', + 'M3' : '\ensuremath{M_3}', 'M_bl' : r'$M(b,l)$', 'angle_bl' : r'$\alpha$', } typical_systematics_latex = {"typical_systematics_electron": "Electron trigger efficiency \& electron selection", "typical_systematics_muon": "Muon trigger efficiency \& muon selection", - "typical_systematics_btagging": "btagging", - "typical_systematics_JES": "Jet Energy Scale", - "typical_systematics_JER": "Jet Energy Resolution", - "typical_systematics_PU": "pileup", - "typical_systematics_hadronisation": "hadronisation", + "typical_systematics_btagging": "b tagging", + "typical_systematics_JES": "Jet energy scale", + "typical_systematics_JER": "Jet energy resolution", + "typical_systematics_PU": "Pileup", + "typical_systematics_hadronisation": "Hadronization", "typical_systematics_QCD_shape": "QCD shape", - "typical_systematics_PDF": "PDF uncertainties", - "typical_systematics_top_mass": "top mass", - "typical_systematics_theoretical": "Theoretical systematics", - 'typical_systematics_background_other': 'Background (other)', - 'typical_systematics_MET': '$E_{T}^{miss}$ uncertainties', - 'typical_systematics_pt_reweight': '$p_\mathrm{T}$ reweighting' + "typical_systematics_PDF": "PDF", + "typical_systematics_top_mass": "Top quark mass", + "typical_systematics_theoretical": "Fact./Renorm. scales and matching threshold", + 'typical_systematics_background_other': 'Background Normalization', + 'typical_systematics_MET': '\met', + 'typical_systematics_pt_reweight': 'Top quark $p_\mathrm{T}$ reweighting' } channel_latex = { 'electron' : r"e + jets", - 'muon' : r"$\mu$ + jets", - 'combined' : r"e, $\mu$ + jets combined", + 'muon' : r"\boldmath$\mu$ + jets", + 'combined' : r"e, \boldmath$\mu$ + jets combined", + } \ No newline at end of file diff --git a/config/plots/HT_7TeV_electron_channel_QCD_shape_comparison.json b/config/plots/HT_7TeV_electron_channel_QCD_shape_comparison.json new file mode 100644 index 00000000..5451587b --- /dev/null +++ b/config/plots/HT_7TeV_electron_channel_QCD_shape_comparison.json @@ -0,0 +1,51 @@ +{ + "command": "compare-hists", + "files": [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_8th_draft/7TeV/central/ElectronHad_5050pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ], + "histograms": [ + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/HT_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCD non iso e+jets/MET/HT_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/HT_2btags" + ], + "labels": [ + "Conversion 0 btag", + "Non Iso 0 btag", + "Conversion 2 btags" + ], + "output_folder": "plots/QCD_shape_e", + "output_format": ["pdf"], + "name_prefix": "compare_QCD_btag_regions_in_HT_electron_channel_", + "plot_type": "shape_comparison", + "ratio_y_limits": [ + [0.95, 1.05] + ], + "title": [ + "Comparison of $H_T$ \\\\(QCD control regions: 0,1 and 2 btags) $\\sqrt{s}$ = 8 TeV" + ], + "x_axis_title": [ + "$H_T$ [GeV]" + ], + "x_limits": [ + [120, 1000] + ], + "y_axis_title": [ + "normalised to unit area" + ], + "y_limits": [ + [0, 0.25] + ], + "rebin" : [ + [120.0, 185.0, 215.0, 247.0, 283.0, 323.0, 365.0, 409.0, 458.0, 512.0, 570.0, 629.0, 691.0, 769.0, 1000.0] + ], + "colours": ["green", "blue", "red"], + "normalise":1, + "fill_area":false, + + "subtractFiles" : [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/central/VJets_5050pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/central/SingleTop_5050pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/central/TTJet_5050pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ] + +} \ No newline at end of file diff --git a/config/plots/HT_8TeV_electron_channel.json b/config/plots/HT_8TeV_electron_channel.json new file mode 100644 index 00000000..5fc29c30 --- /dev/null +++ b/config/plots/HT_8TeV_electron_channel.json @@ -0,0 +1,49 @@ +{ + "command": "compare-files", + "files": [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/SingleElectron_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/QCD_Electron_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/VJets_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/SingleTop_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ], + "histograms": [ + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/HT_2orMoreBtags" + ], + "labels": [ + "Data", + "QCD", + "W/Z + jets", + "single-top", + "$t\\bar{t}$" + ], + "output_folder": "plots/8TeV", + "output_format": ["png"], + "plot_type": "data_mc_comparison", + "ratio_y_limits": [ + [ + 0.8, + 1.3 + ] + ], + "title": [ + "$H_T$ at $\\sqrt{s}$ = 8 TeV, e+jets" + ], + "x_axis_title": [ + "$H_T [GeV]$" + ], + "x_limits": [ + [ + 120, + 1000 + ] + ], + "rebin" : [ + [120.0, 185.0, 215.0, 247.0, 283.0, 323.0, 365.0, 409.0, 458.0, 512.0, 570.0, 629.0, 691.0, 769.0, 1000.0] + ], + "y_axis_title": [ + "events" + ], + "colours": ["black", "yellow", "green", "magenta", "red"], + "mc_error": [0.15] +} \ No newline at end of file diff --git a/config/plots/HT_8TeV_electron_channel_QCD_shape_comparison.json b/config/plots/HT_8TeV_electron_channel_QCD_shape_comparison.json new file mode 100644 index 00000000..673c8236 --- /dev/null +++ b/config/plots/HT_8TeV_electron_channel_QCD_shape_comparison.json @@ -0,0 +1,51 @@ +{ + "command": "compare-hists", + "files": [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_8th_draft/8TeV/central/SingleElectron_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ], + "histograms": [ + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/HT_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCD non iso e+jets/MET/HT_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/HT_2btags" + ], + "labels": [ + "Conversion 0 btag", + "Non Iso 0 btag", + "Conversion 2 btags" + ], + "output_folder": "plots/QCD_shape_e", + "output_format": ["pdf"], + "name_prefix": "compare_QCD_btag_regions_in_HT_electron_channel_", + "plot_type": "shape_comparison", + "ratio_y_limits": [ + [0.95, 1.05] + ], + "title": [ + "Comparison of $H_T$ \\\\(QCD control regions: 0,1 and 2 btags) $\\sqrt{s}$ = 8 TeV" + ], + "x_axis_title": [ + "$H_T$ [GeV]" + ], + "x_limits": [ + [120, 1000] + ], + "y_axis_title": [ + "normalised to unit area" + ], + "y_limits": [ + [0, 0.25] + ], + "rebin" : [ + [120.0, 185.0, 215.0, 247.0, 283.0, 323.0, 365.0, 409.0, 458.0, 512.0, 570.0, 629.0, 691.0, 769.0, 1000.0] + ], + "colours": ["green", "blue", "red"], + "normalise":1, + "fill_area":false, + + "subtractFiles" : [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/VJets_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/SingleTop_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ] + +} \ No newline at end of file diff --git a/config/plots/MET_8TeV_electron_channel.json b/config/plots/MET_8TeV_electron_channel.json index 2d0b1cfe..d9d539dc 100644 --- a/config/plots/MET_8TeV_electron_channel.json +++ b/config/plots/MET_8TeV_electron_channel.json @@ -8,10 +8,10 @@ "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" ], "histograms": [ - "TTbar_plus_X_analysis/EPlusJets/Ref selection/MET/patType1CorrectedPFMet/MET_2orMoreBtags" + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/MET_2orMoreBtags" ], "labels": [ - "data", + "Data", "QCD", "W/Z + jets", "single-top", @@ -38,6 +38,9 @@ 300 ] ], + "rebin" : [ + [0.0, 27.0, 52.0, 87.0, 130.0, 172.0, 300] + ], "y_axis_title": [ "events/5 GeV" ], diff --git a/config/plots/MET_8TeV_electron_channel_QCD_shape_comparison.json b/config/plots/MET_8TeV_electron_channel_QCD_shape_comparison.json index 047202df..db39fe9f 100644 --- a/config/plots/MET_8TeV_electron_channel_QCD_shape_comparison.json +++ b/config/plots/MET_8TeV_electron_channel_QCD_shape_comparison.json @@ -4,14 +4,14 @@ "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_8th_draft/8TeV/central/SingleElectron_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" ], "histograms": [ + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/MET_0btag", "TTbar_plus_X_analysis/EPlusJets/QCD non iso e+jets/MET/patType1CorrectedPFMet/MET_0btag", - "TTbar_plus_X_analysis/EPlusJets/QCD non iso e+jets/MET/patType1CorrectedPFMet/MET_1btag", - "TTbar_plus_X_analysis/EPlusJets/QCD non iso e+jets/MET/patType1CorrectedPFMet/MET_2btags" + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/MET_2btags" ], "labels": [ - "0 btag", - "1 btag", - "2 btags" + "Conversion 0 btag", + "Non Iso 0 btag", + "Conversion 2 btags" ], "output_folder": "plots/QCD_shape_e", "output_format": ["pdf"], @@ -33,11 +33,19 @@ "normalised to unit area" ], "y_limits": [ - [0, 0.6] + [0, 1.0] ], "rebin" : [ [0.0, 27.0, 52.0, 87.0, 130.0, 172.0, 300] ], - "colours": ["green", "yellow", "red"], - "normalise":1 + "colours": ["green", "blue", "red"], + "normalise":1, + "fill_area":false, + + "subtractFiles" : [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/VJets_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/SingleTop_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ] + } \ No newline at end of file diff --git a/config/plots/ST_7TeV_electron_channel_QCD_shape_comparison.json b/config/plots/ST_7TeV_electron_channel_QCD_shape_comparison.json new file mode 100644 index 00000000..4a42e764 --- /dev/null +++ b/config/plots/ST_7TeV_electron_channel_QCD_shape_comparison.json @@ -0,0 +1,51 @@ +{ + "command": "compare-hists", + "files": [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_8th_draft/7TeV/central/ElectronHad_5050pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ], + "histograms": [ + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/ST_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCD non iso e+jets/MET/patType1CorrectedPFMet/ST_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/ST_2btags" + ], + "labels": [ + "Conversion 0 btag", + "Non Iso 0 btag", + "Conversion 2 btags" + ], + "output_folder": "plots/QCD_shape_e", + "output_format": ["pdf"], + "name_prefix": "compare_QCD_btag_regions_in_ST_electron_channel_", + "plot_type": "shape_comparison", + "ratio_y_limits": [ + [0.95, 1.05] + ], + "title": [ + "Comparison of $S_T$ \\\\(QCD control regions: 0,1 and 2 btags) $\\sqrt{s}$ = 8 TeV" + ], + "x_axis_title": [ + "$S_T$ [GeV]" + ], + "x_limits": [ + [146, 1200] + ], + "y_axis_title": [ + "normalised to unit area" + ], + "y_limits": [ + [0, 0.5] + ], + "rebin" : [ + [146.0, 277.0, 319.0, 361.0, 408.0, 459.0, 514.0, 573.0, 637.0, 705.0, 774.0, 854.0, 940.0, 1200.0] + ], + "colours": ["green", "blue", "red"], + "normalise":1, + "fill_area":false, + + "subtractFiles" : [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/central/VJets_5050pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/central/SingleTop_5050pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/7TeV/central/TTJet_5050pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ] + +} \ No newline at end of file diff --git a/config/plots/ST_8TeV_electron_channel.json b/config/plots/ST_8TeV_electron_channel.json new file mode 100644 index 00000000..f4e50a99 --- /dev/null +++ b/config/plots/ST_8TeV_electron_channel.json @@ -0,0 +1,49 @@ +{ + "command": "compare-files", + "files": [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/SingleElectron_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/QCD_Electron_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/VJets_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/SingleTop_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ], + "histograms": [ + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/ST_2orMoreBtags" + ], + "labels": [ + "Data", + "QCD", + "W/Z + jets", + "single-top", + "$t\\bar{t}$" + ], + "output_folder": "plots/8TeV", + "output_format": ["png"], + "plot_type": "data_mc_comparison", + "ratio_y_limits": [ + [ + 0.8, + 1.3 + ] + ], + "title": [ + "$S_T$ at $\\sqrt{s}$ = 8 TeV, e+jets" + ], + "x_axis_title": [ + "$S_T [GeV]$" + ], + "x_limits": [ + [ + 146, + 2000 + ] + ], + "rebin" : [ + [146.0, 277.0, 319.0, 361.0, 408.0, 459.0, 514.0, 573.0, 637.0, 705.0, 774.0, 854.0, 940.0, 1200.0] + ], + "y_axis_title": [ + "events" + ], + "colours": ["black", "yellow", "green", "magenta", "red"], + "mc_error": [0.15] +} \ No newline at end of file diff --git a/config/plots/ST_8TeV_electron_channel_QCD_shape_comparison.json b/config/plots/ST_8TeV_electron_channel_QCD_shape_comparison.json new file mode 100644 index 00000000..f5176925 --- /dev/null +++ b/config/plots/ST_8TeV_electron_channel_QCD_shape_comparison.json @@ -0,0 +1,51 @@ +{ + "command": "compare-hists", + "files": [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_8th_draft/8TeV/central/SingleElectron_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ], + "histograms": [ + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/ST_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCD non iso e+jets/MET/patType1CorrectedPFMet/ST_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/ST_2btags" + ], + "labels": [ + "Conversion 0 btag", + "Non Iso 0 btag", + "Conversion 2 btags" + ], + "output_folder": "plots/QCD_shape_e", + "output_format": ["pdf"], + "name_prefix": "compare_QCD_btag_regions_in_ST_electron_channel_", + "plot_type": "shape_comparison", + "ratio_y_limits": [ + [0.95, 1.05] + ], + "title": [ + "Comparison of $S_T$ \\\\(QCD control regions: 0,1 and 2 btags) $\\sqrt{s}$ = 8 TeV" + ], + "x_axis_title": [ + "$S_T$ [GeV]" + ], + "x_limits": [ + [146, 1200] + ], + "y_axis_title": [ + "normalised to unit area" + ], + "y_limits": [ + [0, 0.25] + ], + "rebin" : [ + [146.0, 277.0, 319.0, 361.0, 408.0, 459.0, 514.0, 573.0, 637.0, 705.0, 774.0, 854.0, 940.0, 1200.0] + ], + "colours": ["green", "blue", "red"], + "normalise":1, + "fill_area":false, + + "subtractFiles" : [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/VJets_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/SingleTop_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ] + +} \ No newline at end of file diff --git a/config/plots/WPT_8TeV_electron_channel_QCD_shape_comparison.json b/config/plots/WPT_8TeV_electron_channel_QCD_shape_comparison.json new file mode 100644 index 00000000..c90bd423 --- /dev/null +++ b/config/plots/WPT_8TeV_electron_channel_QCD_shape_comparison.json @@ -0,0 +1,51 @@ +{ + "command": "compare-hists", + "files": [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_8th_draft/8TeV/central/SingleElectron_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ], + "histograms": [ + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/WPT_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCD non iso e+jets/MET/patType1CorrectedPFMet/WPT_0btag", + "TTbar_plus_X_analysis/EPlusJets/QCDConversions/MET/patType1CorrectedPFMet/WPT_2btags" + ], + "labels": [ + "Conversion 0 btag", + "Non Iso 0 btag", + "Conversion 2 btags" + ], + "output_folder": "plots/QCD_shape_e", + "output_format": ["pdf"], + "name_prefix": "compare_QCD_btag_regions_in_WPT_electron_channel_", + "plot_type": "shape_comparison", + "ratio_y_limits": [ + [0.95, 1.05] + ], + "title": [ + "" + ], + "x_axis_title": [ + "$p_T^W$ [GeV]" + ], + "x_limits": [ + [0, 300] + ], + "y_axis_title": [ + "Arbitrary Units" + ], + "y_limits": [ + [0, 0.45] + ], + "rebin" : [ + [0.0, 27.0, 52.0, 78.0, 105.0, 134.0, 166.0, 200.0, 237.0, 300.0] + ], + "colours": ["green", "blue", "red"], + "normalise":1, + "fill_area":false, + + "subtractFiles" : [ + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/VJets_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/SingleTop_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root", + "/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_6th_draft/8TeV/central/TTJet_19584pb_PFElectron_PFMuon_PF2PATJets_PFMET.root" + ] + +} \ No newline at end of file diff --git a/config/variable_binning.py b/config/variable_binning.py index 0b618401..e21f5cb6 100644 --- a/config/variable_binning.py +++ b/config/variable_binning.py @@ -2,7 +2,7 @@ 'absolute_eta' : [round( i * 0.2, 2 ) for i in range ( int( 3 / 0.2 ) + 1 )], 'M3' : [i * 25 for i in range ( int( 1000 / 25 ) + 1 )], 'M_bl' : [i * 10 for i in range ( int( 1000 / 20 ) + 1 )], - 'angle_bl' : [round( i * 0.2, 2 ) for i in range ( int( 4 / 0.2 ) + 1 )], + 'angle_bl' : [round( i * 0.2, 2 ) for i in range ( int( 3.2 / 0.2 ) + 1 )], } bin_edges = { 'MET':[0.0, 27.0, 52.0, 87.0, 130.0, 172.0, 300.0], @@ -31,9 +31,9 @@ upper_edge = bin_edges[variable][i + 1] bin_widths[variable].append( upper_edge - lower_edge ) bin_name = '%d-%d' % ( int( lower_edge ), int( upper_edge ) ) - bin_name_latex = '%d--%d~\GeV' % ( int( lower_edge ), int( upper_edge ) ) - if ( i + 1 ) == number_of_edges - 1: - bin_name = '%d-inf' % int( lower_edge ) - bin_name_latex = '$\\geq %d$~\GeV' % int( lower_edge ) + bin_name_latex = '%d--%d' % ( int( lower_edge ), int( upper_edge ) ) + # if ( i + 1 ) == number_of_edges - 1: + # bin_name = '%d-inf' % int( lower_edge ) + # bin_name_latex = '$\\geq %d$' % int( lower_edge ) variable_bins_ROOT[variable].append( bin_name ) variable_bins_latex[bin_name] = bin_name_latex diff --git a/cp b/cp new file mode 100644 index 00000000..e69de29b diff --git a/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/produceUnfoldingHistograms.py b/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/produceUnfoldingHistograms.py index 0ea9c605..74643d7b 100644 --- a/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/produceUnfoldingHistograms.py +++ b/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/produceUnfoldingHistograms.py @@ -8,6 +8,8 @@ from scaleFactors import * +import sys + class channel: def __init__(self, channelName, treeName, outputDirName): self.channelName = channelName @@ -76,16 +78,28 @@ def copyEventFilterHist( inputFile, outputFile ): fileNames = { '8TeV' : { - 'central' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_central_8TeV.root', + 'central':'/storage/ec6821/NTupleProd/CMSSW_5_3_23/src/TTJets_8TeV_central.root', + # 'central' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_central_8TeV.root', 'scaleup' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_scaleup_8TeV.root', 'scaledown' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_scaledown_8TeV.root', 'matchingup' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_matchingup_8TeV.root', 'matchingdown' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_matchingdown_8TeV.root', 'powheg' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_powhegpythia_8TeV.root', + 'powhegV2pythia' : '/storage/ec6821/NTupleProd/CMSSW_5_3_23/src/TTJets_PowhegPythia_new__8TeV.root', 'powhegherwig' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_powhegherwig_8TeV.root', + # 'powhegherwig_new' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_powhegherwig_NEW_8TeV.root', + 'powhegherwig_new' : '/storage/ec6821/NTupleProd/CMSSW_5_3_23/src/TTJets_nTuple_53X_mc.root', 'mcatnlo' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mcatnlo_8TeV.root', - 'massdown' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_169_5_8TeV.root', - 'massup' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_173_5_8TeV.root', + # 'massdown' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_169_5_8TeV.root', + # 'massup' : '/hdfs/TopQuarkGroup/mc/8TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_173_5_8TeV.root', + + # 'massdown' : '/storage/ec6821/DailyPythonScripts/legacy/CMSSW_7_4_7/src/DailyPythonScripts/TTJets_mass_169_5_0.root', + # 'massup' : '/storage/ec6821/DailyPythonScripts/legacy/CMSSW_7_4_7/src/DailyPythonScripts/unfolding_TTJets_mass_173_5_8TeV.root', + + 'massdown' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_169_5_8TeV.root', + 'massup' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_173_5_8TeV.root', + + }, '7TeV' : { 'central' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_central_7TeV.root', @@ -93,8 +107,10 @@ def copyEventFilterHist( inputFile, outputFile ): 'scaleup' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_scaleup_7TeV.root', 'matchingdown' :'/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_matchingdown_7TeV.root', 'matchingup' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_matchingup_7TeV.root', - 'massdown' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_169_5_7TeV.root', - 'massup' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_173_5_7TeV.root', + # 'massdown' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_169_5_7TeV.root', + # 'massup' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_173_5_7TeV.root', + 'massdown' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_166_5_7TeV.root', + 'massup' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_mass_178_5_7TeV.root', 'powheg' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_powhegpythia_7TeV.root', 'powhegherwig' : '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/unfolding_TTJets_powhegherwig_7TeV.root', } @@ -116,9 +132,14 @@ def main(): parser.add_option('-n', action='store_true', dest='donothing', default=False) parser.add_option('-e', action='store_true', dest='extraHists', default=False) parser.add_option('-f',action='store_true', dest='fineBinned', default=False) + parser.add_option('--eightToSeven', action='store_true', dest='eightToSeven', default=False) (options, _) = parser.parse_args() + if options.eightToSeven and not int(options.centreOfMassEnergy) == 8: + print "Error : Reweighting from 8 TeV to 7 TeV, but input sample is not 8 TeV" + sys.exit() + # Input file name file_name = 'crap.root' if int(options.centreOfMassEnergy) == 7: @@ -127,7 +148,7 @@ def main(): file_name = fileNames['8TeV'][options.sample] else: print "Error: Unrecognised centre of mass energy." - + print file_name # Output file name outputFileName = 'crap.root' outputFileDir = 'unfolding/' @@ -139,11 +160,17 @@ def main(): elif options.pdfWeight != 0: outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_pdfWeight_%i.root' % ( energySuffix, options.pdfWeight ) elif options.sample != 'central': - outputFileName = outputFileDir+'/unfolding_TTJets_%s_%s_asymmetric.root' % ( energySuffix, options.sample ) + if options.eightToSeven : + outputFileName = outputFileDir+'/unfolding_TTJets_%s_%s_8To7_asymmetric.root' % ( energySuffix, options.sample ) + else : + outputFileName = outputFileDir+'/unfolding_TTJets_%s_%s_asymmetric.root' % ( energySuffix, options.sample ) elif options.fineBinned: outputFileName = outputFileDir+'/unfolding_TTJets_%s.root' % ( energySuffix ) else: - outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric.root' % energySuffix + if options.eightToSeven : + outputFileName = outputFileDir+'/unfolding_TTJets_%s_8To7_asymmetric.root' % energySuffix + else : + outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric.root' % energySuffix with root_open( file_name, 'read' ) as f, root_open( outputFileName, 'recreate') as out: @@ -174,6 +201,7 @@ def main(): for variable in bin_edges: if options.debug and variable != 'HT' : continue + if variable == 'MT' : continue print '--->Doing variable :',variable @@ -207,6 +235,10 @@ def main(): offlineWeight += ' * '+pdfWeight genWeight += ' * '+pdfWeight pass + + if options.eightToSeven: + genWeight += ' * ( unfolding.comWeight )' + offlineWeight += ' * ( unfolding.comWeight )' # Scale factors # scaleFactor = getScaleFactor( options.centreOfMassEnergy, channel.channelName ) @@ -218,7 +250,24 @@ def main(): truth = Hist( bin_edges[variable], name='truth') measured = Hist( bin_edges[variable], name='measured') fake = Hist( bin_edges[variable], name='fake') + + h_nJets = Hist([-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5], name='nJets') + h_nJetsG20 = Hist([-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5], name='nJetsG20') + h_nJetsG30 = Hist([-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5], name='nJetsG30') + h_htG20 = Hist([i*20 for i in range(5,26)], name='htG20') + h_htG30 = Hist([i*20 for i in range(5,26)], name='htG30') + + jetPt = Hist([i*2.5 for i in range(8,20)], name='jetPt') + + jetPt1 = Hist([i*5 for i in range(3,40)], name='jetPt1') + jetPt2 = Hist([i*5 for i in range(3,40)], name='jetPt2') + jetPt3 = Hist([i*5 for i in range(3,40)], name='jetPt3') + jetPt4 = Hist([i*5 for i in range(3,40)], name='jetPt4') + jetPt5 = Hist([i*5 for i in range(3,40)], name='jetPt5') + jetPt6 = Hist([i*5 for i in range(3,40)], name='jetPt6') + jetPt7 = Hist([i*5 for i in range(3,40)], name='jetPt7') + # 2D histograms response = Hist2D( bin_edges[variable], bin_edges[variable], name='response') response_without_fakes = Hist2D( bin_edges[variable], bin_edges[variable], name='response_without_fakes') @@ -244,6 +293,10 @@ def main(): tree.Draw(genVariable,genWeight+'*'+genSelection,hist=truth) tree.Draw(recoVariable,offlineWeight+'*'+offlineSelection,hist=measured) tree.Draw(recoVariable,offlineWeight+'*'+fakeSelection,hist=fake) + # tree.Draw("@(unfolding.jetPt).size()",genWeight+'*'+genSelection,hist=nJets) + # tree.Draw('unfolding.jetPt',genWeight+'*'+genSelection,hist=jetPt) + # tree.Draw('unfolding.recoNJets',offlineSelection+'&&'+genSelection,hist=h_nJets) + # 2D tree.Draw(recoVariable+':'+genVariable,offlineWeight+'*'+offlineSelection,hist=response) tree.Draw(recoVariable+':'+genVariable,offlineWeight+'* ('+offlineSelection+'&&'+genSelection +')',hist=response_without_fakes) @@ -252,11 +305,85 @@ def main(): if options.extraHists: tree.Draw( 'unfolding.puWeight','unfolding.OfflineSelection',hist=puOffline) pass - + + + # tree.SetBranchStatus('*',0) + # tree.SetBranchStatus('unfolding.jetPt',1) + # tree.SetBranchStatus('unfolding.GenSelection',1) + # tree.SetBranchStatus('unfolding.OfflineSelection',1) + # print tree.GetEntries() + # nEvents = 0 + # for event in tree: + # nEvents += 1 + # genSelection = event.__getattr__('unfolding.GenSelection') + # offlineSelection = event.__getattr__('unfolding.OfflineSelection') + + # if genSelection == 1 and offlineSelection == 1 : + # jets = event.__getattr__('unfolding.jetPt') + # nJetsG20 = 0 + # nJetsG30 = 0 + # htG20 = 0 + # htG30 = 0 + + # if len(jets) > 0 : + # jetPt1.Fill(jets[0]) + + # if len(jets) > 1: + # jetPt2.Fill(jets[1]) + + # if len(jets) > 2: + # jetPt3.Fill(jets[2]) + + # if len(jets) > 3: + # jetPt4.Fill(jets[3]) + + # if len(jets) > 4: + # jetPt5.Fill(jets[4]) + + # if len(jets) > 5: + # jetPt6.Fill(jets[5]) + + # if len(jets) > 6: + # jetPt7.Fill(jets[6]) + # for jetPt in jets: + # if jetPt >= 20: + # nJetsG20 += 1 + # htG20 += jetPt + # if jetPt >= 30: + # nJetsG30 += 1 + # htG30 += jetPt + # h_nJetsG20.Fill(nJetsG20) + # h_nJetsG30.Fill(nJetsG30) + # h_htG20.Fill(htG20) + # h_htG30.Fill(htG30) + # if nEvents > 100000 : break # Output histgorams to file outputDir.cd() + # truth.Scale(1/truth.Integral()) truth.Write() measured.Write() + + # h_nJetsG20.Scale(1/h_nJetsG20.Integral()) + # h_nJetsG30.Scale(1/h_nJetsG30.Integral()) + # h_htG20.Scale(1/h_htG20.Integral()) + # h_htG30.Scale(1/h_htG30.Integral()) + + h_nJets.Write() + h_nJetsG20.Write() + h_nJetsG30.Write() + h_htG20.Write() + h_htG30.Write() + + jetPt1.Write() + jetPt2.Write() + jetPt3.Write() + jetPt4.Write() + jetPt5.Write() + jetPt6.Write() + jetPt7.Write() + + # jetPt.Scale( 1/jetPt.Integral() ) + # jetPt.Write() fake.Write() response.Write() response_without_fakes.Write() diff --git a/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runCondor.sh b/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runCondor.sh index 72923475..4416ec20 100755 --- a/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runCondor.sh +++ b/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runCondor.sh @@ -17,13 +17,13 @@ git submodule init && git submodule update ./setup.sh eval `scramv1 runtime -sh` . environment.sh -rm -r bltUnfoldAsymm -mkdir bltUnfoldAsymm +rm -r unfolding +mkdir unfolding time python experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runJobsCrab.py -j $1 echo "bltUnfoldAsymm folder contents:" -ls bltUnfoldAsymm -tar -cf bltUnfoldAsymm.$2.$1.tar bltUnfoldAsymm +ls unfolding +tar -cf bltUnfoldAsymm.$2.$1.tar unfolding mv bltUnfoldAsymm.$2.$1.tar ../../../ echo "DailyPythonScripts folder contents:" ls diff --git a/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runJobsCrab.py b/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runJobsCrab.py index 0b3ab0a4..69aa51e7 100644 --- a/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runJobsCrab.py +++ b/experimental/BLTUnfold/createAsymmetricBinningUnfoldingFiles/runJobsCrab.py @@ -4,57 +4,60 @@ jobs = [ # 8 TeV # Central - '--centreOfMassEnergy 8 -s central', + # '--centreOfMassEnergy 8 -s central', - # Scale up/down - '--centreOfMassEnergy 8 -s scaleup', - '--centreOfMassEnergy 8 -s scaledown', + # # Scale up/down + # '--centreOfMassEnergy 8 -s scaleup', + # '--centreOfMassEnergy 8 -s scaledown', - # Matching up/down - '--centreOfMassEnergy 8 -s matchingup', - '--centreOfMassEnergy 8 -s matchingdown', + # # Matching up/down + # '--centreOfMassEnergy 8 -s matchingup', + # '--centreOfMassEnergy 8 -s matchingdown', - # # Other generators - '--centreOfMassEnergy 8 -s powheg', - '--centreOfMassEnergy 8 -s powhegherwig', - '--centreOfMassEnergy 8 -s mcatnlo', + # # # Other generators + # '--centreOfMassEnergy 8 -s powheg', + # '--centreOfMassEnergy 8 -s powhegherwig', + # '--centreOfMassEnergy 8 -s mcatnlo', - # Mass up/down - '--centreOfMassEnergy 8 -s massup', - '--centreOfMassEnergy 8 -s massdown', + # # Mass up/down + # '--centreOfMassEnergy 8 -s massup', + # '--centreOfMassEnergy 8 -s massdown', - # Top pt reweighting - '--centreOfMassEnergy 8 --topPtReweighting', + # # Top pt reweighting + # '--centreOfMassEnergy 8 --topPtReweighting', # 7 TeV # Central - '--centreOfMassEnergy 7 -s central', + # '--centreOfMassEnergy 7 -s central', # Scale up/down - '--centreOfMassEnergy 7 -s scaleup', - '--centreOfMassEnergy 7 -s scaledown', + # '--centreOfMassEnergy 7 -s scaleup', + # '--centreOfMassEnergy 7 -s scaledown', # Matching up/down - '--centreOfMassEnergy 7 -s matchingup', - '--centreOfMassEnergy 7 -s matchingdown', + # '--centreOfMassEnergy 7 -s matchingup', + # '--centreOfMassEnergy 7 -s matchingdown', # # Other generators - '--centreOfMassEnergy 7 -s powheg', - '--centreOfMassEnergy 7 -s powhegherwig', + # '--centreOfMassEnergy 7 -s powheg', + # '--centreOfMassEnergy 7 -s powhegherwig', - # Mass up/down - '--centreOfMassEnergy 7 -s massup', - '--centreOfMassEnergy 7 -s massdown', - - # Top pt reweighting - '--centreOfMassEnergy 7 --topPtReweighting', + # # Mass up/down + # '--centreOfMassEnergy 7 -s massup', + # '--centreOfMassEnergy 7 -s massdown', + + '--centreOfMassEnergy 8 -s massup --eightToSeven', + '--centreOfMassEnergy 8 -s massdown --eightToSeven', + + # # Top pt reweighting + # '--centreOfMassEnergy 7 --topPtReweighting', ] # Add pdf variations to list of jobs -for variation in range(1,45+1): - jobs.append('--centreOfMassEnergy 8 -p %i' % variation) - jobs.append('--centreOfMassEnergy 7 -p %i' % variation) - pass +# for variation in range(1,45+1): +# # jobs.append('--centreOfMassEnergy 8 -p %i' % variation) +# jobs.append('--centreOfMassEnergy 7 -p %i' % variation) +# pass # print len(jobs) parser = OptionParser() diff --git a/experimental/BLTUnfold/mergeUnfoldingBLT/merge_unfolding_BLT_files_on_DICE.py b/experimental/BLTUnfold/mergeUnfoldingBLT/merge_unfolding_BLT_files_on_DICE.py index 1ea813d8..d0f1e4c9 100644 --- a/experimental/BLTUnfold/mergeUnfoldingBLT/merge_unfolding_BLT_files_on_DICE.py +++ b/experimental/BLTUnfold/mergeUnfoldingBLT/merge_unfolding_BLT_files_on_DICE.py @@ -32,7 +32,8 @@ # first set up lists of jobs to be run -sample_to_BLT_filepath_dictionary_7TeV = {"central":"TTJets_MSDecays_central_TuneZ2_7TeV-madgraph-tauola/crab_TTJets_central_7TeV_madgraph_BLTUnfold_NoSkim", +sample_to_BLT_filepath_dictionary_7TeV = { + "central":"TTJets_MSDecays_central_TuneZ2_7TeV-madgraph-tauola/crab_TTJets_central_7TeV_madgraph_BLTUnfold_NoSkim", "scaleup": "TTJets_MSDecays_scaleup_TuneZ2star_7TeV-madgraph-tauola/crab_TTJets_scaleup_7TeV_madgraph_BLTUnfold_NoSkim", "scaledown": "TTJets_MSDecays_scaledown_TuneZ2star_7TeV-madgraph-tauola/crab_TTJets_scaledown_7TeV_madgraph_BLTUnfold_NoSkim", "matchingup": "TTJets_MSDecays_matchingup_7TeV-madgraph-tauola/crab_TTJets_matchingup_7TeV_madgraph_BLTUnfold_NoSkim", @@ -43,16 +44,17 @@ "powhegherwig": "TT_weights_CT10_AUET2_7TeV-powheg-herwig/crab_TT_CT10_AUET2_7TeV_powheg_herwig_BLTUnfold_NoSkim_2" } -sample_to_BLT_filepath_dictionary_8TeV = {"central":"TTJets_MassiveBinDECAY_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_central_8TeV_madgraph_BLTUnfold_NoSkim", - "scaleup": "TTJets_scaleup_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_scaleup_8TeV_madgraph_BLTUnfold_NoSkim", - "scaledown": "TTJets_scaledown_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_scaledown_8TeV_madgraph_BLTUnfold_NoSkim", - "matchingup": "TTJets_matchingup_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_matchingup_8TeV_madgraph_BLTUnfold_NoSkim", - "matchingdown": "TTJets_matchingdown_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_matchingdown_8TeV_madgraph_BLTUnfold_NoSkim", +sample_to_BLT_filepath_dictionary_8TeV = { + # "central":"TTJets_MassiveBinDECAY_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_central_8TeV_madgraph_BLTUnfold_NoSkim", + # "scaleup": "TTJets_scaleup_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_scaleup_8TeV_madgraph_BLTUnfold_NoSkim", + # "scaledown": "TTJets_scaledown_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_scaledown_8TeV_madgraph_BLTUnfold_NoSkim", + # "matchingup": "TTJets_matchingup_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_matchingup_8TeV_madgraph_BLTUnfold_NoSkim", + # "matchingdown": "TTJets_matchingdown_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_matchingdown_8TeV_madgraph_BLTUnfold_NoSkim", "mass_173_5":"TTJets_MSDecays_mass173_5_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_mass_173_5_8TeV_madgraph_BLTUnfold_NoSkim", "mass_169_5": "TTJets_MSDecays_mass169_5_TuneZ2star_8TeV-madgraph-tauola/crab_TTJets_mass_169_5_8TeV_madgraph_BLTUnfold_NoSkim", - "powhegpythia": "TT_CT10_TuneZ2star_8TeV-powheg-tauola/crab_TT_CT10_8TeV_powheg_tauola_BLTUnfold_NoSkim", - "powhegherwig": "TT_CT10_AUET2_8TeV-powheg-herwig/crab_TT_CT10_AUET2_8TeV_powheg_herwig_BLTUnfold_NoSkim", - "mcatnlo": "TT_8TeV-mcatnlo/crab_TT_8TeV_mcatnlo_BLTUnfold_NoSkim" + # "powhegpythia": "TT_CT10_TuneZ2star_8TeV-powheg-tauola/crab_TT_CT10_8TeV_powheg_tauola_BLTUnfold_NoSkim", + # "powhegherwig": "TT_CT10_AUET2_8TeV-powheg-herwig/crab_TT_CT10_AUET2_8TeV_powheg_herwig_BLTUnfold_NoSkim", + # "mcatnlo": "TT_8TeV-mcatnlo/crab_TT_8TeV_mcatnlo_BLTUnfold_NoSkim" } if options.com == 7: @@ -90,9 +92,10 @@ # get current working directory current_working_directory = os.getcwd() -#make_folder_if_not_exists(current_working_directory + "/mergeBLTUnfold/") +# make_folder_if_not_exists(current_working_directory + "/mergeBLTUnfold/") # if the output file doesn't already exist, merge! +print current_working_directory #if not os.path.exists( current_working_directory + "/mergeBLTUnfold/" + output_file): if not os.path.exists( current_working_directory + "/" + output_file): merge_ROOT_files(input_files, output_file, compression = 7, waitToFinish=True) @@ -103,17 +106,18 @@ # HAVE NOT CHECKED IF THE FOLLOWING COMMENTED CODE WORKS. # Now move output file to hdfs # Check if file already exists on hdfs -# if os.path.exists( output_file_hdfs ): -# print "Output file on hdfs already exists. Removing and replacing with new version." -# command = 'hadoop fs -rm -skipTrash ' + output_file_hdfs.split('/hdfs')[-1] -# p = subprocess.Popen(command, shell=True) -# p.wait() - -# print '\nStarting rsync' -# output_log_file = output_file.replace(".root", ".log") -# command = 'rsync --verbose --progress --stats --compress --recursive --times --update %s %s >> %s' % (output_file,output_file_hdfs,output_log_file) -# print command -# p = subprocess.Popen(command, shell=True) -# p.wait() +output_file_hdfs = '/hdfs/TopQuarkGroup/mc/7TeV/v11/NoSkimUnfolding/BLT/temp/' + output_file +if os.path.exists( output_file_hdfs ): + print "Output file on hdfs already exists. Removing and replacing with new version." + command = 'hadoop fs -rm -skipTrash ' + output_file_hdfs.split('/hdfs')[-1] + p = subprocess.Popen(command, shell=True) + p.wait() + +print '\nStarting rsync' +output_log_file = output_file.replace(".root", ".log") +command = 'rsync --verbose --progress --stats --compress --recursive --times --update %s %s >> %s' % (output_file,output_file_hdfs,output_log_file) +print command +p = subprocess.Popen(command, shell=True) +p.wait() print 'ALL DONE!' diff --git a/experimental/checkHT.py b/experimental/checkHT.py new file mode 100644 index 00000000..0fb8f480 --- /dev/null +++ b/experimental/checkHT.py @@ -0,0 +1,210 @@ +from rootpy.tree import Tree +from rootpy.plotting import Hist, Hist2D, Canvas +from rootpy.io import root_open, File +from ROOT import gPad + +powhegHerwigFile = '/storage/ec6821/NTupleProd/CMSSW_5_3_23/src/TTJets_PowhegHerwig_8TeV.root' +powhegV2HerwigppFile = '/storage/ec6821/NTupleProd/CMSSW_5_3_23/src/TTJets_newPowhegHerwig_8TeV.root' +madgraphFile = '/storage/ec6821/NTupleProd/CMSSW_5_3_23/src/TTJets_Madgraph_8TeV.root' + +powhegHerwigHists = {} +powhegV2HerwigppHists = {} +madgraphHists = {} + +htBins = [120.0, 185.0, 215.0, 247.0, 283.0, 323.0, 365.0, 409.0, 458.0, 512.0, 570.0, 629.0, 691.0, 769.0, 1000.0] +for file in [powhegHerwigFile,powhegV2HerwigppFile,madgraphFile]: + + histStorage = None + if file == powhegHerwigFile: + histStorage = powhegHerwigHists + elif file == powhegV2HerwigppFile: + histStorage = powhegV2HerwigppHists + elif file == madgraphFile: + histStorage = madgraphHists + + h_genHT20 = Hist(htBins, name='genHT20') + h_genNJets20 = Hist(10,0.5,10.5, name='genNJets20') + h_genHT30 = Hist(htBins, name='genHT30') + h_genNJets30 = Hist(10,0.5,10.5, name='genNJets30') + + h_HT20 = Hist(htBins, name='HT20') + h_nJets20 = Hist(10,0.5,10.5, name='nJets20') + h_HT30 = Hist(htBins, name='HT30') + h_nJets30 = Hist(10,0.5,10.5, name='nJets30') + + h_jetPt = [ + Hist(16,20,100,name='jet1Pt'), + Hist(16,20,100,name='jet2Pt'), + Hist(16,20,100,name='jet3Pt'), + Hist(16,20,100,name='jet4Pt'), + Hist(16,20,100,name='jet5Pt'), + Hist(16,20,100,name='jet6Pt'), + Hist(16,20,100,name='jet7Pt'), + ] + + h_recoJetPt = Hist(16,20,100,name='recoJetPt') + + with root_open( file, 'read' ) as f: + tree = f.Get('rootTupleTreeEPlusJets/ePlusJetsTree') + + tree.SetBranchStatus('*',0) + tree.SetBranchStatus('unfolding.genHT',1) + tree.SetBranchStatus('unfolding.genNJets',1) + tree.SetBranchStatus('unfolding.GenSelection',1) + tree.SetBranchStatus('unfolding.OfflineSelection',1) + tree.SetBranchStatus('unfolding.genJetPt',1) + tree.SetBranchStatus('unfolding.jetPt',1) + tree.SetBranchStatus('unfolding.jetPtSmear',1) + + print tree.GetEntries() + nEvents = 0 + for event in tree: + + passesOffline = event.__getattr__('unfolding.OfflineSelection') + passesGen = event.__getattr__('unfolding.GenSelection') + if not ( passesOffline and passesGen ) : continue + + nEvents += 1 + genHT = event.__getattr__('unfolding.genHT') + genNJets = event.__getattr__('unfolding.genNJets') + genJetPt = event.__getattr__('unfolding.genJetPt') + jetPt = event.__getattr__('unfolding.jetPt') + jetPtSmear = event.__getattr__('unfolding.jetPtSmear') + + genHT20 = 0 + genNJet20 = 0 + genHT30 = 0 + genNJet30 = 0 + genJetIndex = 0 + for pt in genJetPt: + if pt > 20 : + genNJet20 += 1 + genHT20 += pt + if genJetIndex < 7: + h_jetPt[genJetIndex].Fill(pt) + genJetIndex += 1 + + + if pt > 30 : + genNJet30 += 1 + genHT30 += pt + + HT20 = 0 + nJet20 = 0 + HT30 = 0 + nJet30 = 0 + for pt, ptSmear in zip( jetPt, jetPtSmear ): + smearedPt = pt * ptSmear + if smearedPt > 20 : + nJet20 += 1 + HT20 += smearedPt + h_recoJetPt.Fill( smearedPt ) + + if smearedPt > 30 : + nJet30 += 1 + HT30 += smearedPt + + if nJet30 < 4 : + print "Fewer than 4 reco jets above 30 : ",nJet30 + h_genHT20.Fill(genHT20) + h_genNJets20.Fill(genNJet20) + h_genHT30.Fill( genHT30 ) + h_genNJets30.Fill( genNJet30 ) + + h_HT20.Fill( HT20 ) + h_nJets20.Fill( nJet20 ) + h_HT30.Fill( HT30 ) + h_nJets30.Fill( nJet30 ) + + if nEvents > 5939: + break + histStorage['genHT20'] = h_genHT20.Clone() + histStorage['genNJets20'] = h_genNJets20.Clone() + histStorage['genHT30'] = h_genHT30.Clone() + histStorage['genNJets30'] = h_genNJets30.Clone() + + histStorage['HT20'] = h_HT20.Clone() + histStorage['nJets20'] = h_nJets20.Clone() + histStorage['HT30'] = h_HT30.Clone() + histStorage['nJets30'] = h_nJets30.Clone() + + histStorage['recoJetPt'] = h_recoJetPt.Clone() + + for jetIndex in range(0,7): + histStorage['jetPt%i' % jetIndex ] = h_jetPt[jetIndex].Clone() + +canvas = Canvas(width=700,height=500) +for name,h in powhegHerwigHists.iteritems(): + # h.Scale(1/h.Integral()) + h.linecolor = 'red' + h.linewidth = 1 + +for name,h in powhegV2HerwigppHists.iteritems(): + # h.Scale(1/h.Integral()) + h.linecolor = 'blue' + h.linewidth = 1 + +for name,h in madgraphHists.iteritems(): + # h.Scale(1/h.Integral()) + h.linecolor = 'green' + h.linewidth = 1 + +print powhegHerwigHists['genHT20'].Integral() +print powhegV2HerwigppHists['genHT20'].Integral() +print madgraphHists['genHT20'].Integral() +powhegHerwigHists['genHT20'].Draw('HIST X0') +powhegV2HerwigppHists['genHT20'].Draw('SAME HIST X0') +madgraphHists['genHT20'].Draw('SAME HIST X0') +raw_input('...') + +powhegHerwigHists['genNJets20'].Draw('HIST X0') +powhegV2HerwigppHists['genNJets20'].Draw('SAME HIST X0') +madgraphHists['genNJets20'].Draw('SAME HIST X0') +raw_input('...') + +# powhegHerwigHists['genHT30'].Draw('HIST X0') +# powhegV2HerwigppHists['genHT30'].Draw('SAME HIST X0') +# madgraphHists['genHT30'].Draw('SAME HIST X0') +# raw_input('...') + +# powhegHerwigHists['genNJets30'].Draw('HIST X0') +# powhegV2HerwigppHists['genNJets30'].Draw('SAME HIST X0') +# madgraphHists['genNJets30'].Draw('SAME HIST X0') +# raw_input('...') + +powhegV2HerwigppHists['HT20'].Draw('HIST X0') +powhegHerwigHists['HT20'].Draw('SAME HIST X0') +madgraphHists['HT20'].Draw('SAME HIST X0') +raw_input('...') + +# powhegHerwigHists['HT30'].Draw('HIST X0') +# powhegV2HerwigppHists['HT30'].Draw('SAME HIST X0') +# madgraphHists['HT30'].Draw('SAME HIST X0') +# raw_input('...') + +# powhegHerwigHists['nJets20'].Draw('HIST X0') +# powhegV2HerwigppHists['nJets20'].Draw('SAME HIST X0') +# madgraphHists['nJets20'].Draw('SAME HIST X0') +# raw_input('...') + +# powhegHerwigHists['nJets30'].Draw('HIST X0') +# powhegV2HerwigppHists['nJets30'].Draw('SAME HIST X0') +# madgraphHists['nJets30'].Draw('SAME HIST X0') + +# raw_input('...') + +# powhegHerwigHists['recoJetPt'].Draw('HIST X0') +# powhegV2HerwigppHists['recoJetPt'].Draw('SAME HIST X0') +# madgraphHists['recoJetPt'].Draw('SAME HIST X0') +# raw_input('...') + +# gPad.SetLogy() +for jetIndex in range(0,7): + print jetIndex + powhegV2HerwigppHists['jetPt%i' % jetIndex ].Draw('HIST X0') + powhegHerwigHists['jetPt%i' % jetIndex ].Draw('SAME HIST X0') + madgraphHists['jetPt%i' % jetIndex ].Draw('SAME HIST X0') + print powhegHerwigHists['jetPt%i' % jetIndex ].Integral() + print powhegV2HerwigppHists['jetPt%i' % jetIndex ].Integral() + print madgraphHists['jetPt%i' % jetIndex ].Integral() + raw_input('...') \ No newline at end of file diff --git a/experimental/compareCOMReweight.py b/experimental/compareCOMReweight.py new file mode 100644 index 00000000..c549b7cd --- /dev/null +++ b/experimental/compareCOMReweight.py @@ -0,0 +1,118 @@ +# imports +from config.variable_binning import bin_edges +from rootpy.io import root_open, File +from rootpy.plotting import Canvas, Pad, Legend +from config.latex_labels import variables_latex +from ROOT import gStyle, gROOT +gStyle.SetOptStat(0) +gROOT.SetBatch() + +def normaliseHistogramToOne( histogram ): + normalisationFactor = histogram.Integral() + if normalisationFactor > 0: + histogram.Scale( 1 / normalisationFactor) + +def setHistogramColour( histogram, colour ): + histogram.linecolor = colour + histogram.markercolor = colour + +def getRatioOfFirstToSecond( first_histogram, second_histogram ): + return first_histogram / second_histogram + +def setupPads(): + canvas = Canvas(width=700,height=500) + plot_pad = Pad( xlow = 0., ylow = 0.3, xup = 1., yup = 1.) + plot_pad.SetBottomMargin(0.01) + ratio_pad = Pad( xlow = 0., ylow = 0., xup = 1., yup = 0.3) + ratio_pad.SetTopMargin(0.04) + ratio_pad.SetBottomMargin(0.2) + + return canvas, plot_pad, ratio_pad + +# files +seven_tev_file_name = '/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_8th_draft/7TeV/unfolding/unfolding_TTJets_7TeV_asymmetric.root' +eight_tev_file_name = '/hdfs/TopQuarkGroup/results/histogramfiles/AN-14-071_8th_draft/8TeV/unfolding/unfolding_TTJets_8TeV_asymmetric.root' +reweighted_eight_tev_file_name = '/storage/ec6821/DailyPythonScripts/legacy/CMSSW_7_4_7/src/DailyPythonScripts/unfolding/unfolding_TTJets_8TeV_8To7_asymmetric.root' +output_file_directory = 'plots/compareCOMReweight/' +# channels +channels = ['electron', 'muon'] + +with root_open( seven_tev_file_name, 'read' ) as seven_tev_file, \ + root_open( eight_tev_file_name, 'read' ) as eight_tev_file, \ + root_open( reweighted_eight_tev_file_name, 'read' ) as reweighted_eight_tev_file: + + for channel in channels: + for variable in bin_edges: + if variable == 'MT' : continue + print variable + + directory_template = 'unfolding_{variable}_analyser_{channel}_channel' + if variable != 'HT': + directory_template += '_patType1CorrectedPFMet' + + directory_name = directory_template.format( + variable = variable, + channel = channel, + ) + + for distribution in ['truth', 'measured']: + + truth_seven_tev = seven_tev_file.Get( directory_name + '/' + distribution ) + truth_eight_tev = eight_tev_file.Get( directory_name + '/' + distribution ) + truth_reweighted_eight_tev = reweighted_eight_tev_file.Get( directory_name + '/' + distribution ) + + normaliseHistogramToOne( truth_seven_tev ) + normaliseHistogramToOne( truth_eight_tev ) + normaliseHistogramToOne( truth_reweighted_eight_tev ) + + ratio_reweighted_eight_to_seven = getRatioOfFirstToSecond( truth_reweighted_eight_tev, truth_seven_tev ) + ratio_eight_to_seven = getRatioOfFirstToSecond( truth_eight_tev, truth_seven_tev ) + + canvas, plot_pad, ratio_pad = setupPads() + plot_pad.Draw() + ratio_pad.Draw() + setHistogramColour( truth_seven_tev, 'red' ) + setHistogramColour( truth_eight_tev, 'green' ) + setHistogramColour( truth_reweighted_eight_tev, 'blue' ) + + plot_pad.cd() + truth_seven_tev.Draw() + truth_seven_tev.yaxis.SetTitle('Arbitrary Units') + truth_seven_tev.xaxis.SetLabelSize(0) + truth_eight_tev.Draw('SAME') + truth_reweighted_eight_tev.Draw('SAME') + + setHistogramColour( ratio_reweighted_eight_to_seven, 'blue') + setHistogramColour( ratio_eight_to_seven, 'green') + + legend = Legend(3) + legend.AddEntry(truth_seven_tev, label='7 TeV', style='LEP') + legend.AddEntry(truth_eight_tev, label='8 TeV', style='LEP') + legend.AddEntry(truth_reweighted_eight_tev, label='Reweighted 8 TeV', style='LEP') + legend.SetBorderSize(0) + legend.Draw() + + ratio_pad.cd() + ratio_reweighted_eight_to_seven.SetLabelSize(0.1,'X') + ratio_reweighted_eight_to_seven.SetTitleSize(0.1, 'X') + ratio_reweighted_eight_to_seven.SetTitleOffset(0, 'X') + ratio_reweighted_eight_to_seven.xaxis.SetTitle(variables_latex[variable] + '/[GeV]') + + ratio_reweighted_eight_to_seven.yaxis.SetRangeUser(0.9,1.1) + # ratio_reweighted_eight_to_seven.yaxis.SetNdivisions(400,1) + ratio_reweighted_eight_to_seven.SetLabelSize(0.1,'Y') + ratio_reweighted_eight_to_seven.yaxis.SetNdivisions(402); + ratio_reweighted_eight_to_seven.yaxis.SetTitleSize(20); + ratio_pad.SetGridy() + ratio_reweighted_eight_to_seven.Draw() + ratio_eight_to_seven.Draw('SAME') + # plot_pad.Update() + # ratio_pad.Update() + # canvas.Update() + output_pdf_template = '{directory}/{variable}_{channel}_{distribution}.pdf' + output_pdf_filename = output_pdf_template.format( directory = output_file_directory, variable = variable, channel = channel, distribution = distribution) + canvas.Update() + canvas.SaveAs(output_pdf_filename) + # raw_input('...') + + diff --git a/experimental/mergeBATOutputFilesOnDICE/merge_samples_onDICE.py b/experimental/mergeBATOutputFilesOnDICE/merge_samples_onDICE.py index 6f9d3dac..ec8b5ddb 100755 --- a/experimental/mergeBATOutputFilesOnDICE/merge_samples_onDICE.py +++ b/experimental/mergeBATOutputFilesOnDICE/merge_samples_onDICE.py @@ -56,13 +56,20 @@ # Make list of all samples to be merged allJobs = [] for category in config.categories_and_prefixes.keys(): + # if category != 'central' : + # continue for sample, input_samples in sample_summations.iteritems(): # Only consider these samples - if not sample in ['WJets', 'DYJets', 'VJets-matchingup', - 'VJets-matchingdown', 'VJets-scaleup', - 'VJets-scaledown','QCD_Electron', - 'QCD_Muon', 'VJets', - 'SingleTop']: # + if not sample in [ + 'WJets', + 'DYJets', + # 'VJets-matchingup', + # 'VJets-matchingdown', 'VJets-scaleup', + # 'VJets-scaledown','QCD_Electron', + 'QCD_Muon', + 'VJets', + 'SingleTop' + ]: # continue # Only consider these samples for central if sample in ['WJets', 'DYJets', 'VJets-matchingup', @@ -100,6 +107,8 @@ output_file = output_file_hdfs.replace("/hdfs/TopQuarkGroup/results/histogramfiles", current_working_directory) input_files = [config.general_category_templates[category] % input_sample for input_sample in input_samples] +print input_files +print output_file if not os.path.exists( output_file ): merge_ROOT_files( input_files, output_file, compression = 7, waitToFinish=True ) print "merging ", sample @@ -116,6 +125,7 @@ print '\nStarting rsync' output_log_file = output_file.replace(".root", ".log") +p = subprocess.Popen('touch %s' % output_log_file, shell=True) command = 'rsync --verbose --progress --stats --compress --recursive --times --update %s %s >> %s' % (output_file,output_file_hdfs,output_log_file) print command p = subprocess.Popen(command, shell=True) diff --git a/experimental/mergeBATOutputFilesOnDICE/move_BAT_output_files_to_hdfs.py b/experimental/mergeBATOutputFilesOnDICE/move_BAT_output_files_to_hdfs.py index 19ca0461..af29839e 100644 --- a/experimental/mergeBATOutputFilesOnDICE/move_BAT_output_files_to_hdfs.py +++ b/experimental/mergeBATOutputFilesOnDICE/move_BAT_output_files_to_hdfs.py @@ -67,9 +67,9 @@ make_folder_if_not_exists( path_to_AN_folder + "/" + category ) #check if folder is empty and exit if not - if len(os.listdir( path_to_AN_folder + category )) != 0: - print "Folder " + path_to_AN_folder + "/" + category + " on hdfs is not empty. Exiting because histogram files with same names will be overwritten." - sys.exit() + # if len(os.listdir( path_to_AN_folder + category )) != 0: + # print "Folder " + path_to_AN_folder + "/" + category + " on hdfs is not empty. Exiting because histogram files with same names will be overwritten." + # sys.exit() print "running command: ", command p = subprocess.Popen(command, shell=True) diff --git a/experimental/mergeBATOutputFilesOnDICE/submitMerge.description b/experimental/mergeBATOutputFilesOnDICE/submitMerge.description index d07480e6..3850724c 100644 --- a/experimental/mergeBATOutputFilesOnDICE/submitMerge.description +++ b/experimental/mergeBATOutputFilesOnDICE/submitMerge.description @@ -3,7 +3,7 @@ Universe = vanilla Output = mergeBAToutput.job.$(cluster).$(process).out Error = mergeBAToutput.job.$(cluster).$(process).err Log = mergeBAToutput.job.$(cluster).$(process).log -arguments = $(process) 8 +arguments = $(process) 7 transfer_input_files = dps.tar should_transfer_files = YES @@ -15,4 +15,4 @@ request_memory=2000 # use the ENV that is provided getenv = true -queue 2 +queue 46 diff --git a/src/cross_section_measurement/01_get_ttjet_normalisation.py b/src/cross_section_measurement/01_get_ttjet_normalisation.py index 9de82756..26c5c730 100644 --- a/src/cross_section_measurement/01_get_ttjet_normalisation.py +++ b/src/cross_section_measurement/01_get_ttjet_normalisation.py @@ -98,14 +98,19 @@ def calculate_normalisation(self): if self.measurement.__class__ == tools.measurement.Systematic: self.measurement.scale_histograms() histograms = self.measurement.histograms + for sample, hist in histograms.items(): bin_edges = variable_bin_edges[self.variable] histograms[sample] = rebin_asymmetric(hist, bin_edges) for sample, hist in histograms.items(): # TODO: this should be a list of bin-contents + # self.initial_normalisation[ + # sample] = hist_to_value_error_tuplelist(fix_overflow(hist)) self.initial_normalisation[ - sample] = hist_to_value_error_tuplelist(fix_overflow(hist)) + sample] = hist_to_value_error_tuplelist(hist) + nBins = hist.GetNbinsX() + if self.method == self.BACKGROUND_SUBTRACTION and sample != 'TTJet': self.normalisation[sample] = self.initial_normalisation[sample] @@ -218,6 +223,7 @@ def main(): } measurement_files = glob.glob(input_template.format(**inputs)) for f in measurement_files: + print (f) measurement = tools.measurement.Measurement.fromJSON(f) # for each measurement norm = TTJetNormalisation( diff --git a/src/cross_section_measurement/02_unfold_and_measure.py b/src/cross_section_measurement/02_unfold_and_measure.py index 3403b5ab..7a7da490 100644 --- a/src/cross_section_measurement/02_unfold_and_measure.py +++ b/src/cross_section_measurement/02_unfold_and_measure.py @@ -6,12 +6,17 @@ # rootpy from rootpy.io import File from rootpy.plotting import Hist2D +from rootpy.matrix import Matrix +from rootpy import asrootpy +import root_numpy as rnp +import numpy as np +from ROOT import TH2F # DailyPythonScripts import config.RooUnfold as unfoldCfg from config.variable_binning import bin_widths, bin_edges from config import XSectionConfig from tools.Calculation import calculate_xsection, calculate_normalised_xsection, \ -combine_complex_results +combine_complex_results, calculate_covariance_for_normalised_xsection from tools.hist_utilities import hist_to_value_error_tuplelist, \ value_error_tuplelist_to_hist from tools.Unfolding import Unfolding, get_unfold_histogram_tuple @@ -30,9 +35,38 @@ def unfold_results( results, category, channel, k_value, h_truth, h_measured, h_ unfoldCfg.Hreco = 0 else: unfoldCfg.Hreco = options.Hreco - + h_unfolded_data = unfolding.unfold( h_data ) - + + # unfolding.unfoldObject.GetCov() + + covariance_matrix = None + if category == 'central': + covariance_matrix = asrootpy( unfolding.unfoldObject.Ereco(options.Hreco) ).to_numpy() + print covariance_matrix + + corr = np.array(np.zeros((covariance_matrix.shape[0], covariance_matrix.shape[1]) )) + + for i in range(0,covariance_matrix.shape[0]): + for j in range(0,covariance_matrix.shape[1]): + corr[i,j] = covariance_matrix[i,j] / np.sqrt( covariance_matrix[i,i] * covariance_matrix[j,j] ) + print 'Correlation' + print corr + # for i_row in range(0,covariance_matrix.shape[0]): + # for i_col in range(0,covariance_matrix.shape[0]): + # print i_row,i_col + # print covariance_matrix[i_row,i_col] + # # print h_unfolded_data.GetBinContent(i_row+1),h_unfolded_data.GetBinContent(i_col+1) + # print h_unfolded_data.GetBinError(i_row+1)*h_unfolded_data.GetBinError(i_col+1) + + # print list( h_unfolded_data.y() ) + # unfolding.unfoldObject.ErecoV(options.Hreco).Draw() + # raw_input() + # cov_matrix = + # cov_matrix.Draw('COLZ text') + # print cov_matrix + # raw_input('...') + if options.write_unfolding_objects: # export the D and SV distributions SVD_path = path_to_JSON + '/unfolding_objects/' + channel + '/kv_' + str( k_value ) + '/' @@ -73,7 +107,7 @@ def unfold_results( results, category, channel, k_value, h_truth, h_measured, h_ unfoldingObjectFile.Close() del unfolding - return hist_to_value_error_tuplelist( h_unfolded_data ) + return hist_to_value_error_tuplelist( h_unfolded_data ), covariance_matrix def data_covariance_matrix( data ): values = list( data ) @@ -85,10 +119,10 @@ def data_covariance_matrix( data ): return cov_matrix def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ): - global variable, met_type, path_to_JSON, file_for_unfolding, file_for_powheg_pythia, file_for_powheg_herwig, file_for_ptreweight, files_for_pdfs + global variable, met_type, path_to_JSON, file_for_unfolding, file_for_powheg_v2_pythia, file_for_powheg_v2_herwig, file_for_ptreweight, files_for_pdfs global centre_of_mass, luminosity, ttbar_xsection, load_fakes, method if centre_of_mass == 8: - global file_for_mcatnlo + global file_for_mcatnlo, file_for_powheg_v1_herwig, file_for_powheg_v1_pythia global file_for_matchingdown, file_for_matchingup, file_for_scaledown, file_for_scaleup global file_for_massdown, file_for_massup global ttbar_generator_systematics, ttbar_theory_systematics, pdf_uncertainties @@ -101,8 +135,10 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ): ttbar_theory_systematic_prefix + 'scaleup':file_for_scaleup, ttbar_theory_systematic_prefix + 'massdown':file_for_massdown, ttbar_theory_systematic_prefix + 'massup':file_for_massup, - ttbar_theory_systematic_prefix + 'powheg_pythia':file_for_powheg_pythia, - ttbar_theory_systematic_prefix + 'powheg_herwig':file_for_powheg_herwig, + ttbar_theory_systematic_prefix + 'powheg_v1_pythia':file_for_powheg_v1_pythia, + ttbar_theory_systematic_prefix + 'powheg_v1_herwig':file_for_powheg_v1_herwig, + ttbar_theory_systematic_prefix + 'powheg_v2_pythia':file_for_powheg_v2_pythia, + ttbar_theory_systematic_prefix + 'powheg_v2_herwig':file_for_powheg_v2_herwig, ttbar_theory_systematic_prefix + 'ptreweight':file_for_ptreweight, } @@ -148,7 +184,7 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ): load_fakes = load_fakes ) - h_truth_POWHEG_PYTHIA, _, _, _ = get_unfold_histogram_tuple( inputfile = file_for_powheg_pythia, + h_truth_powheg_v2_pythia, _, _, _ = get_unfold_histogram_tuple( inputfile = file_for_powheg_v2_pythia, variable = variable, channel = channel, met_type = met_type, @@ -157,7 +193,7 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ): luminosity = luminosity, load_fakes = load_fakes ) - h_truth_POWHEG_HERWIG, _, _, _ = get_unfold_histogram_tuple( inputfile = file_for_powheg_herwig, + h_truth_powheg_v2_herwig, _, _, _ = get_unfold_histogram_tuple( inputfile = file_for_powheg_v2_herwig, variable = variable, channel = channel, met_type = met_type, @@ -167,6 +203,8 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ): load_fakes = load_fakes ) h_truth_MCATNLO = None + h_truth_powheg_v1_herwig = None + h_truth_powheg_v1_pythia = None if centre_of_mass == 8: h_truth_MCATNLO, _, _, _ = get_unfold_histogram_tuple( inputfile = file_for_mcatnlo, variable = variable, @@ -177,6 +215,26 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ): luminosity = luminosity, load_fakes = load_fakes ) + h_truth_powheg_v1_herwig, _, _, _ = get_unfold_histogram_tuple( inputfile = file_for_powheg_v1_herwig, + variable = variable, + channel = channel, + met_type = met_type, + centre_of_mass = centre_of_mass, + ttbar_xsection = ttbar_xsection, + luminosity = luminosity, + load_fakes = load_fakes + ) + + h_truth_powheg_v1_pythia, _, _, _ = get_unfold_histogram_tuple( inputfile = file_for_powheg_v1_pythia, + variable = variable, + channel = channel, + met_type = met_type, + centre_of_mass = centre_of_mass, + ttbar_xsection = ttbar_xsection, + luminosity = luminosity, + load_fakes = load_fakes + ) + h_truth_matchingdown, _, _, _ = get_unfold_histogram_tuple( inputfile = file_for_matchingdown, variable = variable, channel = channel, @@ -223,30 +281,34 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ): luminosity = luminosity, load_fakes = load_fakes ) - + print 'Overflow in madgraph :',h_truth.GetBinContent(h_truth.GetNbinsX()),h_truth.GetBinContent(h_truth.GetNbinsX() + 2) MADGRAPH_results = hist_to_value_error_tuplelist( h_truth ) MADGRAPH_ptreweight_results = hist_to_value_error_tuplelist( h_truth_ptreweight ) - POWHEG_PYTHIA_results = hist_to_value_error_tuplelist( h_truth_POWHEG_PYTHIA ) - POWHEG_HERWIG_results = hist_to_value_error_tuplelist( h_truth_POWHEG_HERWIG ) + powheg_v2_pythia_results = hist_to_value_error_tuplelist( h_truth_powheg_v2_pythia ) + powheg_v2_herwig_results = hist_to_value_error_tuplelist( h_truth_powheg_v2_herwig ) MCATNLO_results = None + powheg_v1_herwig_results = None + powheg_v1_pythia_results = None if centre_of_mass == 8: MCATNLO_results = hist_to_value_error_tuplelist( h_truth_MCATNLO ) - + powheg_v1_herwig_results = hist_to_value_error_tuplelist( h_truth_powheg_v1_herwig ) + powheg_v1_pythia_results = hist_to_value_error_tuplelist( h_truth_powheg_v1_pythia ) + matchingdown_results = hist_to_value_error_tuplelist( h_truth_matchingdown ) matchingup_results = hist_to_value_error_tuplelist( h_truth_matchingup ) scaledown_results = hist_to_value_error_tuplelist( h_truth_scaledown ) scaleup_results = hist_to_value_error_tuplelist( h_truth_scaleup ) - TTJet_fit_results_unfolded = unfold_results( TTJet_fit_results, - category, - channel, - k_value, - h_truth, - h_measured, - h_response, - h_fakes, - method - ) + TTJet_fit_results_unfolded, covariance_matrix = unfold_results( TTJet_fit_results, + category, + channel, + k_value, + h_truth, + h_measured, + h_response, + h_fakes, + method + ) normalisation_unfolded = { 'TTJet_measured' : TTJet_fit_results, @@ -254,18 +316,23 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, k_value ): 'MADGRAPH': MADGRAPH_results, 'MADGRAPH_ptreweight': MADGRAPH_ptreweight_results, # other generators - 'POWHEG_PYTHIA': POWHEG_PYTHIA_results, - 'POWHEG_HERWIG': POWHEG_HERWIG_results, + 'powheg_v2_pythia': powheg_v2_pythia_results, + 'powheg_v2_herwig': powheg_v2_herwig_results, # systematics 'matchingdown': matchingdown_results, 'matchingup': matchingup_results, 'scaledown': scaledown_results, 'scaleup': scaleup_results } + if centre_of_mass == 8: normalisation_unfolded['MCATNLO'] = MCATNLO_results + normalisation_unfolded['powheg_v1_herwig'] = powheg_v1_herwig_results + normalisation_unfolded['powheg_v1_pythia'] = powheg_v1_pythia_results + + return normalisation_unfolded, covariance_matrix + - return normalisation_unfolded def calculate_xsections( normalisation, category, channel, k_value = None ): global variable, met_type, path_to_JSON @@ -277,11 +344,16 @@ def calculate_xsections( normalisation, category, channel, k_value = None ): TTJet_xsection_unfolded = calculate_xsection( normalisation['TTJet_unfolded'], luminosity, branching_ratio ) # L in pb1 MADGRAPH_xsection = calculate_xsection( normalisation['MADGRAPH'], luminosity, branching_ratio ) # L in pb1 MADGRAPH_ptreweight_xsection = calculate_xsection( normalisation['MADGRAPH_ptreweight'], luminosity, branching_ratio ) # L in pb1 - POWHEG_PYTHIA_xsection = calculate_xsection( normalisation['POWHEG_PYTHIA'], luminosity, branching_ratio ) # L in pb1 - POWHEG_HERWIG_xsection = calculate_xsection( normalisation['POWHEG_HERWIG'], luminosity, branching_ratio ) # L in pb1 + powheg_v2_pythia_xsection = calculate_xsection( normalisation['powheg_v2_pythia'], luminosity, branching_ratio ) # L in pb1 + powheg_v2_herwig_xsection = calculate_xsection( normalisation['powheg_v2_herwig'], luminosity, branching_ratio ) # L in pb1 MCATNLO_xsection = None + powheg_v1_herwig_xsection = None + powheg_v1_pythia_xsection = None if centre_of_mass == 8: MCATNLO_xsection = calculate_xsection( normalisation['MCATNLO'], luminosity, branching_ratio ) # L in pb1 + powheg_v1_herwig_xsection = calculate_xsection( normalisation['powheg_v1_herwig'], luminosity, branching_ratio ) # L in pb1 + powheg_v1_pythia_xsection = calculate_xsection( normalisation['powheg_v1_pythia'], luminosity, branching_ratio ) # L in pb1 + matchingdown_xsection = calculate_xsection( normalisation['matchingdown'], luminosity, branching_ratio ) # L in pb1 matchingup_xsection = calculate_xsection( normalisation['matchingup'], luminosity, branching_ratio ) # L in pb1 scaledown_xsection = calculate_xsection( normalisation['scaledown'], luminosity, branching_ratio ) # L in pb1 @@ -291,8 +363,8 @@ def calculate_xsections( normalisation, category, channel, k_value = None ): 'TTJet_unfolded' : TTJet_xsection_unfolded, 'MADGRAPH': MADGRAPH_xsection, 'MADGRAPH_ptreweight': MADGRAPH_ptreweight_xsection, - 'POWHEG_PYTHIA': POWHEG_PYTHIA_xsection, - 'POWHEG_HERWIG': POWHEG_HERWIG_xsection, + 'powheg_v2_pythia': powheg_v2_pythia_xsection, + 'powheg_v2_herwig': powheg_v2_herwig_xsection, # systematics 'matchingdown': matchingdown_xsection, 'matchingup': matchingup_xsection, @@ -301,7 +373,9 @@ def calculate_xsections( normalisation, category, channel, k_value = None ): } if centre_of_mass == 8: xsection_unfolded['MCATNLO'] = MCATNLO_xsection - + xsection_unfolded['powheg_v1_herwig'] = powheg_v1_herwig_xsection + xsection_unfolded['powheg_v1_pythia'] = powheg_v1_pythia_xsection + if k_value: filename = path_to_JSON + '/xsection_measurement_results/%s/kv%d/%s/xsection_%s.txt' % ( channel, k_value, category, met_type ) elif not channel == 'combined': @@ -311,17 +385,22 @@ def calculate_xsections( normalisation, category, channel, k_value = None ): write_data_to_JSON( xsection_unfolded, filename ) -def calculate_normalised_xsections( normalisation, category, channel, k_value = None, normalise_to_one = False ): +def calculate_normalised_xsections( normalisation, category, channel, k_value = None, normalise_to_one = False, debug=False ): global variable, met_type, path_to_JSON TTJet_normalised_xsection = calculate_normalised_xsection( normalisation['TTJet_measured'], bin_widths[variable], normalise_to_one ) - TTJet_normalised_xsection_unfolded = calculate_normalised_xsection( normalisation['TTJet_unfolded'], bin_widths[variable], normalise_to_one ) + TTJet_normalised_xsection_unfolded = calculate_normalised_xsection( normalisation['TTJet_unfolded'], bin_widths[variable], normalise_to_one, debug ) MADGRAPH_normalised_xsection = calculate_normalised_xsection( normalisation['MADGRAPH'], bin_widths[variable], normalise_to_one ) MADGRAPH_ptreweight_normalised_xsection = calculate_normalised_xsection( normalisation['MADGRAPH_ptreweight'], bin_widths[variable], normalise_to_one ) - POWHEG_PYTHIA_normalised_xsection = calculate_normalised_xsection( normalisation['POWHEG_PYTHIA'], bin_widths[variable], normalise_to_one ) - POWHEG_HERWIG_normalised_xsection = calculate_normalised_xsection( normalisation['POWHEG_HERWIG'], bin_widths[variable], normalise_to_one ) + powheg_v2_pythia_normalised_xsection = calculate_normalised_xsection( normalisation['powheg_v2_pythia'], bin_widths[variable], normalise_to_one ) + powheg_v2_herwig_normalised_xsection = calculate_normalised_xsection( normalisation['powheg_v2_herwig'], bin_widths[variable], normalise_to_one ) MCATNLO_normalised_xsection = None + powheg_v1_herwig_normalised_xsection = None + powheg_v1_pythia_normalised_xsection = None if centre_of_mass == 8: MCATNLO_normalised_xsection = calculate_normalised_xsection( normalisation['MCATNLO'], bin_widths[variable], normalise_to_one ) + powheg_v1_herwig_normalised_xsection = calculate_normalised_xsection( normalisation['powheg_v1_herwig'], bin_widths[variable], normalise_to_one ) + powheg_v1_pythia_normalised_xsection = calculate_normalised_xsection( normalisation['powheg_v1_pythia'], bin_widths[variable], normalise_to_one ) + matchingdown_normalised_xsection = calculate_normalised_xsection( normalisation['matchingdown'], bin_widths[variable], normalise_to_one ) matchingup_normalised_xsection = calculate_normalised_xsection( normalisation['matchingup'], bin_widths[variable], normalise_to_one ) scaledown_normalised_xsection = calculate_normalised_xsection( normalisation['scaledown'], bin_widths[variable], normalise_to_one ) @@ -331,8 +410,8 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value = 'TTJet_unfolded' : TTJet_normalised_xsection_unfolded, 'MADGRAPH': MADGRAPH_normalised_xsection, 'MADGRAPH_ptreweight': MADGRAPH_ptreweight_normalised_xsection, - 'POWHEG_PYTHIA': POWHEG_PYTHIA_normalised_xsection, - 'POWHEG_HERWIG': POWHEG_HERWIG_normalised_xsection, + 'powheg_v2_pythia': powheg_v2_pythia_normalised_xsection, + 'powheg_v2_herwig': powheg_v2_herwig_normalised_xsection, # systematics 'matchingdown': matchingdown_normalised_xsection, 'matchingup': matchingup_normalised_xsection, @@ -341,18 +420,19 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value = } if centre_of_mass == 8: normalised_xsection['MCATNLO'] = MCATNLO_normalised_xsection + normalised_xsection['powheg_v1_herwig'] = powheg_v1_herwig_normalised_xsection + normalised_xsection['powheg_v1_pythia'] = powheg_v1_pythia_normalised_xsection if not channel == 'combined': filename = path_to_JSON + '/xsection_measurement_results/%s/kv%d/%s/normalised_xsection_%s.txt' % ( channel, k_value, category, met_type ) else: filename = path_to_JSON + '/xsection_measurement_results/%s/%s/normalised_xsection_%s.txt' % ( channel, category, met_type ) - if normalise_to_one: filename = filename.replace( 'normalised_xsection', 'normalised_to_one_xsection' ) write_data_to_JSON( normalised_xsection, filename ) if __name__ == '__main__': - set_root_defaults( msg_ignore_level = 3001 ) + set_root_defaults( set_batch=False, msg_ignore_level = 3001 ) # setup parser = OptionParser() parser.add_option( "-p", "--path", dest = "path", default = 'data/', @@ -398,11 +478,15 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value = path_to_files = measurement_config.path_to_files file_for_unfolding = File( measurement_config.unfolding_madgraph, 'read' ) - file_for_powheg_pythia = File( measurement_config.unfolding_powheg_pythia, 'read' ) - file_for_powheg_herwig = File( measurement_config.unfolding_powheg_herwig, 'read' ) + file_for_powheg_v2_pythia = File( measurement_config.unfolding_powheg_v2_pythia, 'read' ) + file_for_powheg_v2_herwig = File( measurement_config.unfolding_powheg_v2_herwig, 'read' ) file_for_mcatnlo = None + file_for_powheg_v1_herwig = None + file_for_powheg_v1_pythia = None if centre_of_mass == 8: file_for_mcatnlo = File( measurement_config.unfolding_mcatnlo, 'read' ) + file_for_powheg_v1_herwig = File( measurement_config.unfolding_powheg_v1_herwig, 'read' ) + file_for_powheg_v1_pythia = File( measurement_config.unfolding_powheg_v1_pythia, 'read' ) file_for_ptreweight = File ( measurement_config.unfolding_ptreweight, 'read' ) files_for_pdfs = { 'PDFWeights_%d' % index : File ( measurement_config.unfolding_pdfweights[index] ) for index in range( 1, 45 ) } @@ -434,7 +518,9 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value = # ttbar theory systematics, including pt reweighting and hadronisation systematic ttbar_theory_systematics = [ ttbar_theory_systematic_prefix + 'ptreweight' ] - ttbar_theory_systematics.extend( [ttbar_theory_systematic_prefix + 'powheg_pythia', ttbar_theory_systematic_prefix + 'powheg_herwig'] ) + ttbar_theory_systematics.extend( [ttbar_theory_systematic_prefix + 'powheg_v2_pythia', ttbar_theory_systematic_prefix + 'powheg_v2_herwig'] ) + if centre_of_mass == 8: + ttbar_theory_systematics.extend( [ttbar_theory_systematic_prefix + 'powheg_v1_pythia', ttbar_theory_systematic_prefix + 'powheg_v1_herwig'] ) categories.extend( ttbar_theory_systematics ) # Add mass systematics @@ -529,11 +615,16 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value = filename = '' # get unfolded normalisation - unfolded_normalisation_electron = get_unfolded_normalisation( TTJet_fit_results_electron, category, 'electron', k_value_electron ) - unfolded_normalisation_muon = get_unfolded_normalisation( TTJet_fit_results_muon, category, 'muon', k_value_muon ) + unfolded_normalisation_electron, covariance_electron = get_unfolded_normalisation( TTJet_fit_results_electron, category, 'electron', k_value_electron ) + unfolded_normalisation_muon, covariance_muon = get_unfolded_normalisation( TTJet_fit_results_muon, category, 'muon', k_value_muon ) + covariance_combined = None if combine_before_unfolding: - unfolded_normalisation_combined = get_unfolded_normalisation( TTJet_fit_results_combined, category, 'combined', k_value_combined ) + unfolded_normalisation_combined, covariance_combined = get_unfolded_normalisation( TTJet_fit_results_combined, category, 'combined', k_value_combined ) else: + if category == 'central': + covariance_combined = covariance_electron + covariance_combined += covariance_muon + # covariance_combined = asrootpy( covariance_combined ) unfolded_normalisation_combined = combine_complex_results( unfolded_normalisation_electron, unfolded_normalisation_muon ) filename = path_to_JSON + '/xsection_measurement_results/electron/kv%d/%s/normalisation_%s.txt' % ( k_value_electron_central, category, met_type ) @@ -542,6 +633,9 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value = write_data_to_JSON( unfolded_normalisation_muon, filename ) filename = path_to_JSON + '/xsection_measurement_results/combined/%s/normalisation_%s.txt' % ( category, met_type ) write_data_to_JSON( unfolded_normalisation_combined, filename ) + + + # if measurement_config.include_higgs: # # now the same for the Higgs @@ -555,23 +649,55 @@ def calculate_normalised_xsections( normalisation, category, channel, k_value = # else: # unfolded_normalisation_combined_higgs = combine_complex_results( unfolded_normalisation_electron_higgs, unfolded_normalisation_muon_higgs ) # -# filename = path_to_JSON + '/xsection_measurement_results/electron/kv%d/%s/normalisation_%s_Higgs.txt' % ( k_value_electron_central, category, met_type ) -# write_data_to_JSON( unfolded_normalisation_electron_higgs, filename ) -# filename = path_to_JSON + '/xsection_measurement_results/muon/kv%d/%s/normalisation_%s_Higgs.txt' % ( k_value_muon_central, category, met_type ) -# write_data_to_JSON( unfolded_normalisation_muon_higgs, filename ) -# filename = path_to_JSON + '/xsection_measurement_results/combined/%s/normalisation_%s_Higgs.txt' % ( category, met_type ) -# write_data_to_JSON( unfolded_normalisation_combined_higgs, filename ) + # filename = path_to_JSON + '/xsection_measurement_results/electron/kv%d/%s/normalisation_%s_Higgs.txt' % ( k_value_electron_central, category, met_type ) + # write_data_to_JSON( unfolded_normalisation_electron_higgs, filename ) + # filename = path_to_JSON + '/xsection_measurement_results/muon/kv%d/%s/normalisation_%s_Higgs.txt' % ( k_value_muon_central, category, met_type ) + # write_data_to_JSON( unfolded_normalisation_muon_higgs, filename ) + # filename = path_to_JSON + '/xsection_measurement_results/combined/%s/normalisation_%s_Higgs.txt' % ( category, met_type ) + # write_data_to_JSON( unfolded_normalisation_combined_higgs, filename ) # measure xsection calculate_xsections( unfolded_normalisation_electron, category, 'electron', k_value_electron_central ) calculate_xsections( unfolded_normalisation_muon, category, 'muon', k_value_muon_central ) calculate_xsections( unfolded_normalisation_combined, category, 'combined' ) - + calculate_normalised_xsections( unfolded_normalisation_electron, category, 'electron', k_value_electron_central ) calculate_normalised_xsections( unfolded_normalisation_muon, category, 'muon', k_value_muon_central ) - calculate_normalised_xsections( unfolded_normalisation_combined, category, 'combined' ) + calculate_normalised_xsections( unfolded_normalisation_combined, category, 'combined', debug=True ) normalise_to_one = True calculate_normalised_xsections( unfolded_normalisation_electron, category, 'electron', k_value_electron_central, normalise_to_one ) calculate_normalised_xsections( unfolded_normalisation_muon, category, 'muon', k_value_muon_central, normalise_to_one ) calculate_normalised_xsections( unfolded_normalisation_combined, category, 'combined', normalise_to_one ) + + if category == 'central': + normalised_xsection_covariance_electron = calculate_covariance_for_normalised_xsection( covariance_electron, unfolded_normalisation_electron['TTJet_unfolded'], bin_widths[variable] ) + normalised_xsection_covariance_muon = calculate_covariance_for_normalised_xsection( covariance_muon, unfolded_normalisation_muon['TTJet_unfolded'], bin_widths[variable] ) + normalised_xsection_covariance_combined = calculate_covariance_for_normalised_xsection( covariance_combined, unfolded_normalisation_combined['TTJet_unfolded'], bin_widths[variable] ) + + # Write covariance matrices to file + filename = path_to_JSON + '/xsection_measurement_results/electron/kv%d/%s/covariance.txt' % ( k_value_electron_central, category ) + np.savetxt( filename, normalised_xsection_covariance_electron, delimiter = ',' ) + + # # cov_file = File( filename, 'recreate' ) + # # covariance_electron.Write() + # # normalised_xsection_covariance_electron.(filename,',') + # # cov_file.Close() + + filename = path_to_JSON + '/xsection_measurement_results/muon/kv%d/%s/covariance.txt' % ( k_value_muon_central, category ) + np.savetxt( filename, normalised_xsection_covariance_muon, delimiter = ',' ) + + # # cov_file = File( filename, 'recreate' ) + # # covariance_muon.Write() + # # normalised_xsection_covariance_muon.Write() + # # cov_file.Close() + + filename = path_to_JSON + '/xsection_measurement_results/combined/%s/covariance.txt' % ( category ) + np.savetxt( filename, normalised_xsection_covariance_combined, delimiter = ',' ) + for i in range(0,normalised_xsection_covariance_combined.shape[0]): + print np.sqrt( normalised_xsection_covariance_combined[i,i] ) + + # cov_file = File( filename, 'recreate' ) + # covariance_combined.Write() + # normalised_xsection_covariance_combined.Write() + # cov_file.Close() diff --git a/src/cross_section_measurement/03_calculate_systematics.py b/src/cross_section_measurement/03_calculate_systematics.py index 6ca05ca9..792fd5e7 100644 --- a/src/cross_section_measurement/03_calculate_systematics.py +++ b/src/cross_section_measurement/03_calculate_systematics.py @@ -20,16 +20,34 @@ from config import XSectionConfig from tools.file_utilities import read_data_from_JSON, write_data_to_JSON from tools.Calculation import calculate_lower_and_upper_PDFuncertainty, \ -calculate_lower_and_upper_systematics, combine_errors_in_quadrature - +calculate_lower_and_upper_systematics, combine_errors_in_quadrature, \ +calculate_lower_and_upper_systematics_properly, calculate_covariance_of_systematics_03, calculate_covariance_of_systematics_properly +import numpy as np +from decimal import * def read_normalised_xsection_measurement( category, channel ): global path_to_JSON, met_type, met_uncertainties_list, k_values filename = '' - + path_to_use = path_to_JSON + k_value_to_use = k_values[channel] + # Use 8 TeV samples for hadronisation systematic + if 'powheg_v1_pythia' in category or 'powheg_v1_herwig' in category: + if '7TeV' in path_to_JSON: + path_to_use = path_to_JSON.replace('7TeV','8TeV') + config8TeV = XSectionConfig( 8 ) + if channel == 'electron': + k_value_to_use = config8TeV.k_values_electron[variable] + elif channel == 'muon': + k_value_to_use = config8TeV.k_values_muon[variable] + if category in met_uncertainties_list and variable == 'HT': - filename = path_to_JSON + '/' + channel + '/kv' + str( k_values[channel] ) + '/central/normalised_xsection_' + met_type + '.txt' + filename = path_to_use + '/' + channel + '/kv' + str( k_value_to_use ) + '/central/normalised_xsection_' + met_type + '.txt' else: - filename = path_to_JSON + '/' + channel + '/kv' + str( k_values[channel] ) + '/' + category + '/normalised_xsection_' + met_type + '.txt' + filename = path_to_use + '/' + channel + '/kv' + str( k_value_to_use ) + '/' + category + '/normalised_xsection_' + met_type + '.txt' + + # if category in met_uncertainties_list and variable == 'HT': + # filename = path_to_use + '/' + channel + '/kv' + str( k_value_to_use ) + '/central/xsection_' + met_type + '.txt' + # else: + # filename = path_to_use + '/' + channel + '/kv' + str( k_value_to_use ) + '/' + category + '/xsection_' + met_type + '.txt' if channel == 'combined': filename = filename.replace( 'kv' + str( k_values[channel] ), '' ) @@ -66,62 +84,223 @@ def read_normalised_xsection_systematics( list_of_systematics, channel ): return systematics, systematics_unfolded -def summarise_systematics( list_of_central_measurements, dictionary_of_systematics, pdf_calculation = False, hadronisation_systematic = False, mass_systematic = False, kValueSystematic = False ): +def summarise_systematics( list_of_central_measurements, dictionary_of_systematics, pdf_calculation = False, hadronisation_systematic = False, mass_systematic = False, kValueSystematic = False , oneway = False, debug = False ): global symmetrise_errors # number of bins number_of_bins = len( list_of_central_measurements ) down_errors = [0] * number_of_bins up_errors = [0] * number_of_bins + errors_for_covariance = [{}] * number_of_bins + mass_systematic_signs = [] + + # if debug: + # print 'Central:',list_of_central_measurements for bin_i in range( number_of_bins ): central_value = list_of_central_measurements[bin_i][0] # 0 = value, 1 = error error_down, error_up = 0, 0 + e_cov = {} if pdf_calculation: pdf_uncertainty_values = {systematic:measurement[bin_i][0] for systematic, measurement in dictionary_of_systematics.iteritems()} - error_down, error_up = calculate_lower_and_upper_PDFuncertainty( central_value, pdf_uncertainty_values ) + error_down, error_up, dict_of_unc = calculate_lower_and_upper_PDFuncertainty( central_value, pdf_uncertainty_values ) + + # sign = 1 + # if max(error_down, error_up) == error_down: + # sign = -1 + # newsign = np.sign( error_up - error_down ) + + # print sign,newsign,error_up, error_down + if symmetrise_errors: error_down = max( error_down, error_up ) error_up = max( error_down, error_up ) + # e_cov = {'PDF' : max( error_down, error_up ) * sign } + e_cov = dict_of_unc + if debug: + print dict_of_unc + print e_cov + raw_input('...') elif hadronisation_systematic: # always symmetric: absolute value of the difference between powheg_herwig and powheg_pythia - powheg_herwig = dictionary_of_systematics['TTJets_powheg_herwig'][bin_i][0] - powheg_pythia = dictionary_of_systematics['TTJets_powheg_pythia'][bin_i][0] + powheg_herwig = dictionary_of_systematics['TTJets_powheg_v1_herwig'][bin_i][0] + powheg_pythia = dictionary_of_systematics['TTJets_powheg_v1_pythia'][bin_i][0] difference = powheg_herwig - powheg_pythia mean = (powheg_herwig + powheg_pythia)/2.0 + sign = np.sign(difference) + # print central_value,powheg_herwig, powheg_pythia,difference,mean,sign + # print np.sign( central_value - powheg_herwig), np.sign( central_value - powheg_pythia), sign, difference/mean*100 difference = abs(difference) # now scale the error to the central value relative_error = difference/mean error_down = relative_error * central_value error_up = error_down + e_cov = { 'hadronisation' : error_up * sign } elif mass_systematic: - list_of_systematics = [systematic[bin_i][0] for systematic in dictionary_of_systematics.values()] - error_down, error_up = calculate_lower_and_upper_systematics( central_value, list_of_systematics, False ) + list_of_systematics = {systematicName : systematic[bin_i][0] for systematicName, systematic in dictionary_of_systematics.iteritems()} + error_down, error_up, e_cov = calculate_lower_and_upper_systematics( central_value, list_of_systematics, symmetrise_errors=False, mass = True, debug = True ) + if debug: + print 'Before scale :',error_down, error_up, e_cov + + + # Scale errors calculated using very different top masses - error_down, error_up = scaleTopMassSystematicErrors( [error_down], [error_up] ) + error_down, error_up, e_cov = scaleTopMassSystematicErrors( [error_down], [error_up], e_cov ) error_down = error_down[0] error_up = error_up[0] + + sign = 1 + if e_cov['TTJets_massup'] > 0: + if e_cov['TTJets_massup'] > e_cov['TTJets_massdown'] : sign = 1 + else : sign = -1 + elif e_cov['TTJets_massup'] < 0: + if e_cov['TTJets_massup'] > e_cov['TTJets_massdown'] : sign = 1 + else : sign = -1 + + if debug: + print 'After scale :',error_down, error_up, e_cov + + # if debug: + # print error_down, error_up, e_cov + mass_systematic_signs.append(sign) + # print 'MASS' + # print error_down, error_up,e_cov + # e_cov = { 'TTJets_massdown' : error_down } + # e_cov = { 'TTJets_massup' : error_up } + # e_cov = { 'TTJets_massdown' : error_down, 'TTJets_massup' : error_up} elif kValueSystematic: - list_of_systematics = [systematic[bin_i][0] for systematic in dictionary_of_systematics.values()] - error_down, error_up = calculate_lower_and_upper_systematics( central_value, list_of_systematics, True ) + list_of_systematics = {systematicName : systematic[bin_i][0] for systematicName, systematic in dictionary_of_systematics.iteritems()} + error_down, error_up, e_cov = calculate_lower_and_upper_systematics( central_value, list_of_systematics, True ) else: - list_of_systematics = [systematic[bin_i][0] for systematic in dictionary_of_systematics.values()] - error_down, error_up = calculate_lower_and_upper_systematics( central_value, list_of_systematics, symmetrise_errors ) + list_of_systematics = {systematicName : systematic[bin_i][0] for systematicName, systematic in dictionary_of_systematics.iteritems()} + error_down, error_up, e_cov = calculate_lower_and_upper_systematics( central_value, list_of_systematics, symmetrise_errors, debug=debug ) + if debug: + print 'After calculate_lower...' + sys_in_this_bin = {systematicName : systematic[bin_i] for systematicName, systematic in dictionary_of_systematics.iteritems()} + print list_of_central_measurements[bin_i] + print sys_in_this_bin + print error_down, error_up + # print e_cov + # print 'DOING A PROPER JOB' + # dictionary_of_systematics_for_this_bin = { category : systematic[bin_i][0] for category, systematic in dictionary_of_systematics.iteritems()} + # error_down, error_up = calculate_lower_and_upper_systematics_properly( central_value, dictionary_of_systematics_for_this_bin ) down_errors[bin_i] = error_down up_errors[bin_i] = error_up - - return down_errors, up_errors -def scaleTopMassSystematicErrors( error_down, error_up ): + errors_for_covariance[bin_i] = e_cov + if debug: + print 'Calculating covariance' + print 'Errors' + print down_errors + print up_errors + print errors_for_covariance + # print errors_for_covariance + + covariance_matrix = calculate_covariance_of_systematics_03(errors_for_covariance, list_of_central_measurements, mass_systematic, hadronisation_systematic, pdf_calculation, oneway, debug = debug ) + + if mass_systematic: + if debug: + print 'Original' + print np.sign( covariance_matrix ) + + for i_row in range(0, covariance_matrix.shape[0]): + for i_col in range(0, covariance_matrix.shape[1]): + cov = covariance_matrix[i_row,i_col] + sign = mass_systematic_signs[i_row] * mass_systematic_signs[i_col] + if np.sign(cov) != sign: + covariance_matrix[i_row,i_col] *= -1.0 + # print mass_systematic_signs + if debug: + print 'New' + print np.sign( covariance_matrix ) + + for i_row in range(0, covariance_matrix.shape[0]): + print np.sqrt( covariance_matrix[i_row, i_row] ) / list_of_central_measurements[bin_i][0] + + # proper_cov = calculate_covariance_of_systematics_properly( dictionary_of_systematics, list_of_central_measurements ) + + # cov = proper_cov + # if hadronisation_systematic or mass_systematic or pdf_calculation: + # cov = covariance_matrix + + # if debug: + # print 'Covariance' + # print covariance_matrix + # print 'Proper covariance' + # print proper_cov + + # print 'Errors from covariance' + # for i in range(0,number_of_bins): + # print np.sqrt( covariance_matrix[i,i] ) + + return down_errors, up_errors, covariance_matrix + # return down_errors, up_errors, cov + +def scaleTopMassSystematicErrors( error_down, error_up, e_cov=None, debug = False ): error_down_new, error_up_new = [], [] - for down,up in zip( error_down,error_up ): + + # print error_down + # print error_up + for down,up in zip( error_down,error_up ): upMassDifference = measurement_config.topMasses[2] - measurement_config.topMasses[1] downMassDifference = measurement_config.topMasses[1] - measurement_config.topMasses[0] error_down_new.append( down * measurement_config.topMassUncertainty / downMassDifference ) error_up_new.append( up * measurement_config.topMassUncertainty / upMassDifference ) - return error_down_new, error_up_new + + if e_cov != None: + + + print down, up, e_cov, upMassDifference, downMassDifference, measurement_config.topMassUncertainty + + if np.sign( e_cov['TTJets_massdown'] ) == np.sign( e_cov['TTJets_massup'] ): + print 'Signs are the same' + sign = np.sign( e_cov['TTJets_massdown'] ) + down_new = down * measurement_config.topMassUncertainty / downMassDifference + up_new = up * measurement_config.topMassUncertainty / upMassDifference + + new_error = down_new + if down_new == 0: + new_error = up_new + + if abs( e_cov['TTJets_massdown'] ) > abs( e_cov['TTJets_massup'] ): + e_cov['TTJets_massdown'] = new_error * sign + e_cov['TTJets_massup'] = 0 + else : + e_cov['TTJets_massup'] = new_error * sign + e_cov['TTJets_massdown'] = 0 + # e_cov['TTJets_massdown'] = 1.0 * e_cov['TTJets_massdown'] * measurement_config.topMassUncertainty / downMassDifference + # e_cov['TTJets_massup'] = 1.0 * e_cov['TTJets_massup'] * measurement_config.topMassUncertainty / upMassDifference + else: + + if abs( down ) == abs( e_cov['TTJets_massdown'] ): + e_cov['TTJets_massdown'] = down * np.sign( e_cov['TTJets_massdown'] ) * measurement_config.topMassUncertainty / downMassDifference + elif abs( down ) == abs( e_cov['TTJets_massup'] ): + print 'Down is up :',down, measurement_config.topMassUncertainty, upMassDifference + e_cov['TTJets_massup'] = down * np.sign( e_cov['TTJets_massup'] ) * measurement_config.topMassUncertainty / downMassDifference + + if abs( up ) == abs( e_cov['TTJets_massup'] ): + e_cov['TTJets_massup'] = up * np.sign ( e_cov['TTJets_massup'] ) * measurement_config.topMassUncertainty / upMassDifference + elif abs( up ) == abs( e_cov['TTJets_massdown'] ): + e_cov['TTJets_massdown'] = up * np.sign ( e_cov['TTJets_massdown'] ) * measurement_config.topMassUncertainty / upMassDifference + + print e_cov + # print 'Scaled :', down,up, e_cov + + # down_sign = np.sign(e_cov['TTJets_massdown']) + # up_sign = np.sign(e_cov['TTJets_massup']) + + # print e_cov + # print 'Down scale :',measurement_config.topMassUncertainty / downMassDifference, down_sign, down + # print 'Up scale :',measurement_config.topMassUncertainty / upMassDifference, up_sign, up + # e_cov['TTJets_massdown'] = down_sign * down * measurement_config.topMassUncertainty / downMassDifference + # e_cov['TTJets_massup'] = up_sign * up * measurement_config.topMassUncertainty / upMassDifference + # print e_cov + + + if e_cov != None: return error_down_new, error_up_new, e_cov + else: return error_down_new, error_up_new def get_measurement_with_lower_and_upper_errors( list_of_central_measurements, lists_of_lower_systematic_errors, lists_of_upper_systematic_errors ): ''' @@ -134,19 +313,27 @@ def get_measurement_with_lower_and_upper_errors( list_of_central_measurements, l n_entries = len( list_of_central_measurements ) complete_measurement = [( 0, 0, 0 )] * n_entries + print 'In table (?)' + # print lists_of_lower_systematic_errors + # print lists_of_upper_systematic_errors for index in range( n_entries ): central_value, central_error = list_of_central_measurements[index] # 0 = value, 1 = error lower_errors = [error[index] for error in lists_of_lower_systematic_errors] upper_errors = [error[index] for error in lists_of_upper_systematic_errors] - # add central error to the list - lower_errors.append( central_error ) - upper_errors.append( central_error ) + + # # add central error to the list + # lower_errors.append( central_error ) + # upper_errors.append( central_error ) + # calculate total errors total_lower_error = combine_errors_in_quadrature( lower_errors ) total_upper_error = combine_errors_in_quadrature( upper_errors ) if symmetrise_errors: total_lower_error = max( total_lower_error, total_upper_error ) total_upper_error = max( total_lower_error, total_upper_error ) + + print total_lower_error / central_value * 100 + complete_measurement[index] = ( central_value, total_lower_error, total_upper_error ) return complete_measurement @@ -199,7 +386,7 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio b_tag_bin = translate_options[options.bjetbin] path_to_JSON = options.path + '/' + str( options.CoM ) + 'TeV/' + variable + '/xsection_measurement_results/' symmetrise_errors = options.symmetrise_errors - + # set up lists for systematics ttbar_generator_systematics_list = [ttbar_theory_systematic_prefix + systematic for systematic in measurement_config.generator_systematics] vjets_generator_systematics_list = [vjets_theory_systematic_prefix + systematic for systematic in measurement_config.generator_systematics] @@ -212,7 +399,7 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio # ttbar theory systematics: ptreweighting, hadronisation systematic (powheg_pythia - powheg_herwig) ttbar_ptreweight_systematic_list = [ttbar_theory_systematic_prefix + 'ptreweight'] - ttbar_hadronisation_systematic_list = [ttbar_theory_systematic_prefix + 'powheg_pythia', ttbar_theory_systematic_prefix + 'powheg_herwig'] + ttbar_hadronisation_systematic_list = [ttbar_theory_systematic_prefix + 'powheg_v1_pythia', ttbar_theory_systematic_prefix + 'powheg_v1_herwig'] # 45 PDF uncertainties pdf_uncertainties = ['PDFWeights_%d' % index for index in range( 1, 45 )] @@ -232,11 +419,15 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio other_uncertainties_list.append( 'QCD_shape' ) other_uncertainties_list.extend( rate_changing_systematics_list ) - for channel in ['electron', 'muon', 'combined']: -# print "channel = ", channel + # for channel in ['electron', 'muon', 'combined']: + for channel in ['combined']: # read central measurement central_measurement, central_measurement_unfolded = read_normalised_xsection_measurement( 'central', channel ) + nBins = len(central_measurement) + systematic_covariance = np.array( np.zeros((nBins,nBins )) ) + systematic_correlation = np.array( np.zeros((nBins,nBins )) ) + # read groups of systematics ttbar_generator_systematics, ttbar_generator_systematics_unfolded = read_normalised_xsection_systematics( list_of_systematics = ttbar_generator_systematics_list, channel = channel ) ttbar_ptreweight_systematic, ttbar_ptreweight_systematic_unfolded = read_normalised_xsection_systematics( list_of_systematics = ttbar_ptreweight_systematic_list, channel = channel ) @@ -250,32 +441,97 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio other_systematics, other_systematics_unfolded = read_normalised_xsection_systematics( list_of_systematics = other_uncertainties_list, channel = channel ) # get the minimal and maximal deviation for each group of systematics # ttbar generator systematics (factorisation scale and matching threshold) - ttbar_generator_min, ttbar_generator_max = summarise_systematics( central_measurement, ttbar_generator_systematics ) - ttbar_generator_min_unfolded, ttbar_generator_max_unfolded = summarise_systematics( central_measurement_unfolded, ttbar_generator_systematics_unfolded ) + ttbar_generator_min, ttbar_generator_max, covariance = summarise_systematics( central_measurement, ttbar_generator_systematics ) + ttbar_generator_min_unfolded, ttbar_generator_max_unfolded, covariance = summarise_systematics( central_measurement_unfolded, ttbar_generator_systematics_unfolded) + systematic_covariance += covariance + + + + # for i in range(0,nBins): + # for j in range(0,nBins): + # # systematic_covariance[i,j] = Decimal( systematic_covariance[i,j] ).quantize(Decimal('0.00000000000001')) + # systematic_correlation[i,j] = covariance[i,j] / np.sqrt(covariance[i,i] * covariance[j,j]) + # print systematic_correlation + # print covariance + # raw_input('...') + # ttbar theory systematics (pt reweighting and hadronisation) - ttbar_ptreweight_min, ttbar_ptreweight_max = summarise_systematics( central_measurement, ttbar_ptreweight_systematic ) - ttbar_ptreweight_min_unfolded, ttbar_ptreweight_max_unfolded = summarise_systematics( central_measurement_unfolded, ttbar_ptreweight_systematic_unfolded ) - ttbar_hadronisation_min, ttbar_hadronisation_max = summarise_systematics( central_measurement, ttbar_hadronisation_systematic, hadronisation_systematic = True ) - ttbar_hadronisation_min_unfolded, ttbar_hadronisation_max_unfolded = summarise_systematics( central_measurement_unfolded, ttbar_hadronisation_systematic_unfolded, hadronisation_systematic = True ) + ttbar_ptreweight_min, ttbar_ptreweight_max, covariance = summarise_systematics( central_measurement, ttbar_ptreweight_systematic ) + ttbar_ptreweight_min_unfolded, ttbar_ptreweight_max_unfolded, covariance = summarise_systematics( central_measurement_unfolded, ttbar_ptreweight_systematic_unfolded, oneway = True ) + systematic_covariance += covariance + + # for i in range(0,nBins): + # for j in range(0,nBins): + # # systematic_covariance[i,j] = Decimal( systematic_covariance[i,j] ).quantize(Decimal('0.00000000000001')) + # systematic_correlation[i,j] = covariance[i,j] / np.sqrt(covariance[i,i] * covariance[j,j]) + # print systematic_correlation + # print covariance + # raw_input('...') + + ttbar_hadronisation_min, ttbar_hadronisation_max, covariance = summarise_systematics( central_measurement, ttbar_hadronisation_systematic, hadronisation_systematic = True ) + ttbar_hadronisation_min_unfolded, ttbar_hadronisation_max_unfolded, covariance = summarise_systematics( central_measurement_unfolded, ttbar_hadronisation_systematic_unfolded, hadronisation_systematic = True ) + systematic_covariance += covariance + + # Top mass systematic - ttbar_mass_min, ttbar_mass_max = summarise_systematics( central_measurement, ttbar_mass_systematic, mass_systematic = True ) - ttbar_mass_min_unfolded, ttbar_mass_max_unfolded = summarise_systematics( central_measurement_unfolded, ttbar_mass_systematic_unfolded, mass_systematic = True ) + ttbar_mass_min, ttbar_mass_max, covariance = summarise_systematics( central_measurement, ttbar_mass_systematic, mass_systematic = True ) + ttbar_mass_min_unfolded, ttbar_mass_max_unfolded, covariance = summarise_systematics( central_measurement_unfolded, ttbar_mass_systematic_unfolded, mass_systematic = True ) + systematic_covariance += covariance + # k Value systematic - kValue_min, kValue_max = summarise_systematics( central_measurement, kValue_systematic, kValueSystematic = True) - kValue_min_unfolded, kValue_max_unfolded = summarise_systematics( central_measurement_unfolded, kValue_systematic_unfolded, kValueSystematic = True) + kValue_min, kValue_max, covariance = summarise_systematics( central_measurement, kValue_systematic, kValueSystematic = True) + kValue_min_unfolded, kValue_max_unfolded, covariance = summarise_systematics( central_measurement_unfolded, kValue_systematic_unfolded, kValueSystematic = True) # Take up variation as the down variation also. kValue_systematic['kValue_down'] = kValue_systematic['kValue_up'] kValue_systematic_unfolded['kValue_down'] = kValue_systematic_unfolded['kValue_up'] # 45 PDFs - pdf_min, pdf_max = summarise_systematics( central_measurement, pdf_systematics, pdf_calculation = True ) - pdf_min_unfolded, pdf_max_unfolded = summarise_systematics( central_measurement_unfolded, pdf_systematics_unfolded, pdf_calculation = True ) + pdf_min, pdf_max, covariance = summarise_systematics( central_measurement, pdf_systematics, pdf_calculation = True ) + pdf_min_unfolded, pdf_max_unfolded, covariance = summarise_systematics( central_measurement_unfolded, pdf_systematics_unfolded, pdf_calculation = True ) + systematic_covariance += covariance + + # MET - met_min, met_max = summarise_systematics( central_measurement, met_systematics ) - met_min_unfolded, met_max_unfolded = summarise_systematics( central_measurement_unfolded, met_systematics_unfolded ) + met_min, met_max, covariance = summarise_systematics( central_measurement, met_systematics ) + met_min_unfolded, met_max_unfolded, covariance = summarise_systematics( central_measurement_unfolded, met_systematics_unfolded ) + if variable != 'HT': + systematic_covariance += covariance + + + # other - other_min, other_max = summarise_systematics( central_measurement, other_systematics ) - other_min_unfolded, other_max_unfolded = summarise_systematics( central_measurement_unfolded, other_systematics_unfolded ) + other_min, other_max, covariance = summarise_systematics( central_measurement, other_systematics ) + print 'OTHER SYSTEMATICS' + other_min_unfolded, other_max_unfolded, covariance = summarise_systematics( central_measurement_unfolded, other_systematics_unfolded ) + systematic_covariance += covariance + + + + # for i in range(0,nBins): + # for j in range(0,nBins): + # # systematic_covariance[i,j] = Decimal( systematic_covariance[i,j] ).quantize(Decimal('0.00000000000001')) + # systematic_correlation[i,j] = covariance[i,j] / np.sqrt(covariance[i,i] * covariance[j,j]) + # print systematic_correlation + # print covariance + + # raw_input('...') + + print 'DONE' + filename = path_to_JSON + '/combined/central/covariance_systematic.txt' + np.savetxt( filename, systematic_covariance, delimiter = ',' ) + print channel + # # print systematic_covariance + for i in range(0,nBins): + print np.sqrt(systematic_covariance[i,i]) / central_measurement_unfolded[i][0] * 100 + for i in range(0,nBins): + for j in range(0,nBins): + # systematic_covariance[i,j] = Decimal( systematic_covariance[i,j] ).quantize(Decimal('0.00000000000001')) + systematic_correlation[i,j] = systematic_covariance[i,j] / np.sqrt(systematic_covariance[i,i] * systematic_covariance[j,j]) + + # print 'covariance' + # print systematic_covariance + # print 'Correlation' + # print systematic_correlation # get the central measurement with fit, unfolding and systematic errors combined central_measurement_with_systematics = get_measurement_with_lower_and_upper_errors( central_measurement, [ttbar_generator_min, ttbar_ptreweight_min, @@ -308,16 +564,26 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio pdf_max, met_max, other_max] ) central_measurement_unfolded_with_systematics = get_measurement_with_lower_and_upper_errors( central_measurement_unfolded, - [ttbar_generator_min_unfolded, ttbar_ptreweight_min_unfolded, + [ + ttbar_generator_min_unfolded, + ttbar_ptreweight_min_unfolded, ttbar_hadronisation_min_unfolded, ttbar_mass_min_unfolded, -# kValue_min_unfolded, - pdf_min_unfolded, met_min_unfolded, other_min_unfolded], - [ttbar_generator_max_unfolded, ttbar_ptreweight_max_unfolded, +# # # # kValue_min_unfolded, + pdf_min_unfolded, + met_min_unfolded, + other_min_unfolded + ], + [ + ttbar_generator_max_unfolded, + ttbar_ptreweight_max_unfolded, ttbar_hadronisation_max_unfolded, ttbar_mass_max_unfolded, -# kValue_max_unfolded, - pdf_max_unfolded, met_max_unfolded, other_max_unfolded] ) +# # # # kValue_max_unfolded, + pdf_max_unfolded, + met_max_unfolded, + other_max_unfolded + ] ) central_measurement_unfolded_with_systematics_but_without_ttbar_theory = get_measurement_with_lower_and_upper_errors( central_measurement_unfolded, [pdf_min_unfolded, met_min_unfolded, other_min_unfolded, ttbar_mass_min_unfolded, diff --git a/src/cross_section_measurement/04_make_plots_matplotlib.py b/src/cross_section_measurement/04_make_plots_matplotlib.py index ea7cd8a4..b5f84991 100644 --- a/src/cross_section_measurement/04_make_plots_matplotlib.py +++ b/src/cross_section_measurement/04_make_plots_matplotlib.py @@ -10,9 +10,10 @@ from tools.file_utilities import read_data_from_JSON, make_folder_if_not_exists from tools.hist_utilities import value_error_tuplelist_to_hist, spread_x,\ value_tuplelist_to_hist, value_errors_tuplelist_to_graph, graph_to_value_errors_tuplelist +from tools.plotting import adjust_ratio_ticks from math import sqrt # rootpy & matplotlib -from ROOT import kRed, kGreen, kMagenta, kBlue, kBlack +from ROOT import kRed, kGreen, kMagenta, kBlue, kBlack, kCyan, kOrange from tools.ROOT_utils import set_root_defaults import matplotlib as mpl from tools.plotting import get_best_max_y @@ -64,11 +65,16 @@ def read_xsection_measurement_results( category, channel ): # true distributions h_normalised_xsection_MADGRAPH = value_error_tuplelist_to_hist( normalised_xsection_unfolded['MADGRAPH'], bin_edges[variable] ) h_normalised_xsection_MADGRAPH_ptreweight = value_error_tuplelist_to_hist( normalised_xsection_unfolded['MADGRAPH_ptreweight'], bin_edges[variable] ) - h_normalised_xsection_POWHEG_PYTHIA = value_error_tuplelist_to_hist( normalised_xsection_unfolded['POWHEG_PYTHIA'], bin_edges[variable] ) - h_normalised_xsection_POWHEG_HERWIG = value_error_tuplelist_to_hist( normalised_xsection_unfolded['POWHEG_HERWIG'], bin_edges[variable] ) + h_normalised_xsection_POWHEG_V2_PYTHIA = value_error_tuplelist_to_hist( normalised_xsection_unfolded['powheg_v2_pythia'], bin_edges[variable] ) + h_normalised_xsection_POWHEG_V2_HERWIG = value_error_tuplelist_to_hist( normalised_xsection_unfolded['powheg_v2_herwig'], bin_edges[variable] ) h_normalised_xsection_MCATNLO = None + h_normalised_xsection_POWHEG_V1_HERWIG = None + h_normalised_xsection_POWHEG_V1_PYTHIA = None if measurement_config.centre_of_mass_energy == 8: h_normalised_xsection_MCATNLO = value_error_tuplelist_to_hist( normalised_xsection_unfolded['MCATNLO'], bin_edges[variable] ) + h_normalised_xsection_POWHEG_V1_HERWIG = value_error_tuplelist_to_hist( normalised_xsection_unfolded['powheg_v1_herwig'], bin_edges[variable] ) + h_normalised_xsection_POWHEG_V1_PYTHIA = value_error_tuplelist_to_hist( normalised_xsection_unfolded['powheg_v1_pythia'], bin_edges[variable] ) + h_normalised_xsection_mathchingup = value_error_tuplelist_to_hist( normalised_xsection_unfolded['matchingup'], bin_edges[variable] ) h_normalised_xsection_mathchingdown = value_error_tuplelist_to_hist( normalised_xsection_unfolded['matchingdown'], bin_edges[variable] ) h_normalised_xsection_scaleup = value_error_tuplelist_to_hist( normalised_xsection_unfolded['scaleup'], bin_edges[variable] ) @@ -76,10 +82,14 @@ def read_xsection_measurement_results( category, channel ): histograms_normalised_xsection_different_generators.update( {'MADGRAPH':h_normalised_xsection_MADGRAPH, 'MADGRAPH_ptreweight':h_normalised_xsection_MADGRAPH_ptreweight, - 'POWHEG_PYTHIA':h_normalised_xsection_POWHEG_PYTHIA, - 'POWHEG_HERWIG':h_normalised_xsection_POWHEG_HERWIG} ) + 'powheg_v2_pythia':h_normalised_xsection_POWHEG_V2_PYTHIA, + 'powheg_v2_herwig':h_normalised_xsection_POWHEG_V2_HERWIG, + } ) if measurement_config.centre_of_mass_energy == 8: histograms_normalised_xsection_different_generators.update( {'MCATNLO':h_normalised_xsection_MCATNLO} ) + histograms_normalised_xsection_different_generators.update( {'powheg_v1_herwig':h_normalised_xsection_POWHEG_V1_HERWIG} ) + histograms_normalised_xsection_different_generators.update( {'powheg_v1_pythia':h_normalised_xsection_POWHEG_V1_PYTHIA} ) + histograms_normalised_xsection_systematics_shifts.update( {'MADGRAPH':h_normalised_xsection_MADGRAPH, 'MADGRAPH_ptreweight':h_normalised_xsection_MADGRAPH_ptreweight, @@ -299,8 +309,11 @@ def get_cms_labels( channel ): lepton = 'e + jets' elif channel == 'muon': lepton = '$\mu$ + jets' + lepton = r'\boldmath$\mu$ + jets' else: - lepton = 'e, $\mu$ + jets combined' + # lepton = 'e, $\mu$ + jets combined' + lepton = r'e, \boldmath$\mu$ + jets combined' + # channel_label = '%s, $\geq$ 4 jets, %s' % ( lepton, b_tag_bins_latex[b_tag_bin] ) channel_label = lepton template = '%.1f fb$^{-1}$ (%d TeV)' @@ -364,8 +377,8 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True plt.xlabel( '$%s$ [GeV]' % variables_latex[variable], CMS.x_axis_title ) axes.minorticks_on() - - plt.ylabel( r'$\frac{1}{\sigma} \frac{d\sigma}{d' + variables_latex[variable] + '} \left[\mathrm{GeV}^{-1}\\right]$', CMS.y_axis_title ) + + plt.ylabel( r'$\displaystyle\frac{1}{\sigma} \frac{d\sigma}{d' + variables_latex[variable] + '} (\mathrm{GeV}^{-1})$', CMS.y_axis_title ) plt.tick_params( **CMS.axis_label_major ) plt.tick_params( **CMS.axis_label_minor ) @@ -388,7 +401,7 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True # now draw data points with only statistical error bars with caps if not draw_background_subtraction: - rplt.errorbar( hist_data, axes = axes, label = 'do_not_show', xerr = None, capsize = 15, capthick = 3, elinewidth = 2, zorder = len( histograms ) + 2 ) + rplt.errorbar( hist_data, axes = axes, label = 'do_not_show', xerr = None, capsize = 8, capthick = 3, elinewidth = 2, zorder = len( histograms ) + 2 ) elif draw_background_subtraction: data_histograms_only_stat_error_to_plot = [hist_data, hist_data_background_subtraction] graphs = spread_x(data_histograms_only_stat_error_to_plot, bin_edges[variable]) @@ -411,37 +424,56 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True if show_before_unfolding: rplt.errorbar( hist_measured, axes = axes, label = 'data (before unfolding)', xerr = None, zorder = len( histograms ) ) + dashes = {} for key, hist in sorted( histograms.iteritems() ): + zorder = sorted( histograms, reverse = False ).index( key ) + if key == 'MADGRAPH' and zorder != len(histograms) - 3: + zorder = len(histograms) - 3 + elif key != 'MADGRAPH' and not 'unfolded' in key: + while zorder >= len(histograms) - 3: + zorder = zorder - 1 + if not 'unfolded' in key and not 'measured' in key: hist.linewidth = 4 + linestyle = None # setting colours - if 'POWHEG_PYTHIA' in key or 'matchingdown' in key: - hist.linestyle = 'longdashdot' + if 'powheg_v2_pythia' in key or 'matchingdown' in key: + dashes[key] = [25,5,5,5,5,5,5,5] + hist.SetLineColor( kCyan + 2 ) + elif 'powheg_v1_pythia' in key: + dashes[key] = [10,5,5,10,5,5] hist.SetLineColor( kBlue ) - elif 'POWHEG_HERWIG' in key or 'scaledown' in key: - hist.linestyle = 'dashed' + elif key == 'powheg_v1_herwig' or 'scaledown' in key: + dashes[key] = [5,5] hist.SetLineColor( kGreen ) elif 'MADGRAPH_ptreweight' in key: - hist.linestyle = 'dashed' - hist.SetLineColor( kBlack ) + dashes[key] = [20,5] + hist.SetLineColor( kRed + 3 ) elif 'MADGRAPH' in key: - hist.linestyle = 'solid' + linestyle = 'solid' + dashes[key] = None hist.SetLineColor( kRed + 1 ) - elif 'matchingup' in key: - hist.linestyle = 'verylongdashdot' - hist.linecolor = 'orange' + elif 'matchingup' in key or 'powheg_v2_herwig' in key: + dashes[key] = [25,5,5,5,25,5,5,5] + hist.SetLineColor( kMagenta + 1 ) elif 'MCATNLO' in key or 'scaleup' in key: - hist.linestyle = 'dotted' - hist.SetLineColor( kMagenta + 3 ) - rplt.hist( hist, axes = axes, label = measurements_latex[key], zorder = sorted( histograms, reverse = True ).index( key ) ) - + dashes[key] = [5,5,10,5] + hist.SetLineColor( kOrange + 5 ) + # hist.linecolor = 'orange' + + if linestyle != None: + hist.linestyle = linestyle + line, h = rplt.hist( hist, axes = axes, label = measurements_latex[key], zorder = zorder ) + if dashes[key] != None: + line.set_dashes(dashes[key]) + h.set_dashes(dashes[key]) handles, labels = axes.get_legend_handles_labels() # making data first in the list data_label_index = labels.index( 'data' ) data_handle = handles[data_label_index] labels.remove( 'data' ) handles.remove( data_handle ) - labels.insert( 0, 'unfolded data' ) + labels.insert( 0, 'Unfolded data' ) handles.insert( 0, data_handle ) new_handles, new_labels = [], [] @@ -450,10 +482,12 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True new_handles.append( handle ) new_labels.append( label ) - legend_location = (0.98, 0.88) + legend_location = (0.95, 0.88) if variable == 'MT': legend_location = (0.9, 0.6) # 0.98, 0.88 CMS.legend_properties['size'] = 27 + if variable == 'WPT' and measurement_config.centre_of_mass_energy == 8: + legend_location = (0.97, 0.9) # 0.98, 0.88 # elif variable == 'ST': # legend_location = (0.95, 0.88) plt.legend( new_handles, new_labels, numpoints = 1, prop = CMS.legend_properties, frameon = False, bbox_to_anchor=legend_location, @@ -463,10 +497,10 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True plt.title( label,loc='right', **CMS.title ) # CMS text # note: fontweight/weight does not change anything as we use Latex text!!! - plt.text(0.05, 0.98, r"\textbf{CMS}", transform=axes.transAxes, fontsize=42, + plt.text(0.05, 0.965, r"\textbf{CMS}", transform=axes.transAxes, fontsize=42, verticalalignment='top',horizontalalignment='left') # channel text plt.text(0.95, 0.98 - axes.text(0.95, 0.98, r"\emph{%s}" %channel_label, transform=axes.transAxes, fontsize=42, + axes.text(0.95, 0.965, r"%s" %channel_label, transform=axes.transAxes, fontsize=42, verticalalignment='top',horizontalalignment='right') if show_ratio: @@ -480,10 +514,11 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) ) ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) - plt.xlabel( '$%s$ [GeV]' % variables_latex[variable], CMS.x_axis_title ) + plt.xlabel( '$%s$ (GeV)' % variables_latex[variable], CMS.x_axis_title ) + ax1.xaxis.labelpad = 20 plt.tick_params( **CMS.axis_label_major ) plt.tick_params( **CMS.axis_label_minor ) - plt.ylabel( '$\\frac{\\textrm{pred.}}{\\textrm{data}}$', CMS.y_axis_title ) + plt.ylabel( '$\\displaystyle\\frac{\\mathrm{pred.}}{\\mathrm{data}}$', CMS.y_axis_title ) ax1.yaxis.set_label_coords(-0.115, 0.8) #draw a horizontal line at y=1 for data plt.axhline(y = 1, color = 'black', linewidth = 2) @@ -492,7 +527,10 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True if not 'unfolded' in key and not 'measured' in key: ratio = hist.Clone() ratio.Divide( hist_data ) #divide by data - rplt.hist( ratio, axes = ax1, label = 'do_not_show' ) + line, h = rplt.hist( ratio, axes = ax1, label = 'do_not_show' ) + if dashes[key] != None: + line.set_dashes(dashes[key]) + h.set_dashes(dashes[key]) stat_lower = hist_data.Clone() stat_upper = hist_data.Clone() @@ -510,36 +548,59 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True syst_upper.SetBinContent( bin_i, 1 + syst_errors[bin_i-1][2]/syst_errors[bin_i-1][0] ) if category == 'central': rplt.fill_between( syst_lower, syst_upper, ax1, - color = 'yellow', alpha = 0.5 ) + color = 'gold' ) - rplt.fill_between( stat_upper, stat_lower, ax1, color = '0.75', - alpha = 0.5 ) + rplt.fill_between( stat_upper, stat_lower, ax1, color = '0.6' ) # legend for ratio plot - p_stat = mpatches.Patch(facecolor='0.75', label='Stat.',alpha = 0.5, edgecolor='black' ) - p_stat_and_syst = mpatches.Patch(facecolor='yellow', label=r'Stat. $\oplus$ Syst.', alpha = 0.5, edgecolor='black' ) - l1 = ax1.legend(handles = [p_stat], loc = 'upper left', - frameon = False, prop = {'size':26}) - - ax1.legend(handles = [p_stat_and_syst], loc = 'lower left', - frameon = False, prop = {'size':26}) - ax1.add_artist(l1) - + p_stat = mpatches.Patch(facecolor='0.6', label='Stat.', edgecolor='black' ) + p_stat_and_syst = mpatches.Patch(facecolor='gold', label=r'Stat. $\oplus$ Syst.', edgecolor='black' ) + + ratioLegendLocation = 'best' if variable == 'MET': - ax1.set_ylim( ymin = 0.7, ymax = 1.3 ) - ax1.yaxis.set_major_locator( MultipleLocator( 0.2 ) ) -# ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + ax1.set_ylim( ymin = 0.6, ymax = 1.4 ) + # ax1.yaxis.set_major_locator( MultipleLocator( 0.3 ) ) + # ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + if measurement_config.centre_of_mass_energy == 7: + ratioLegendLocation = 'upper left' + elif measurement_config.centre_of_mass_energy == 8: + ratioLegendLocation = 'upper left' + if variable == 'MT': ax1.set_ylim( ymin = 0.8, ymax = 1.2 ) - ax1.yaxis.set_major_locator( MultipleLocator( 0.2 ) ) - ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) - elif variable == 'HT' or variable == 'ST': - ax1.set_ylim( ymin = 0.5, ymax = 1.5 ) - ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) ) - ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + # ax1.yaxis.set_major_locator( MultipleLocator( 0.2 ) ) + # ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + elif variable == 'ST': + ax1.set_ylim( ymin = 0.6, ymax = 1.4 ) + # ax1.yaxis.set_major_locator( MultipleLocator( 0.3 ) ) + # ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + if measurement_config.centre_of_mass_energy == 7: + ratioLegendLocation = 'lower center' + elif measurement_config.centre_of_mass_energy == 8: + ratioLegendLocation = 'lower center' + elif variable == 'HT': + ax1.set_ylim( ymin = 0.6, ymax = 1.4 ) + # ax1.yaxis.set_major_locator( MultipleLocator( 0.3 ) ) + # ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + if measurement_config.centre_of_mass_energy == 7: + ratioLegendLocation = 'lower center' + elif measurement_config.centre_of_mass_energy == 8: + ratioLegendLocation = 'lower center' elif variable == 'WPT': - ax1.set_ylim( ymin = 0.75, ymax = 1.5 ) - ax1.yaxis.set_major_locator( MultipleLocator( 0.25 ) ) - ax1.yaxis.set_minor_locator( MultipleLocator( 0.05 ) ) + ax1.set_ylim( ymin = 0.6, ymax = 1.4 ) + # ax1.yaxis.set_major_locator( MultipleLocator( 0.2 ) ) + # ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + if measurement_config.centre_of_mass_energy == 7: + ratioLegendLocation = 'upper left' + elif measurement_config.centre_of_mass_energy == 8: + ratioLegendLocation = 'upper left' + adjust_ratio_ticks(ax1.yaxis, n_ticks = 3, y_limits = ax1.get_ylim() ) + + l1 = ax1.legend(handles = [p_stat, p_stat_and_syst], loc = ratioLegendLocation, + frameon = False, prop = {'size':26}, ncol = 2) + ax1.add_artist(l1) + + if variable == 'WPT' and measurement_config.centre_of_mass_energy == 8: + axes.set_ylim( axes.get_ylim()[0], axes.get_ylim()[1] * 1.12 ) if CMS.tight_layout: plt.tight_layout() @@ -709,7 +770,8 @@ def get_unit_string(fit_variable): all_measurements.extend( met_uncertainties ) all_measurements.extend( new_uncertainties ) all_measurements.extend( rate_changing_systematics ) - for channel in ['electron', 'muon', 'combined']: + # for channel in ['electron', 'muon', 'combined']: + for channel in ['combined']: for category in all_measurements: if not category == 'central' and not options.draw_systematics: continue @@ -726,14 +788,14 @@ def get_unit_string(fit_variable): if met_type == 'PFMETJetEnDown': met_type = 'patPFMetJetEnDown' - if not channel == 'combined' and options.additional_plots: - #Don't make additional plots for e.g. generator systematics, mass systematics, k value systematics and pdf systematics because they are now done \ - #in the unfolding process with BLT unfolding files. - if category in ttbar_generator_systematics or category in ttbar_mass_systematics or category in kValue_systematics or category in pdf_uncertainties: - continue - fit_templates, fit_results = read_fit_templates_and_results_as_histograms( category, channel ) - make_template_plots( fit_templates, category, channel ) - plot_fit_results( fit_results, category, channel ) + # if not channel == 'combined' and options.additional_plots: + # #Don't make additional plots for e.g. generator systematics, mass systematics, k value systematics and pdf systematics because they are now done \ + # #in the unfolding process with BLT unfolding files. + # if category in ttbar_generator_systematics or category in ttbar_mass_systematics or category in kValue_systematics or category in pdf_uncertainties: + # continue + # fit_templates, fit_results = read_fit_templates_and_results_as_histograms( category, channel ) + # make_template_plots( fit_templates, category, channel ) + # plot_fit_results( fit_results, category, channel ) # change back to original MET type met_type = translate_options[options.metType] @@ -741,7 +803,6 @@ def get_unit_string(fit_variable): met_type = 'patMETsPFlow' histograms_normalised_xsection_different_generators, histograms_normalised_xsection_systematics_shifts = read_xsection_measurement_results( category, channel ) - make_plots( histograms_normalised_xsection_different_generators, category, output_folder, 'normalised_xsection_' + channel + '_different_generators' ) make_plots( histograms_normalised_xsection_systematics_shifts, category, output_folder, 'normalised_xsection_' + channel + '_systematics_shifts' ) diff --git a/src/cross_section_measurement/05_make_tables.py b/src/cross_section_measurement/05_make_tables.py index d96b8bbb..ad8093a8 100644 --- a/src/cross_section_measurement/05_make_tables.py +++ b/src/cross_section_measurement/05_make_tables.py @@ -1,15 +1,28 @@ from __future__ import division # the result of the division will be always a float from optparse import OptionParser from copy import deepcopy -from config.latex_labels import variables_latex, measurements_latex, met_systematics_latex, samples_latex, typical_systematics_latex +from config.latex_labels import variables_latex, measurements_latex, met_systematics_latex, samples_latex, typical_systematics_latex, variables_latex_macros, variables_hepdata from config.variable_binning import variable_bins_latex, variable_bins_ROOT from config import XSectionConfig -from tools.Calculation import getRelativeError +from tools.Calculation import getRelativeError, calculate_covariance_of_systematics, get_correlation_matrix from tools.file_utilities import read_data_from_JSON, make_folder_if_not_exists from lib import read_normalisation, read_initial_normalisation import math import os.path -from numpy import median +import numpy as np +from ROOT import TMath +from decimal import * + +def fix_trailing_zeroes( number ): + nMissingZeros = 2 - len( number.replace('.','') ) + if not '.' in number and len(number) < 2: number += '.' + if number[0] == '0': + while len( number.split('.')[-1] ) < 2 : number += '0' + else : + while len( number.replace('.','') ) < 2 : number += '0' + + + return number def read_xsection_measurement_results_with_errors(channel): global path_to_JSON, variable, k_values, met_type @@ -24,7 +37,7 @@ def read_xsection_measurement_results_with_errors(channel): normalised_xsection_measured_unfolded = {'measured':normalised_xsection_unfolded['TTJet_measured'], 'unfolded':normalised_xsection_unfolded['TTJet_unfolded']} - + file_name = file_template.replace('.txt', '_with_errors.txt') normalised_xsection_unfolded_with_errors = read_data_from_JSON( file_name ) @@ -50,7 +63,20 @@ def read_xsection_measurement_results_with_errors(channel): normalised_xsection_new_errors = read_data_from_JSON( file_name ) normalised_xsection_measured_unfolded.update({'measured_with_systematics':normalised_xsection_unfolded_with_errors['TTJet_measured'], - 'unfolded_with_systematics':normalised_xsection_unfolded_with_errors['TTJet_unfolded']}) + 'unfolded_with_systematics':normalised_xsection_unfolded_with_errors['TTJet_unfolded']},) + + normalised_xsection_measured_unfolded.update({'madgraph':normalised_xsection_unfolded['MADGRAPH'], + 'madgraph_ptreweight':normalised_xsection_unfolded['MADGRAPH_ptreweight'], + 'powheg_v2_pythia':normalised_xsection_unfolded['powheg_v2_pythia'], + 'powheg_v2_herwig':normalised_xsection_unfolded['powheg_v2_herwig'], + }) + if measurement_config.centre_of_mass_energy == 8: + normalised_xsection_measured_unfolded.update({ + 'mcatnlo':normalised_xsection_unfolded['MCATNLO'], + 'powheg_v1_pythia':normalised_xsection_unfolded['powheg_v1_pythia'], + 'powheg_v1_herwig':normalised_xsection_unfolded['powheg_v1_herwig'], + }, + ) normalised_xsection_measured_errors = normalised_xsection_ttbar_generator_errors['TTJet_measured'] normalised_xsection_measured_errors.update(normalised_xsection_PDF_errors['TTJet_measured']) @@ -244,26 +270,37 @@ def print_xsections(xsections, channel, toFile = True, print_before_unfolding = printout += '\n' printout += '\\begin{table}[htbp]\n' - printout += '\\setlength{\\tabcolsep}{2pt}\n' + printout += '\\setlength{\\tabcolsep}{8pt}\n' printout += '\\centering\n' - printout += '\\caption{Normalised \\ttbar cross section measurement with respect to \\%s variable\n' % variable - printout += 'at a centre-of-mass energy of %d TeV ' % measurement_config.centre_of_mass_energy + printout += '\\caption{Normalized \\ttbar differential cross section measurements with respect to the %s variable\n' % variables_latex_macros[variable] + printout += 'at a center-of-mass energy of %d TeV ' % measurement_config.centre_of_mass_energy if channel == 'combined': printout += '(combination of electron and muon channels).' else: printout += '(%s channel).' % channel - printout += ' The errors shown are combined statistical, fit and unfolding errors ($^\dagger$) and systematic uncertainty ($^\star$).}\n' + + printout += ' The rightmost three columns show the relative uncertainties on the measured values, in percent. The statistical and systematic uncertainties are listed separately, and are combined in quadrature to give the overall relative uncertainty.}\n' + + # printout += ' The first uncertainty is the statistical uncertainty, and the second is the systematic uncertainty. The overall relative uncertainty is also given in percent.}\n' + # printout += ' The first uncertainty is the statistical uncertainty, and the second is the systematic uncertainty. The overall uncertainty is also given. All uncertainties are relative uncertainties in percent.}\n' printout += '\\label{tab:%s_xsections_%dTeV_%s}\n' % (variable, measurement_config.centre_of_mass_energy, channel) #printout += '\\resizebox{\\columnwidth}{!} {\n' - printout += '\\begin{tabular}{lrrrr}\n' + printout += '\\begin{tabular}{ccccc}\n' printout += '\\hline\n' - printout += '$%s$ bin [\\GeV] & \\multicolumn{4}{c}{$\sigma_{meas} \\left(\\times 10^{3}\\right)$}' % variables_latex[variable] + # printout += '%s & $\sigma_{\mathrm{meas}}^{\mathrm{norm}}$ & $\pm \\textrm{ stat.}$ & $\pm \\textrm{ syst.}$ & Relative \\\\ \n' % variables_latex_macros[variable] + printout += '%s & $1/\sigma\ \mathrm{d}\sigma/\mathrm{d}%s$ & $\pm \\textrm{ stat.}$ & $\pm \\textrm{ syst.}$ & Relative \\\\ \n' % ( variables_latex_macros[variable], variables_latex_macros[variable]) + +# r'$\displaystyle\frac{1}{\sigma} \frac{d\sigma}{d' + variables_latex[variable] + '} \left[\mathrm{GeV}^{-1}\\right]$', CMS.y_axis_title + + # printout += '(\GeV) & \multicolumn{3}{c}{($\\times 10^3 \\,\\mathrm{GeV}^{-1}$)} & uncertainty (\%)' + printout += '(\GeV) & ($\\mathrm{GeV}^{-1}$) & (\%) & (\%) & uncertainty (\%)' printout += '\\\\ \n\hline\n' - scale = 1000 - + bins = variable_bins_ROOT[variable] assert(len(bins) == len(xsections['unfolded_with_systematics'])) - + # print 'Output table' + outputRelErrors = [] + # print 'Table' for bin_i, variable_bin in enumerate(bins): if print_before_unfolding: value, stat_error = xsections['measured'][bin_i] @@ -271,20 +308,52 @@ def print_xsections(xsections, channel, toFile = True, print_before_unfolding = else: value, stat_error = xsections['unfolded'][bin_i] _, total_error_up, total_error_down = xsections['unfolded_with_systematics'][bin_i] + + + # extracting the systematic error from the total in quadrature - syst_error_up = math.sqrt(total_error_up**2 - stat_error**2) - syst_error_down = math.sqrt(total_error_down**2 - stat_error**2) + # syst_error_up = math.sqrt(total_error_up**2 - stat_error**2) + # syst_error_down = math.sqrt(total_error_down**2 - stat_error**2) + syst_error_up = getRelativeError( value, math.sqrt(total_error_up**2 - stat_error**2) ) + syst_error_down = getRelativeError( value, math.sqrt(total_error_down**2 - stat_error**2) ) + + stat_error = getRelativeError( value, stat_error ) + #relative errors for percentages total_relativeError_up = getRelativeError(value, total_error_up) total_relativeError_down = getRelativeError(value, total_error_down) + + scale = 10 + exponent = 1 + value *= scale + while value < 1: + value *= 10 + exponent += 1 + if total_error_up == total_error_down: - printout += '%s & ' % variable_bins_latex[variable_bin] + ' $%.2f$ & $ \pm~ %.2f^\\dagger$ & $ \pm~ %.2f^\\star$ & ' % (value * scale, stat_error * scale, syst_error_up * scale) +\ - '$(%.2f' % (total_relativeError_up * 100) + '\%)$' + + relErrorToPrint = fix_trailing_zeroes( '%.2g' % ( total_relativeError_up * 100 ) ) + stat_error_to_print = fix_trailing_zeroes( '%.2g' % ( stat_error * 100 ) ) + syst_error_to_print = fix_trailing_zeroes( '%.2g' % ( syst_error_up * 100 ) ) + + + + # relErrorToPrint = '%.2g' % (total_relativeError_up * 100) + + # printout += '%s & ' % variable_bins_latex[variable_bin] + ' $%.2f$ & $%.2f$ & $%.2f$ & ' % (value * scale, stat_error * scale, syst_error_up * scale) +\ + # '$%s' % (relErrorToPrint) + '$' + # # '$%.f' % (total_relativeError_up) * 100) + '$' + + printout += '%s & ' % variable_bins_latex[variable_bin] + ' $%.2f \\times 10^{-%i}$ & $%s$ & $%s$ & ' % (value, exponent, stat_error_to_print, syst_error_to_print) +\ + '$%s' % (relErrorToPrint) + '$' + # '$%.f' % (total_relativeError_up) * 100) + '$' + # print stat_error_to_print, syst_error_to_print,relErrorToPrint + outputRelErrors.append( relErrorToPrint ) else: printout += '%s & ' % variable_bins_latex[variable_bin] + ' $%.2f$ & $ \pm~ %.2f^\\dagger$ & $ ~^{+%.2f}_{-%.2f}^\\star$ & ' % (value * scale, stat_error * scale, syst_error_up * scale, syst_error_down * scale) +\ '$(^{+%.2f}_{-%.2f}' % (total_relativeError_up * 100, total_relativeError_down * 100) + '\%)$' printout += '\\\\ \n' - + print 'Done' printout += '\\hline \n' printout += '\\end{tabular}\n' #printout += '}\n' #for resizebox @@ -304,6 +373,8 @@ def print_xsections(xsections, channel, toFile = True, print_before_unfolding = else: print printout + return outputRelErrors + def print_error_table(central_values, errors, channel, toFile = True, print_before_unfolding = False): global output_folder, variable, k_values, met_type, b_tag_bin, all_measurements bins = variable_bins_ROOT[variable] @@ -418,13 +489,16 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr else: assert(len(bins) == len(central_values['unfolded'])) + if variable == "HT": + typical_systematics_order.remove("typical_systematics_MET") + for bin_i, variable_bin in enumerate(bins): if print_before_unfolding: central_value = central_values['measured'][bin_i][0] else: central_value = central_values['unfolded'][bin_i][0] - for systematic_group in typical_systematics.keys(): + for systematic_group in typical_systematics_order: for source in all_measurements: if source in typical_systematics[systematic_group]: if met_type in source: @@ -441,11 +515,10 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr rows_for_typical_systematics_table = {} - if variable == "HT": - del(typical_systematics["typical_systematics_MET"]) - for systematic_group in typical_systematics.keys(): + allErrors = [] + for systematic_group in typical_systematics_order: typical_systematics_row = [] - typical_systematics_row.append(typical_systematics_latex[systematic_group] + ' (\%)') + typical_systematics_row.append(typical_systematics_latex[systematic_group] ) for bin_i, variable_bin in enumerate(bins): sum = 0. @@ -475,11 +548,20 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr typical_systematics_row.append(sum) label = typical_systematics_row.pop(0) - #for each systematic group, take the median of the values across all bins - value = median(typical_systematics_row) - text = '%.2f' % (value*100) + #for each systematic group, take the median of the values across all bins + allErrors.append( typical_systematics_row ) + value = np.median(typical_systematics_row) + text = None + if value * 100 < 0.1 : + text = '$<0.1$' + else : + text = '%.1f' % (value*100) rows_for_typical_systematics_table[systematic_group] = [label, text] + allErrors2 = np.power( np.array(allErrors), 2) + sumAllErrors = np.sqrt( np.sum(allErrors2,axis=0) ) + medianAllErrors = np.median( sumAllErrors ) + printout = '%% ' + '=' * 60 printout += '\n' printout += '%% Typical systematics table for %s channel, k-value %s, met type %s, %s b-tag region\n' % (channel, str(k_values[channel]), met_type, b_tag_bin) @@ -490,8 +572,8 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr printout += '\\begin{table}[htbp]\n' printout += '\\centering\n' - printout += '\\caption{Typical systematic uncertainties in percent (median values) for the normalised \\ttbar\n' - printout += '\\differential cross section measurement at \ensuremath{\roots=%d\TeV} ' + printout += '\\caption{Typical systematic uncertainties in percent (median values) in the normalized \\ttbar cross\n' + printout += 'section measurement at a center-of-mass energy of 8 TeV' if channel == 'combined': printout += '(combination of electron and muon channels). ' else: @@ -500,25 +582,27 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr printout += '\\label{tab:typical_systematics_%dTeV_%s}\n' % (measurement_config.centre_of_mass_energy, channel) printout += '\\resizebox{\\columnwidth}{!} {\n' - printout += '\\begin{tabular}{l' + 'r'*len(variable_bins_ROOT) + '}\n' + printout += '\\tiny\n' + printout += '\\begin{tabular}{c' + 'c'*len(variable_bins_ROOT) + '}\n' printout += '\\hline\n' - + printout += ' & \\multicolumn{4}{c}{Relative uncertainty (\\%) } \\\\\n' header = 'Uncertainty source ' - header += '& %s' % (variables_latex[variable]) + header += '& %s' % (variables_latex_macros[variable]) header += ' ' printout += header printout += '\n\\hline\n' - for systematic_group in sorted(rows_for_typical_systematics_table.keys()): + for systematic_group in typical_systematics_order: if systematic_group == 'central': continue for item in rows_for_typical_systematics_table[systematic_group]: printout += item + ' & ' printout = printout.rstrip('& ') - printout += ' \n' + printout += '\n' printout += '\\hline \n' + printout += 'Total & %.1f \n' % (medianAllErrors*100) printout += '\\hline \n' printout += '\\end{tabular}\n' printout += '}\n' @@ -537,9 +621,11 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr lines = output_file.readlines() for line_number, line in enumerate (lines): if line.startswith("Uncertainty source"): - lines[line_number] = lines[line_number].strip() + " & " + variables_latex[variable] + "\n" - elif variable == "HT" and line.startswith("$E_{T}^{miss}$ uncertainties"): + lines[line_number] = lines[line_number].strip() + " & " + variables_latex_macros[variable] + "\n" + elif variable == "HT" and line.startswith(typical_systematics_latex["typical_systematics_MET"]): lines[line_number] = lines[line_number].strip() + " & - \n" + elif line.startswith("Total"): + lines[line_number] = lines[line_number].strip() + " & %.1f \n" % (medianAllErrors*100) else: for table_entry in enumerate(typical_systematics_latex): if line.startswith(typical_systematics_latex[table_entry[1]]): @@ -554,6 +640,159 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr else: print printout +def calculate_chi2(xsections,covariance): + interestingModels = [] + if measurement_config.centre_of_mass_energy == 8: + interestingModels = [ 'madgraph', 'madgraph_ptreweight', 'mcatnlo', 'powheg_v1_pythia', 'powheg_v1_herwig', 'powheg_v2_pythia', 'powheg_v2_herwig' ] + else: + interestingModels = [ 'madgraph', 'madgraph_ptreweight', 'powheg_v2_pythia', 'powheg_v2_herwig' ] + + # interestingModels = [ 'madgraph', 'madgraph_ptreweight' ] + + # nDigits = int( abs( round( np.log10( maxCovariance / 1e7 ), 0 ) ) ) + + # nDigitsXsec = -3 + unfolded_data = [ float('%.3g' % i[0]) for i in xsections['unfolded'] ] + # print nDigitsXsec + # print unfolded_data + # print unfolded_data + # print covariance + # print np.linalg.inv(covariance) + # diag_cov = np.array(np.zeros((len(unfolded_data),len(unfolded_data)))) + # for i in range (0,len(unfolded_data)): + # diag_cov[i,i] = covariance[i,i] + # covariance = diag_cov + # print diag_cov + # covariance = np.abs(covariance) + # covariance[2,0] *= 0.5 + # covariance[0,2] *= 0.5 + # covariance[2,4] *= 0.5 + # covariance[4,2] *= 0.5 + + # print covariance + unc = [] + for i in range (0,covariance.shape[0]): + unc.append( np.sqrt( covariance[i,i] ) ) + correlation = get_correlation_matrix( unc, covariance) + + # maxCovariance = np.max(np.abs(covariance) ) + # print maxCovariance + # eps = np.finfo(float).eps + # print 'NDigits :',nDigits + # for i in range (0,covariance.shape[0]): + # for j in range (0,covariance.shape[1]): + # covariance[i,j] = round( covariance[i,j], nDigits ) + # print covariance[i,j] + # if abs( maxCovariance - covariance[i,j] ) < 1e-2 * maxCovariance : + # print 'Small covariance : ',maxCovariance,covariance[i,j],maxCovariance - covariance[i,j] + # covariance[i,j] = Decimal(0) + # else: + # covariance[i,j] = Decimal( covariance[i,j] ) + + # print correlation + # print np.max(correlation) + # inverse = np.linalg.inv(covariance) + + # print inverse + # # maxInverse = np.max(inverse) + # # for i in range (0,inverse.shape[0]): + # # for j in range (0,inverse.shape[1]): + # # if abs( maxInverse - inverse[i,j]) < eps: + # # inverse[i,j] == 0. + # # print covariance[i,j] print 'FLOAT INFO :',np.finfo(float).eps + + # print 'CHECK :',inverse.dot( covariance ) + + # print unfolded_data + # print covariance + + for model in interestingModels: + xsecs = np.array( [ float('%.3g' % i[0]) for i in xsections[model] ] ) + + # for i,j in zip( xsecs, unfolded_data ): + # print i,j + # xsecs = np.array( [ i*0.98 for i in unfolded_data ] ) + # print xsecs + # print model + # print unfolded_data,xsecs + # print unfolded_data - xsecs + # # print unc + # # c = 0 + # # for i,u in zip( unfolded_data - xsecs, unc): + # # c += i**2/u**2 + # print np.linalg.inv(covariance) + # print c + # print np.linalg.inv(covariance)[0,:] + # print (unfolded_data - xsecs).dot ( np.linalg.inv(covariance)[:,0] ) + # n = 0 + # for i,j in zip( (unfolded_data - xsecs), np.linalg.inv(covariance)[:,0]): + # print i,j,i*j + # n += i * j + # print n + # print (unfolded_data - xsecs).dot( np.linalg.inv(covariance) ) + # diff = (unfolded_data - xsecs) + # variance = [diag_cov[i,i] for i in range(0,diag_cov.shape[0])] + # diffOverError = [ d**2/v for d,v in zip( diff, variance) ] + # print diff + # print variance + # print diffOverError + # simple_chi2 = sum( diffOverError ) + # print simple_chi2 / diag_cov.shape[0] + # print ( ( unfolded_data - xsecs ).dot( np.linalg.inv(covariance) ) ) + # print (unfolded_data - xsecs) + chi2 = ( np.transpose( unfolded_data - xsecs ).dot( np.linalg.inv(covariance) ) ).dot( unfolded_data - xsecs ) + ndf = covariance.shape[0] + prob = TMath.Prob( chi2, ndf ) + print model,chi2,ndf,prob + # print model,chi2, chi2// ( diag_cov.shape[0] - 1 ) + +def print_covariance( full_covariance, variable, com ): + nBins = full_covariance.shape[0] + + printout = '' + printout += '*dscomment: Covariance matrix for the normalized tt differential cross section measurements with respect to the $%s$ variable at a center-of-mass energy of %s TeV. Both statistical and systematic effects are considered.\n' % ( variables_hepdata[variable], str(com) ) + printout += '*reackey: P P --> TOP TOPBAR X\n' + printout += '*obskey: DSIG/D%s\n' % ( variable ) + printout += '*qual: Cross section : $d\sigma/d%s$\n' % ( variables_hepdata[variable] ) + printout += '*qual: SQRT(S) IN GEV : %s000.0\n' % ( com ) + + printout += '*yheader' + i = 1 + while i <= nBins: + printout += ': ' + printout += str(i) + i += 1 + printout += '\n' + + printout += '*xheader: Bin of $%s$ \n' % ( variables_hepdata[variable] ) + + printout += '*data: x ' + + i = 1 + while i <= nBins: + printout += ': y ' + i += 1 + printout += '\n' + # printout += '*data: x : y : y : y \n' + + + + maxCovariance = np.max( np.abs(full_covariance) ) + nDigits = int( abs( round( np.log10( maxCovariance / 1e9 ), 0 ) ) ) + i = 1 + for row in full_covariance: + printout += str(i) + i += 1 + printout += ';\t' + for element in row: + printout += '{0:.7g};\t'.format( round( element, nDigits ) ) + printout += '\n' + + printout += '*dataend\n' + + print '\n' + print printout + if __name__ == '__main__': parser = OptionParser() parser.add_option("-p", "--path", dest="path", default='data/', @@ -579,6 +818,7 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr vjets_theory_systematic_prefix = measurement_config.vjets_theory_systematic_prefix met_systematics_suffixes = measurement_config.met_systematics_suffixes typical_systematics = measurement_config.typical_systematics + typical_systematics_order = measurement_config.typical_systematics_order variable = options.variable output_folder = options.output_folder @@ -624,11 +864,71 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr all_measurements.extend(new_uncertainties) all_measurements.extend(rate_changing_systematics) - for channel in ['electron', 'muon', 'combined']: + other_uncertainties_list = deepcopy( measurement_config.categories_and_prefixes.keys() ) + other_uncertainties_list.extend( vjets_generator_systematics ) + other_uncertainties_list.append( 'QCD_shape' ) + other_uncertainties_list.extend( rate_changing_systematics ) + + # for channel in ['electron', 'muon', 'combined']: + for channel in ['combined']: normalised_xsection_measured_unfolded, normalised_xsection_measured_errors, normalised_xsection_unfolded_errors = read_xsection_measurement_results_with_errors(channel) - print_xsections(normalised_xsection_measured_unfolded, channel, toFile = True, print_before_unfolding = False) - print_xsections(normalised_xsection_measured_unfolded, channel, toFile = True, print_before_unfolding = True) + stat_covariance = np.loadtxt(path_to_JSON + '/' + variable + '/xsection_measurement_results/combined/central/covariance.txt',delimiter=',') + syst_covariance = np.loadtxt(path_to_JSON + '/' + variable + '/xsection_measurement_results/combined/central/covariance_systematic.txt',delimiter=',') + + # print stat_covariance + # print syst_covariance + full_covariance = stat_covariance + syst_covariance + # full_covariance = stat_covariance + + # for i in range(0,full_covariance.shape[0]): + # for j in range(0,full_covariance.shape[1]): + # if i != j : + # full_covariance[i,j] = 0 + # print full_covariance + + xsec = normalised_xsection_measured_unfolded['unfolded_with_systematics'] + + nBins = full_covariance.shape[0] + + for i in range(0,nBins): + stat = '%.2g' % ( np.sqrt( stat_covariance[i,i] ) / xsec[i][0] * 100 ) + sys = '%.2g' % ( np.sqrt( syst_covariance[i,i] ) / xsec[i][0] * 100 ) + full = '%.2g' % ( np.sqrt( full_covariance[i,i] ) / xsec[i][0] * 100 ) + print i, stat, sys, full + + + tableErrors = print_xsections(normalised_xsection_measured_unfolded, channel, toFile = True, print_before_unfolding = False) + # print_xsections(normalised_xsection_measured_unfolded, channel, toFile = True, print_before_unfolding = True) + + # for i in range(0,nBins): + # print i,stat_covariance[i,i],np.sqrt(stat_covariance[i,i]) / xsec[i][0] * 100 + # print i,np.sqrt(syst_covariance[i,i]) / xsec[i][0] * 100 + # print i,'%.2g' % ( np.sqrt(stat_covariance[i,i]) / xsec[i][0] * 100),'%.2g' % (np.sqrt(syst_covariance[i,i]) / xsec[i][0] * 100),'%.2g' % ( np.sqrt(full_covariance[i,i]) / xsec[i][0] * 100 ) + # print '%s %.2g' % ( tableErrors[i], np.sqrt(full_covariance[i,i]) / xsec[i][0] * 100 ) + + # print 'Stat' + # calculate_chi2(normalised_xsection_measured_unfolded,stat_covariance) + # print 'Sys' + # calculate_chi2(normalised_xsection_measured_unfolded,syst_covariance) + # print 'Full' + # maxCovariance = np.max(np.abs(full_covariance) ) + # # print maxCovariance + # # eps = np.finfo(float).eps + # nDigits = int( abs( round( np.log10( maxCovariance / 1e7 ), 0 ) ) ) + # print nDigits + # # print unfolded_data + # for i in range(0,full_covariance.shape[0]): + # for j in range(0,full_covariance.shape[1]): + # full_covariance[i,j] = round( full_covariance[i,j], nDigits ) + # # print full_covariance + calculate_chi2(normalised_xsection_measured_unfolded,full_covariance) + + path = output_folder + '/' + str(measurement_config.centre_of_mass_energy) + 'TeV/' + variable + filename = path + '/fullCovariance.txt' + np.savetxt( filename, full_covariance, delimiter = ',' ) + + print_covariance( full_covariance, variable, measurement_config.centre_of_mass_energy ) print_error_table(normalised_xsection_measured_unfolded, normalised_xsection_unfolded_errors, channel, toFile = True, print_before_unfolding = False) print_error_table(normalised_xsection_measured_unfolded, normalised_xsection_measured_errors, channel, toFile = True, print_before_unfolding = True) @@ -637,9 +937,9 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr print_typical_systematics_table(normalised_xsection_measured_unfolded, normalised_xsection_unfolded_errors, channel, toFile = True, print_before_unfolding = False) # print_typical_systematics_table(normalised_xsection_measured_unfolded, normalised_xsection_measured_errors, channel, toFile = True, print_before_unfolding = True) - if not channel == 'combined': - fit_input = read_initial_normalisation(path_to_JSON, variable, 'central', channel, met_type) - fit_results = read_normalisation(path_to_JSON, variable, 'central', channel, met_type) - print_fit_results_table(fit_input, fit_results, channel, toFile = True) + # if not channel == 'combined': + # fit_input = read_initial_normalisation(path_to_JSON, variable, 'central', channel, met_type) + # fit_results = read_normalisation(path_to_JSON, variable, 'central', channel, met_type) + # print_fit_results_table(fit_input, fit_results, channel, toFile = True) diff --git a/src/cross_section_measurement/98b_fit_cross_checks.py b/src/cross_section_measurement/98b_fit_cross_checks.py index ae8be4c7..6eb324dc 100644 --- a/src/cross_section_measurement/98b_fit_cross_checks.py +++ b/src/cross_section_measurement/98b_fit_cross_checks.py @@ -125,6 +125,7 @@ def create_latex_string( fit_var_input ): initial_values_electron = {} initial_values_muon = {} for fit_variable in fit_variables: + path = path_to_JSON + fit_variable + '/' + str( come ) + 'TeV/' fit_results_electron[fit_variable] = read_normalisation( path, variable, @@ -149,16 +150,16 @@ def create_latex_string( fit_var_input ): 'muon', met_type ) - if not 'closure' in path_to_JSON: - fit_results_electron['before'] = read_normalisation( 'data_single_var_fit/8TeV/', - variable, - category, - 'electron', - met_type ) - fit_results_muon['before'] = read_normalisation( 'data_single_var_fit/8TeV/', - variable, - category, - 'muon', - met_type ) + # if not 'closure' in path_to_JSON: + # fit_results_electron['before'] = read_normalisation( 'data_single_var_fit/8TeV/', + # variable, + # category, + # 'electron', + # met_type ) + # fit_results_muon['before'] = read_normalisation( 'data_single_var_fit/8TeV/', + # variable, + # category, + # 'muon', + # met_type ) plot_fit_results( fit_results_electron, initial_values_electron, 'electron' ) plot_fit_results( fit_results_muon, initial_values_muon, 'muon' ) diff --git a/src/cross_section_measurement/compare_fit_background_subtraction.py b/src/cross_section_measurement/compare_fit_background_subtraction.py new file mode 100644 index 00000000..54e46ea5 --- /dev/null +++ b/src/cross_section_measurement/compare_fit_background_subtraction.py @@ -0,0 +1,166 @@ +''' +Approval conditions for TOP-15-013 +''' +from __future__ import division +from tools.plotting import Histogram_properties, compare_histograms, Plot +from tools.file_utilities import read_data_from_JSON +from tools.hist_utilities import value_error_tuplelist_to_hist, value_tuplelist_to_hist +from config.variable_binning import bin_edges +from config.latex_labels import variables_latex +from math import sqrt +def getHistograms( fileName, variable, label, plotErrors=False ): + + data = read_data_from_JSON(fileName)[label] + for bin in data: + while len(bin) > 2 : + del bin[-1] + + # Interested in errors, rather than central value + if plotErrors: + for bin in data: + binContent = bin[0] + del bin[0] + bin.append(0) + bin[0] /= binContent + bin[0] *= 100 + # h = value_tuplelist_to_hist( data, bin_edges[variable] ) + # else: + h=value_error_tuplelist_to_hist( data, bin_edges[variable] ) + return h + +def useOnlyNormalisationErrors(h, fileName, variable, label): + bkgHists = {} + + bkgErr = { + 'QCD' : 1., + 'SingleTop' : 0.3, + 'V+Jets' : 0.3 + } + + print list(h.y()) + var = None + firstBkg = True + for bkg in ['QCD', 'SingleTop', 'V+Jets']: + bkgHist = getHistograms(fileName, variable, bkg) + + if firstBkg: + var = ( bkgHist * ( bkgErr[bkg] ) / h ) ** 2 + firstBkg = False + else : + var += ( bkgHist * ( bkgErr[bkg] ) / h ) ** 2 + + # Add luminosity + var += 0.022**2 + + for bin in range(0,h.GetNbinsX()): + newBinError2 = var.GetBinContent(bin) + var.SetBinContent(bin, sqrt(newBinError2)) + h.SetBinError(bin, sqrt(newBinError2)) + + + + +def compare_yield_methods( measurement='fit_results', + com=8, + channel='electron', + compareErrors=False, + statOnly=False): + + file_template = None + label = None + fileSuffix = '' + if compareErrors and not statOnly: + fileSuffix = '_with_errors' + if measurement == 'fit_results': + file_template = 'data/{method}/{com}TeV/{variable}/central/normalisation_{channel}_patType1CorrectedPFMet{suffix}.txt' + label = 'TTJet' + elif measurement == 'normalised_xsection': + file_template = 'data/{method}/{com}TeV/{variable}/xsection_measurement_results/{channel}/central/normalised_xsection_patType1CorrectedPFMet{suffix}.txt' + label = 'TTJet_unfolded' + + variables = ['MET', 'HT', 'ST','WPT'] + # variables = ['HT'] + + for variable in variables: + fit = file_template.format( + method='absolute_eta_M3_angle_bl', + com=com, + variable=variable, + channel=channel, + suffix=fileSuffix + ) + background = file_template.format( + method='normalisation/background_subtraction', + com=com, + variable=variable, + channel=channel, + suffix=fileSuffix + ) + + h_fit = getHistograms( fit, variable, label, compareErrors ) + h_bkg = getHistograms( background, variable, label, compareErrors ) + + # if measurement is 'fit_results': + # useOnlyNormalisationErrors(h_bkg, background, variable, label) + # for bin in range(0,h_fit.GetNbinsX()): + # h_fit.SetBinError(bin,0) + + + + properties = Histogram_properties() + properties.name = '{0}_compare_yield_methods_{1}_{2}_{3}'.format(measurement, variable, com, channel) + if compareErrors: + if statOnly: + properties.name = properties.name.replace('yield','StatErrors') + else: + properties.name = properties.name.replace('yield','StatPlusSysErrors') + properties.title = 'Comparison of ttbar yield methods' + properties.path = 'plots/compareFitBackground/{0}TeV/{1}/'.format(com,variable) + properties.has_ratio = True + properties.xerr = True + properties.x_limits = ( + bin_edges[variable][0], bin_edges[variable][-1]) + properties.x_axis_title = variables_latex[variable] + + if measurement == 'normalised_xsection': + if compareErrors: + if statOnly: + properties.y_axis_title = r'Statistical error in $\frac{1}{\sigma} \frac{d\sigma}{d' + \ + variables_latex[variable] + '}$' + properties.ratio_y_limits = [0.5,1.5] + else: + properties.y_axis_title = r'Total error in $\frac{1}{\sigma} \frac{d\sigma}{d' + \ + variables_latex[variable] + '}$' + properties.ratio_y_limits = [0.5,1.5] + else : + properties.y_axis_title = r'$\frac{1}{\sigma} \frac{d\sigma}{d' + \ + variables_latex[variable] + '}$' + properties.ratio_y_limits = [0.9,1.1] + else : + properties.y_axis_title = r'$t\bar{t}$ normalisation' + properties.ratio_y_limits = [0.7,1.3] + + histograms = {'Fit': h_fit, 'Background subtraction': h_bkg} + plot = Plot(histograms, properties) + plot.draw_method = 'errorbar' + compare_histograms(plot) + +if __name__ == '__main__': + import sys + if '-d' in sys.argv: + from tools.logger import log + log.setLevel(log.DEBUG) + + # compare_yield_methods( measurement='fit_results', com=7, channel='electron') + # compare_yield_methods( measurement='fit_results', com=7, channel='muon') + # compare_yield_methods( measurement='fit_results', com=8, channel='electron') + # compare_yield_methods( measurement='fit_results', com=8, channel='muon') + + # compare_yield_methods( measurement='normalised_xsection', com=8, channel='combined') + # compare_yield_methods( measurement='normalised_xsection', com=7, channel='combined') + + compare_yield_methods( measurement='normalised_xsection', com=8, channel='combined', compareErrors=True) + compare_yield_methods( measurement='normalised_xsection', com=7, channel='combined', compareErrors=True) + + compare_yield_methods( measurement='normalised_xsection', com=8, channel='combined', compareErrors=True, statOnly=True) + compare_yield_methods( measurement='normalised_xsection', com=7, channel='combined', compareErrors=True, statOnly=True) \ No newline at end of file diff --git a/src/cross_section_measurement/create_measurement.py b/src/cross_section_measurement/create_measurement.py index cb1ea875..2c968d6a 100644 --- a/src/cross_section_measurement/create_measurement.py +++ b/src/cross_section_measurement/create_measurement.py @@ -128,13 +128,13 @@ def create_measurement(com, category, variable, channel): 'Ref selection', config.electron_control_region) if category == 'QCD_shape': qcd_template = qcd_template.replace( - 'Ref selection', config.electron_control_region_systematic) + config.electron_control_region, config.electron_control_region_systematic) else: qcd_template = qcd_template.replace( 'Ref selection', config.muon_control_region) if category == 'QCD_shape': qcd_template = qcd_template.replace( - 'Ref selection', config.muon_control_region_systematic) + config.muon_control_region, config.muon_control_region_systematic) m_qcd.addSample('TTJet', config.ttbar_category_templates[ template_category], qcd_template, False) diff --git a/src/cross_section_measurement/make_control_plots.py b/src/cross_section_measurement/make_control_plots.py index d368e8e4..b5042b0c 100644 --- a/src/cross_section_measurement/make_control_plots.py +++ b/src/cross_section_measurement/make_control_plots.py @@ -9,6 +9,8 @@ from tools.hist_utilities import prepare_histograms, clean_control_region from tools.ROOT_utils import get_histograms_from_files, set_root_defaults from tools.latex import setup_matplotlib +from math import sqrt +# from numpy import sqrt # latex, font, etc setup_matplotlib() @@ -42,6 +44,22 @@ def get_normalisation_error( normalisation ): total_error += number[1] return total_error / total_normalisation +def get_total_error_on_mc( h_ttjet, h_others, ttbar_uncertainty ): + totalError2 = [] + for bin in range(1,h_ttjet.GetNbinsX() ): + statError = h_ttjet.GetBinError( bin ) + normError = h_ttjet.GetBinContent( bin ) * ttbar_uncertainty + ttjetError2 = statError ** 2 + normError ** 2 + totalError2.append( ttjetError2 ) + + for h in h_others: + for bin in range(1,h.GetNbinsX() ): + statError = h.GetBinError( bin ) + totalError2[ bin - 1 ] += statError ** 2 + + totalError = [ sqrt(i) for i in totalError2 ] + return totalError + def compare_shapes( channel, x_axis_title, y_axis_title, control_region_1, control_region_2, name_region_1, name_region_2, @@ -105,12 +123,11 @@ def make_plot( channel, x_axis_title, y_axis_title, legend_location = ( 0.98, 0.78 ), cms_logo_location = 'right', log_y = False, legend_color = False, - ratio_y_limits = [0.3, 1.7], + ratio_y_limits = [0.9, 1.1], normalise = False, ): global output_folder, measurement_config, category, normalise_to_fit global preliminary, norm_variable, sum_bins, b_tag_bin, histogram_files - qcd_data_region = '' title = title_template % ( measurement_config.new_luminosity / 1000., measurement_config.centre_of_mass_energy ) normalisation = None @@ -189,8 +206,28 @@ def make_plot( channel, x_axis_title, y_axis_title, signal_region_hists['V+Jets'], signal_region_hists['SingleTop'], signal_region_hists['TTJet']] - histogram_lables = ['data', 'QCD', 'V+Jets', 'Single-Top', samples_latex['TTJet']] - histogram_colors = ['black', 'yellow', 'green', 'magenta', 'red'] + histogram_lables = ['Data', 'QCD', 'W/Z+Jets', 'Single t', samples_latex['TTJet']] + # histogram_colors = ['black', 'yellow', 'green', 'magenta', 'red'] + # histogram_colors = ['black', 390, 413, 616, 633] + histogram_colors = ['black', 390, 414, 606, 633] + + totalErrorList = get_total_error_on_mc( signal_region_hists['TTJet'], + [ signal_region_hists['V+Jets'], signal_region_hists['SingleTop'], qcd_from_data ], + mc_uncertainty ) + + # print signal_region_hists['data'].Integral() + # print signal_region_hists['TTJet'].Integral() + # print signal_region_hists['SingleTop'].Integral() + # print signal_region_hists['V+Jets'].Integral() + # print qcd_from_data.Integral() + + # sumMC = signal_region_hists['TTJet'].Integral() + signal_region_hists['SingleTop'].Integral() + signal_region_hists['V+Jets'].Integral() + qcd_from_data.Integral() + # print 'Total MC :',sumMC + # print 'Percents' + # print 'TTJet :',signal_region_hists['TTJet'].Integral() / sumMC * 100 + # print 'Single Top :',signal_region_hists['SingleTop'].Integral() / sumMC * 100 + # print 'V+Jets :',signal_region_hists['V+Jets'].Integral() / sumMC * 100 + # print 'QCD :',qcd_from_data.Integral() / sumMC * 100 histogram_properties = Histogram_properties() histogram_properties.name = name_prefix + b_tag_bin @@ -205,10 +242,7 @@ def make_plot( channel, x_axis_title, y_axis_title, histogram_properties.xerr = None # workaround for rootpy issue #638 histogram_properties.emptybins = True - if b_tag_bin: - histogram_properties.additional_text = channel_latex[channel] + ', ' + b_tag_bins_latex[b_tag_bin] - else: - histogram_properties.additional_text = channel_latex[channel] + histogram_properties.additional_text = channel_latex[channel] histogram_properties.legend_location = legend_location histogram_properties.cms_logo_location = cms_logo_location histogram_properties.preliminary = preliminary @@ -219,15 +253,16 @@ def make_plot( channel, x_axis_title, y_axis_title, if normalise_to_fit: histogram_properties.mc_error = get_normalisation_error( normalisation ) - histogram_properties.mc_errors_label = 'fit uncertainty' + histogram_properties.mc_errors_label = 'Fit uncertainty' else: - histogram_properties.mc_error = mc_uncertainty + # histogram_properties.mc_error = mc_uncertainty + histogram_properties.mc_error = totalErrorList histogram_properties.mc_errors_label = 'MC unc.' - make_data_mc_comparison_plot( histograms_to_draw, histogram_lables, histogram_colors, - histogram_properties, save_folder = output_folder, - show_ratio = False, normalise = normalise, - ) + # make_data_mc_comparison_plot( histograms_to_draw, histogram_lables, histogram_colors, + # histogram_properties, save_folder = output_folder, + # show_ratio = False, normalise = normalise, + # ) histogram_properties.name += '_with_ratio' loc = histogram_properties.legend_location # adjust legend location as it is relative to canvas! @@ -274,8 +309,11 @@ def make_plot( channel, x_axis_title, y_axis_title, # this is shown as \ttbar (or MC) uncertainty on the plots # in fact, it takes the uncertainty on the whole MC stack # although unimportant, this needs revision - mc_uncertainty = 0.10 - + # mc_uncertainty = 0.10 + mc_uncertainty = 0.0643 # 8 TeV + # mc_uncertainty = 0.0683 # 7 TeV + # mc_uncertainty = 0 + histogram_files = { 'TTJet': measurement_config.ttbar_category_templates[category], 'V+Jets': measurement_config.VJets_category_templates[category], @@ -303,37 +341,37 @@ def make_plot( channel, x_axis_title, y_axis_title, preliminary = False b_tag_bin = '2orMoreBtags' - norm_variable = 'MET' + norm_variable = 'HT' # comment out plots you don't want include_plots = [ - 'eta', - 'pT', + # 'eta', + # 'pT', 'MET', - 'MET log', - 'MET phi', + # # 'MET log', + # # 'MET phi', 'HT', 'ST', 'WPT', - 'MT', - 'M3', - 'angle_bl', - 'bjet invariant mass', - 'b-tag multiplicity', - 'b-tag multiplicity reweighted', - 'jet multiplicity', - 'n vertex', - 'n vertex reweighted', + # # 'MT', + # 'M3', + # 'angle_bl', + # 'bjet invariant mass', + # 'b-tag multiplicity', + # 'b-tag multiplicity reweighted', + # 'jet multiplicity', + # 'n vertex', + # 'n vertex reweighted', ] additional_qcd_plots = [ - 'eta in MET bins', - 'eta in HT bins', - 'eta in ST bins', - 'eta in MT bins', - 'eta in WPT bins', - 'QCD PFReliso non-iso', - 'QCD PFReliso', - 'QCD eta', - 'QCD eta shapes', + # 'eta in MET bins', + # 'eta in HT bins', + # 'eta in ST bins', + # 'eta in MT bins', + # 'eta in WPT bins', + # 'QCD PFReliso non-iso', + # 'QCD PFReliso', + # 'QCD eta', + # 'QCD eta shapes', ] if make_additional_QCD_plots: include_plots.extend( additional_qcd_plots ) @@ -341,13 +379,14 @@ def make_plot( channel, x_axis_title, y_axis_title, # lepton |eta| ################################################### if 'eta' in include_plots: + print '------> eta' make_plot( 'electron', x_axis_title = '$\left|\eta(\mathrm{e})\\right|$', y_axis_title = 'Events/(0.1)', signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/Electron/electron_AbsEta_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'electron_AbsEta_', - x_limits = [0, 2.6], + x_limits = [0, 2.5], rebin = 10, legend_location = ( 0.98, 0.78 ), cms_logo_location = 'right', @@ -358,7 +397,7 @@ def make_plot( channel, x_axis_title, y_axis_title, signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/Muon/muon_AbsEta_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'muon_AbsEta_', - x_limits = [0, 2.6], + x_limits = [0, 2.5], rebin = 10, legend_location = ( 0.98, 0.78 ), cms_logo_location = 'right', @@ -368,7 +407,7 @@ def make_plot( channel, x_axis_title, y_axis_title, ################################################### if 'pT' in include_plots: make_plot( 'electron', - x_axis_title = '$p_\mathrm{T}(\mathrm{e})$ [GeV]', + x_axis_title = '$p_\mathrm{T}(\mathrm{e})$ (GeV)', y_axis_title = 'Events/(10 GeV)', signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/Electron/electron_pT_' + b_tag_bin, qcd_data_region_btag = '0btag', @@ -379,7 +418,7 @@ def make_plot( channel, x_axis_title, y_axis_title, cms_logo_location = 'right', ) make_plot( 'muon', - x_axis_title = '$p_\mathrm{T}(\mu)$ [GeV]', + x_axis_title = '$p_\mathrm{T}(\mu)$ (GeV)', y_axis_title = 'Events/(10 GeV)', signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/Muon/muon_pT_' + b_tag_bin, qcd_data_region_btag = '0btag', @@ -393,34 +432,38 @@ def make_plot( channel, x_axis_title, y_axis_title, # MET ################################################### if 'MET' in include_plots: + print '------> MET' + make_plot( 'electron', - x_axis_title = '$%s$ [GeV]' % variables_latex['MET'], - y_axis_title = 'Events/(5 GeV)', + x_axis_title = '$%s$ (GeV)' % variables_latex['MET'], + y_axis_title = 'Events/(10 GeV)', signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/MET/patType1CorrectedPFMet/MET_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'EPlusJets_patType1CorrectedPFMet_', x_limits = [0, 200], - rebin = 1, - legend_location = ( 0.98, 0.78 ), - cms_logo_location = 'right', + ratio_y_limits = [0.7, 1.3], + rebin = 2, + legend_location = ( 0.89, 0.78 ), + cms_logo_location = 'left', ) make_plot( 'muon', - x_axis_title = '$%s$ [GeV]' % variables_latex['MET'], - y_axis_title = 'Events/(5 GeV)', + x_axis_title = '$%s$ (GeV)' % variables_latex['MET'], + y_axis_title = 'Events/(10 GeV)', signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/MET/patType1CorrectedPFMet/MET_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'MuPlusJets_patType1CorrectedPFMet_', x_limits = [0, 200], - rebin = 1, - legend_location = ( 0.98, 0.78 ), - cms_logo_location = 'right', + ratio_y_limits = [0.7, 1.3], + rebin = 2, + legend_location = ( 0.89, 0.78 ), + cms_logo_location = 'left', ) ################################################### # MET log ################################################### if 'MET log' in include_plots: make_plot( 'electron', - x_axis_title = '$%s$ [GeV]' % variables_latex['MET'], + x_axis_title = '$%s$ (GeV)' % variables_latex['MET'], y_axis_title = 'Events/(20 GeV)', signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/MET/patType1CorrectedPFMet/MET_' + b_tag_bin, qcd_data_region_btag = '0btag', @@ -432,7 +475,7 @@ def make_plot( channel, x_axis_title, y_axis_title, log_y = True, ) make_plot( 'muon', - x_axis_title = '$%s$ [GeV]' % variables_latex['MET'], + x_axis_title = '$%s$ (GeV)' % variables_latex['MET'], y_axis_title = 'Events/(20 GeV)', signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/MET/patType1CorrectedPFMet/MET_' + b_tag_bin, qcd_data_region_btag = '0btag', @@ -476,81 +519,90 @@ def make_plot( channel, x_axis_title, y_axis_title, ################################################### norm_variable = 'HT' if 'HT' in include_plots: + print '------> HT' make_plot( 'electron', - x_axis_title = '$%s$ [GeV]' % variables_latex['HT'], - y_axis_title = 'Events/(20 GeV)', + x_axis_title = '$%s$ (GeV)' % variables_latex['HT'], + y_axis_title = 'Events/(40 GeV)', signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/MET/HT_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'EPlusJets_HT_', - x_limits = [100, 1000], - rebin = 4, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', + x_limits = [120, 1000], + ratio_y_limits = [0.7, 1.3], + rebin = 8, + legend_location = ( 0.89, 0.78 ), + cms_logo_location = 'left', ) make_plot( 'muon', - x_axis_title = '$%s$ [GeV]' % variables_latex['HT'], - y_axis_title = 'Events/(20 GeV)', + x_axis_title = '$%s$ (GeV)' % variables_latex['HT'], + y_axis_title = 'Events/(40 GeV)', signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/MET/HT_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'MuPlusJets_HT_', - x_limits = [100, 1000], - rebin = 4, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', + x_limits = [120, 1000], + ratio_y_limits = [0.7, 1.3], + rebin = 8, + legend_location = ( 0.89, 0.78 ), + cms_logo_location = 'left', ) ################################################### # ST ################################################### norm_variable = 'ST' if 'ST' in include_plots: + print '------> ST' make_plot( 'electron', - x_axis_title = '$%s$ [GeV]' % variables_latex['ST'], - y_axis_title = 'Events/(20 GeV)', + x_axis_title = '$%s$ (GeV)' % variables_latex['ST'], + y_axis_title = 'Events/(40 GeV)', signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/MET/patType1CorrectedPFMet/ST_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'EPlusJets_patType1CorrectedPFMet_ST_', x_limits = [150, 1200], - rebin = 4, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', + ratio_y_limits = [0.7, 1.3], + rebin = 8, + legend_location = ( 0.89, 0.78 ), + cms_logo_location = 'left', ) make_plot( 'muon', - x_axis_title = '$%s$ [GeV]' % variables_latex['ST'], - y_axis_title = 'Events/(20 GeV)', + x_axis_title = '$%s$ (GeV)' % variables_latex['ST'], + y_axis_title = 'Events/(40 GeV)', signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/MET/patType1CorrectedPFMet/ST_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'MuPlusJets_patType1CorrectedPFMet_ST_', x_limits = [150, 1200], - rebin = 4, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', + ratio_y_limits = [0.7, 1.3], + rebin = 8, + legend_location = ( 0.89, 0.78 ), + cms_logo_location = 'left', ) ################################################### # WPT ################################################### norm_variable = 'WPT' if 'WPT' in include_plots: + print '------> WPT' make_plot( 'electron', - x_axis_title = '$%s$ [GeV]' % variables_latex['WPT'], - y_axis_title = 'Events/(10 GeV)', + x_axis_title = '$%s$ (GeV)' % variables_latex['WPT'], + y_axis_title = 'Events/(20 GeV)', signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/MET/patType1CorrectedPFMet/WPT_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'EPlusJets_patType1CorrectedPFMet_WPT_', - x_limits = [0, 500], - rebin = 10, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', + x_limits = [0, 400], + ratio_y_limits = [0.7, 1.3], + rebin = 20, + legend_location = ( 0.89, 0.78 ), + cms_logo_location = 'left', ) make_plot( 'muon', - x_axis_title = '$%s$ [GeV]' % variables_latex['WPT'], - y_axis_title = 'Events/(10 GeV)', + x_axis_title = '$%s$ (GeV)' % variables_latex['WPT'], + y_axis_title = 'Events/(20 GeV)', signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/MET/patType1CorrectedPFMet/WPT_' + b_tag_bin, qcd_data_region_btag = '0btag', name_prefix = 'MuPlusJets_patType1CorrectedPFMet_WPT_', - x_limits = [0, 500], - rebin = 10, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', + x_limits = [0, 400], + ratio_y_limits = [0.7, 1.3], + rebin = 20, + legend_location = ( 0.89, 0.78 ), + cms_logo_location = 'left', ) ################################################### # MT @@ -558,7 +610,7 @@ def make_plot( channel, x_axis_title, y_axis_title, norm_variable = 'MT' if 'MT' in include_plots: make_plot( 'electron', - x_axis_title = '$%s$ [GeV]' % variables_latex['MT'], + x_axis_title = '$%s$ (GeV)' % variables_latex['MT'], y_axis_title = 'Events/(5 GeV)', signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/MET/patType1CorrectedPFMet/Transverse_Mass_' + b_tag_bin, qcd_data_region_btag = '0btag', @@ -569,7 +621,7 @@ def make_plot( channel, x_axis_title, y_axis_title, cms_logo_location = 'right', ) make_plot( 'muon', - x_axis_title = '$%s$ [GeV]' % variables_latex['MT'], + x_axis_title = '$%s$ (GeV)' % variables_latex['MT'], y_axis_title = 'Events/(5 GeV)', signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/MET/patType1CorrectedPFMet/Transverse_Mass_' + b_tag_bin, qcd_data_region_btag = '0btag', @@ -583,7 +635,7 @@ def make_plot( channel, x_axis_title, y_axis_title, ################################################### # M3 ################################################### - norm_variable = 'MT' + norm_variable = 'HT' if 'M3' in include_plots: # M3 histograms are not plotted in the histogram output files from analysis software # so sum the M3 histograms from the BAT output histogram file and sum over all bins @@ -591,27 +643,29 @@ def make_plot( channel, x_axis_title, y_axis_title, tmp = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/Binned_MT_Analysis/MT_with_patType1CorrectedPFMet_bin_%s/M3_' + b_tag_bin regions = [tmp % bin_i for bin_i in variable_bins_ROOT['MT']] make_plot( 'electron', - x_axis_title = '$M3$ [GeV]', + x_axis_title = fit_variables_latex['M3'] + '(GeV)', y_axis_title = 'Events/(20 GeV)', signal_region = regions, qcd_data_region_btag = '0btag', name_prefix = 'EPlusJets_M3_', x_limits = [0, 1000], + ratio_y_limits = [0.8, 1.2], rebin = 4, - legend_location = ( 0.95, 0.78 ), + legend_location = ( 0.94, 0.78 ), cms_logo_location = 'right', ) tmp = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/Binned_MT_Analysis/MT_with_patType1CorrectedPFMet_bin_%s/M3_' + b_tag_bin regions = [tmp % bin_i for bin_i in variable_bins_ROOT['MT']] make_plot( 'muon', - x_axis_title = '$M3$ [GeV]', + x_axis_title = fit_variables_latex['M3'] + '(GeV)', y_axis_title = 'Events/(20 GeV)', signal_region = regions, qcd_data_region_btag = '0btag', name_prefix = 'MuPlusJets_M3_', x_limits = [0, 1000], + ratio_y_limits = [0.8, 1.2], rebin = 4, - legend_location = ( 0.95, 0.78 ), + legend_location = ( 0.94, 0.78 ), cms_logo_location = 'right', ) ################################################### @@ -624,12 +678,12 @@ def make_plot( channel, x_axis_title, y_axis_title, tmp = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/Binned_MT_Analysis/MT_with_patType1CorrectedPFMet_bin_%s/angle_bl_' + b_tag_bin regions = [tmp % bin_i for bin_i in variable_bins_ROOT['MT']] make_plot( 'electron', - x_axis_title = fit_variables_latex['angle_bl'], + x_axis_title = fit_variables_latex['angle_bl'] + ' [radians]', y_axis_title = 'Events/(0.2)', signal_region = regions, qcd_data_region_btag = '1btag', name_prefix = 'EPlusJets_angle_bl_', - x_limits = [0, 4], + x_limits = [0, 3.2], rebin = 2, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', @@ -637,12 +691,12 @@ def make_plot( channel, x_axis_title, y_axis_title, tmp = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/Binned_MT_Analysis/MT_with_patType1CorrectedPFMet_bin_%s/angle_bl_' + b_tag_bin regions = [tmp % bin_i for bin_i in variable_bins_ROOT['MT']] make_plot( 'muon', - x_axis_title = fit_variables_latex['angle_bl'], + x_axis_title = fit_variables_latex['angle_bl'] + ' [radians]', y_axis_title = 'Events/(0.2)', signal_region = regions, qcd_data_region_btag = '1btag', name_prefix = 'MuPlusJets_angle_bl_', - x_limits = [0, 4], + x_limits = [0, 3.2], rebin = 2, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', @@ -689,7 +743,7 @@ def make_plot( channel, x_axis_title, y_axis_title, qcd_data_region_btag = '', use_qcd_data_region = False, name_prefix = 'EPlusJets_N_BJets', - x_limits = [1.5, 7.5], + x_limits = [-0.5, 7.5], rebin = 1, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', @@ -702,7 +756,7 @@ def make_plot( channel, x_axis_title, y_axis_title, qcd_data_region_btag = '', use_qcd_data_region = False, name_prefix = 'MuPlusJets_N_BJets', - x_limits = [1.5, 7.5], + x_limits = [-0.5, 7.5], rebin = 1, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', @@ -717,7 +771,7 @@ def make_plot( channel, x_axis_title, y_axis_title, qcd_data_region_btag = '', use_qcd_data_region = False, name_prefix = 'EPlusJets_N_BJets_reweighted', - x_limits = [1.5, 7.5], + x_limits = [-0.5, 7.5], rebin = 1, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', @@ -730,68 +784,68 @@ def make_plot( channel, x_axis_title, y_axis_title, qcd_data_region_btag = '', use_qcd_data_region = False, name_prefix = 'MuPlusJets_N_BJets_reweighted', - x_limits = [1.5, 7.5], + x_limits = [-0.5, 7.5], rebin = 1, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', log_y = False, ) - if 'b-tag multiplicity' in include_plots: - b_tag_bin = '' - make_plot( 'electron', - x_axis_title = 'B-tag multiplicity', - y_axis_title = 'Events', - signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/N_BJets', - qcd_data_region_btag = '', - use_qcd_data_region = False, - name_prefix = 'EPlusJets_N_BJets_logy', - x_limits = [1.5, 7.5], - rebin = 1, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', - log_y = True, - ) - make_plot( 'muon', - x_axis_title = 'B-tag multiplicity', - y_axis_title = 'Events', - signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/N_BJets', - qcd_data_region_btag = '', - use_qcd_data_region = False, - name_prefix = 'MuPlusJets_N_BJets_logy', - x_limits = [1.5, 7.5], - rebin = 1, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', - log_y = True, - ) - if 'b-tag multiplicity reweighted' in include_plots: - b_tag_bin = '' - make_plot( 'electron', - x_axis_title = 'B-tag multiplicity', - y_axis_title = 'Events', - signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/N_BJets_reweighted', - qcd_data_region_btag = '', - use_qcd_data_region = False, - name_prefix = 'EPlusJets_N_BJets_logy_reweighted', - x_limits = [1.5, 7.5], - rebin = 1, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', - log_y = True, - ) - make_plot( 'muon', - x_axis_title = 'B-tag multiplicity', - y_axis_title = 'Events', - signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/N_BJets', - qcd_data_region_btag = '', - use_qcd_data_region = False, - name_prefix = 'MuPlusJets_N_BJets_logy_reweighted', - x_limits = [1.5, 7.5], - rebin = 1, - legend_location = ( 0.95, 0.78 ), - cms_logo_location = 'right', - log_y = True, - ) + # if 'b-tag multiplicity' in include_plots: + # b_tag_bin = '' + # make_plot( 'electron', + # x_axis_title = 'B-tag multiplicity', + # y_axis_title = 'Events', + # signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/N_BJets', + # qcd_data_region_btag = '', + # use_qcd_data_region = False, + # name_prefix = 'EPlusJets_N_BJets_logy', + # x_limits = [1.5, 7.5], + # rebin = 1, + # legend_location = ( 0.95, 0.78 ), + # cms_logo_location = 'right', + # log_y = True, + # ) + # make_plot( 'muon', + # x_axis_title = 'B-tag multiplicity', + # y_axis_title = 'Events', + # signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/N_BJets', + # qcd_data_region_btag = '', + # use_qcd_data_region = False, + # name_prefix = 'MuPlusJets_N_BJets_logy', + # x_limits = [1.5, 7.5], + # rebin = 1, + # legend_location = ( 0.95, 0.78 ), + # cms_logo_location = 'right', + # log_y = True, + # ) + # if 'b-tag multiplicity reweighted' in include_plots: + # b_tag_bin = '' + # make_plot( 'electron', + # x_axis_title = 'B-tag multiplicity', + # y_axis_title = 'Events', + # signal_region = 'TTbar_plus_X_analysis/EPlusJets/Ref selection/N_BJets_reweighted', + # qcd_data_region_btag = '', + # use_qcd_data_region = False, + # name_prefix = 'EPlusJets_N_BJets_logy_reweighted', + # x_limits = [1.5, 7.5], + # rebin = 1, + # legend_location = ( 0.95, 0.78 ), + # cms_logo_location = 'right', + # log_y = True, + # ) + # make_plot( 'muon', + # x_axis_title = 'B-tag multiplicity', + # y_axis_title = 'Events', + # signal_region = 'TTbar_plus_X_analysis/MuPlusJets/Ref selection/N_BJets', + # qcd_data_region_btag = '', + # use_qcd_data_region = False, + # name_prefix = 'MuPlusJets_N_BJets_logy_reweighted', + # x_limits = [1.5, 7.5], + # rebin = 1, + # legend_location = ( 0.95, 0.78 ), + # cms_logo_location = 'right', + # log_y = True, + # ) ################################################### # jet multiplicity ################################################### @@ -1204,6 +1258,7 @@ def make_plot( channel, x_axis_title, y_axis_title, # QCD lepton |eta| ################################################### if 'QCD eta' in include_plots: + print '----> QCD eta' b_tag_bin = '0btag' make_plot( 'electron', x_axis_title = '$\left|\eta(\mathrm{e})\\right|$', diff --git a/src/cross_section_measurement/outputYodaResults.py b/src/cross_section_measurement/outputYodaResults.py new file mode 100644 index 00000000..971d7b40 --- /dev/null +++ b/src/cross_section_measurement/outputYodaResults.py @@ -0,0 +1,84 @@ +from tools.file_utilities import make_folder_if_not_exists, read_data_from_JSON +from config.variable_binning import bin_edges + + +# measurement_config = XSectionConfig(13) + +path = 'data/normalisation/background_subtraction/' +variables = bin_edges.keys() + +number = { + 'MET' : '5', + 'HT' : '6', + 'ST' : '7', + 'WPT' : '8' +} + +kvalues = { + 'MET' : '3', + 'HT' : '3', + 'ST' : '4', + 'WPT' : '3' +} + +for variable in variables: + if variable == 'MT': continue + # if variable != 'HT' : continue + print variable + # path_to_JSON = '{path}/{com}TeV/{variable}/xsection_measurement_results/combined/central/normalised_xsection_patType1CorrectedPFMet.txt' + path_to_JSON = '{path}/{com}TeV/{variable}/xsection_measurement_results/muon/kv{kValue}/central/normalised_xsection_patType1CorrectedPFMet.txt' + filename = path_to_JSON.format(path = path, com = 8, + variable = variable, + kValue = kvalues[variable] + ) + + normalised_xsection_unfolded = read_data_from_JSON( filename ) + + xsections = normalised_xsection_unfolded['MADGRAPH'] + + # h_normalised_xsection_MADGRAPH = value_error_tuplelist_to_hist( normalised_xsection_unfolded['MADGRAPH'], bin_edges[variable] ) + + # file_template = '{path}/{category}/{name}_{channel}_{method}{suffix}.txt' + # filename = file_template.format( + # path = path_to_JSON, + # category = 'central', + # name = 'normalised_xsection', + # channel = 'combined', + # method = 'RooUnfoldSvd', + # suffix = '', + # ) + # normalised_xsection_unfolded = read_data_from_JSON( filename ) + + # # xsections = normalised_xsection_measured_unfolded['unfolded_with_systematics'] + # xsections = normalised_xsection_unfolded['powhegPythia8'] + + + edges = bin_edges[variable] + + + + print "# BEGIN YODA_SCATTER2D /CMS_2016_I1473674/d0%s-x01-y01" % number[variable] + print "Path=/CMS_2016_I1473674/d0%s-x01-y01" % number[variable] + print "Type=Scatter2D" + print "# xval xerr- xerr+ yval yerr- yerr+" + + for i_xsec in range ( 0, len( xsections ) ): + xsec = xsections[i_xsec][0] + xsec_error = xsections[i_xsec][1] + xlow = edges[i_xsec] + xup = edges[i_xsec+1] + xwidth = xup - xlow + xcentre = xlow + xwidth / 2 + + line = '{xcentre} {xerr_down} {xerr_up} {y} {yerr_down} {yerr_up}'.format( + xcentre = xcentre, + xerr_down = xwidth / 2, + xerr_up = xwidth / 2, + y = xsec, + yerr_down = xsec_error, + yerr_up = xsec_error + ) + + print line + print "# END YODA_SCATTER2D" + print "\n" diff --git a/src/cross_section_measurement/plot_covariance_matrix.py b/src/cross_section_measurement/plot_covariance_matrix.py new file mode 100644 index 00000000..a3602685 --- /dev/null +++ b/src/cross_section_measurement/plot_covariance_matrix.py @@ -0,0 +1,104 @@ +from optparse import OptionParser +from config import XSectionConfig +import numpy as np +from rootpy.plotting import Hist2D +from root_numpy import array2hist +from tools.Calculation import get_correlation_matrix + + + + + +if __name__ == '__main__': + parser = OptionParser() + 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='tables/', + help="set path to save tables") + parser.add_option("-v", "--variable", dest="variable", default='MET', + help="set variable to plot (MET, HT, ST, MT, WPT)") + parser.add_option("-m", "--metType", dest="metType", default='type1', + help="set MET type used in the analysis of MET-dependent variables") + parser.add_option("-c", "--centre-of-mass-energy", dest="CoM", default=8, type=int, + help="set the centre of mass energy for analysis. Default = 8 [TeV]") + + (options, args) = parser.parse_args() + measurement_config = XSectionConfig(options.CoM) + # caching of variables for shorter access + translate_options = measurement_config.translate_options + + + variable = options.variable + output_folder = options.output_folder + if not output_folder.endswith('/'): + output_folder += '/' + + met_type = translate_options[options.metType] + path_to_JSON = options.path + '/' + str(measurement_config.centre_of_mass_energy) + 'TeV/' + + for channel in ['combined']: + + stat_covariance = np.loadtxt(path_to_JSON + '/' + variable + '/xsection_measurement_results/combined/central/covariance.txt',delimiter=',') + syst_covariance = np.loadtxt(path_to_JSON + '/' + variable + '/xsection_measurement_results/combined/central/covariance_systematic.txt',delimiter=',') + + full_covariance = stat_covariance + syst_covariance + # full_covariance = syst_covariance + + print full_covariance + nBins = full_covariance.shape[0] + + uncertainties = [ np.sqrt(full_covariance[i,i]) for i in range(0,nBins) ] + + full_correlation = get_correlation_matrix( uncertainties, full_covariance ) + + hist = Hist2D(nBins, 0, nBins, nBins, 0, nBins, type='F') + array2hist(full_correlation, hist) + + hist.Draw('COLZ TEXT') + raw_input('...') + + +# 1. +# madgraph 260.200698419 14 2.22235586282e-47 +# madgraph_ptreweight 206.092857526 14 3.1199844713e-36 +# mcatnlo 452.430154807 14 1.09007959737e-87 +# powheg_v1_pythia 1136.08556935 14 9.45606334174e-234 +# powheg_v1_herwig 767.50906492 14 9.80253562557e-155 +# powheg_v2_pythia 172.213828364 14 2.4440527553e-29 +# powheg_v2_herwig 2595.64265928 14 0.0 + +# 0.995 +# madgraph 44.952092151 14 4.15303289684e-05 +# madgraph_ptreweight 8.97303169932 14 0.832774277818 +# mcatnlo 23.4724573835 14 0.0530032732625 +# powheg_v1_pythia 39.072729903 14 0.000355412330874 +# powheg_v1_herwig 25.1522914114 14 0.0330933635701 +# powheg_v2_pythia 20.3633234788 14 0.119096482135 +# powheg_v2_herwig 102.494382619 14 1.57464239767e-15 + +# 0.99 +# madgraph 30.9649954797 14 0.0056063902049 +# madgraph_ptreweight 6.52232220208 14 0.951563093841 +# mcatnlo 16.1987455777 14 0.301388411893 +# powheg_v1_pythia 29.7711155148 14 0.00820516299424 +# powheg_v1_herwig 18.6963603902 14 0.176878931589 +# powheg_v2_pythia 12.6196539729 14 0.556666472407 +# powheg_v2_herwig 65.2314532301 14 1.39097221014e-08 + +# 0.9 +# madgraph 10.7635419935 14 0.704503726933 +# madgraph_ptreweight 2.51432745404 14 0.999668848356 +# mcatnlo 5.15094707419 14 0.983599020582 +# powheg_v1_pythia 13.1771312971 14 0.512626005209 +# powheg_v1_herwig 7.06757553801 14 0.932075478849 +# powheg_v2_pythia 2.9027501205 14 0.999233557149 +# powheg_v2_herwig 21.0768834128 14 0.0996784514207 + +# 0 +# madgraph 9.67054408874 14 0.785852080078 +# madgraph_ptreweight 6.17360385996 14 0.961924159492 +# mcatnlo 9.50466821867 14 0.797427986823 +# powheg_v1_pythia 12.1628479675 14 0.593226598367 +# powheg_v1_herwig 5.64351092282 14 0.974692894781 +# powheg_v2_pythia 10.217000073 14 0.746152402623 +# powheg_v2_herwig 53.3771850831 14 1.64091512776e-06 diff --git a/tools/Calculation.py b/tools/Calculation.py index 1d3cd8dd..b1b8d35b 100644 --- a/tools/Calculation.py +++ b/tools/Calculation.py @@ -4,12 +4,14 @@ @author: kreczko ''' from __future__ import division -from uncertainties import ufloat +import uncertainties as u import numpy from math import sqrt from config.met_systematics import metsystematics_sources from rootpy import asrootpy +from rootpy.plotting import Hist2D from config.variable_binning import bin_edges +from decimal import * def calculate_xsection(inputs, luminosity, efficiency=1.): ''' @@ -25,7 +27,124 @@ def calculate_xsection(inputs, luminosity, efficiency=1.): add_result((xsection, xsection_error)) return result -def calculate_normalised_xsection(inputs, bin_widths, normalise_to_one=False): +def get_correlation_matrix(uncertainties,covariance): + nBins = len(uncertainties) + correlation = numpy.array( numpy.zeros( (nBins, nBins) ) ) + for i in range(0,nBins): + for j in range(0,nBins): + cov = covariance[i,j] + unc_i = uncertainties[i] + unc_j = uncertainties[j] + correlation[i,j] = cov / ( unc_i * unc_j) + return correlation + +def calculate_covariance_for_normalised_xsection(covariance, inputs, bin_widths,): + new_covariance = covariance.copy() + + values = [u.ufloat( i[0], i[1] ) for i in inputs] + normalisation = sum( values ) + + nominal_values = [v.nominal_value for v in values] + + values_correlated = u.correlated_values( nominal_values, covariance.tolist() ) + print 'Original values :',values_correlated + print 'Original correlation :',u.correlation_matrix(values_correlated) + norm = sum(values_correlated) + + norm_values_correlated = [] + + for v,width in zip( values_correlated, bin_widths ): + norm_values_correlated.append( v / width / norm ) + + print 'New values :',norm_values_correlated + print 'New correlation :',u.correlation_matrix(norm_values_correlated) + new_covariance = numpy.array( u.covariance_matrix(norm_values_correlated) ) + print 'New covariance :',u.covariance_matrix(norm_values_correlated) + + # n_rows = covariance.shape[0] + # n_cols = covariance.shape[1] + + # uncertainties = [numpy.sqrt(covariance[i,i]) for i in range(0,n_rows)] + + + # # for i in range(0,n_rows): + # # print numpy.sqrt(covariance[i,i]),values[i].std_dev + # correlation = get_correlation_matrix( uncertainties, covariance ) + # print 'Original correlation' + # print correlation + + # # print correlation + # for i_row in range(0,n_rows): + # for i_col in range(0,n_cols): + + # cor_ij = correlation[i_row,i_col] + + # xsec_i = ( values[i_row] / bin_widths[i_row] / normalisation ) + # xsec_j = ( values[i_col] / bin_widths[i_col] / normalisation ) + + # new_element = xsec_i.std_dev * xsec_j.std_dev * cor_ij + # # value_row = u.ufloat( values[i_row].nominal_value, numpy.sqrt( abs(covariance[i_row,i_col]) ) ) + # # value_column = u.ufloat( values[i_col].nominal_value, numpy.sqrt( abs(covariance[i_row,i_col]) ) ) + + # # xsec_i_xsec_j = u.ufloat( (values[i_row] * values[i_col]).nominal_value, abs( covariance[i_row,i_col] ) ) + # # value_i_value_j = value_row * value_column + + + + # # normalisation = u.ufloat(0,0) + # # for i in range(0,n_rows): + # # if i == n_rows: + # # # normalisation += value_row + # # elif i == n_cols: + # # # normalisation += value_column + # # else: + # # # normalisation += values[i] + + # # if i_row == i_col: + # # print i_row, i_col, new_element, numpy.sqrt(new_element) + # # print value_row,value_column,value_i_value_j + # # print normalisation + # # new_element = value_i_value_j / bin_widths[i_row] / bin_widths[i_col] / normalisation / normalisation + # # new_element = new_element.std_dev * numpy.sign( covariance[i_row,i_col] ) + + # # simple = covariance[i_row,i_col] / bin_widths[i_row] / bin_widths[i_col] / normalisation ** 2 + # # new_element = covariance[i_row, i_col] / values[i_row].nominal_value ** 2 + norm2Error / normalisation.nominal_value ** 2 + # # print covariance[i_row, i_col] / values[i_row].nominal_value ** 2 + # # print norm2Error + # # print normalisation.nominal_value + # # print normalisation.nominal_value ** 2 + # # print norm2Error / normalisation.nominal_value ** 2 + # # print 'Sqrt this :',new_element + # # new_element = numpy.sqrt( new_element ) + # # new_element *= values[i_row].nominal_value * values[i_col].nominal_value + # # new_element /= ( bin_widths[i_col] * bin_widths[i_row] ) + # # new_element = covariance[i_row,i_col] / bin_widths[i_row] / bin_widths[i_col] / normalisation.nominal_value ** 2 + # # / ( bin_widths[i_col] * bin_widths[i_row] * normalisation * normalisation ) + # # if i_row == i_col: + # # # # print values[i_row] + # # # # # print numpy.sqrt( values[i_row].nominal_value / values[i_row].std_dev ) ** 2 + (norm2Error / normalisation ** 2) * values[i_row].nominal_value / bin_widths[i_row] / normalisation + + # # # # print 'Original : ',covariance[i_row, i_col] + # # # # print bin_widths[i_col],bin_widths[i_row],normalisation + # # # print xsec_i_xsec_j + # # print values[i_row].nominal_value,values[i_row].std_dev + # # print (values[i_row] * values[i_col]).std_dev, covariance[i_row,i_col] + # # print values[i_row] * values[i_col], covariance[i_row,i_col] + + # # print 'New :',new_element, numpy.sqrt(new_element) + # # print 'Simple : ',simple, numpy.sqrt(simple.nominal_value) + # new_covariance[i_row, i_col] = new_element + + # new_uncertainties = [numpy.sqrt(new_covariance[i,i]) for i in range(0,n_rows)] + # print 'New uncertainties :',new_uncertainties + + # new_correlation = get_correlation_matrix( new_uncertainties, new_covariance ) + # print 'New correlation' + # print new_correlation + + return new_covariance + +def calculate_normalised_xsection(inputs, bin_widths, normalise_to_one=False, debug=False): """ Calculates normalised average x-section for each bin: 1/N *1/bin_width sigma_i There are two ways to calculate this @@ -35,14 +154,33 @@ def calculate_normalised_xsection(inputs, bin_widths, normalise_to_one=False): @param inputs: list of value-error pairs @param bin_widths: bin widths of the inputs """ - values = [ufloat( i[0], i[1] ) for i in inputs] + values = [u.ufloat( i[0], i[1] ) for i in inputs] normalisation = 0 if normalise_to_one: normalisation = sum( [value / bin_width for value, bin_width in zip( values, bin_widths )] ) else: normalisation = sum( values ) + xsections = [( 1 / bin_width ) * value / normalisation for value, bin_width in zip( values, bin_widths )] result = [(xsection.nominal_value, xsection.std_dev) for xsection in xsections] + + # if debug: + + # # for value, bin_width in zip( values, bin_widths ): + # # print value.nominal_value,value.std_dev + + # for xsec in xsections: + # print xsec.nominal_value, xsec.std_dev + # # print 'Here' + # # xsec = ( 1 / bin_width ) * value / normalisation.nominal_value + # # print value.nominal_value, bin_width, normalisation.nominal_value, xsec.nominal_value + # # print value.std_dev, bin_width, normalisation.std_dev, xsec.std_dev + + # # print ( 1 / bin_width ) * value / normalisation + # # print ( 1 / bin_width ) * value / normalisation.nominal_value + # # print value, bin_width, normalisation, xsec + # # print xsec.nominal_value,xsec.std_dev, numpy.sqrt( (value.std_dev / value.nominal_value) ** 2 + ( normalisation.std_dev / normalisation.nominal_value ) ** 2 ) * ( 1 / bin_width ) * value.nominal_value / normalisation.nominal_value + return result def decombine_result(combined_result, original_ratio): @@ -58,7 +196,7 @@ def decombine_result(combined_result, original_ratio): if original_ratio == 0: return combined_result, (0,0) else: - combined_result = ufloat(combined_result[0], combined_result[1]) + combined_result = u.ufloat(combined_result[0], combined_result[1]) sample_1 = combined_result * original_ratio / (1 + original_ratio) sample_2 = combined_result - sample_1 @@ -120,13 +258,47 @@ def calculate_lower_and_upper_PDFuncertainty(central_measurement, pdf_uncertaint negative.append(pdf_uncertainty) else: positive.append(pdf_uncertainty) - + # print 'Negative :',negative + # print 'Positive :',positive + # print 'Max sign :',numpy.sign( sum(max(x - central_measurement, y - central_measurement, 0) for x, y in zip(negative, positive)) ) + # print 'Min sign :',numpy.sign( sum(max(central_measurement - x, central_measurement - y, 0) for x, y in zip(negative, positive)) ) + dict_of_unc = {} + # for i in range(1,45): + # list_of_unc.append([]) + pdf_max = numpy.sqrt(sum(max(x - central_measurement, y - central_measurement, 0) ** 2 for x, y in zip(negative, positive))) pdf_min = numpy.sqrt(sum(max(central_measurement - x, central_measurement - y, 0) ** 2 for x, y in zip(negative, positive))) + + + for x,y,i in zip(negative, positive, range(0,44)): + maxUnc = 0. + if pdf_max > pdf_min: + maxUnc = max( x - central_measurement, y - central_measurement, 0 ) + else: + maxUnc = max( central_measurement - x, central_measurement - y, 0 ) + + # maxUnc = max( abs(pMax), abs(pMin) ) + + sign = 0. + if maxUnc == abs( x - central_measurement ): + sign = -1.0 + elif maxUnc == abs( y - central_measurement ): + sign = 1.0 + elif maxUnc == abs( central_measurement - x ): + sign = 1.0 + elif maxUnc == abs( central_measurement - y ): + sign = -1.0 + + maxUnc = maxUnc * sign + # print x, y, pMax, pMin, maxUnc + dict_of_unc[str(i)] = maxUnc + + # print dict_of_unc + # raw_input('...') - return pdf_min, pdf_max + return pdf_min, pdf_max, dict_of_unc -def calculate_lower_and_upper_systematics(central_measurement, list_of_systematics, symmetrise_errors = False): +def calculate_lower_and_upper_systematics(central_measurement, list_of_systematics, symmetrise_errors = False, debug=False, mass = False): ''' More generic replacement for calculateTotalUncertainty. Calculates the total negative and positve systematics. @param central_measurement: measurement from the central sample @@ -135,23 +307,450 @@ def calculate_lower_and_upper_systematics(central_measurement, list_of_systemati ''' negative_error = 0 positive_error = 0 - for systematic in list_of_systematics: + + positive_error_dictionary = {} + negative_error_dictionary = {} + + for name,systematic in list_of_systematics.iteritems(): deviation = abs(systematic) - abs(central_measurement) - + if debug: + print name, systematic, central_measurement, deviation, (systematic/central_measurement - 1)*100 + if deviation > 0: positive_error += deviation**2 + positive_error_dictionary[name] = deviation + negative_error_dictionary[name] = 0 else: negative_error += deviation**2 - + positive_error_dictionary[name] = 0 + negative_error_dictionary[name] = deviation + negative_error = sqrt(negative_error) positive_error = sqrt(positive_error) - + + # print negative_error_dictionary + # print positive_error_dictionary + dictionary_of_errors_to_use = {} + + if debug: + print 'In calculate' + print symmetrise_errors + print negative_error, positive_error + if symmetrise_errors: + if negative_error > positive_error: + dictionary_of_errors_to_use = negative_error_dictionary + else: + dictionary_of_errors_to_use = positive_error_dictionary + if debug: + print dictionary_of_errors_to_use + negative_error = max(negative_error, positive_error) positive_error = max(negative_error, positive_error) + else : + for source,value in negative_error_dictionary.iteritems(): + if value == 0: + dictionary_of_errors_to_use[source] = positive_error_dictionary[source] + else: + dictionary_of_errors_to_use[source] = negative_error_dictionary[source] + if debug : + print 'Debug :',negative_error,positive_error,dictionary_of_errors_to_use + # print negative_error_dictionary + # print positive_error_dictionary + return negative_error, positive_error, dictionary_of_errors_to_use + +def calculate_covariance_of_systematics_03(errors, central, mass_systematic=False, hadronisation=False, pdf=False, oneway=False, debug = False): + + all_systematic_labels = errors[0].keys() + # print all_systematic_labels + systematic_labels = [] + # print all_systematic_labels + for systematic in all_systematic_labels: + if 'patType1CorrectedPFMet' in systematic: + systematic = systematic.split('patType1CorrectedPFMet')[-1] + + if 'down' in systematic or '-' in systematic or ( 'min' in systematic and systematic != 'luminosity+') or 'Down' in systematic : + systematic_labels.append( systematic ) + elif not ( 'up' in systematic or '+' in systematic or 'max' in systematic or 'Up' in systematic ): + systematic_labels.append( systematic ) + + # Now have list of down systematics (or the only systematic e.g. QCD_shape) + if debug : + print 'In calculate covariance' + print systematic_labels + + # For each systematic, need a list of errors in each bin + errors_for_each_syst = {} + for syst in systematic_labels: + errors_for_each_syst[syst] = [] + + + nBins = len(errors) + for bin, c in zip(errors,central): + # for e in bin.values(): + # totalE2 += e*e + nSystematics = 0 + for systName, error in bin.iteritems(): + if systName in systematic_labels: + nSystematics += 1 + down_error = error + sign = numpy.sign( down_error ) + + # Get up variation + upSource = None + if 'down' in systName: + upSource = systName.replace('down', 'up') + elif '-' in systName: + upSource = systName.replace('-', '+') + elif 'min' in systName and systName != 'luminosity+': + upSource = systName.replace('min', 'max') + elif 'Down' in systName: + upSource = systName.replace('Down', 'Up') + + up_error = 0 + if upSource in bin.keys(): + nSystematics += 1 + up_error = bin[upSource] + if sign == 0: + sign = numpy.sign( up_error ) + + # if debug: + # print systName, upSource, down_error, up_error + + sign = 0 + if hadronisation or oneway: + sign = numpy.sign(down_error) + else: + if down_error == 0: + if up_error > 0: + sign = 1 + elif up_error < 0: + sign = -1 + elif up_error == 0: + if down_error > 0: + sign = -1 + elif down_error < 0: + sign = 1 + else: + sign = numpy.sign( up_error - down_error ) + # if debug: + # print 'Up down have same sign :',upSource,up_error,systName,down_error + # print up_error / c[0] * 100, down_error / c[0] * 100 + # print sign + # t = numpy.sqrt(up_error**2 + down_error**2) + # print t/c[0] * 100 + + if mass_systematic: + + sign = numpy.sign( up_error - down_error ) + # total_error = numpy.sqrt(up_error**2 + down_error**2) * sign + total_error = 0 + + if numpy.sign( up_error ) == numpy.sign( down_error ): + total_error = numpy.sqrt(up_error**2 + down_error**2) * sign + elif abs(up_error) > abs(down_error): + total_error = abs( up_error ) * sign + else: + total_error = abs( down_error ) * sign + + if debug: + print 'Mass',up_error,down_error,sign,total_error + + else: + total_error = numpy.sqrt(up_error**2 + down_error**2) * sign + + errors_for_each_syst[systName].append(total_error) + + if debug: + print errors_for_each_syst + + # Have 1D hist of errors for each systematic source + # Construct covariance matrix for each source + total_covariance_matrix = numpy.array( numpy.zeros((nBins,nBins )) ) + for source,e in errors_for_each_syst.iteritems(): + covariance_matrix = numpy.array( numpy.zeros((nBins,nBins )) ) + correlation_matrix = numpy.array( numpy.zeros((nBins,nBins )) ) + + maxError = max( [abs(error/c[0]) for error,c in zip( e, central ) ] ) + # print e + # print central + # for error, c in zip( e,central ): + # print error/c[0] + + + for i_row in range(0,nBins): + for i_col in range(0,nBins): + cor = 1.0 + # if pdf: + # cor = 0.0 + + # if 'JES' in source or 'JER' in source or 'BJet' in source: + # cor = 0.9 + + # if 'hadronisation' in source: + # cor = 1 + + # if 'TTJets_scaledown' in source: + # cor = 0.9 + + # if 'hadronisation' in source or 'VJets' in source or 'TTJets_scale' in source or 'TTJets_matching' in source or mass_systematic: + # if abs(e[i_row]) < central[i_row][1] or abs(e[i_col]) < central[i_col][1] : + # if debug: + # print central[i_row],central[i_col],e[i_row],e[i_col] + # cor = 0. + # else: + # if abs(e[i_row]) < central[i_row][1] * 0.2 or abs(e[i_col]) < central[i_col][1] * 0.2 : + # if debug: + # print central[i_row],central[i_col],e[i_row],e[i_col] + # cor = 0. + + # if 'TTJets_ptreweight' in source: + # cor = 0 + # if 'JES' in source or 'BJet' in source or 'PU' in source or 'SingleTop' in source or 'PDF' in source: + # cor = 0 + + # if maxError < 0.01: cor = 0. + + + if i_row == i_col: cor = 1. + + + + # if debug: + # print e[i_row],e[i_col] + + + + cov = e[i_row] * e[i_col] * cor + covariance_matrix[i_row, i_col] = cov + if cov != 0: + correlation_matrix[i_row, i_col] = cov / abs( e[i_row] * e[i_col] ) + else: + correlation_matrix[i_row, i_col] = 0 + + if debug : + print 'Correlation for',source + print correlation_matrix + raw_input('...') + total_covariance_matrix += covariance_matrix + + # for i_row in range(0,nBins): + # print i_row, total_covariance_matrix[i_row,i_row],sqrt(total_covariance_matrix[i_row,i_row]) + + + return total_covariance_matrix + +def calculate_covariance_of_systematics(all_categories, all_variations): + nBins = len(all_variations['central']) + total_covariance_matrix = numpy.array( numpy.zeros((nBins,nBins )) ) + + all_categories_errors = {} + for systematic in all_categories: + if 'patType1CorrectedPFMet' in systematic: + systematic = systematic.split('patType1CorrectedPFMet')[-1] + + if 'down' in systematic or '-' in systematic or ( 'min' in systematic and systematic != 'luminosity+') or 'Down' in systematic : + all_categories_errors[systematic] = [] + elif not ( 'up' in systematic or '+' in systematic or 'max' in systematic or 'Up' in systematic ): + all_categories_errors[systematic] = [] + + errors = [] + for bin in range(0,nBins): + + # print all_categories + # print all_variations.keys() + error_for_bin = 0 + negative_error = 0 + positive_error = 0 + + sources_for_this_bin = [] + negative_sources = [] + positive_sources = [] + + for systematic in all_categories: + if 'patType1CorrectedPFMet' in systematic: + systematic = systematic.split('patType1CorrectedPFMet')[-1] + + deviation = all_variations[systematic][bin] + if deviation > 0: + positive_error += deviation**2 + positive_sources.append( systematic ) + else: + negative_error += deviation**2 + negative_sources.append( systematic ) + + negative_error = sqrt(negative_error) + positive_error = sqrt(positive_error) + + + if negative_error > positive_error: + sources_for_this_bin = negative_sources + error_for_bin = negative_error + else : + sources_for_this_bin = positive_sources + error_for_bin = positive_error + + print sources_for_this_bin + for source in sources_for_this_bin: + if 'down' in systematic or '-' in systematic or ( 'min' in systematic and systematic != 'luminosity+') or 'Down' in systematic : + # Check to see if 'up' version is also present + upSource = None + if 'down' in systematic: + upSource = source.replace('down', 'up') + elif '-' in systematic: + upSource = source.replace('-', '+') + elif 'min' in systematic and systematic != 'luminosity+': + upSource = source.replace('min', 'max') + elif 'Down' in systematic: + upSource = source.replace('Down', 'Up') + + if upSource in sources_for_this_bin: + print 'Up and down in here!',source,upSource + + print source + # all_categories_errors[source].append(all_variations[source][bin]) + + print sources_for_this_bin + errors.append(error_for_bin) + + for i in range(0,nBins): + total_covariance_matrix[i,i] = errors[i] * errors[i] + + return total_covariance_matrix + +def calculate_covariance_of_systematics_properly(dictionary_of_systematics, central_measurement, debug=False): + # print all_categories + # print all_variations.keys() + if debug: + print central_measurement + print dictionary_of_systematics + all_categories = dictionary_of_systematics.keys() + + all_categories_errors = {} + for systematic in all_categories: + if 'patType1CorrectedPFMet' in systematic: + systematic = systematic.split('patType1CorrectedPFMet')[-1] + + if 'down' in systematic or '-' in systematic or ( 'min' in systematic and systematic != 'luminosity+') or 'Down' in systematic : + all_categories_errors[systematic] = [] + elif not ( 'up' in systematic or '+' in systematic or 'max' in systematic or 'Up' in systematic ): + all_categories_errors[systematic] = [] + + nBins = len(dictionary_of_systematics[all_categories[0]]) + + total_covariance_matrix = numpy.array( numpy.zeros((nBins,nBins )) ) + + for down_source in all_categories_errors: + max_variation = [] + + down_variation = dictionary_of_systematics[down_source] + + up_variation = down_variation + if 'down' in down_source: + up_variation = dictionary_of_systematics[down_source.replace('down','up')] + elif '-' in down_source: + up_variation = dictionary_of_systematics[down_source.replace('-','+')] + elif 'min' in systematic and systematic != 'luminosity+': + up_variation = dictionary_of_systematics[down_source.replace('min','max')] + elif 'Down' in down_source: + up_variation = dictionary_of_systematics[down_source.replace('Down','Up')] + + # print down_source + + for u,d,c in zip( up_variation, down_variation, central_measurement): + up_deviation = abs(u[0]) - abs(c[0]) + down_deviation = abs(d[0]) - abs(c[0]) + if debug: + print 'Differences' + print c[0],u[0],d[0] + print up_deviation,down_deviation + m = max( abs(up_deviation), abs(down_deviation) ) + sign = 0 + # if m == abs(up_deviation): + # sign = numpy.sign(up_deviation) + # elif m == abs(down_deviation): + # sign = numpy.sign(down_deviation) + # if debug: + # print 'Signs :',sign,numpy.sign( up_deviation - down_deviation ) + sign = numpy.sign( up_deviation - down_deviation ) + max_variation.append( m * sign ) + + covariance_matrix = numpy.array( numpy.zeros((nBins,nBins )) ) + + for i_row in range(0,nBins): + for i_col in range(0,nBins): + covariance = max_variation[i_row] * max_variation[i_col] + covariance_matrix[i_row,i_col] = covariance + + if debug: + print down_source + print max_variation + uncertainties = [numpy.sqrt(covariance_matrix[i,i]) for i in range(0,nBins)] + cor = get_correlation_matrix( uncertainties, covariance_matrix) + print cor + raw_input('...') + total_covariance_matrix += covariance_matrix + + if debug: + uncertainties = [numpy.sqrt(total_covariance_matrix[i,i]) for i in range(0,nBins)] + cor = get_correlation_matrix( uncertainties, total_covariance_matrix) + print 'Total correlation matrix' + print cor + raw_input('...') + # for category in all_categories: + # if 'down' in category or '-' in category or 'min' in category: + return total_covariance_matrix + + +def calculate_lower_and_upper_systematics_properly(central_measurement, dictionary_of_systematics): + + all_systematics = dictionary_of_systematics.keys() + systematic_categories = [] + print all_systematics + for systematic in all_systematics: + if 'down' in systematic or '-' in systematic : + systematic_categories.append( systematic ) + elif not ( 'up' in systematic or '+' in systematic ): + systematic_categories.append( systematic ) + print systematic_categories + negative_error = 0 + positive_error = 0 + + for category in systematic_categories: + down_variation = dictionary_of_systematics[category] - central_measurement + + up_variation = down_variation + if 'down' in category: + up_variation = dictionary_of_systematics[category.replace('down','up')] - central_measurement + elif '-' in category: + up_variation = dictionary_of_systematics[category.replace('-','+')] - central_measurement + + max_variation = max( abs(up_variation), abs(down_variation) ) + sign = 0 + if max_variation == abs(up_variation): + sign = numpy.sign(up_variation) + elif max_variation == abs(down_variation): + sign = numpy.sign(down_variation) + + print category,down_variation,up_variation, max_variation, sign + + # for systematic in list_of_systematics: + # deviation = abs(systematic) - abs(central_measurement) + + # if deviation > 0: + # positive_error += deviation**2 + # else: + # negative_error += deviation**2 + + # negative_error = sqrt(negative_error) + # positive_error = sqrt(positive_error) - return negative_error, positive_error + # if symmetrise_errors: + # negative_error = max(negative_error, positive_error) + # positive_error = max(negative_error, positive_error) + return negative_error, positive_error + def combine_errors_in_quadrature(list_of_errors): list_of_errors_squared = [error**2 for error in list_of_errors] sum_of_errors_squared = sum(list_of_errors_squared) @@ -286,3 +885,4 @@ def which_variable_bin(variable, value): else: break return variable_bin + diff --git a/tools/HistSet.py b/tools/HistSet.py index 57bda69d..2eb757db 100644 --- a/tools/HistSet.py +++ b/tools/HistSet.py @@ -69,7 +69,7 @@ def plot( self, plot_options = {} ): colours = colours, histogram_properties = histogram_properties, fill_area = fill_area, - make_ratio = True, + make_ratio = False, alpha = alpha, save_folder = output_folder, save_as = output_format, diff --git a/tools/Unfolding.py b/tools/Unfolding.py index fb4ad5c7..8a9d4f90 100644 --- a/tools/Unfolding.py +++ b/tools/Unfolding.py @@ -261,9 +261,9 @@ def get_unfold_histogram_tuple( h_measured.Scale( lumiweight ) h_response.Scale( lumiweight ) - h_truth, h_measured, h_response = [ fix_overflow( hist ) for hist in [h_truth, h_measured, h_response] ] - if load_fakes: - h_fakes = fix_overflow( h_fakes ) + # h_truth, h_measured, h_response = [ fix_overflow( hist ) for hist in [h_truth, h_measured, h_response] ] + # if load_fakes: + # h_fakes = fix_overflow( h_fakes ) return h_truth, h_measured, h_response, h_fakes diff --git a/tools/latex.py b/tools/latex.py index d794f883..57db084f 100644 --- a/tools/latex.py +++ b/tools/latex.py @@ -10,7 +10,6 @@ import os from distutils.spawn import find_executable - def setup_matplotlib(): ''' Seup matplotlib with all the latex fancyness we have @@ -18,6 +17,7 @@ def setup_matplotlib(): rc( 'font', **CMS.font ) rc( 'text', usetex = True ) rcParams['text.latex.preamble'] = compile_package_list() + rcParams['mathtext.default'] = 'regular' def compile_package_list(): ''' diff --git a/tools/plotting.py b/tools/plotting.py index 0c1e169e..f4815fee 100644 --- a/tools/plotting.py +++ b/tools/plotting.py @@ -16,7 +16,7 @@ from copy import deepcopy import matplotlib.gridspec as gridspec from matplotlib.ticker import MultipleLocator, FixedLocator -from itertools import cycle +from itertools import cycle, combinations from tools.latex import setup_matplotlib setup_matplotlib() @@ -47,7 +47,7 @@ class Histogram_properties: #If True (the default) then plot bins with zero content otherwise only # show bins with nonzero content. emptybins = False - + formats = ['pdf'] def __init__( self, dictionary = {} ): for name, value in dictionary.iteritems(): @@ -74,11 +74,90 @@ class PlotConfig: ''' general_options = ['files', 'histograms', 'labels', 'plot_type', 'output_folder', 'output_format', 'command', 'data_index', 'normalise', - 'show_ratio', 'show_stat_errors_on_mc', 'colours', 'name_prefix'] + 'show_ratio', 'show_stat_errors_on_mc', 'colours', 'name_prefix', 'subtract-files'] def __init__( self, config_file, **kwargs ): self.config_file = config_file +class Plot(object): + ''' + A class to define a plot + ''' + delegate_attr = attribute_names=['name', 'title', 'formats'] + def __init__(self, histograms, properties, **kwargs): + self.__draw_method = 'errorbar' + self.__properties = properties + self._path = properties.path + self.__histograms = histograms + if self._path != '' and not self._path.endswith('/'): + self._path += '/' + self.__errorbands = [] + if kwargs.has_key('errorbands'): + self.__errorbands = kwargs.pop('errorbands') + + def add_error_band(self, errorband): + self.__errorbands.append(errorband) + + @property + def errorbands(self): + return self.__errorbands + + @property + def properties(self): + return self.__properties + + def save(self): + make_folder_if_not_exists(self._path) + for f in self.__properties.formats: + file_name = '{path}{name}.{format}' + file_name = file_name.format( + path = self._path, + name = self.__properties.name, + format = f, + ) + plt.savefig(file_name) + + def cleanup(self): + # clear current figure and axes + plt.clf() + plt.cla() + plt.close() + + def __enter__(self): + self.cleanup() + return self + + def __exit__(self, exc_type, exc_value, traceback): + self.cleanup() + + @property + def draw_method(self): + return self.__draw_method + + @draw_method.setter + def draw_method(self, method): + if method in rplt.__dict__.keys(): + self.__draw_method = method + else: + raise ValueError('Invalid draw method') + + @property + def histograms(self): + return self.__histograms + + @property + def show_ratios(self): + return self.__properties.has_ratio and len(self.__histograms) > 1 + + @property + def fix_to_zero(self): + return self.__properties.fix_to_zero + +# def __getattr__(self, name): +# print name +# if name in Plot.delegate_attr: +# return getattr(self.__properties, name) + def make_data_mc_comparison_plot( histograms = [], histogram_lables = [], histogram_colors = [], @@ -94,6 +173,7 @@ def make_data_mc_comparison_plot( histograms = [], save_folder = check_save_folder(save_folder) # make copies in order not to mess with existing histograms histograms_ = deepcopy(histograms) + stack = HistStack() add_mc = stack.Add for index, histogram in enumerate( histograms_ ): @@ -122,7 +202,10 @@ def make_data_mc_comparison_plot( histograms = [], plt.figure( figsize = CMS.figsize, dpi = CMS.dpi, facecolor = CMS.facecolor ) if show_ratio: ratio = data.Clone( 'ratio' ) - ratio.Divide( sum( stack.GetHists() ) ) + mcHist = sum( stack.GetHists() ) + for bin_i in range(1,mcHist.GetNbinsX()): + mcHist.SetBinError( bin_i, 0 ) + ratio.Divide( mcHist ) ratio.SetMarkerSize( 3 ) gs = gridspec.GridSpec( 2, 1, height_ratios = [5, 1] ) axes = plt.subplot( gs[0] ) @@ -134,14 +217,64 @@ def make_data_mc_comparison_plot( histograms = [], axes.set_ylim( ymin = 1e-2 ) mc_error = histogram_properties.mc_error - if mc_error > 0: + if isinstance( mc_error, list ): + stack_lower = sum( stack.GetHists() ) + stack_upper = stack_lower.Clone( 'upper' ) + + for bin_i in range( 1, stack_lower.GetNbinsX() ): + stack_lower.SetBinContent( bin_i, stack_lower.GetBinContent( bin_i ) - mc_error[bin_i - 1] ) + stack_upper.SetBinContent( bin_i, stack_upper.GetBinContent( bin_i ) + mc_error[bin_i - 1] ) + # stack_lower -= mc_error + # stack_upper += mc_error + # print 'Stack : ', list(stack.y()) + # print 'Stack upper : ',list(stack_upper.y()) + # print 'Stack lower : ',list(stack_lower.y()) + +# plt.fill(x,np.sin(x),color='blue',alpha=0.5) +# plt.fill(x,np.sin(x),color='None',alpha=0.5,edgecolor='blue',hatch='/') + + # rplt.fill_between( stack_upper, + # stack_lower, axes, + # color = 'black', + # # edgecolor='black', + # # hatch='//', + # linewidth='0.0', + # # facecolor = '0.75', + # # alpha = 0., + # zorder = len(histograms_) + 1 + # ) + + # rplt.fill_between( stack_upper, + # stack_lower, axes, + # # color = 'black', + # edgecolor='black', + # hatch='//', + # linewidth='0.0', + # # facecolor = '0.75', + # # alpha = 1., + # zorder = len(histograms_) + 1 + # ) + + rplt.fill_between( stack_upper, stack_lower, axes, + # facecolor = '0.75', + facecolor = 'none', + edgecolor='black', + hatch='//', + linewidth='0.0', + # alpha = '0.0000001' + # facecolor = '0.75', + # alpha = 0.1, + # alpha = '1.0', + zorder = len(histograms_) + 1 + ) + elif mc_error > 0: stack_lower = sum( stack.GetHists() ) stack_upper = stack_lower.Clone( 'upper' ) stack_lower.Scale( 1 - mc_error ) stack_upper.Scale( 1 + mc_error ) rplt.fill_between( stack_upper, stack_lower, axes, facecolor = '0.75', - alpha = 0.5, hatch = '/', + alpha = 0.5, zorder = len(histograms_) + 1 ) if not mc_error > 0 and show_stat_errors_on_mc: stack_lower = sum( stack.GetHists() ) @@ -150,8 +283,9 @@ def make_data_mc_comparison_plot( histograms = [], for bin_i in range( 1, stack_lower.GetNbinsX() ): stack_lower.SetBinContent( bin_i, stack_lower.GetBinContent( bin_i ) - mc_errors[bin_i - 1] ) stack_upper.SetBinContent( bin_i, stack_upper.GetBinContent( bin_i ) + mc_errors[bin_i - 1] ) + rplt.fill_between( stack_upper, stack_lower, axes, facecolor = '0.75', - alpha = 0.5, hatch = '/', + alpha = 0.5, zorder = len(histograms_) + 1 ) # a comment on zorder: the MC stack should be always at the very back (z = 1), @@ -165,14 +299,20 @@ def make_data_mc_comparison_plot( histograms = [], # put legend into the correct order (data is always first!) handles, labels = axes.get_legend_handles_labels() - data_label_index = labels.index( 'data' ) + data_label_index = labels.index( 'Data' ) data_handle = handles[data_label_index] - labels.remove( 'data' ) + labels.remove( 'Data' ) handles.remove( data_handle ) - labels.insert( 0, 'data' ) + labels.insert( 0, 'Data' ) handles.insert( 0, data_handle ) if mc_error > 0 or ( not mc_error > 0 and show_stat_errors_on_mc ): - p1 = Rectangle( ( 0, 0 ), 1, 1, fc = "0.75", alpha = 0.5, hatch = '/' ) + p1 = Rectangle( ( 0, 0 ), 1, 1, + color = 'none', + edgecolor='black', + hatch='//', + linewidth='0.0', + # fc = "0.75", alpha = 0.5, + ) handles.append( p1 ) labels.append( histogram_properties.mc_errors_label ) @@ -188,16 +328,18 @@ def make_data_mc_comparison_plot( histograms = [], x_limits = histogram_properties.x_limits y_limits = histogram_properties.y_limits - if len( x_limits ) == 2: - axes.set_xlim( xmin = x_limits[0], xmax = x_limits[1] ) if len( y_limits ) == 2: axes.set_ylim( ymin = y_limits[0], ymax = y_limits[1] ) else: y_max = get_best_max_y(histograms_, x_limits=x_limits) * histogram_properties.y_max_scale - axes.set_ylim( ymin = 0, ymax = y_max ) - if histogram_properties.set_log_y: - if not len( y_limits ) == 2: # if not user set y-limits, set default - axes.set_ylim( ymin = 1e-1 ) + y_limits = [ 0, y_max ] + if histogram_properties.set_log_y: + y_limits = [ 0.1, y_max * 10 ] + + axes.set_ylim( ymin = y_limits[0], ymax = y_limits[1] ) + + if len( x_limits ) == 2: + axes.set_xlim( xmin = x_limits[0], xmax = x_limits[1] ) #draw a red vertical line if needed: if draw_vertical_line != 0: @@ -211,22 +353,52 @@ def make_data_mc_comparison_plot( histograms = [], # Add horizontal line at y=1 on ratio plot ax1.axhline(y=1, linewidth = 1) set_labels( plt, histogram_properties, show_x_label = True, show_title = False ) - plt.ylabel( r'$\frac{\mathrm{data}}{\mathrm{pred.}}$', CMS.y_axis_title ) + # plt.ylabel( r'$\frac{\mathrm{data}}{\mathrm{pred.}}$', CMS.y_axis_title ) + plt.ylabel( '$\\displaystyle\\frac{\\mathrm{data}}{\\mathrm{pred.}}$', CMS.y_axis_title ) + ax1.yaxis.set_label_coords(-0.115, 0.8) rplt.errorbar( ratio, emptybins = histogram_properties.emptybins, axes = ax1, xerr = histogram_properties.xerr, ) if len( x_limits ) == 2: ax1.set_xlim( xmin = x_limits[0], xmax = x_limits[1] ) if len( histogram_properties.ratio_y_limits ) == 2: - ax1.set_ylim( ymin = histogram_properties.ratio_y_limits[0], - ymax = histogram_properties.ratio_y_limits[1] ) + ax1.set_ylim( ymin = histogram_properties.ratio_y_limits[0]-0.02, + ymax = histogram_properties.ratio_y_limits[1]+0.02 ) + ax1.xaxis.labelpad = round(ax1.xaxis.labelpad * 3,0) # dynamic tick placement - adjust_ratio_ticks(ax1.yaxis, n_ticks = 3) + adjust_ratio_ticks(ax1.yaxis, n_ticks = 3, y_limits = histogram_properties.ratio_y_limits) + + if isinstance( mc_error, list ): + ratio_lower = ratio.Clone('lower') + ratio_upper = ratio.Clone('upper') + mc_hist = sum( stack.GetHists() ) + for bin_i in range( 1, ratio_lower.GetNbinsX() ): + if ratio_lower.GetBinContent( bin_i ) != 0 : + ratio_lower.SetBinContent( bin_i, 1 - mc_error[bin_i - 1] / mc_hist.GetBinContent( bin_i ) ) + else : + ratio_lower.SetBinContent(bin_i, 1) + + if ratio_upper.GetBinContent( bin_i ) != 0 : + ratio_upper.SetBinContent( bin_i, 1 + mc_error[bin_i - 1] / mc_hist.GetBinContent( bin_i ) ) + else : + ratio_upper.SetBinContent(bin_i, 1) + + rplt.fill_between( ratio_upper, ratio_lower, axes, + # facecolor = '0.75', + color = 'none', + edgecolor='black', + hatch='//', + linewidth='0.0', + # facecolor = '0.75', + # alpha = 0.5, + # alpha = 0.5, + # zorder = len(histograms_) + 1 + ) if CMS.tight_layout: plt.tight_layout() - + for save in save_as: plt.savefig( save_folder + histogram_properties.name + '.' + save ) @@ -350,7 +522,7 @@ def make_shape_comparison_plot( shapes = [], shape.fillstyle = 'solid' else: shape.linewidth = 5 - + print shapes_ if not histogram_properties.y_limits: histogram_properties.y_limits = [0, get_best_max_y(shapes_, False)] # plot with matplotlib @@ -385,10 +557,14 @@ def make_shape_comparison_plot( shapes = [], ncol = histogram_properties.legend_columns ) l1.set_zorder(102) #add error bars - graphs = spread_x(shapes_, list(shapes_[0].xedges())) - for graph in graphs: - rplt.errorbar( graph, axes = axes ) + # graphs = spread_x(shapes_, list(shapes_[0].xedges())) + graphs = shapes_ + for graph in graphs: + graph.elinewidth = 2 + (_, caps, _) =rplt.errorbar( graph, axes = axes, xerr=0, elinewidth = 2, capsize=5 ) + for c in caps: + c.set_markeredgewidth(3) adjust_axis_limits(axes, histogram_properties, shapes_) if make_ratio: plt.setp( axes.get_xticklabels(), visible = False ) @@ -413,7 +589,7 @@ def make_shape_comparison_plot( shapes = [], ax1.set_ylim( ymin = histogram_properties.ratio_y_limits[0], ymax = histogram_properties.ratio_y_limits[1] ) # dynamic tick placement - adjust_ratio_ticks(ax1.yaxis, n_ticks = 3) + adjust_ratio_ticks(ax1.yaxis, n_ticks = 3, y_limits = histogram_properties.ratio_y_limits) if CMS.tight_layout: plt.tight_layout() @@ -480,7 +656,7 @@ def make_plot( histogram, histogram_label, histogram_properties = Histogram_prop if CMS.tight_layout: plt.tight_layout() - + for save in save_as: plt.savefig( save_folder + histogram_properties.name + '.' + save ) plt.close() @@ -578,6 +754,7 @@ def set_labels( plt, histogram_properties, show_x_label = True, if show_x_label: plt.xlabel( histogram_properties.x_axis_title, CMS.x_axis_title ) plt.ylabel( histogram_properties.y_axis_title, CMS.y_axis_title ) + plt.tick_params( **CMS.axis_label_major ) plt.tick_params( **CMS.axis_label_minor ) if show_title: @@ -585,16 +762,20 @@ def set_labels( plt, histogram_properties, show_x_label = True, if not axes: return + # axes.xaxis.labelpad = 20 + # CMS text # note: fontweight/weight does not change anything as we use Latex text!!! - logo_location = (0.05, 0.98) + logo_location = (0.05, 0.96) prelim_location = (0.05, 0.92) - additional_location = (0.95, 0.98) + # additional_location = (0.94, 0.98) + additional_location = (0.75, 0.86) + loc = histogram_properties.cms_logo_location if loc == 'right': - logo_location = (0.95, 0.98) + logo_location = (0.95, 0.96) prelim_location = (0.95, 0.92) - additional_location = (0.95, 0.86) + # additional_location = (0.75, 0.86) plt.text(logo_location[0], logo_location[1], r"\textbf{CMS}", transform=axes.transAxes, fontsize=42, @@ -605,7 +786,7 @@ def set_labels( plt, histogram_properties, show_x_label = True, verticalalignment='top',horizontalalignment=loc) # channel text axes.text(additional_location[0], additional_location[1], - r"\emph{%s}" %histogram_properties.additional_text, + r"%s" %histogram_properties.additional_text, transform=axes.transAxes, fontsize=40, verticalalignment='top', horizontalalignment='right') @@ -659,16 +840,106 @@ def check_save_folder(save_folder): return save_folder -def adjust_ratio_ticks( axis, n_ticks = 3 ): +def adjust_ratio_ticks( axis, n_ticks = 3, y_limits = None ): # dynamic tick placement ticks = axis.get_ticklocs() tick_min, tick_max = ticks[0], ticks[-1] + + # Check if these are outside of the y_limits. Use those instead if so. + if y_limits != None: + tick_min = y_limits[0] + tick_max = y_limits[-1] + # limit to 3 ticks tick_distance = abs( tick_max - tick_min ) / ( n_ticks + 1 ) includes_one = tick_max > 1 and tick_min < 1 if includes_one: - axis.set_major_locator( FixedLocator( [tick_min + tick_distance/2, 1, tick_max - tick_distance/2] ) ) + if y_limits != None: + axis.set_major_locator( FixedLocator( [round(tick_min,1), 1, round(tick_max,1)] ) ) + else : + axis.set_major_locator( FixedLocator( [round(tick_min + tick_distance/2,1), 1, round(tick_max - tick_distance/2,1)] ) ) else: axis.set_major_locator( MultipleLocator( tick_distance ) ) axis.set_minor_locator( MultipleLocator( tick_distance / 2 ) ) +def compare_histograms(plot): + histograms = plot.histograms + properties = plot.properties + # the cycling needs to be better + colors = ['green', 'red', 'blue', 'magenta'] + colorcycler = cycle( colors ) + markers = [20, 23, 22, 33, 21, 29] + markercycler = cycle( markers ) + + plt.figure( figsize = CMS.figsize, dpi = CMS.dpi, facecolor = CMS.facecolor ) + axes = plt.axes() + if plot.show_ratios: + gs = gridspec.GridSpec( 2, 1, height_ratios = [5, 1] ) + axes = plt.subplot( gs[0] ) + plot_function = rplt.__dict__[plot.draw_method] + for label, histogram in histograms.items(): + histogram.markersize = 2 + histogram.markerstyle = next( markercycler ) + histogram.color = next( colorcycler ) + plot_function( histogram, axes = axes, label = label, + emptybins = properties.emptybins, + yerr = True, + xerr = properties.xerr, elinewidth = 2 ) + + set_labels( plt, properties, show_x_label = not plot.show_ratios, axes = axes ) + errorbands = plot.errorbands + handles, labels = axes.get_legend_handles_labels() + for band in errorbands: + band.draw(axes) + p1 = Rectangle( ( 0, 0 ), 1, 1, fc = "0.75", alpha = 0.5, hatch = '/' , + label = band.name) + handles.append( p1 ) + labels.append( band.name ) + adjust_axis_limits( axes, properties, histograms.values() ) + + # or sort legend by labels + import operator + hl = sorted(zip(handles, labels), + key=operator.itemgetter(1)) + handles2, labels2 = zip(*hl) + l1 = axes.legend( + handles2, labels2, numpoints = 1, + frameon = properties.legend_color, + bbox_to_anchor = properties.legend_location, + bbox_transform=plt.gcf().transFigure, + prop = CMS.legend_properties, + ncol = properties.legend_columns, + ) + l1.set_zorder(102) + + ratios = {} + if plot.show_ratios: + for (l1, l2) in combinations(histograms.keys(), 2): + label = '{0}/{1}'.format(l1, l2) + h = histograms[l1].clone() + h.Divide(histograms[l2]) + h.SetName(label) + ratios[label] = h + plt.setp( axes.get_xticklabels(), visible = False ) + axes_ratio = plt.subplot( gs[1] ) + axes_ratio.set_ylim(properties.ratio_y_limits[0],properties.ratio_y_limits[-1]) + axes_ratio.minorticks_on() + axes_ratio.grid( True, 'major', linewidth = 1) + axes_ratio.axhline(y=1, linewidth = 1, linestyle = 'dashed', color = 'black') + set_labels( plt, properties, show_x_label = True, show_title = False ) + plt.ylabel('ratio') + axes_ratio.yaxis.set_label_coords(-0.115, 0.8) + for label, ratio in ratios.items(): + plot_function( ratio, xerr = properties.xerr, + emptybins = properties.emptybins, + axes = axes_ratio, elinewidth = 2 ) + + adjust_ratio_ticks(axes_ratio.yaxis, n_ticks = 3, y_limits = properties.ratio_y_limits ) + axes_ratio.set_ylim(properties.ratio_y_limits[0],properties.ratio_y_limits[-1]) + # adjust_axis_limits( axes_ratio, properties, ratios.values() ) + axes_ratio.set_xlim(properties.x_limits[0], properties.x_limits[-1]) + if CMS.tight_layout: + plt.tight_layout() + + plot.save() + plt.close() \ No newline at end of file