From 4e373688ec369956c1e860fffde08658006bc6b8 Mon Sep 17 00:00:00 2001 From: Emyr Clement Date: Tue, 3 Oct 2017 17:07:09 +0100 Subject: [PATCH 1/3] Zoom in on lepton eta systematic plots. --- dps/analysis/xsection/05_make_systematic_plots.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/dps/analysis/xsection/05_make_systematic_plots.py b/dps/analysis/xsection/05_make_systematic_plots.py index 0de7bdf0..de8472da 100644 --- a/dps/analysis/xsection/05_make_systematic_plots.py +++ b/dps/analysis/xsection/05_make_systematic_plots.py @@ -28,6 +28,8 @@ def plot_systematic_uncertainties(systematic_uncertainties, bin_edges, variable, x_limits = [bin_edges[0], bin_edges[-1]] # y_limits = [-0.6,0.6] y_limits = [0,0.4] + if 'abs_lepton_eta' in variable and 'normalised' in output_folder: + y_limits = [0,0.05] fig_syst = plt.figure( figsize = ( 20, 16 ), dpi = 400, facecolor = 'white' ) ax_syst = fig_syst.add_subplot(1, 1, 1) From 684bdfa27682aa83f3b71f8f5b29bd720301a122 Mon Sep 17 00:00:00 2001 From: Emyr Clement Date: Tue, 3 Oct 2017 17:12:30 +0100 Subject: [PATCH 2/3] Calculate uncertainties and covariance by averaging up/down variations. --- dps/utils/systematic.py | 163 ++++++++++++++++++++++++++++++++-------- 1 file changed, 131 insertions(+), 32 deletions(-) diff --git a/dps/utils/systematic.py b/dps/utils/systematic.py index 985ebbc2..b2d07eb1 100644 --- a/dps/utils/systematic.py +++ b/dps/utils/systematic.py @@ -14,6 +14,7 @@ from dps.config.latex_labels import variables_latex, variables_NonLatex from dps.config.xsection import XSectionConfig from copy import deepcopy +from itertools import product measurement_config = XSectionConfig( 13 ) template = '%.1f fb$^{-1}$ (%d TeV)' title = template % ( measurement_config.new_luminosity/1000, measurement_config.centre_of_mass_energy) @@ -381,6 +382,7 @@ def get_scale_envelope(d_scale_syst, central): scale.append(down[index][i]) down['envelope_down']=scale + up_diff = up.subtract(central, axis='index') up_diff = up_diff.abs() up_diff['max'] = up_diff.max(axis = 1) @@ -459,7 +461,9 @@ def get_symmetrised_systematic_uncertainty( options, syst_unc_x_secs ): Combine PDFs and alphaS systematics ''' xsections_with_symmetrised_systematics = deepcopy(syst_unc_x_secs) + xsections_with_uncertainties_and_covariance = {} central_measurement = syst_unc_x_secs['central'] + d_Cov_Cor = {} for systematic, variation in syst_unc_x_secs.iteritems(): if (systematic in ['PDF', 'CT14', 'MMHT14']): # Replace all PDF weights with full PDF combination @@ -471,58 +475,110 @@ def get_symmetrised_systematic_uncertainty( options, syst_unc_x_secs ): xsections_with_symmetrised_systematics[systematic] = [ pdf_sym, pdf_sign - ] + ] elif systematic == 'central': xsections_with_symmetrised_systematics['central'] = central_measurement + xsections_with_uncertainties_and_covariance['central'] = central_measurement else: lower_measurement = variation[0] upper_measurement = variation[1] isTopMassSystematic = True if systematic == 'TTJets_mass' else False - symmetrised_uncertainties, signed_uncertainties = get_symmetrised_errors( - central_measurement, - upper_measurement, - lower_measurement, - options, + # The commented code in the remainder of this function + # is used to calcualte the uncertainties and covariance + # by taking the maximum uncertainty in each bin + # symmetrised_uncertainties, signed_uncertainties = get_symmetrised_errors( + # central_measurement, + # upper_measurement, + # lower_measurement, + # options, + # isTopMassSystematic, + # systematic + # ) + + # xsections_with_symmetrised_systematics[systematic] = [ + # symmetrised_uncertainties, + # signed_uncertainties, + # ] + + # This code replace the previous lines (similar for the rest of this function) + # and calculates the uncertainty and covariance by averaging the up/down + # variations of each source + # If only one variation is available (e.g. a colour reconnection sample) + # then it will give identical results to the previous method + symmetrised_uncertainties, covariance_matrix = get_average_covariance_matrix( + central_measurement, + upper_measurement, + lower_measurement, + options, isTopMassSystematic, - ) + systematic + ) - xsections_with_symmetrised_systematics[systematic] = [ + xsections_with_uncertainties_and_covariance[systematic] = [ symmetrised_uncertainties, - signed_uncertainties, - ] + covariance_matrix, + ] + d_Cov_Cor[systematic] = [ covariance_matrix, correlation_from_covariance(covariance_matrix) ] # Add the inputMC statistics - xsections_with_symmetrised_systematics['inputMC'], _ = add_inputMC_systematic(options) + # xsections_with_symmetrised_systematics['inputMC'], _ = add_inputMC_systematic(options) + inputMCUncertainties, inputMCCovariance = add_inputMC_systematic(options) + xsections_with_uncertainties_and_covariance['inputMC'] = [ inputMCUncertainties[0], inputMCCovariance ] + d_Cov_Cor['inputMC'] = inputMCCovariance - # Generate the Covariance Matrices - d_Cov_Cor = generate_covariance_matrices( - options, - xsections_with_symmetrised_systematics - ) + # # Generate the Covariance Matrices + # d_Cov_Cor = generate_covariance_matrices( + # options, + # xsections_with_symmetrised_systematics + # ) # Combine Covariance Matrices - if 'BJet' in xsections_with_symmetrised_systematics and 'LightJet' in xsections_with_symmetrised_systematics: + # if 'BJet' in xsections_with_symmetrised_systematics and 'LightJet' in xsections_with_symmetrised_systematics: + # # Combine LightJet and BJet Systematics + # bJet = xsections_with_symmetrised_systematics['BJet'][0] + # lightJet = xsections_with_symmetrised_systematics['LightJet'][0] + # bJet_tot = [combine_errors_in_quadrature([e1, e2]) for e1, e2 in zip(bJet, lightJet)] + # xsections_with_symmetrised_systematics['BJet'][0] = bJet_tot + # del xsections_with_symmetrised_systematics['LightJet'] + + # # Combine the covariance matrices + # l_Cov = [d_Cov_Cor['BJet'][0], d_Cov_Cor['LightJet'][0]] + # d_Cov_Cor['BJet'] = combine_covariance_matrices(l_Cov) + # del d_Cov_Cor['LightJet'] + if 'BJet' in xsections_with_uncertainties_and_covariance and 'LightJet' in xsections_with_uncertainties_and_covariance: # Combine LightJet and BJet Systematics - bJet = xsections_with_symmetrised_systematics['BJet'][0] - lightJet = xsections_with_symmetrised_systematics['LightJet'][0] + bJet = xsections_with_uncertainties_and_covariance['BJet'][0] + lightJet = xsections_with_uncertainties_and_covariance['LightJet'][0] bJet_tot = [combine_errors_in_quadrature([e1, e2]) for e1, e2 in zip(bJet, lightJet)] - xsections_with_symmetrised_systematics['BJet'][0] = bJet_tot - del xsections_with_symmetrised_systematics['LightJet'] + xsections_with_uncertainties_and_covariance['BJet'][0] = bJet_tot + del xsections_with_uncertainties_and_covariance['LightJet'] # Combine the covariance matrices l_Cov = [d_Cov_Cor['BJet'][0], d_Cov_Cor['LightJet'][0]] d_Cov_Cor['BJet'] = combine_covariance_matrices(l_Cov) del d_Cov_Cor['LightJet'] - if 'PDF' in xsections_with_symmetrised_systematics and 'TTJets_alphaS' in xsections_with_symmetrised_systematics: + # if 'PDF' in xsections_with_symmetrised_systematics and 'TTJets_alphaS' in xsections_with_symmetrised_systematics: + # # Combine PDF with alphaS variations + # pdf = xsections_with_symmetrised_systematics['PDF'][0] + # alphaS = xsections_with_symmetrised_systematics['TTJets_alphaS'][0] + # pdf_tot = [combine_errors_in_quadrature([e1, e2]) for e1, e2 in zip(alphaS, pdf)] + # xsections_with_symmetrised_systematics['PDF'][0] = pdf_tot + # del xsections_with_symmetrised_systematics['TTJets_alphaS'] + + # # Combine the covariance matrices + # l_Cov = [d_Cov_Cor['PDF'][0], d_Cov_Cor['TTJets_alphaS'][0]] + # d_Cov_Cor['PDF'] = combine_covariance_matrices(l_Cov) + # del d_Cov_Cor['TTJets_alphaS'] + if 'PDF' in xsections_with_uncertainties_and_covariance and 'TTJets_alphaS' in xsections_with_uncertainties_and_covariance: # Combine PDF with alphaS variations - pdf = xsections_with_symmetrised_systematics['PDF'][0] - alphaS = xsections_with_symmetrised_systematics['TTJets_alphaS'][0] + pdf = xsections_with_uncertainties_and_covariance['PDF'][0] + alphaS = xsections_with_uncertainties_and_covariance['TTJets_alphaS'][0] pdf_tot = [combine_errors_in_quadrature([e1, e2]) for e1, e2 in zip(alphaS, pdf)] - xsections_with_symmetrised_systematics['PDF'][0] = pdf_tot - del xsections_with_symmetrised_systematics['TTJets_alphaS'] + xsections_with_uncertainties_and_covariance['PDF'][0] = pdf_tot + del xsections_with_uncertainties_and_covariance['TTJets_alphaS'] # Combine the covariance matrices l_Cov = [d_Cov_Cor['PDF'][0], d_Cov_Cor['TTJets_alphaS'][0]] @@ -534,11 +590,49 @@ def get_symmetrised_systematic_uncertainty( options, syst_unc_x_secs ): d_Cov_Cor['inputMC'] = add_covariance( options, d_Cov_Cor, addInputMCStat = True ) d_Cov_Cor['Total'] = add_covariance( options, d_Cov_Cor, addTotal = True ) save_covariance_matrices( options, d_Cov_Cor ) - return xsections_with_symmetrised_systematics + return xsections_with_uncertainties_and_covariance + + +def get_average_covariance_matrix( central_measurement, upper_measurement, lower_measurement, options, isTopMassSystematic, systematic ) : + ''' + Calculates the covariance matrix for each source, by averaging the up/down variations in each bin: + C(i,j) = 1 / 2 * ( up(i) * up(j) + down(i) * down(j) ) + + Returns the covariance matrix, and the uncertainties in each bin (which are just the sqrt of the diagonal elements) + + ''' + + number_of_bins = len(central_measurement) + bin_indices = range(number_of_bins) + covariance_matrix = np.matrix( np.zeros( ( number_of_bins, number_of_bins ) ) ) + for i, j in product( bin_indices, bin_indices ): + central_i = central_measurement[i][0] + central_j = central_measurement[j][0] + + upper_variation_i = central_i - upper_measurement[i] + upper_variation_j = central_j - upper_measurement[j] + + lower_variation_i = central_i - lower_measurement[i] + lower_variation_j = central_j - lower_measurement[j] + + if isTopMassSystematic: + upper_variation_i, lower_variation_i = scaleTopMassSystematic( upper_variation_i, lower_variation_i, options['topMasses'], options['topMassUncertainty'] ) + upper_variation_j, lower_variation_j = scaleTopMassSystematic( upper_variation_j, lower_variation_j, options['topMasses'], options['topMassUncertainty'] ) + covariance_ij = 0.5 * ( upper_variation_i * upper_variation_j + lower_variation_i * lower_variation_j ) + + covariance_matrix[i,j] = covariance_ij + + # Extract uncertainties (for compatability later on) + uncertainties = [] + for i in bin_indices: + uncertainties.append( np.sqrt( covariance_matrix[i,i] ) ) + + return uncertainties, covariance_matrix + # @profile(stream=fp) -def get_symmetrised_errors(central_measurement, upper_measurement, lower_measurement, options, isTopMassSystematic=False ): +def get_symmetrised_errors(central_measurement, upper_measurement, lower_measurement, options, isTopMassSystematic=False, systematic=None ): ''' Returns the symmetric error in each bin for a specific systematic and also the sign of the systematic. Sign is used for calculating the covariance matrices. @@ -553,6 +647,7 @@ def get_symmetrised_errors(central_measurement, upper_measurement, lower_measure number_of_bins = len(central_measurement) symm_uncerts = [] sign_uncerts = [] + for c, u, l in zip(central_measurement, upper_measurement, lower_measurement): central = c[0] # Getting the measurement, not the stat unc [xsec, unc] upper = u @@ -624,7 +719,8 @@ def get_measurement_with_total_systematic_uncertainty(options, xsec_with_symmetr for systematic, measurement in xsec_with_symmetrised_systematics.iteritems(): if systematic in ['central','CT14', 'MMHT14']: continue sys_unc += measurement[0][bin_i]**2 - measurement_with_total_uncertainty.append( [central[0], sqrt(sys_unc), sqrt(sys_unc)] ) + + measurement_with_total_uncertainty.append( [central[0], np.sqrt(sys_unc), np.sqrt(sys_unc)] ) return measurement_with_total_uncertainty @@ -647,11 +743,11 @@ def generate_covariance_matrices(options, xsec_with_symmetrised_systematics): sign = measurement[1] # Create the matrices in numpy.matrix format - covariance_matrix, correlation_matrix = generate_covariance_matrix(number_of_bins, systematic, sign) + covariance_matrix, correlation_matrix = generate_covariance_matrix(number_of_bins, systematic, sign, syst) d_Cov_Cor[syst] = [covariance_matrix, correlation_matrix] return d_Cov_Cor -def generate_covariance_matrix(number_of_bins, systematic, sign): +def generate_covariance_matrix(number_of_bins, systematic, sign, syst): ''' Variance_ii = (Unc_i) * (Unc_i) Covariance_ij = (Sign_i*Unc_i) * (Sign_j*Unc_j) @@ -679,7 +775,10 @@ def correlation_from_covariance(covariance_matrix): correlation_matrix = np.matrix( 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] ): - correlation_matrix[i,j] = covariance_matrix[i,j] / sqrt( covariance_matrix[i,i] * covariance_matrix[j,j] ) + if ( covariance_matrix[i,i] * covariance_matrix[j,j] == 0 ): + correlation_matrix[i,j] = 0 + else: + correlation_matrix[i,j] = covariance_matrix[i,j] / sqrt( covariance_matrix[i,i] * covariance_matrix[j,j] ) return correlation_matrix ### Add Covariances From ad0e284813fceaff2b7e20c1f12413982c99c5a0 Mon Sep 17 00:00:00 2001 From: Emyr Clement Date: Tue, 3 Oct 2017 17:15:01 +0100 Subject: [PATCH 3/3] Treat QCD uncertainties in different channels as uncorrelated. --- .../xsection/02_unfold_and_measure.py | 16 +++++++- dps/config/xsection.py | 40 ++++++++++++++----- 2 files changed, 44 insertions(+), 12 deletions(-) diff --git a/dps/analysis/xsection/02_unfold_and_measure.py b/dps/analysis/xsection/02_unfold_and_measure.py index 35f1f13e..44d5767f 100644 --- a/dps/analysis/xsection/02_unfold_and_measure.py +++ b/dps/analysis/xsection/02_unfold_and_measure.py @@ -1249,7 +1249,6 @@ def parse_arguments(): for category in all_measurements: if run_just_central and not category == 'central': continue - if ( variable in measurement_config.variables_no_met ) and (category in measurement_config.met_specific_systematics): continue print 'Unfolding category {}'.format(category) @@ -1280,13 +1279,26 @@ def parse_arguments(): elif category == 'Electron_up' or category == 'Electron_down': normalisation_results_electron = read_tuple_from_file( electron_file ) normalisation_results_muon = read_tuple_from_file( path_to_DF + '/central/normalisation_muon.txt' ) + elif category == 'QCD_shape_electron' or category == 'QCD_cross_section_electron': + if category == 'QCD_shape_electron': + normalisation_results_electron = read_tuple_from_file( electron_file.replace('QCD_shape_electron','QCD_shape') ) + else: + normalisation_results_electron = read_tuple_from_file( electron_file.replace('QCD_cross_section_electron','QCD_cross_section') ) + + normalisation_results_muon = read_tuple_from_file( path_to_DF + '/central/normalisation_muon.txt' ) + elif category == 'QCD_shape_muon' or category == 'QCD_cross_section_muon': + normalisation_results_electron = read_tuple_from_file( path_to_DF + '/central/normalisation_electron.txt' ) + if category == 'QCD_shape_muon': + normalisation_results_muon = read_tuple_from_file( muon_file.replace('QCD_shape_muon','QCD_shape') ) + else: + normalisation_results_muon = read_tuple_from_file( muon_file.replace('QCD_cross_section_muon','QCD_cross_section') ) else: normalisation_results_electron = read_tuple_from_file( electron_file ) normalisation_results_muon = read_tuple_from_file( muon_file ) # Add additional 50% uncertainty here instead of 03. Can do it before combining channels. # I think this is the correct way of doing it. - if category == 'QCD_cross_section': + if category == 'QCD_cross_section' or category == 'QCD_cross_section_electron': for i in range(0, len(normalisation_results_electron['QCD']) ): normalisation_results_electron['QCD'][i] = tuple([1.5*x for x in normalisation_results_electron['QCD'][i]]) diff --git a/dps/config/xsection.py b/dps/config/xsection.py index d32c63b3..f945efd4 100644 --- a/dps/config/xsection.py +++ b/dps/config/xsection.py @@ -86,7 +86,7 @@ def __fill_defaults__( self ): # self.path_to_unfolding_histograms = '/hdfs/TopQuarkGroup/run2/unfolding/13TeV/Moriond2017/' # self.path_to_unfolding_histograms = '/hdfs/TopQuarkGroup/run2/unfolding/13TeV/EPS2017/' - self.path_to_unfolding_histograms = 'unfolding/13TeV/' + self.path_to_unfolding_histograms = 'unfolding/13TeV_debug/' path_to_unfolding_histograms = self.path_to_unfolding_histograms @@ -203,6 +203,12 @@ def __fill_defaults__( self ): 'QCD_cross_section', 'QCD_shape', + 'QCD_cross_section_electron', + 'QCD_shape_electron', + + 'QCD_cross_section_muon', + 'QCD_shape_muon', + 'QCD_other_control_region', 'QCD_signal_MC' ] @@ -243,14 +249,24 @@ def __fill_defaults__( self ): self.measurements = self.normalisation_systematics + self.generator_systematics self.list_of_systematics = { - # Theoretical Uncertainties (Rate Changing) + # # Theoretical Uncertainties (Rate Changing) 'V+Jets_cross_section' : ['V+Jets_cross_section+', 'V+Jets_cross_section-'], - 'QCD_cross_section' : ['QCD_cross_section', 'QCD_cross_section'], + # 'QCD_cross_section' : ['QCD_cross_section', 'QCD_cross_section'], 'SingleTop_cross_section' : ['SingleTop_cross_section+', 'SingleTop_cross_section-'], 'luminosity' : ['luminosity+', 'luminosity-'], # QCD Shape - 'QCD_shape' : ['QCD_shape', 'QCD_shape'], - # Generator Uncertainties + # 'QCD_shape' : ['QCD_shape', 'QCD_shape'], + + 'QCD_shape_electron' : ['QCD_shape_electron', 'QCD_shape_electron'], + 'QCD_cross_section_electron' : ['QCD_cross_section_electron', 'QCD_cross_section_electron'], + + 'QCD_shape_muon' : ['QCD_shape_muon', 'QCD_shape_muon'], + 'QCD_cross_section_muon' : ['QCD_cross_section_muon', 'QCD_cross_section_muon'], + + # 'QCD_other_control_region' : ['QCD_other_control_region', 'QCD_other_control_region'], + + + # # Generator Uncertainties 'TTJets_mass' : ['TTJets_massup', 'TTJets_massdown'], 'TTJets_ue' : ['TTJets_ueup', 'TTJets_uedown'], 'TTJets_topPt' : ['TTJets_topPt', 'TTJets_topPt'], @@ -260,10 +276,10 @@ def __fill_defaults__( self ): 'TTJets_fsrup', 'TTJets_fsrdown', 'TTJets_isrup', 'TTJets_isrdown' ], - # 'TTJets_fsr' : ['TTJets_fsrup', 'TTJets_fsrdown'], - # 'TTJets_isr' : ['TTJets_isrup', 'TTJets_isrdown'], + # # 'TTJets_fsr' : ['TTJets_fsrup', 'TTJets_fsrdown'], + # # 'TTJets_isr' : ['TTJets_isrup', 'TTJets_isrdown'], - 'TTJets_alphaS' : ['TTJets_alphaSup', 'TTJets_alphaSdown'], + # 'TTJets_alphaS' : ['TTJets_alphaSup', 'TTJets_alphaSdown'], 'TTJets_hdamp' : ['TTJets_hdampup', 'TTJets_hdampdown'], 'TTJets_semiLepBr' : ['TTJets_semiLepBrup', 'TTJets_semiLepBrdown'], 'TTJets_frag' : ['TTJets_fragup', 'TTJets_fragdown'], @@ -274,7 +290,7 @@ def __fill_defaults__( self ): 'TTJets_CR_GluonMove' : ['TTJets_GluonMove', 'TTJets_GluonMove'], # 'TTJets_CR_GluonMove_erdOn' : ['TTJets_GluonMove_erdOn', 'TTJets_GluonMove_erdOn'], - # Event Reweighting + # # Event Reweighting 'PileUp' : ['PileUp_up', 'PileUp_down'], 'JES' : ['JES_up', 'JES_down'], 'JER' : ['JER_up', 'JER_down'], @@ -284,7 +300,7 @@ def __fill_defaults__( self ): 'Electron' : ['Electron_up', 'Electron_down'], 'Muon' : ['Muon_up', 'Muon_down'], # PDF Uncertainties - 'PDF' : ['PDF', 'PDF'], + # 'PDF' : ['PDF', 'PDF'], # MET Uncertainties 'ElectronEn' : ['ElectronEnUp', 'ElectronEnDown'], 'MuonEn' : ['MuonEnUp', 'MuonEnDown'], @@ -322,6 +338,10 @@ def __fill_defaults__( self ): 'SingleTop_cross_section', 'QCD_cross_section', 'QCD_shape', + 'QCD_cross_section_electron', + 'QCD_shape_electron', + 'QCD_cross_section_muon', + 'QCD_shape_muon', ] self.systematic_group_met = [ 'ElectronEn',