diff --git a/bin/x_01b_all_vars b/bin/x_01b_all_vars index 45cda92c..9f1b3df6 100755 --- a/bin/x_01b_all_vars +++ b/bin/x_01b_all_vars @@ -5,7 +5,7 @@ N_JOBS=4 i=0 -for var in MET HT ST WPT; do +for var in MET HT ST WPT lepton_pt abs_lepton_eta NJets; do echo "Fitting distribution: $var" nohup time python src/cross_section_measurement/01_get_ttjet_normalisation.py -v $var -i config/measurements/background_subtraction &> logs/01_${var}_bgs_13TeV_fullPS.log & let i+=1 @@ -17,7 +17,7 @@ for var in MET HT ST WPT; do done wait; -for var in MET HT ST WPT; do +for var in MET HT ST WPT lepton_pt abs_lepton_eta NJets; do echo "Fitting distribution: $var" nohup time python src/cross_section_measurement/01_get_ttjet_normalisation.py -v $var -i config/measurements/background_subtraction --visiblePS &> logs/01_${var}_bgs_13TeV_visiblePS.log & let i+=1 diff --git a/bin/x_02b_all_vars b/bin/x_02b_all_vars index 0ddf41bc..95349025 100755 --- a/bin/x_02b_all_vars +++ b/bin/x_02b_all_vars @@ -1,7 +1,7 @@ #!/bin/bash echo "This will take a while ... grab a coffee/tea/water" mkdir -p logs -N_JOBS=3 +N_JOBS=4 i=0 diff --git a/bin/x_05_all_vars b/bin/x_05_all_vars index e7ba6081..ab4fe72e 100755 --- a/bin/x_05_all_vars +++ b/bin/x_05_all_vars @@ -4,7 +4,7 @@ mkdir -p logs mkdir -p plots fit_var="M3,angle_bl" nice_fit_var=`echo $fit_var | sed 's/,/_/g'` -N_JOBS=4 +N_JOBS=1 echo "Using the fit variable(s): $fit_var" # input_folder=/hdfs/TopQuarkGroup/run2/dpsData/data/normalisation/background_subtraction/ @@ -12,23 +12,24 @@ input_folder=data/normalisation/background_subtraction/ output_folder=tables/background_subtraction/ com=13 typical_systematics="${output_folder}/${com}TeV/FullPS/typical_systematics_${com}TeV_combined.tex" -typical_systematics_vis="${output_folder}/${com}TeV/VisiblePS/typical_systematics_${com}TeV_combined_vis.tex" -if [ -f $typical_systematics ]; +typical_systematics_vis="${output_folder}/${com}TeV/VisiblePS/typical_systematics_${com}TeV_combined.tex" +if [ -e $typical_systematics ]; then echo "Cleaning up old typical systematics file for full phase space" rm $typical_systematics fi -if [ -f $typical_systematics_vis ]; +if [ -e $typical_systematics_vis ]; then echo "Cleaning up old typical systematics file for visible phase space" rm $typical_systematics_vis fi +echo "Visible phase space ..." i=0 for var in MET HT ST WPT lepton_pt abs_lepton_eta NJets; do echo "Tabulating diff. x-section for distribution: $var" - nohup time python src/cross_section_measurement/05_make_tables.py -v $var -c $com -p $input_folder -o $output_folder -a &> logs/05_${var}_table_${com}TeV.log & + nohup time python src/cross_section_measurement/05_make_tables.py -v $var -c 13 -p $input_folder -o $output_folder -a --visiblePS &> logs/05_${var}_table_13TeV_vis.log & let i+=1 if (( $i % N_JOBS == 0 )) then @@ -37,18 +38,17 @@ for var in MET HT ST WPT lepton_pt abs_lepton_eta NJets; do fi done -wait; - -echo "Now visible phase space ..." i=0 -for var in MET HT ST WPT lepton_pt abs_lepton_eta NJets; do - echo "Tabulating diff. x-section for distribution: $var" - nohup time python src/cross_section_measurement/05_make_tables.py -v $var -c 13 -p $input_folder -o $output_folder -a --visiblePS &> logs/05_${var}_table_13TeV_vis.log & - let i+=1 - if (( $i % N_JOBS == 0 )) - then - echo "Waiting on the above to finish." - wait; - fi -done +# for var in MET HT ST WPT lepton_pt abs_lepton_eta NJets; do +# echo "Tabulating diff. x-section for distribution: $var" +# nohup time python src/cross_section_measurement/05_make_tables.py -v $var -c $com -p $input_folder -o $output_folder -a &> logs/05_${var}_table_${com}TeV.log & +# let i+=1 +# if (( $i % N_JOBS == 0 )) +# then +# echo "Waiting on the above to finish." +# wait; +# fi +# done + +wait; echo "All done!" diff --git a/config/CMS.py b/config/CMS.py index 6c0aa8cf..84651f4a 100644 --- a/config/CMS.py +++ b/config/CMS.py @@ -22,7 +22,7 @@ } y_axis_title_small = { - 'fontsize':40, + 'fontsize':50, 'position' : (0, 1.), 'verticalalignment': 'bottom', 'horizontalalignment': 'right' diff --git a/config/cross_section_config.py b/config/cross_section_config.py index 8d3f41cd..a4526fd9 100644 --- a/config/cross_section_config.py +++ b/config/cross_section_config.py @@ -5,7 +5,7 @@ class XSectionConfig(): current_analysis_path = '/hdfs/TopQuarkGroup/run2/atOutput/' known_centre_of_mass_energies = [7,8,13] # has to be separate as many variables depend on it - luminosities = {7:5050, 8:19584, 13:40.028} + luminosities = {7:5050, 8:19584, 13:41.62912} parameters = ['SingleTop_category_templates', 'SingleTop_category_templates_trees', 'SingleTop_file', 'VJets_category_templates', 'VJets_category_templates_trees', 'analysis_types', 'categories_and_prefixes', 'central_general_template', @@ -185,6 +185,8 @@ def __fill_defaults__( self ): # 'LightJet_down':'_minusLightJet', # 'LightJet_up':'_plusLightJet', + 'PileUpSystematic' : '', + # Other MET uncertainties not already included 'ElectronEnUp' : '', @@ -292,6 +294,10 @@ def __fill_defaults__( self ): self.ttbar_scaledown_category_templates_trees = path_to_files + '/TTJets_PowhegPythia8_scaledown_tree.root' self.ttbar_mtop1695_category_templates_trees = path_to_files + '/TTJets_PowhegPythia8_mtop1695_tree.root' self.ttbar_mtop1755_category_templates_trees = path_to_files + '/TTJets_PowhegPythia8_mtop1755_tree.root' + self.ttbar_jesup_category_templates_trees = path_to_files + '/TTJets_PowhegPythia8_plusJES_tree.root' + self.ttbar_jesdown_category_templates_trees = path_to_files + '/TTJets_PowhegPythia8_minusJES_tree.root' + self.ttbar_jerup_category_templates_trees = path_to_files + '/TTJets_PowhegPythia8_plusJER_tree.root' + self.ttbar_jerdown_category_templates_trees = path_to_files + '/TTJets_PowhegPythia8_minusJER_tree.root' self.data_muon_category_templates = { 'central': self.data_file_muon, @@ -347,6 +353,26 @@ def __fill_defaults__( self ): self.unfolding_matching_up = self.unfolding_matching_up_raw.replace( '.root', '_asymmetric.root' ) 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 + self.unfolding_Lepton_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_leptondown_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_Lepton_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_leptonup_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_jes_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_jesdown_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_jes_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_jesup_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_jer_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_jerdown_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_jer_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_jerup_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_bjet_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_bjetdown_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_bjet_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_bjetup_asymmetric.root' % self.centre_of_mass_energy + + self.unfolding_ElectronEn_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_ElectronEnDown_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_ElectronEn_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_ElectronEnUp_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_MuonEn_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_MuonEnDown_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_MuonEn_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_MuonEnUp_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_TauEn_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_TauEnDown_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_TauEn_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_TauEnUp_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_UnclusteredEn_down = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_UnclusteredEnDown_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_UnclusteredEn_up = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_UnclusteredEnUp_asymmetric.root' % self.centre_of_mass_energy + + self.unfolding_PUSystematic = path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_pileupSystematic_asymmetric.root' % self.centre_of_mass_energy + self.unfolding_pdfweights = {index : path_to_unfolding_histograms + 'unfolding_TTJets_%dTeV_asymmetric_generatorWeight_%d.root' % (self.centre_of_mass_energy, index) for index in range( 9, 109 )} @@ -393,8 +419,8 @@ def __fill_defaults__( self ): 'Lepton selection efficiency': [('Electron_down', 'Electron_up'), ('Muon_down', 'Muon_up')], 'Jet energy \& resolution': [('JES_down', 'JES_up', 'JER_down', 'JER_up')], 'B-tagging' : [('BJet_down', 'BJet_up')], - 'MET' : [('ElectronEnDown', 'ElectronEnUp'), ('MuonEnDown','MuonEnUp'),('TauEnDown','TauEnUp'),('UnclusteredEnDown','UnclusteredEnUp')], - 'Normalisation': [ + '\met' : [('ElectronEnDown', 'ElectronEnUp'), ('MuonEnDown','MuonEnUp'),('TauEnDown','TauEnUp'),('UnclusteredEnDown','UnclusteredEnUp')], + 'Background Normalisation': [ # ('TTJet_cross_section-', 'TTJet_cross_section+'), ('SingleTop_cross_section-', 'SingleTop_cross_section+'), ('luminosity-', 'luminosity+'), @@ -404,6 +430,7 @@ def __fill_defaults__( self ): 'Hadronisation': [('TTJets_hadronisation', 'TTJets_hadronisation')], 'NLO generator': [('TTJets_NLOgenerator', 'TTJets_NLOgenerator')], 'PDF': [('PDF_total_lower', 'PDF_total_upper')], + 'Pileup' : [('PileUpSystematic','PileUpSystematic')], 'QCD shape': [('QCD_shape', 'QCD_shape')] } self.typical_systematics = [] @@ -524,7 +551,8 @@ def __fill_defaults_13TeV__( self ): middle = self.middle path_to_files = self.path_to_files - self.new_luminosity = 40.028 # pb^-1 + # self.new_luminosity = 40.028 # pb^-1 + self.new_luminosity = 41.62912 # pb^-1 self.ttbar_xsection = 831.76 # pb self.rate_changing_systematics = {#TODO check where this is used diff --git a/config/latex_labels.py b/config/latex_labels.py index 73222c95..f0d81ad1 100644 --- a/config/latex_labels.py +++ b/config/latex_labels.py @@ -28,12 +28,22 @@ 'ttbarM': '\ensuremath{M_\mathrm{t\\bar{t}}}', 'ttbarRap': '\ensuremath{y_{\mathrm{t\\bar{t}}}}', 'NJets': '\ensuremath{N_{\mathrm{Jets}}}', - 'lepton_pt': '\ensuremath{ \mathrm{lepton} p_{\mathrm{T}} }', - 'lepton_eta': '\ensuremath{ \mathrm{lepton} \eta }', - 'abs_lepton_eta': '\ensuremath{ \mathrm{lepton} |\eta| }', + 'lepton_pt': '\ensuremath{ p_{\mathrm{T}}^\mathrm{l} }', + 'lepton_eta': '\ensuremath{ \eta^\mathrm{l} }', + 'abs_lepton_eta': '\ensuremath{ |\eta^\mathrm{l}| }', 'bjets_pt': '\ensuremath{ \mathrm{b-jet} p_{\mathrm{T}} }', 'bjets_eta': '\ensuremath{ \mathrm{b-jet} \eta }', + 'sigmaietaieta' : '\ensuremath{\sigma_{i\eta i \eta}}', +} +variables_NonLatex = { + 'MET': 'MET', + 'HT': 'HT', + 'ST': 'ST', + 'WPT': 'WPT', + 'NJets': 'N Jets', + 'lepton_pt': 'lepton pt', + 'abs_lepton_eta': 'lepton eta', } control_plots_latex = { @@ -53,19 +63,23 @@ 'MADGRAPH': '$t\\bar{t}$ (MADGRAPH+Pythia)', 'MADGRAPH_ptreweight': '$t\\bar{t}$ (MADGRAPH+$p_\mathrm{T}^\mathrm{reweight}$)', 'amcatnlo': 'aMC@NLO', - 'madgraphMLM': 'MADGRAPHMLM', + 'madgraphMLM': 'Madgraph', 'POWHEG_PYTHIA': '$t\\bar{t}$ (POWHEG+Pythia)', - 'POWHEG_HERWIG': '$t\\bar{t}$ (POWHEG+Herwig)', - 'powhegPythia8': 'POWHEG+Pythia8', + 'POWHEG_HERWIG': 'Powheg Herwig++', + 'powhegPythia8': 'Powheg Pythia8', 'pythia8': '$t\\bar{t}$ (Pythia8)', '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)', + # 'scaledown': '$t\\bar{t}$ ($Q^{2}$ down)', + # 'scaleup': '$t\\bar{t}$ ($Q^{2}$ up)', + 'scaledown': '$Q^{2}$ down', + 'scaleup': '$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_scaledown': '$t\\bar{t}$ ($Q^{2}$ down)', + # 'TTJets_scaleup': '$t\\bar{t}$ ($Q^{2}$ up)', + 'TTJets_scaledown': '$Q^{2}$ down', + 'TTJets_scaleup': '$Q^{2}$ up', 'TTJets_massdown': '$t\\bar{t}$ (top mass down)', 'TTJets_massup': '$t\\bar{t}$ (top mass up)', 'VJets_matchingdown': 'V+jets (matching down)', @@ -84,6 +98,7 @@ 'LightJet_up': 'b-tagging mis-tag rate $+1\sigma$', 'PU_down': 'Pile-up $-1\sigma$', 'PU_up': 'Pile-up $+1\sigma$', + 'PileUpSystematic' : 'Pile-up', 'central': 'central', #'ptreweight_max': '$p_\mathrm{T}(t,\\bar{t})$ reweighting', 'TTJets_hadronisation': 'Hadronisation uncertainty', @@ -113,8 +128,10 @@ 'SingleTop_cross_section+': 'Single top cross section \ensuremath{+1\sigma}', 'luminosity-': 'Luminosity \ensuremath{-1\sigma}', 'luminosity+': 'Luminosity \ensuremath{+1\sigma}', - 'massdown': '\ensuremath{t\\bar{t}} (top mass down)', - 'massup': '\ensuremath{t\\bar{t}} (top mass up)', + # 'massdown': '\ensuremath{t\\bar{t}} (top mass down)', + # 'massup': '\ensuremath{t\\bar{t}} (top mass up)', + 'massdown': 'Top mass down', + 'massup': 'Top mass up', } met_systematics_latex = { diff --git a/config/variable_binning.py b/config/variable_binning.py index f3ab141d..4705d3e4 100644 --- a/config/variable_binning.py +++ b/config/variable_binning.py @@ -88,6 +88,8 @@ 'WPT' : [i * 25 for i in range ( 0, 17 )], 'HT' : [i * 50 for i in range ( 0, 21 )], 'ST' : [i * 50 for i in range ( 2, 25 )], + # 'sigmaietaieta' : [i * 0.001 for i in range ( 12, 40 )], + 'sigmaietaieta' : [i * 0.002 for i in range ( 0, 20 )], } diff --git a/experimental/BLTUnfold/produceUnfoldingHistograms.py b/experimental/BLTUnfold/produceUnfoldingHistograms.py index a438b355..1e47d26c 100644 --- a/experimental/BLTUnfold/produceUnfoldingHistograms.py +++ b/experimental/BLTUnfold/produceUnfoldingHistograms.py @@ -93,6 +93,24 @@ def getFileName( com, sample, measurementConfig ) : 'scaledown' : measurementConfig.ttbar_scaledown_category_templates_trees, 'massdown' : measurementConfig.ttbar_mtop1695_category_templates_trees, 'massup' : measurementConfig.ttbar_mtop1755_category_templates_trees, + 'jesdown' : measurementConfig.ttbar_jesdown_category_templates_trees, + 'jesup' : measurementConfig.ttbar_jesup_category_templates_trees, + 'jerdown' : measurementConfig.ttbar_jerdown_category_templates_trees, + 'jerup' : measurementConfig.ttbar_jerup_category_templates_trees, + 'bjetdown' : measurementConfig.ttbar_category_templates_trees['central'], + 'bjetup' : measurementConfig.ttbar_category_templates_trees['central'], + 'leptondown' : measurementConfig.ttbar_category_templates_trees['central'], + 'leptonup' : measurementConfig.ttbar_category_templates_trees['central'], + 'pileupSystematic' : measurementConfig.ttbar_category_templates_trees['central'], + + 'ElectronEnUp' : measurementConfig.ttbar_category_templates_trees['central'], + 'ElectronEnDown' : measurementConfig.ttbar_category_templates_trees['central'], + 'MuonEnUp' : measurementConfig.ttbar_category_templates_trees['central'], + 'MuonEnDown' : measurementConfig.ttbar_category_templates_trees['central'], + 'TauEnUp' : measurementConfig.ttbar_category_templates_trees['central'], + 'TauEnDown' : measurementConfig.ttbar_category_templates_trees['central'], + 'UnclusteredEnUp' : measurementConfig.ttbar_category_templates_trees['central'], + 'UnclusteredEnDown' : measurementConfig.ttbar_category_templates_trees['central'], }, } @@ -109,6 +127,7 @@ def main(): parser.add_option('--topPtReweighting', action='store_true', dest='applyTopPtReweighting', default=False ) parser.add_option('-c', '--centreOfMassEnergy', dest='centreOfMassEnergy', type='int', default=13 ) parser.add_option('--generatorWeight', type='int', dest='generatorWeight', default=-1 ) + parser.add_option('--nGeneratorWeights', type='int', dest='nGeneratorWeights', default=1 ) parser.add_option('-s', '--sample', dest='sample', default='central') parser.add_option('-d', '--debug', action='store_true', dest='debug', default=False) parser.add_option('-n', action='store_true', dest='donothing', default=False) @@ -124,276 +143,337 @@ def main(): if int(options.centreOfMassEnergy) == 13: # file_name = fileNames['13TeV'][options.sample] file_name = getFileName('13TeV', options.sample, measurement_config) - if options.generatorWeight >= 0: - file_name = 'localInputFile.root' + # if options.generatorWeight >= 0: + # file_name = 'localInputFile.root' else: print "Error: Unrecognised centre of mass energy." + generatorWeightsToRun = [] + if options.nGeneratorWeights > 1 : + for i in range (0, options.nGeneratorWeights): + generatorWeightsToRun.append( options.generatorWeight + i ) + elif options.generatorWeight >= 0 : + generatorWeightsToRun.append(options.generatorWeight) + else: generatorWeightsToRun.append( -1 ) + # Output file name outputFileName = 'crap.root' outputFileDir = 'unfolding/%sTeV/' % options.centreOfMassEnergy make_folder_if_not_exists(outputFileDir) energySuffix = '%sTeV' % ( options.centreOfMassEnergy ) - - if options.applyTopPtReweighting: - outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_withTopPtReweighting.root' % energySuffix - elif options.generatorWeight == 4: - outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_scaleUpWeight.root' % ( energySuffix ) - elif options.generatorWeight == 8: - outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_scaleDownWeight.root' % ( energySuffix ) - elif options.generatorWeight >= 0: - outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_generatorWeight_%i.root' % ( energySuffix, options.generatorWeight ) - elif options.sample != 'central': - 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 - with root_open( file_name, 'read' ) as f, root_open( outputFileName, 'recreate') as out: - - for channel in channels: - if options.debug and channel.channelName != 'muPlusJets' : continue + for meWeight in generatorWeightsToRun : + if options.applyTopPtReweighting: + outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_withTopPtReweighting.root' % energySuffix + elif meWeight == 4: + outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_scaleUpWeight.root' % ( energySuffix ) + elif meWeight == 8: + outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_scaleDownWeight.root' % ( energySuffix ) + elif meWeight >= 9 and meWeight <= 108: + outputFileName = outputFileDir+'/unfolding_TTJets_%s_asymmetric_generatorWeight_%i.root' % ( energySuffix, meWeight ) + elif options.sample != 'central': + 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 + + with root_open( file_name, 'read' ) as f, root_open( outputFileName, 'recreate') as out: - print 'Channel : ',channel.channelName - - # Get the tree - tree = f.Get("TTbar_plus_X_analysis/Unfolding/Unfolding") - nEntries = tree.GetEntries() - print "Number of entries in tree : ", nEntries - - if options.generatorWeight >= 0 : - tree.AddFriend('TTbar_plus_X_analysis/Unfolding/GeneratorSystematicWeights') - tree.SetBranchStatus('genWeight_*',0) - tree.SetBranchStatus('genWeight_%i' % options.generatorWeight, 1) - - # Keep record of generator weight - if options.generatorWeight >= 0: - generatorWeight = '( genWeight_%i )' % options.generatorWeight - generatorWeightHist = Hist( 50, 0.4, 1.6, name='generatorWeights_'+channel.channelName ) - if not options.donothing: - tree.Draw( generatorWeight, hist=generatorWeightHist) - outputDir = 0 - if not ( out.FindObject('generatorWeights') ): - outputDir = out.mkdir('generatorWeights') - else : - outputDir = out.Get('generatorWeights') - outputDir.cd() - generatorWeightHist.Write() - pass + for channel in channels: + if options.debug and channel.channelName != 'muPlusJets' : continue + + print 'Channel : ',channel.channelName + + # Get the tree + treeName = "TTbar_plus_X_analysis/Unfolding/Unfolding" + if options.sample == "jesup": + treeName += "_JESUp" + elif options.sample == "jesdown": + treeName += "_JESDown" + elif options.sample == "jerup": + treeName += "_JERUp" + elif options.sample == "jerdown": + treeName += "_JERDown" + + tree = f.Get(treeName) + nEntries = tree.GetEntries() + print "Number of entries in tree : ", nEntries + + if meWeight >= 0 : + tree.AddFriend('TTbar_plus_X_analysis/Unfolding/GeneratorSystematicWeights') + tree.SetBranchStatus('genWeight_*',0) + tree.SetBranchStatus('genWeight_%i' % meWeight, 1) + + # Keep record of generator weight + if meWeight >= 0: + generatorWeight = '( genWeight_%i )' % meWeight + generatorWeightHist = Hist( 50, 0.4, 1.6, name='generatorWeights_'+channel.channelName ) + if not options.donothing: + tree.Draw( generatorWeight, hist=generatorWeightHist) + outputDir = 0 + if not ( out.FindObject('generatorWeights') ): + outputDir = out.mkdir('generatorWeights') + else : + outputDir = out.Get('generatorWeights') + outputDir.cd() + generatorWeightHist.Write() + pass - # For variables where you want bins to be symmetric about 0, use abs(variable) (but also make plots for signed variable) - allVariablesBins = bin_edges.copy() - for variable in bin_edges: - if 'Rap' in variable: - allVariablesBins['abs_%s' % variable] = [0,bin_edges[variable][-1]] + # For variables where you want bins to be symmetric about 0, use abs(variable) (but also make plots for signed variable) + allVariablesBins = bin_edges.copy() + for variable in bin_edges: + if 'Rap' in variable: + allVariablesBins['abs_%s' % variable] = [0,bin_edges[variable][-1]] - for variable in allVariablesBins: - if options.debug and variable != 'HT' : continue - - print '--->Doing variable :',variable + for variable in allVariablesBins: + if options.debug and variable != 'HT' : continue - # Output dir name - metSuffix='_patType1CorrectedPFMet' - if variable is 'HT': - metSuffix='' - pass - # Make dir in output file - outputDir = out.mkdir('unfolding_'+variable+'_analyser_'+channel.outputDirName+'_channel'+metSuffix) - - # - # Variable names - # - recoVariable = branchNames[variable] - genVariable_particle = genBranchNames_particle[variable] - genVariable_parton = None - if variable in genBranchNames_parton: - genVariable_parton = genBranchNames_parton[variable] - - # - # Weights and selection - # - - # Generator level - genWeight = '( EventWeight * %.4f)' % measurement_config.luminosity_scale - # genWeight = '( unfolding.puWeight )' - genSelection = '' - genSelectionVis = '' - if channel.channelName is 'muPlusJets' : - genSelection = '( isSemiLeptonicMuon == 1 )' - genSelectionVis = '( isSemiLeptonicMuon == 1 && passesGenEventSelection )' - elif channel.channelName is 'ePlusJets' : - genSelection = '( isSemiLeptonicElectron == 1 )' - genSelectionVis = '( isSemiLeptonicElectron == 1 && passesGenEventSelection )' - - # Offline level - # offlineWeight = '( unfolding.bTagWeight * unfolding.puWeight )' - leptonWeight = '1' - - offlineWeight = '( EventWeight * %s * %.4f)' % ( leptonWeight, measurement_config.luminosity_scale ) - offlineSelection = '' - if channel.channelName is 'muPlusJets' : - offlineSelection = '( passSelection == 1 )' - elif channel.channelName is 'ePlusJets' : - offlineSelection = '( passSelection == 2 )' - - # Fake selection - fakeSelection = '( ' + offlineSelection+"&&!"+genSelection +' ) ' - fakeSelectionVis = '( ' + offlineSelection+"&&!"+genSelectionVis +' ) ' - - # Weights derived from variables in tree - if options.applyTopPtReweighting: - ptWeight = topPtWeight( int(options.centreOfMassEnergy) ) - offlineWeight += ' * '+ptWeight - genWeight += ' * '+ptWeight - pass - - # Apply generator weight - if options.generatorWeight >= 0: - generatorWeight = '( genWeight_%i )' % options.generatorWeight - offlineWeight += ' * '+generatorWeight - genWeight += ' * '+generatorWeight - nEntries = 1000000 - print 'Changing nEntries to ',nEntries, "<---- DOES NOT SEEM TO DO ANYTHING" - pass + if options.sample in measurement_config.met_systematics and variable not in ['MET', 'ST', 'WPT']: + continue - # Scale factors - # scaleFactor = getScaleFactor( options.centreOfMassEnergy, channel.channelName ) - # scaleFactor = '( unfolding.leptonWeight )' - # offlineWeight += ' * '+scaleFactor - - # - # Histograms to fill - # - # 1D histograms - truth = Hist( allVariablesBins[variable], name='truth') - truthVis = Hist( bin_edges_vis[variable], name='truthVis') - truth_parton = Hist( allVariablesBins[variable], name='truth_parton') - measured = Hist( allVariablesBins[variable], name='measured') - measuredVis = Hist( bin_edges_vis[variable], name='measuredVis') - fake = Hist( allVariablesBins[variable], name='fake') - - # 2D histograms - response = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response') - response_without_fakes = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_without_fakes') - response_only_fakes = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_only_fakes') - - responseVis_without_fakes = Hist2D( bin_edges_vis[variable], bin_edges_vis[variable], name='responseVis_without_fakes') - responseVis_only_fakes = Hist2D( bin_edges_vis[variable], bin_edges_vis[variable], name='responseVis_only_fakes') - - response_parton = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_parton') - response_without_fakes_parton = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_without_fakes_parton') - response_only_fakes_parton = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_only_fakes_parton') - - if options.fineBinned: - minVar = trunc( allVariablesBins[variable][0] ) - maxVar = trunc( max( tree.GetMaximum(genVariable_particle), tree.GetMaximum( recoVariable ) ) * 1.2 ) - nBins = int(maxVar - minVar) - if variable is 'lepton_eta' or variable is 'bjets_eta': - maxVar = 2.5 - minVar = -2.5 - nBins = 1000 - elif 'abs' in variable and 'eta' in variable: - maxVar = 3.0 - minVar = 0. - nBins = 1000 - elif 'Rap' in variable: - maxVar = 3.0 - minVar = -3.0 - nBins = 1000 - elif 'NJets' in variable: - maxVar = 20.5 - minVar = -0.5 - nBins = 21 - - truth = Hist( nBins, minVar, maxVar, name='truth') - truthVis = Hist( nBins, minVar, maxVar, name='truthVis') - truth_parton = Hist( nBins, minVar, maxVar, name='truth_parton') - measured = Hist( nBins, minVar, maxVar, name='measured') - measuredVis = Hist( nBins, minVar, maxVar, name='measuredVis') - fake = Hist( nBins, minVar, maxVar, name='fake') - response = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response') - response_without_fakes = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_without_fakes') - response_only_fakes = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_only_fakes') - responseVis_without_fakes = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='responseVis_without_fakes') - responseVis_only_fakes = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='responseVis_only_fakes') - - response_parton = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_parton') - response_without_fakes_parton = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_without_fakes_parton') - response_only_fakes_parton = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_only_fakes_parton') - - # Some interesting histograms - puOffline = Hist( 20, 0, 2, name='puWeights_offline') - eventWeight = Hist( 100, -2, 2, name='EventWeight') - - phaseSpaceInfoHist = Hist( 10, 0, 1, name='phaseSpaceInfoHist') - if options.sample == 'central' and options.generatorWeight == -1 : - nVis = ( tree.Draw( '1', genWeight+'*'+genSelectionVis ) ).Integral() - nVisNotOffline = ( tree.Draw( '1', genWeight+'* ( '+genSelectionVis+'&& !'+offlineSelection+')' ) ).Integral() - nOffline = ( tree.Draw( '1', offlineWeight+'*'+offlineSelection ) ).Integral() - nOfflineNotVis = ( tree.Draw( '1', offlineWeight+'* ( '+offlineSelection+'&& !'+genSelectionVis+')' ) ).Integral() - nFull = ( tree.Draw( '1', genWeight+'*'+genSelection ) ).Integral() - phaseSpaceInfoHist.SetBinContent(1, nVisNotOffline / nVis) - phaseSpaceInfoHist.SetBinContent(2, nOfflineNotVis / nOffline) - phaseSpaceInfoHist.SetBinContent(3, nVis / nFull) - - # - # Fill histograms - # - if not options.donothing: - # tree.Draw('(EventWeight * %.4f)' % measurement_config.luminosity_scale,'1',hist=eventWeight) - # 1D - tree.Draw(genVariable_particle,genWeight+'*'+genSelection,hist=truth, nentries=nEntries) - tree.Draw(genVariable_particle,genWeight+'*'+genSelectionVis,hist=truthVis, nentries=nEntries) - if genVariable_parton != None: - tree.Draw(genVariable_parton,genWeight+'*'+genSelection,hist=truth_parton, nentries=nEntries) - tree.Draw(recoVariable,offlineWeight+'*'+offlineSelection,hist=measured, nentries=nEntries) - tree.Draw(recoVariable,offlineWeight+'*'+offlineSelection,hist=measuredVis, nentries=nEntries) - tree.Draw(recoVariable,offlineWeight+'*'+fakeSelection,hist=fake, nentries=nEntries) - # 2D - tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'*'+offlineSelection,hist=response, nentries=nEntries) - tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'* ('+offlineSelection+'&&'+genSelection +')',hist=response_without_fakes, nentries=nEntries) - tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'*'+fakeSelection,hist=response_only_fakes, nentries=nEntries) - - tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'* ('+offlineSelection+'&&'+genSelectionVis +')',hist=responseVis_without_fakes, nentries=nEntries) - tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'*'+fakeSelection,hist=responseVis_only_fakes, nentries=nEntries) - - if genVariable_parton != None: - tree.Draw(recoVariable+':'+genVariable_parton,offlineWeight+'*'+offlineSelection,hist=response_parton) - tree.Draw(recoVariable+':'+genVariable_parton,offlineWeight+'* ('+offlineSelection+'&&'+genSelection +')',hist=response_without_fakes_parton, nentries=nEntries) - tree.Draw(recoVariable+':'+genVariable_parton,offlineWeight+'*'+fakeSelection,hist=response_only_fakes_parton, nentries=nEntries) + + print '--->Doing variable :',variable - if options.extraHists: - tree.Draw( 'unfolding.puWeight','unfolding.OfflineSelection',hist=puOffline) + # Output dir name + metSuffix='_patType1CorrectedPFMet' + if variable is 'HT': + metSuffix='' + pass + # Make dir in output file + outputDir = out.mkdir('unfolding_'+variable+'_analyser_'+channel.outputDirName+'_channel'+metSuffix) + + # + # Variable names + # + recoVariable = branchNames[variable] + if variable in ['MET', 'ST', 'WPT']: + if options.sample == "jesup": + recoVariable += '_METUncertainties[2]' + elif options.sample == "jesdown": + recoVariable += '_METUncertainties[3]' + elif options.sample == "jerup": + recoVariable += '_METUncertainties[0]' + elif options.sample == "jerdown": + recoVariable += '_METUncertainties[1]' + elif options.sample in measurement_config.met_systematics: + recoVariable += '_METUncertainties[%i]' % measurement_config.met_systematics[options.sample] + + genVariable_particle = genBranchNames_particle[variable] + genVariable_parton = None + if variable in genBranchNames_parton: + genVariable_parton = genBranchNames_parton[variable] + + # + # Weights and selection + # + + pileupWeight = "PUWeight" + if options.sample == "pileupSystematic": + pileupWeight = "1" + # Generator level + genWeight = '( EventWeight * %.4f)' % ( measurement_config.luminosity_scale) + # genWeight = '( unfolding.puWeight )' + genSelection = '' + genSelectionVis = '' + if channel.channelName is 'muPlusJets' : + genSelection = '( isSemiLeptonicMuon == 1 )' + genSelectionVis = '( isSemiLeptonicMuon == 1 && passesGenEventSelection )' + elif channel.channelName is 'ePlusJets' : + genSelection = '( isSemiLeptonicElectron == 1 )' + genSelectionVis = '( isSemiLeptonicElectron == 1 && passesGenEventSelection )' + + # Offline level + # offlineWeight = '( unfolding.bTagWeight * unfolding.puWeight )' + leptonWeight = 'LeptonEfficiencyCorrection' + if options.sample == 'leptonup': + leptonWeight = 'LeptonEfficiencyCorrectionUp' + elif options.sample == 'leptondown': + leptonWeight == 'LeptonEfficiencyCorrectionDown' + + bjetWeight = "BJetWeight" + if options.sample == "bjetup": + bjetWeight = "BJetUpWeight" + elif options.sample == "bjetdown": + bjetWeight = "BJetDownWeight" + + offlineWeight = '( EventWeight * %s * %s * %s * %.4f)' % ( pileupWeight, bjetWeight, leptonWeight, measurement_config.luminosity_scale ) + offlineSelection = '' + if channel.channelName is 'muPlusJets' : + offlineSelection = '( passSelection == 1 )' + elif channel.channelName is 'ePlusJets' : + offlineSelection = '( passSelection == 2 )' + + # Fake selection + fakeSelection = '( ' + offlineSelection+"&&!"+genSelection +' ) ' + fakeSelectionVis = '( ' + offlineSelection+"&&!"+genSelectionVis +' ) ' + + # Weights derived from variables in tree + if options.applyTopPtReweighting: + ptWeight = topPtWeight( int(options.centreOfMassEnergy) ) + offlineWeight += ' * '+ptWeight + genWeight += ' * '+ptWeight pass - # - # Output histgorams to file - # - outputDir.cd() - truth.Write() - truthVis.Write() - truth_parton.Write() - measured.Write() - measuredVis.Write() - fake.Write() - response.Write() - response_without_fakes.Write() - response_only_fakes.Write() - response_parton.Write() - response_without_fakes_parton.Write() - response_only_fakes_parton.Write() - responseVis_without_fakes.Write() - responseVis_only_fakes.Write() - - phaseSpaceInfoHist.Write() - eventWeight.Write() - if options.extraHists: - puOffline.Write() + + # Apply generator weight + if meWeight >= 0: + genWeight = '( EventWeight * %.4f * genWeight_%s )' % (measurement_config.luminosity_scale, meWeight) + offlineWeight = '( EventWeight * %s * %s * %s * %.4f * genWeight_%s )' % ( pileupWeight, bjetWeight, leptonWeight, measurement_config.luminosity_scale, meWeight ) + nEntries = 1000000 + print 'Changing nEntries to ',nEntries, "<---- DOES NOT SEEM TO DO ANYTHING" + pass + + # Scale factors + # scaleFactor = getScaleFactor( options.centreOfMassEnergy, channel.channelName ) + # scaleFactor = '( unfolding.leptonWeight )' + # offlineWeight += ' * '+scaleFactor + + # + # Histograms to fill + # + # 1D histograms + truth = Hist( allVariablesBins[variable], name='truth') + truthVis = Hist( bin_edges_vis[variable], name='truthVis') + truth_parton = Hist( allVariablesBins[variable], name='truth_parton') + measured = Hist( allVariablesBins[variable], name='measured') + measuredVis = Hist( bin_edges_vis[variable], name='measuredVis') + fake = Hist( allVariablesBins[variable], name='fake') + + # 2D histograms + response = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response') + response_without_fakes = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_without_fakes') + response_only_fakes = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_only_fakes') + + responseVis_without_fakes = Hist2D( bin_edges_vis[variable], bin_edges_vis[variable], name='responseVis_without_fakes') + responseVis_only_fakes = Hist2D( bin_edges_vis[variable], bin_edges_vis[variable], name='responseVis_only_fakes') + + response_parton = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_parton') + response_without_fakes_parton = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_without_fakes_parton') + response_only_fakes_parton = Hist2D( allVariablesBins[variable], allVariablesBins[variable], name='response_only_fakes_parton') + + if options.fineBinned: + minVar = trunc( allVariablesBins[variable][0] ) + maxVar = trunc( max( tree.GetMaximum(genVariable_particle), tree.GetMaximum( recoVariable ) ) * 1.2 ) + nBins = int(maxVar - minVar) + if variable is 'lepton_eta' or variable is 'bjets_eta': + maxVar = 2.5 + minVar = -2.5 + nBins = 1000 + elif 'abs' in variable and 'eta' in variable: + maxVar = 3.0 + minVar = 0. + nBins = 1000 + elif 'Rap' in variable: + maxVar = 3.0 + minVar = -3.0 + nBins = 1000 + elif 'NJets' in variable: + maxVar = 20.5 + minVar = -0.5 + nBins = 21 + + truth = Hist( nBins, minVar, maxVar, name='truth') + truthVis = Hist( nBins, minVar, maxVar, name='truthVis') + truth_parton = Hist( nBins, minVar, maxVar, name='truth_parton') + measured = Hist( nBins, minVar, maxVar, name='measured') + measuredVis = Hist( nBins, minVar, maxVar, name='measuredVis') + fake = Hist( nBins, minVar, maxVar, name='fake') + response = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response') + response_without_fakes = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_without_fakes') + response_only_fakes = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_only_fakes') + responseVis_without_fakes = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='responseVis_without_fakes') + responseVis_only_fakes = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='responseVis_only_fakes') + + response_parton = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_parton') + response_without_fakes_parton = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_without_fakes_parton') + response_only_fakes_parton = Hist2D( nBins, minVar, maxVar, nBins, minVar, maxVar, name='response_only_fakes_parton') + + # Some interesting histograms + puOffline = Hist( 20, 0, 2, name='puWeights_offline') + eventWeightHist = Hist( 100, -2, 2, name='eventWeightHist') + genWeightHist = Hist( 100, -2, 2, name='genWeightHist') + offlineWeightHist = Hist( 100, -2, 2, name='offlineWeightHist') + + phaseSpaceInfoHist = Hist( 10, 0, 1, name='phaseSpaceInfoHist') + if not options.donothing: + nVis = ( tree.Draw( '1', genWeight+'*'+genSelectionVis ) ).Integral() + nVisNotOffline = ( tree.Draw( '1', genWeight+'* ( '+genSelectionVis+'&& !'+offlineSelection+')' ) ).Integral() + nOffline = ( tree.Draw( '1', offlineWeight+'*'+offlineSelection ) ).Integral() + nOfflineNotVis = ( tree.Draw( '1', offlineWeight+'* ( '+offlineSelection+'&& !'+genSelectionVis+')' ) ).Integral() + nFull = ( tree.Draw( '1', genWeight+'*'+genSelection ) ).Integral() + phaseSpaceInfoHist.SetBinContent(1, nVisNotOffline / nVis) + phaseSpaceInfoHist.SetBinContent(2, nOfflineNotVis / nOffline) + phaseSpaceInfoHist.SetBinContent(3, nVis / nFull) + + nOfflineSL = ( tree.Draw( '1', offlineWeight+'* ( '+offlineSelection+'&& '+genSelection+')' ) ).Integral() + nSL = ( tree.Draw( '1', offlineWeight+'* ( '+genSelection+')' ) ).Integral() + # Selection efficiency for SL ttbar + phaseSpaceInfoHist.SetBinContent(4, nOfflineSL / nSL) + # Fraction of offline that are SL + phaseSpaceInfoHist.SetBinContent(5, nOfflineSL / nOffline) + + + # + # Fill histograms + # + if not options.donothing: + # 1D + + tree.Draw(genVariable_particle,genWeight+'*'+genSelection,hist=truth, nentries=nEntries) + tree.Draw(genVariable_particle,genWeight+'*'+genSelectionVis,hist=truthVis, nentries=nEntries) + # if genVariable_parton != None: + # tree.Draw(genVariable_parton,genWeight+'*'+genSelection,hist=truth_parton, nentries=nEntries) + tree.Draw(recoVariable,offlineWeight+'*'+offlineSelection,hist=measured, nentries=nEntries) + tree.Draw(recoVariable,offlineWeight+'*'+offlineSelection,hist=measuredVis, nentries=nEntries) + tree.Draw(recoVariable,offlineWeight+'*'+fakeSelection,hist=fake, nentries=nEntries) + # 2D + tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'*'+offlineSelection,hist=response, nentries=nEntries) + tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'* ('+offlineSelection+'&&'+genSelection +')',hist=response_without_fakes, nentries=nEntries) + tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'*'+fakeSelection,hist=response_only_fakes, nentries=nEntries) + + tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'* ('+offlineSelection+'&&'+genSelectionVis +')',hist=responseVis_without_fakes, nentries=nEntries) + tree.Draw(recoVariable+':'+genVariable_particle,offlineWeight+'*'+fakeSelection,hist=responseVis_only_fakes, nentries=nEntries) + + # if genVariable_parton != None: + # tree.Draw(recoVariable+':'+genVariable_parton,offlineWeight+'*'+offlineSelection,hist=response_parton) + # tree.Draw(recoVariable+':'+genVariable_parton,offlineWeight+'* ('+offlineSelection+'&&'+genSelection +')',hist=response_without_fakes_parton, nentries=nEntries) + # tree.Draw(recoVariable+':'+genVariable_parton,offlineWeight+'*'+fakeSelection,hist=response_only_fakes_parton, nentries=nEntries) + + if options.extraHists: + tree.Draw('EventWeight',genSelection,hist=eventWeightHist) + tree.Draw(genWeight,genSelection,hist=genWeightHist) + tree.Draw( offlineWeight,offlineSelection,hist=offlineWeightHist) + pass + # + # Output histgorams to file + # + outputDir.cd() + truth.Write() + truthVis.Write() + truth_parton.Write() + measured.Write() + measuredVis.Write() + fake.Write() + response.Write() + response_without_fakes.Write() + response_only_fakes.Write() + response_parton.Write() + response_without_fakes_parton.Write() + response_only_fakes_parton.Write() + responseVis_without_fakes.Write() + responseVis_only_fakes.Write() + + phaseSpaceInfoHist.Write() + genWeightHist.Write() + offlineWeightHist.Write() + eventWeightHist.Write() + if options.extraHists: + puOffline.Write() + pass pass pass pass - pass if __name__ == '__main__': main() diff --git a/experimental/BLTUnfold/runJobsCrab.py b/experimental/BLTUnfold/runJobsCrab.py index b98050c2..3270a1b9 100755 --- a/experimental/BLTUnfold/runJobsCrab.py +++ b/experimental/BLTUnfold/runJobsCrab.py @@ -3,7 +3,7 @@ import os jobs = [ - # 13 TeV + # # # 13 TeV '--centreOfMassEnergy 13 -f', '--centreOfMassEnergy 13 -s central', @@ -11,22 +11,51 @@ '--centreOfMassEnergy 13 -s madgraph', '--centreOfMassEnergy 13 -s herwigpp', - # # # PS scale samples - # # '--centreOfMassEnergy 13 -s scaleup', - # # '--centreOfMassEnergy 13 -s scaledown', - # # ME scale weights + # # # # PS scale samples + # # # '--centreOfMassEnergy 13 -s scaleup', + # # # '--centreOfMassEnergy 13 -s scaledown', + # # # ME scale weights '--centreOfMassEnergy 13 --generatorWeight 4', '--centreOfMassEnergy 13 --generatorWeight 8', '--centreOfMassEnergy 13 -s massup', '--centreOfMassEnergy 13 -s massdown', + # '--centreOfMassEnergy 13 -s jesup', + # '--centreOfMassEnergy 13 -s jesdown', + # '--centreOfMassEnergy 13 -s jerup', + # '--centreOfMassEnergy 13 -s jerdown', + + # '--centreOfMassEnergy 13 -s leptonup', + # '--centreOfMassEnergy 13 -s leptondown', + + # '--centreOfMassEnergy 13 -s bjetup', + # '--centreOfMassEnergy 13 -s bjetdown', + + # '--centreOfMassEnergy 13 -s pileupSystematic', + + # '--centreOfMassEnergy 13 -s ElectronEnUp', + # '--centreOfMassEnergy 13 -s ElectronEnDown' , + # '--centreOfMassEnergy 13 -s MuonEnUp', + # '--centreOfMassEnergy 13 -s MuonEnDown', + # '--centreOfMassEnergy 13 -s TauEnUp', + # '--centreOfMassEnergy 13 -s TauEnDown', + # '--centreOfMassEnergy 13 -s UnclusteredEnUp', + # '--centreOfMassEnergy 13 -s UnclusteredEnDown', ] # # # Add pdf variations to list of jobs -for variation in range(9,109): # <- 9 to 108 makes 100 jobs - jobs.append('--centreOfMassEnergy 13 --generatorWeight %i' % variation) +nPDFPerJob = 6 +minPDF = 9 +maxPDF = 109 +variation = 9 +while variation < maxPDF : + nForThisJob = nPDFPerJob + if variation + nPDFPerJob > maxPDF: + nForThisJob = maxPDF - variation + jobs.append('--centreOfMassEnergy 13 --generatorWeight %i --nGeneratorWeights %i' % (variation, nForThisJob) ) + variation += nPDFPerJob pass def parse_args(parameters = []): diff --git a/experimental/add_control_region.py b/experimental/add_control_region.py index 4d19974b..53d6146b 100644 --- a/experimental/add_control_region.py +++ b/experimental/add_control_region.py @@ -22,9 +22,9 @@ def create_new_trees(input_file, suffix = ''): cr1 = 'TTbar_plus_X_analysis/MuPlusJets/QCD 0.12 < iso <= 0.3' cr2 = 'TTbar_plus_X_analysis/MuPlusJets/QCD iso > 0.3' - with root_open(input_file) as f: - t1 = f.Get(tree1_name) - t2 = f.Get(tree2_name) + with root_open(input_file) as file: + t1 = file.Get(tree1_name) + t2 = file.Get(tree2_name) t1.AddFriend(t2) # h1 = t1.Draw('MET', 'relIso_04_deltaBeta > 0.3') @@ -60,10 +60,22 @@ def create_new_trees(input_file, suffix = ''): new_tree2 = f_in.Get(cr2 + '/FitVariables' + suffix).CloneTree() for f in input_files: - shutil.copy(input_folder + f, f) for suffix in ['', '_JERDown', '_JERUp', '_JESDown', '_JESUp']: + fileToUse = f if 'data' in f and not suffix == '': continue - create_new_trees(f, suffix=suffix) - subprocess.call(['hadoop', 'fs', '-rm', output_folder + f]) - subprocess.call(['hadoop', 'fs', '-copyFromLocal', f, output_folder + f]) + + if suffix == '_JERDown': + fileToUse = f.replace('tree','minusJER_tree') + elif suffix == '_JERUp': + fileToUse = f.replace('tree','plusJER_tree') + elif suffix == '_JESDown': + fileToUse = f.replace('tree','minusJES_tree') + elif suffix == '_JESUp': + fileToUse = f.replace('tree','plusJES_tree') + + shutil.copy(input_folder + fileToUse, fileToUse) + + create_new_trees(fileToUse, suffix=suffix) + subprocess.call(['hadoop', 'fs', '-rm', output_folder + fileToUse]) + subprocess.call(['hadoop', 'fs', '-copyFromLocal', fileToUse, output_folder + fileToUse]) diff --git a/src/cross_section_measurement/01_get_ttjet_normalisation.py b/src/cross_section_measurement/01_get_ttjet_normalisation.py index 00393de1..b857d206 100644 --- a/src/cross_section_measurement/01_get_ttjet_normalisation.py +++ b/src/cross_section_measurement/01_get_ttjet_normalisation.py @@ -24,7 +24,7 @@ from src.cross_section_measurement.lib import closure_tests from tools.file_utilities import write_data_to_JSON from tools.hist_utilities import clean_control_region, \ - hist_to_value_error_tuplelist + hist_to_value_error_tuplelist, fix_overflow import glob import tools.measurement @@ -103,7 +103,7 @@ def calculate_normalisation(self): for sample, hist in histograms.items(): # TODO: this should be a list of bin-contents self.initial_normalisation[ - sample] = hist_to_value_error_tuplelist(hist) + sample] = hist_to_value_error_tuplelist(fix_overflow(hist)) if self.method == self.BACKGROUND_SUBTRACTION and sample != 'TTJet': self.normalisation[sample] = self.initial_normalisation[sample] diff --git a/src/cross_section_measurement/02_unfold_and_measure.py b/src/cross_section_measurement/02_unfold_and_measure.py index 7de0988d..30e80539 100644 --- a/src/cross_section_measurement/02_unfold_and_measure.py +++ b/src/cross_section_measurement/02_unfold_and_measure.py @@ -33,7 +33,6 @@ def unfold_results( results, category, channel, tau_value, h_truth, h_measured, unfoldCfg.Hreco = options.Hreco h_unfolded_data = unfolding.unfold( h_data ) - del unfolding return hist_to_value_error_tuplelist( h_unfolded_data ) @@ -63,13 +62,40 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, tau_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, + + # 'JES_down' : file_for_jesdown, + # 'JES_up' : file_for_jesup, + + # 'JER_down' : file_for_jerdown, + # 'JER_up' : file_for_jerup, + + # 'BJet_up' : file_for_bjetdown, + # 'BJet_down' : file_for_bjetup, + ttbar_theory_systematic_prefix + 'hadronisation' : file_for_powheg_herwig, ttbar_theory_systematic_prefix + 'NLOgenerator' : file_for_amcatnlo, + + # 'ElectronEnUp' : file_for_ElectronEnUp, + # 'ElectronEnDown' : file_for_ElectronEnDown, + # 'MuonEnUp' : file_for_MuonEnUp, + # 'MuonEnDown' : file_for_MuonEnDown, + # 'TauEnUp' : file_for_TauEnUp, + # 'TauEnDown' : file_for_TauEnDown, + # 'UnclusteredEnUp' : file_for_UnclusteredEnUp, + # 'UnclusteredEnDown' : file_for_UnclusteredEnDown, + + # 'Muon_up' : file_for_LeptonUp, + # 'Muon_down' : file_for_LeptonDown, + # 'Electron_up' : file_for_LeptonUp, + # 'Electron_down' : file_for_LeptonDown, + + # 'PileUpSystematic' : file_for_PUSystematic, } h_truth, h_measured, h_response, h_fakes = None, None, None, None # Systematics where you change the response matrix - if category in ttbar_generator_systematics : + if category in ttbar_generator_systematics or category in files_for_systematics : + print 'Doing category',category,'by changing response matrix' h_truth, h_measured, h_response, h_fakes = get_unfold_histogram_tuple( inputfile = files_for_systematics[category], variable = variable, channel = channel, @@ -104,7 +130,7 @@ def get_unfolded_normalisation( TTJet_fit_results, category, channel, tau_value, visiblePS = visiblePS, ) - central_results = hist_to_value_error_tuplelist( h_truth ) + # central_results = hist_to_value_error_tuplelist( h_truth ) TTJet_fit_results_unfolded = unfold_results( TTJet_fit_results, category, channel, @@ -403,7 +429,31 @@ def calculate_normalised_xsections( normalisation, category, channel, normalise_ ### file_for_massdown = File( measurement_config.unfolding_mass_down, 'read' ) file_for_massup = File( measurement_config.unfolding_mass_up, 'read' ) + file_for_jesdown = File( measurement_config.unfolding_jes_down, 'read' ) + file_for_jesup = File( measurement_config.unfolding_jes_up, 'read' ) + ### + file_for_jerdown = File( measurement_config.unfolding_jer_down, 'read' ) + file_for_jerup = File( measurement_config.unfolding_jer_up, 'read' ) ### + file_for_bjetdown = File( measurement_config.unfolding_bjet_down, 'read' ) + file_for_bjetup = File( measurement_config.unfolding_bjet_up, 'read' ) + ### + file_for_LeptonDown = File( measurement_config.unfolding_Lepton_down, 'read' ) + file_for_LeptonUp = File( measurement_config.unfolding_Lepton_up, 'read' ) + ### + file_for_ElectronEnDown = File( measurement_config.unfolding_ElectronEn_down, 'read' ) + file_for_ElectronEnUp = File( measurement_config.unfolding_ElectronEn_up, 'read' ) + ### + file_for_MuonEnDown = File( measurement_config.unfolding_MuonEn_down, 'read' ) + file_for_MuonEnUp = File( measurement_config.unfolding_MuonEn_up, 'read' ) + ### + file_for_TauEnDown = File( measurement_config.unfolding_TauEn_down, 'read' ) + file_for_TauEnUp = File( measurement_config.unfolding_TauEn_up, 'read' ) + ### + file_for_UnclusteredEnDown = File( measurement_config.unfolding_UnclusteredEn_down, 'read' ) + file_for_UnclusteredEnUp = File( measurement_config.unfolding_UnclusteredEn_up, 'read' ) + ### + file_for_PUSystematic = File( measurement_config.unfolding_PUSystematic, 'read') file_for_powhegPythia8 = File( measurement_config.unfolding_powheg_pythia8, 'read') file_for_amcatnlo = File( measurement_config.unfolding_amcatnlo, 'read') @@ -449,14 +499,14 @@ def calculate_normalised_xsections( normalisation, category, channel, normalise_ # all MET uncertainties except JES as this is already included met_uncertainties = [suffix for suffix in measurement_config.met_systematics_suffixes if not 'JetEn' in suffix and not 'JetRes' in suffix] all_measurements = deepcopy( categories ) - # all_measurements.extend( pdf_uncertainties ) + all_measurements.extend( pdf_uncertainties ) all_measurements.extend( ['QCD_shape'] ) all_measurements.extend( rate_changing_systematics ) + all_measurements.extend( ['central_TTJet'] ) print 'Performing unfolding for variable', variable for category in all_measurements: - print 'Doing category ',category - + print 'Doing category ', category if run_just_central and not category == 'central': continue # Don't need to consider MET uncertainties for HT @@ -483,20 +533,38 @@ def calculate_normalised_xsections( normalisation, category, channel, normalise_ # combined_file = path_to_JSON + '/fit_results/' + category + '/fit_results_combined_' + met_type + '.txt' # don't change fit input for ttbar generator/theory systematics and PDF weights - if category in ttbar_generator_systematics or category in ttbar_theory_systematics or category in pdf_uncertainties: + if category in ttbar_generator_systematics or category in pdf_uncertainties: # or category in ttbar_mass_systematics electron_file = path_to_JSON + '/central/normalisation_electron_' + met_type + '.txt' muon_file = path_to_JSON + '/central/normalisation_muon_' + met_type + '.txt' - # combined_file = path_to_JSON + '/central/normalisation_combined_' + met_type + '.txt' + # combined_file = path_to_JSON + '/central/normalisation_combined_' + met_type + '.txt' + elif category in rate_changing_systematics or category == 'QCD_shape': + electron_file = path_to_JSON + '/' + category + '/normalisation_electron_' + met_type + '.txt' + muon_file = path_to_JSON + '/' + category + '/normalisation_muon_' + met_type + '.txt' + elif category == 'central_TTJet': + electron_file = path_to_JSON + '/central/initial_normalisation_electron_' + met_type + '.txt' + muon_file = path_to_JSON + '/central/initial_normalisation_muon_' + met_type + '.txt' + # elif category in met_uncertainties and not 'JES' in category and not 'JER' in category: + # electron_file = path_to_JSON + '/'+category+'/initial_normalisation_electron_' + met_type + '.txt' + # muon_file = path_to_JSON + '/'+category+'/initial_normalisation_muon_' + met_type + '.txt' + elif category != 'central': + # electron_file = path_to_JSON + '/'+category+'/initial_normalisation_electron_' + met_type + '.txt' + # muon_file = path_to_JSON + '/'+category+'/initial_normalisation_muon_' + met_type + '.txt' + # electron_file = path_to_JSON + '/central/normalisation_electron_' + translate_options[options.metType] + '.txt' + # muon_file = path_to_JSON + '/central/normalisation_muon_' + translate_options[options.metType] + '.txt' + electron_file = path_to_JSON + '/' + category + '/normalisation_electron_' + met_type + '.txt' + muon_file = path_to_JSON + '/' + category + '/normalisation_muon_' + met_type + '.txt' fit_results_electron = None fit_results_muon = None if category == 'Muon_up' or category == 'Muon_down': + # fit_results_electron = read_data_from_JSON( path_to_JSON + '/central/initial_normalisation_electron_' + met_type + '.txt' ) fit_results_electron = read_data_from_JSON( path_to_JSON + '/central/normalisation_electron_' + met_type + '.txt' ) fit_results_muon = read_data_from_JSON( muon_file ) elif category == 'Electron_up' or category == 'Electron_down': fit_results_electron = read_data_from_JSON( electron_file ) + # fit_results_muon = read_data_from_JSON( path_to_JSON + '/central/initial_normalisation_muon_' + met_type + '.txt' ) fit_results_muon = read_data_from_JSON( path_to_JSON + '/central/normalisation_muon_' + met_type + '.txt' ) else: fit_results_electron = read_data_from_JSON( electron_file ) diff --git a/src/cross_section_measurement/03_calculate_systematics.py b/src/cross_section_measurement/03_calculate_systematics.py index 45bce67e..b5d061bb 100644 --- a/src/cross_section_measurement/03_calculate_systematics.py +++ b/src/cross_section_measurement/03_calculate_systematics.py @@ -65,7 +65,7 @@ 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, experimentalUncertainty = False, actualCentralMeasurements = [] ): global symmetrise_errors # number of bins number_of_bins = len( list_of_central_measurements ) @@ -103,10 +103,15 @@ def summarise_systematics( list_of_central_measurements, dictionary_of_systemati 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 ) + elif experimentalUncertainty: + 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 ) + actualCentralValue = actualCentralMeasurements[bin_i][0] + error_down = error_down / central_value * actualCentralValue + error_up = error_up / central_value * actualCentralValue 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 ) - down_errors[bin_i] = error_down up_errors[bin_i] = error_up @@ -234,15 +239,17 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio rate_changing_systematics_list = [systematic for systematic in measurement_config.rate_changing_systematics_names] # all other uncertainties (including JES and JER) - other_uncertainties_list = deepcopy( measurement_config.categories_and_prefixes.keys() ) + experimental_uncertainties_list = deepcopy( measurement_config.categories_and_prefixes.keys() ) + ### other_uncertainties_list.extend( vjets_generator_systematics_list ) - other_uncertainties_list.append( 'QCD_shape' ) + other_uncertainties_list = [ 'QCD_shape' ] other_uncertainties_list.extend( rate_changing_systematics_list ) for channel in ['electron', 'muon', 'combined']: # print("channel = ", channel) # read central measurement central_measurement, central_measurement_unfolded = read_normalised_xsection_measurement( 'central', channel ) + # central_TTJet_measurement, central_TTJet_measurement_unfolded = read_normalised_xsection_measurement( 'central_TTJet', channel ) # read groups of systematics ttbar_generator_systematics, ttbar_generator_systematics_unfolded = read_normalised_xsection_systematics( list_of_systematics = ttbar_generator_systematics_list, channel = channel ) @@ -254,12 +261,14 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio # kValue_systematic, kValue_systematic_unfolded = read_normalised_xsection_systematics( list_of_systematics = kValue_systematic_list, channel = channel ) pdf_systematics, pdf_systematics_unfolded = read_normalised_xsection_systematics( list_of_systematics = pdf_uncertainties, channel = channel ) # met_systematics, met_systematics_unfolded = read_normalised_xsection_systematics( list_of_systematics = met_uncertainties_list, channel = channel ) + experimental_systematics, experimental_systematics_unfolded = read_normalised_xsection_systematics( list_of_systematics = experimental_uncertainties_list, channel = channel ) 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 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 ) @@ -281,6 +290,14 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio # 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 ) + + # experimental_min, experimental_max = summarise_systematics( central_TTJet_measurement, experimental_systematics, experimentalUncertainty = True, actualCentralMeasurements = central_measurement) + # experimental_min_unfolded, experimental_max_unfolded = summarise_systematics( central_TTJet_measurement_unfolded, experimental_systematics_unfolded, experimentalUncertainty = True, actualCentralMeasurements = central_measurement_unfolded) + + experimental_min, experimental_max = summarise_systematics( central_measurement, experimental_systematics ) + experimental_min_unfolded, experimental_max_unfolded = summarise_systematics( central_measurement_unfolded, experimental_systematics_unfolded ) + + # 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 ) @@ -292,7 +309,8 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio ### ttbar_mass_min, # ### kValue_min, pdf_min, - ### met_min, + ### met_min, + experimental_min, other_min], [ ttbar_generator_max, #ttbar_ptreweight_max, @@ -301,6 +319,7 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio # ### kValue_max, pdf_max, ### met_max, + experimental_max, other_max] ) ### central_measurement_with_systematics_but_without_ttbar_theory = get_measurement_with_lower_and_upper_errors( central_measurement, ### [pdf_min, met_min, other_min, @@ -319,6 +338,7 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio # ### kValue_min, pdf_min, ### met_min, + experimental_min, other_min], [ # ttbar_hadronisation_max, @@ -326,6 +346,7 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio ###ttbar_mass_max, # ### kValue_max, pdf_max, + experimental_max, ### met_max, other_max] ) @@ -333,10 +354,12 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio [ ttbar_generator_min, pdf_min, - other_min], + experimental_min, + other_min,], [ ttbar_generator_max, pdf_max, + experimental_max, other_max], systematicErrorsOnly = True ) @@ -348,13 +371,15 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio # ### kValue_min_unfolded, pdf_min_unfolded, ### met_min_unfolded, + experimental_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, + pdf_max_unfolded, + experimental_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, @@ -375,7 +400,8 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio ###ttbar_mass_min_unfolded, # ### kValue_min_unfolded, pdf_min_unfolded, - ### met_min_unfolded, + ### met_min_unfolded, + experimental_min_unfolded, other_min_unfolded], [ # ttbar_hadronisation_max_unfolded, @@ -383,6 +409,7 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio ###ttbar_mass_max_unfolded, # ### kValue_max_unfolded, pdf_max_unfolded, + experimental_max_unfolded, ### met_max_unfolded, other_max_unfolded] ) @@ -390,10 +417,12 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio [ ttbar_generator_min_unfolded, pdf_min_unfolded, + experimental_min_unfolded, other_min_unfolded], [ ttbar_generator_max_unfolded, pdf_max_unfolded, + experimental_max_unfolded, other_max_unfolded], systematicErrorsOnly = True ) @@ -422,6 +451,7 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio ###ttbar_mass_systematic = replace_measurement_with_deviation_from_central( central_measurement, ttbar_mass_systematic ) ###kValue_systematic = replace_measurement_with_deviation_from_central( central_measurement, kValue_systematic ) other_systematics = replace_measurement_with_deviation_from_central( central_measurement, other_systematics ) + experimental_systematics = replace_measurement_with_deviation_from_central( central_measurement, experimental_systematics ) ttbar_generator_systematics_unfolded = replace_measurement_with_deviation_from_central( central_measurement_unfolded, ttbar_generator_systematics_unfolded ) pdf_systematics_unfolded = replace_measurement_with_deviation_from_central( central_measurement_unfolded, pdf_systematics_unfolded ) @@ -429,7 +459,8 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio ###ttbar_mass_systematic_unfolded = replace_measurement_with_deviation_from_central( central_measurement_unfolded, ttbar_mass_systematic_unfolded ) ###kValue_systematic_unfolded = replace_measurement_with_deviation_from_central( central_measurement_unfolded, kValue_systematic_unfolded ) other_systematics_unfolded = replace_measurement_with_deviation_from_central( central_measurement_unfolded, other_systematics_unfolded ) - + experimental_systematics_unfolded = replace_measurement_with_deviation_from_central( central_measurement_unfolded, experimental_systematics_unfolded ) + # Scale mass systematic ###ttbar_mass_systematic['TTJets_massdown'], ttbar_mass_systematic['TTJets_massup'] = scaleTopMassSystematicErrors( ttbar_mass_systematic['TTJets_massdown'], ttbar_mass_systematic['TTJets_massup'] ) ###ttbar_mass_systematic_unfolded['TTJets_massdown'], ttbar_mass_systematic_unfolded['TTJets_massup'] = scaleTopMassSystematicErrors( ttbar_mass_systematic_unfolded['TTJets_massdown'], ttbar_mass_systematic_unfolded['TTJets_massup'] ) @@ -454,6 +485,9 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio other_systematics['total_lower'], other_systematics['total_upper'] = other_min, other_max other_systematics_unfolded['total_lower'], other_systematics_unfolded['total_upper'] = other_min_unfolded, other_max_unfolded + experimental_systematics['total_lower'], experimental_systematics['total_upper'] = experimental_min, experimental_max + experimental_systematics_unfolded['total_lower'], experimental_systematics_unfolded['total_upper'] = experimental_min_unfolded, experimental_max_unfolded + #### print 'Generator',ttbar_generator_systematics_unfolded #### print 'k Value',kValue_systematic_unfolded ###new_systematics = {} @@ -471,5 +505,6 @@ def replace_measurement_with_deviation_from_central( central_measurement, dictio ###write_normalised_xsection_measurement( ttbar_mass_systematic, ttbar_mass_systematic_unfolded, channel, summary = 'topMass' ) ###write_normalised_xsection_measurement( kValue_systematic, kValue_systematic_unfolded, channel, summary = 'kValue' ) write_normalised_xsection_measurement( other_systematics, other_systematics_unfolded, channel, summary = 'other' ) + write_normalised_xsection_measurement( experimental_systematics, experimental_systematics_unfolded, channel, summary = 'experimental' ) ###write_normalised_xsection_measurement( new_systematics, new_systematics_unfolded, channel, summary = 'new' ) diff --git a/src/cross_section_measurement/04_make_plots_matplotlib.py b/src/cross_section_measurement/04_make_plots_matplotlib.py index 0d3691ed..87ac587d 100644 --- a/src/cross_section_measurement/04_make_plots_matplotlib.py +++ b/src/cross_section_measurement/04_make_plots_matplotlib.py @@ -250,7 +250,7 @@ def make_template_plots( histograms, category, channel ): plt.text(0.95, 0.95, r"\textbf{CMS}", transform=axes.transAxes, fontsize=42, verticalalignment='top',horizontalalignment='right') # channel text - axes.text(0.95, 0.90, r"\emph{%s}" %channel_label, transform=axes.transAxes, fontsize=40, + axes.text(0.95, 0.95, r"\emph{%s}" %channel_label, transform=axes.transAxes, fontsize=40, verticalalignment='top',horizontalalignment='right') plt.tight_layout() @@ -406,20 +406,28 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True data_handle = handles[data_label_index] labels.remove( 'data' ) handles.remove( data_handle ) - labels.insert( 0, 'unfolded data' ) + labels.insert( 0, 'data' ) handles.insert( 0, data_handle ) new_handles, new_labels = [], [] - for handle, label in zip( handles, labels ): - if not label == 'do_not_show': - new_handles.append( handle ) - new_labels.append( label ) + zipped = dict( zip( labels, handles ) ) + labelOrder = ['data', 'Powheg Pythia8', 'Powheg Herwig++', 'aMC@NLO', 'Madgraph', '$Q^{2}$ up', '$Q^{2}$ down', 'Top mass up', 'Top mass down'] + for label in labelOrder: + if label in labels: + new_handles.append(zipped[label]) + new_labels.append(label) - legend_location = (0.97, 0.88) + print (new_labels) + print (new_handles) + legend_location = (0.97, 0.82) if variable == 'MT': - legend_location = (0.05, 0.88) + legend_location = (0.05, 0.82) elif variable == 'ST': - legend_location = (0.90, 0.88) + legend_location = (0.97, 0.82) + elif variable == 'WPT': + legend_location = (1.0, 0.84) + elif variable == 'abs_lepton_eta': + legend_location = (1.0, 0.94) plt.legend( new_handles, new_labels, numpoints = 1, prop = CMS.legend_properties, frameon = False, bbox_to_anchor=legend_location, bbox_transform=plt.gcf().transFigure ) label, channel_label = get_cms_labels( channel ) @@ -429,6 +437,15 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True # note: fontweight/weight does not change anything as we use Latex text!!! logo_location = (0.05, 0.98) prelim_location = (0.05, 0.92) + channel_location = ( 0.05, 0.86) + if variable == 'WPT': + logo_location = (0.03, 0.98) + prelim_location = (0.03, 0.92) + channel_location = (0.03, 0.86) + elif variable == 'abs_lepton_eta': + logo_location = (0.03, 0.98) + prelim_location = (0.03, 0.92) + channel_location = (0.03, 0.86) plt.text(logo_location[0], logo_location[1], r"\textbf{CMS}", transform=axes.transAxes, fontsize=42, verticalalignment='top',horizontalalignment='left') # preliminary @@ -436,12 +453,17 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True transform=axes.transAxes, fontsize=42, verticalalignment='top',horizontalalignment='left') # channel text - axes.text(0.95, 0.98, r"\emph{%s}" %channel_label, transform=axes.transAxes, fontsize=40, - verticalalignment='top',horizontalalignment='right') + plt.text(channel_location[0], channel_location[1], r"\emph{%s}" %channel_label, transform=axes.transAxes, fontsize=40, + verticalalignment='top',horizontalalignment='left') ylim = axes.get_ylim() if ylim[0] < 0: axes.set_ylim( ymin = 0.) - axes.set_ylim(ymax = ylim[1]*1.2) + if variable == 'WPT': + axes.set_ylim(ymax = ylim[1]*1.3) + elif variable == 'abs_lepton_eta': + axes.set_ylim(ymax = ylim[1]*1.3) + else : + axes.set_ylim(ymax = ylim[1]*1.2) if show_ratio: @@ -493,39 +515,50 @@ def make_plots( histograms, category, output_folder, histname, show_ratio = True if category == 'central': rplt.fill_between( syst_lower, syst_upper, ax1, - color = 'yellow', alpha = 0.5 ) + color = 'yellow' ) rplt.fill_between( stat_upper, stat_lower, ax1, color = '0.75', - alpha = 0.5 ) + ) # 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}) + p_stat = mpatches.Patch(facecolor='0.75', label='Stat.', edgecolor='black' ) + p_stat_and_syst = mpatches.Patch(facecolor='yellow', label=r'Stat. $\oplus$ Syst.', edgecolor='black' ) + l1 = ax1.legend(handles = [p_stat, p_stat_and_syst], loc = 'upper left', + frameon = False, prop = {'size':26}, ncol = 2) - ax1.legend(handles = [p_stat_and_syst], loc = 'lower left', - frameon = False, prop = {'size':30}) + # ax1.legend(handles = [p_stat_and_syst], loc = 'lower left', + # frameon = False, prop = {'size':30}) ax1.add_artist(l1) if variable == 'MET': - ax1.set_ylim( ymin = 0.0, ymax = 2.5 ) + ax1.set_ylim( ymin = 0.5, ymax = 2.4 ) ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) ) # ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) 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., ymax = 2 ) + elif variable == 'HT': + ax1.set_ylim( ymin = 0.5, ymax = 1.8 ) + ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) ) + ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + elif variable == 'ST': + ax1.set_ylim( ymin = 0.5, ymax = 1.8 ) ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) ) ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) elif variable == 'WPT': - ax1.set_ylim( ymin = 0., ymax = 2 ) + ax1.set_ylim( ymin = 0.5, ymax = 1.85 ) ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) ) ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) elif variable == 'NJets': - ax1.set_ylim( ymin = 0.0, ymax = 2.5 ) - + ax1.set_ylim( ymin = 0.5, ymax = 2.5 ) + elif variable == 'abs_lepton_eta': + ax1.set_ylim( ymin = 0.5, ymax = 1.6 ) + ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) ) + ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) + elif variable == 'lepton_pt': + ax1.set_ylim( ymin = 0.5, ymax = 1.8 ) + ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) ) + ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) ) if CMS.tight_layout: plt.tight_layout() @@ -690,6 +723,7 @@ def get_unit_string(fit_variable): # all_measurements.extend( new_uncertainties ) all_measurements.extend( rate_changing_systematics ) for channel in ['electron', 'muon', 'combined']: + # for channel in ['combined']: for category in all_measurements: if not category == 'central' and not options.additional_plots: continue diff --git a/src/cross_section_measurement/05_make_tables.py b/src/cross_section_measurement/05_make_tables.py index 8f91c1e5..dc892c41 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, samples_latex, typical_systematics_latex, met_systematics_latex -from config.variable_binning import variable_bins_latex, variable_bins_ROOT, variable_bins_visiblePS_ROOT, variable_bins_visiblePS_latex +from config.latex_labels import variables_latex, variables_NonLatex, measurements_latex, samples_latex, typical_systematics_latex, met_systematics_latex +from config.variable_binning import variable_bins_latex, variable_bins_ROOT, variable_bins_visiblePS_ROOT, variable_bins_visiblePS_latex, bin_edges_vis, bin_edges from config import XSectionConfig from tools.Calculation import getRelativeError from tools.file_utilities import read_data_from_JSON, make_folder_if_not_exists +from tools.hist_utilities import values_and_errors_to_hist from lib import read_normalisation, read_initial_normalisation import math import os.path from numpy import median +import matplotlib as mpl +mpl.use( 'agg' ) +import matplotlib.pyplot as plt +import rootpy.plotting.root2matplotlib as rplt +from config import CMS +import matplotlib.cm as cm +# use full stpectrum, yet use white for less than vmin=1 events +my_cmap = cm.get_cmap( 'jet' ) +my_cmap.set_under( 'w' ) +from matplotlib import rc +rc( 'font', **CMS.font ) +rc( 'text', usetex = False ) def read_xsection_measurement_results_with_errors(channel): global path_to_JSON, variable, met_type, phase_space @@ -40,6 +53,9 @@ def read_xsection_measurement_results_with_errors(channel): file_name = file_template.replace('.txt', '_PDF_errors.txt') normalised_xsection_PDF_errors = read_data_from_JSON( file_name ) + file_name = file_template.replace('.txt', '_experimental_errors.txt') + normalised_xsection_experimental_errors = read_data_from_JSON( file_name ) + file_name = file_template.replace('.txt', '_other_errors.txt') normalised_xsection_other_errors = read_data_from_JSON( file_name ) @@ -59,7 +75,7 @@ def read_xsection_measurement_results_with_errors(channel): # normalised_xsection_measured_errors.update(normalised_xsection_MET_errors['TTJet_measured']) # normalised_xsection_measured_errors.update(normalised_xsection_topMass_errors['TTJet_measured']) ### normalised_xsection_measured_errors.update(normalised_xsection_kValue_errors['TTJet_measured']) - normalised_xsection_measured_errors.update(normalised_xsection_other_errors['TTJet_measured']) + normalised_xsection_measured_errors.update(normalised_xsection_experimental_errors['TTJet_measured']) # normalised_xsection_measured_errors.update(normalised_xsection_new_errors['TTJet_measured']) normalised_xsection_unfolded_errors = normalised_xsection_other_errors['TTJet_unfolded'] @@ -68,7 +84,7 @@ def read_xsection_measurement_results_with_errors(channel): # normalised_xsection_unfolded_errors.update(normalised_xsection_MET_errors['TTJet_unfolded']) # normalised_xsection_unfolded_errors.update(normalised_xsection_topMass_errors['TTJet_unfolded']) ### normalised_xsection_unfolded_errors.update(normalised_xsection_kValue_errors['TTJet_unfolded']) - normalised_xsection_unfolded_errors.update(normalised_xsection_other_errors['TTJet_unfolded']) + normalised_xsection_unfolded_errors.update(normalised_xsection_experimental_errors['TTJet_unfolded']) # normalised_xsection_unfolded_errors.update(normalised_xsection_new_errors['TTJet_unfolded']) return normalised_xsection_measured_unfolded, normalised_xsection_measured_errors, normalised_xsection_unfolded_errors @@ -322,18 +338,102 @@ def print_xsections(xsections, channel, toFile = True, print_before_unfolding = else: print printout +def make_error_plot( errorHists, bins ): + global output_folder, variable + # For each up/down source, reduce to one set of numbers + symmetricErrorHists = {} + for source, hist in errorHists.iteritems(): + if ( variable == 'HT' or variable == 'NJets' or variable == 'lepton_pt' or variable == 'abs_lepton_eta' ) and source in measurement_config.met_systematics and not 'JES' in source and not 'JER' in source: + continue + + if 'down' in source or '-' in source or 'lower' in source or 'Down' in source: + # Find up version + upHist = None + newSource = '' + if 'down' in source: + upHist = errorHists[source.replace('down','up')] + newSource = source.replace('down','') + elif 'Down' in source: + upHist = errorHists[source.replace('Down','Up')] + newSource = source.replace('Down','') + elif '-' in source: + upHist = errorHists[source.replace('-','+')] + newSource = source.replace('-','') + elif 'lower' in source: + upHist = errorHists[source.replace('lower','upper')] + newSource = source.replace('lower','') + + if newSource[-1] == '_': + newSource = newSource[:-1] + # if '_' in newSource: + # newSource = newSource.replace('_','') + + symmetricErrorHists[newSource] = [] + for errorup, errordown in zip(hist, upHist): + newError = max( abs(errorup), abs(errordown) ) + symmetricErrorHists[newSource].append(newError) + elif 'TTJets_hadronisation' in source or 'QCD_shape' in source or 'TTJets_NLOgenerator' in source: + symmetricErrorHists[source] = [ abs(i) for i in hist ] + + x_limits = [bins[0], bins[-1]] + y_limits = [0,0.6] + plt.figure( figsize = ( 20, 16 ), dpi = 200, facecolor = 'white' ) + + ax0 = plt.axes() + ax0.minorticks_on() + ax0.xaxis.labelpad = 12 + ax0.yaxis.labelpad = 12 + ax0.set_xlim( x_limits ) + plt.tick_params( **CMS.axis_label_major ) + plt.tick_params( **CMS.axis_label_minor ) + + + statisticalErrorHists = values_and_errors_to_hist( errorHists['statistical'], [], bins ) + for source, hist in symmetricErrorHists.iteritems(): + symmetricErrorHists[source] = values_and_errors_to_hist( hist, [], bins ) + + colours = ['silver', 'r', 'tan', 'chartreuse', 'cadetblue', 'dodgerblue', 'pink', 'hotpink', 'coral', 'forestgreen', 'cyan', 'teal', 'crimson', 'darkmagenta', 'olive', 'slateblue', 'deepskyblue', 'orange', 'r' ] + for source, colour in zip( symmetricErrorHists.keys(), colours): + hist = symmetricErrorHists[source] + hist.linewidth = 4 + hist.color = colour + rplt.hist( hist, stacked=False, axes = ax0, cmap = my_cmap, vmin = 1, label = source ) + + statisticalErrorHists.linewidth = 4 + statisticalErrorHists.color = 'black' + statisticalErrorHists.linestyle = 'dashed' + rplt.hist( statisticalErrorHists, stacked=False, axes = ax0, cmap = my_cmap, vmin = 1, label = 'stat.' ) + + ax0.set_ylim( y_limits ) + leg = plt.legend(loc=1,prop={'size':40},ncol=2) + leg.draw_frame(False) + x_title = variables_NonLatex[variable] + if variable in ['HT', 'MET', 'WPT', 'ST', 'lepton_pt']: + x_title += ' [GeV]' + plt.xlabel( x_title, CMS.x_axis_title ) + plt.ylabel( 'Relative Uncertainty', CMS.y_axis_title) + plt.tight_layout() + + path = output_folder + '/' + variable + '/' + make_folder_if_not_exists(path) + file_template = path + '/%s_systematics_%dTeV_%s.pdf' % (variable, measurement_config.centre_of_mass_energy, channel) + plt.savefig(file_template) + pass + def print_error_table(central_values, errors, channel, toFile = True, print_before_unfolding = False): global output_folder, variable, met_type, b_tag_bin, all_measurements, phase_space bins = None bins_latex = None + binEdges = None variable_latex = variables_latex[variable] if phase_space == 'VisiblePS': bins = variable_bins_visiblePS_ROOT[variable] bins_latex = variable_bins_visiblePS_latex[variable] + binEdges = bin_edges_vis[variable] elif phase_space == 'FullPS': bins = variable_bins_ROOT[variable] bins_latex = variable_bins_latex[variable] - + binEdges = bin_edges[variable] printout = '%% ' + '=' * 60 printout += '\n' printout += '%% Systematics table for %s variable, %s channel, met type %s, %s b-tag region\n' % (variable, channel, met_type, b_tag_bin) @@ -367,6 +467,11 @@ def print_error_table(central_values, errors, channel, toFile = True, print_befo else: assert(len(bins) == len(central_values['unfolded'])) + errorHists = {} + errorHists['statistical'] = [] + for source in all_measurements: + errorHists[source] = [] + for bin_i, variable_bin in enumerate(bins): header += '& %s' % (bins_latex[variable_bin]) if print_before_unfolding: @@ -380,6 +485,9 @@ def print_error_table(central_values, errors, channel, toFile = True, print_befo abs_error = errors[source][bin_i] relative_error = getRelativeError(central_value, abs_error) + + errorHists[source].append(relative_error) + text = '%.2f' % (relative_error*100) if rows.has_key(source): rows[source].append(text) @@ -411,9 +519,13 @@ def print_error_table(central_values, errors, channel, toFile = True, print_befo else: value, error = central_values['unfolded'][bin_i] relativeError = getRelativeError(value, error) + errorHists['statistical'].append(relativeError) total_line += ' & %.2f ' % (relativeError * 100) printout += total_line + '\\\\ \n' + if not print_before_unfolding: + make_error_plot( errorHists, binEdges ) + #append the total systematic error to the table total_line = 'Total Sys. (\%)' for bin_i, variable_bin in enumerate(bins): @@ -640,6 +752,7 @@ def print_typical_systematics_table(central_values, errors, channel, toFile = Tr all_measurements.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) diff --git a/src/cross_section_measurement/compareQCDControlRegions.py b/src/cross_section_measurement/compareQCDControlRegions.py new file mode 100644 index 00000000..8d69c255 --- /dev/null +++ b/src/cross_section_measurement/compareQCDControlRegions.py @@ -0,0 +1,577 @@ +from optparse import OptionParser +from config.latex_labels import b_tag_bins_latex, samples_latex, channel_latex, \ + variables_latex, fit_variables_latex, control_plots_latex +from config.variable_binning import bin_edges, fit_variable_bin_edges, control_plots_bins +from config.histogram_colours import histogram_colours as colours +from config import XSectionConfig +from tools.file_utilities import read_data_from_JSON, make_folder_if_not_exists +from tools.plotting import make_data_mc_comparison_plot, Histogram_properties, \ +make_control_region_comparison +from rootpy.plotting import Hist +from tools.hist_utilities import prepare_histograms, clean_control_region, get_normalisation_error, get_fitted_normalisation +from tools.ROOT_utils import get_histograms_from_trees, set_root_defaults +from tools.latex import setup_matplotlib +from uncertainties import ufloat + +# latex, font, etc +setup_matplotlib() + +title_template = '$%.0f$ pb$^{-1}$ (%d TeV)' + +def binWidth(binning): + return ( binning[-1] - binning[0] ) / ( len(binning)-1 ) + + +def getPUWeights(histograms_to_draw, histogram_lables) : + # hists = dict(zip(histogram_lables, histograms_to_draw)) + hists = {} + dataHist = None + mcHist = None + for label, histogram in zip(histogram_lables, histograms_to_draw): + if label == 'data': + dataHist = histogram.Clone() + else : + if mcHist == None: + mcHist = histogram.Clone() + else: + mcHist += histogram + dataValues = list(dataHist.y()) + mcValues = list(mcHist.y()) + + weights = [ data / mc for data, mc in zip(dataValues, mcValues)] + print 'PU weights' + print 'Bin edges :',list(dataHist.xedges()) + print 'Weights : ',weights + +def make_plot( channel, x_axis_title, y_axis_title, + signal_region_tree, + control_region_tree, + branchName, + name_prefix, x_limits, nBins, + use_qcd_data_region = False, + compare_qcd_signal_with_data_control = False, + y_limits = [], + y_max_scale = 1.3, + rebin = 1, + legend_location = ( 0.98, 0.78 ), cms_logo_location = 'right', + log_y = False, + legend_color = False, + ratio_y_limits = [0.3, 2.5], + normalise = False, + ): + global output_folder, measurement_config, category, normalise_to_fit + global preliminary, norm_variable, sum_bins, b_tag_bin, histogram_files + + controlToCompare = [] + if 'electron' in channel : + controlToCompare = ['QCDConversions', 'QCD non iso e+jets'] + elif 'muon' in channel : + controlToCompare = ['QCD iso > 0.3', 'QCD 0.12 < iso <= 0.3'] + + histogramsToCompare = {} + for qcd_data_region in controlToCompare: + print 'Doing ',qcd_data_region + # Input files, normalisations, tree/region names + title = title_template % ( measurement_config.new_luminosity, measurement_config.centre_of_mass_energy ) + normalisation = None + weightBranchSignalRegion = 'EventWeight' + if 'electron' in channel: + histogram_files['data'] = measurement_config.data_file_electron_trees + histogram_files['QCD'] = measurement_config.electron_QCD_MC_category_templates_trees[category] + if normalise_to_fit: + normalisation = normalisations_electron[norm_variable] + # if use_qcd_data_region: + # qcd_data_region = 'QCDConversions' + # # qcd_data_region = 'QCD non iso e+jets' + if not 'QCD' in channel and not 'NPU' in branchName: + weightBranchSignalRegion += ' * ElectronEfficiencyCorrection' + if 'muon' in channel: + histogram_files['data'] = measurement_config.data_file_muon_trees + histogram_files['QCD'] = measurement_config.muon_QCD_MC_category_templates_trees[category] + if normalise_to_fit: + normalisation = normalisations_muon[norm_variable] + # if use_qcd_data_region: + # qcd_data_region = 'QCD iso > 0.3' + if not 'QCD' in channel and not 'NPU' in branchName: + weightBranchSignalRegion += ' * MuonEfficiencyCorrection' + + if not "_NPUNoWeight" in name_prefix: + weightBranchSignalRegion += ' * PUWeight' + + if not "_NBJetsNoWeight" in name_prefix: + weightBranchSignalRegion += ' * BJetWeight' + + selection = '1' + if branchName == 'abs(lepton_eta)' : + selection = 'lepton_eta > -10' + else: + selection = '%s >= 0' % branchName + # if 'QCDConversions' in signal_region_tree: + # selection += '&& isTightElectron' + # print selection + histograms = get_histograms_from_trees( trees = [signal_region_tree, control_region_tree], branch = branchName, weightBranch = weightBranchSignalRegion, files = histogram_files, nBins = nBins, xMin = x_limits[0], xMax = x_limits[-1], selection = selection ) + histograms_QCDControlRegion = None + if use_qcd_data_region: + qcd_control_region = signal_region_tree.replace( 'Ref selection', qcd_data_region ) + histograms_QCDControlRegion = get_histograms_from_trees( trees = [qcd_control_region], branch = branchName, weightBranch = 'EventWeight', files = histogram_files, nBins = nBins, xMin = x_limits[0], xMax = x_limits[-1], selection = selection ) + + # Split histograms up into signal/control (?) + signal_region_hists = {} + control_region_hists = {} + for sample in histograms.keys(): + signal_region_hists[sample] = histograms[sample][signal_region_tree] + + if compare_qcd_signal_with_data_control: + if sample is 'data': + signal_region_hists[sample] = histograms[sample][control_region_tree] + elif sample is 'QCD' : + signal_region_hists[sample] = histograms[sample][signal_region_tree] + else: + del signal_region_hists[sample] + + if use_qcd_data_region: + control_region_hists[sample] = histograms_QCDControlRegion[sample][qcd_control_region] + + # Prepare histograms + if normalise_to_fit: + # only scale signal region to fit (results are invalid for control region) + prepare_histograms( signal_region_hists, rebin = rebin, + scale_factor = measurement_config.luminosity_scale, + normalisation = normalisation ) + elif normalise_to_data: + totalMC = 0 + for sample in signal_region_hists: + if sample is 'data' : continue + totalMC += signal_region_hists[sample].Integral() + newScale = signal_region_hists['data'].Integral() / totalMC + + prepare_histograms( signal_region_hists, rebin = rebin, + scale_factor = newScale, + ) + else: + print measurement_config.luminosity_scale + prepare_histograms( signal_region_hists, rebin = rebin, + scale_factor = measurement_config.luminosity_scale ) + prepare_histograms( control_region_hists, rebin = rebin, + scale_factor = measurement_config.luminosity_scale ) + + # Use qcd from data control region or not + qcd_from_data = None + if use_qcd_data_region: + qcd_from_data = clean_control_region( control_region_hists, + + subtract = ['TTJet', 'V+Jets', 'SingleTop'] ) + # Normalise control region correctly + nBins = signal_region_hists['QCD'].GetNbinsX() + n, error = signal_region_hists['QCD'].integral(0,nBins+1,error=True) + n_qcd_predicted_mc_signal = ufloat( n, error) + + n, error = control_region_hists['QCD'].integral(0,nBins+1,error=True) + n_qcd_predicted_mc_control = ufloat( n, error) + + n, error = qcd_from_data.integral(0,nBins+1,error=True) + n_qcd_control_region = ufloat( n, error) + + if not n_qcd_control_region == 0: + dataDrivenQCDScale = n_qcd_predicted_mc_signal / n_qcd_predicted_mc_control + print 'Overall scale : ',dataDrivenQCDScale + qcd_from_data.Scale( dataDrivenQCDScale.nominal_value ) + signalToControlScale = n_qcd_predicted_mc_signal / n_qcd_control_region + dataToMCscale = n_qcd_control_region / n_qcd_predicted_mc_control + print "Signal to control :",signalToControlScale + print "QCD scale : ",dataToMCscale + else: + qcd_from_data = signal_region_hists['QCD'] + + # Which histograms to draw, and properties + histograms_to_draw = [] + histogram_lables = [] + histogram_colors = [] + + if compare_qcd_signal_with_data_control : + histograms_to_draw = [signal_region_hists['data'], qcd_from_data ] + histogram_lables = ['data', 'QCD'] + histogram_colors = ['black', 'yellow'] + else : + histograms_to_draw = [signal_region_hists['data'], qcd_from_data, + 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 = [colours['data'], colours['QCD'], colours['V+Jets'], colours['Single-Top'], colours['TTJet'] ] + + + print list(qcd_from_data.y()) + histogramsToCompare[qcd_data_region] = qcd_from_data + + print histogramsToCompare + histogram_properties = Histogram_properties() + histogram_properties.name = 'QCD_control_region_comparison_' + channel + '_' + branchName + histogram_properties.title = title + histogram_properties.x_axis_title = x_axis_title + histogram_properties.y_axis_title = y_axis_title + histogram_properties.x_limits = x_limits + histogram_properties.y_limits = y_limits + histogram_properties.mc_error = 0.0 + histogram_properties.legend_location = ( 0.98, 0.78 ) + histogram_properties.ratio_y_limits = ratio_y_limits + if 'electron' in channel: + make_control_region_comparison(histogramsToCompare['QCDConversions'], histogramsToCompare['QCD non iso e+jets'], + name_region_1='Conversions', name_region_2='Non Iso', + histogram_properties=histogram_properties, save_folder=output_folder) + elif 'muon' in channel: + make_control_region_comparison(histogramsToCompare['QCD iso > 0.3'], histogramsToCompare['QCD 0.12 < iso <= 0.3'], + name_region_1='QCD iso > 0.3', name_region_2='QCD 0.12 < iso <= 0.3', + histogram_properties=histogram_properties, save_folder=output_folder) + + # histogram_properties = Histogram_properties() + # histogram_properties.name = name_prefix + b_tag_bin + # if category != 'central': + # histogram_properties.name += '_' + category + # histogram_properties.title = title + # histogram_properties.x_axis_title = x_axis_title + # histogram_properties.y_axis_title = y_axis_title + # histogram_properties.x_limits = x_limits + # histogram_properties.y_limits = y_limits + # histogram_properties.y_max_scale = y_max_scale + # 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.legend_location = legend_location + # histogram_properties.cms_logo_location = cms_logo_location + # histogram_properties.preliminary = preliminary + # histogram_properties.set_log_y = log_y + # histogram_properties.legend_color = legend_color + # if ratio_y_limits: + # histogram_properties.ratio_y_limits = ratio_y_limits + + # # if normalise_to_fit: + # # histogram_properties.mc_error = get_normalisation_error( normalisation ) + # # histogram_properties.mc_errors_label = 'fit uncertainty' + # # else: + # # histogram_properties.mc_error = mc_uncertainty + # # histogram_properties.mc_errors_label = 'MC unc.' + + # if normalise_to_data: + # histogram_properties.name += '_normToData' + # output_folder_to_use = output_folder + # if use_qcd_data_region: + # output_folder_to_use += 'WithQCDFromControl/' + # make_folder_if_not_exists(output_folder_to_use) + + # if branchName == 'NPU': + # getPUWeights(histograms_to_draw, histogram_lables) + + # print output_folder_to_use + # # Actually draw histograms + # make_data_mc_comparison_plot( histograms_to_draw, histogram_lables, histogram_colors, + # histogram_properties, save_folder = output_folder_to_use, + # show_ratio = False, normalise = normalise, + # ) + # histogram_properties.name += '_with_ratio' + # loc = histogram_properties.legend_location + # # adjust legend location as it is relative to canvas! + # histogram_properties.legend_location = ( loc[0], loc[1] + 0.05 ) + # make_data_mc_comparison_plot( histograms_to_draw, histogram_lables, histogram_colors, + # histogram_properties, save_folder = output_folder_to_use, + # show_ratio = True, normalise = normalise, + # ) + + +if __name__ == '__main__': + set_root_defaults() + parser = OptionParser() + parser.add_option( "-p", "--path", dest = "path", default = 'data/M3_angle_bl/', + help = "set path to JSON files" ) + parser.add_option( "-o", "--output_folder", dest = "output_folder", default = 'plots/qcdComparison/', + help = "set path to save plots" ) + 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 = 13, type = int, + help = "set the centre of mass energy for analysis. Default = 13 [TeV]" ) + parser.add_option( "--category", dest = "category", default = 'central', + help = "set the category to take the fit results from (default: central)" ) + parser.add_option( "-n", "--normalise_to_fit", dest = "normalise_to_fit", action = "store_true", + help = "normalise the MC to fit results" ) + parser.add_option( "-d", "--normalise_to_data", dest = "normalise_to_data", action = "store_true", + help = "normalise the MC to data" ) + parser.add_option( "-a", "--additional-plots", action = "store_true", dest = "additional_QCD_plots", + help = "creates a set of QCD plots for exclusive bins for all variables" ) + + ( options, args ) = parser.parse_args() + measurement_config = XSectionConfig( options.CoM ) + # caching of variables for shorter access + translate_options = measurement_config.translate_options + + path_to_JSON = '%s/%dTeV/' % ( options.path, measurement_config.centre_of_mass_energy ) + normalise_to_fit = options.normalise_to_fit + normalise_to_data = options.normalise_to_data + if normalise_to_fit: + output_folder = '%s/after_fit/%dTeV/' % ( options.output_folder, measurement_config.centre_of_mass_energy ) + else: + output_folder = '%s/before_fit/%dTeV/' % ( options.output_folder, measurement_config.centre_of_mass_energy ) + make_folder_if_not_exists( output_folder ) + output_folder_base = output_folder + category = options.category + met_type = translate_options[options.metType] + make_additional_QCD_plots = options.additional_QCD_plots + + # 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 + + histogram_files = { + 'TTJet': measurement_config.ttbar_category_templates_trees[category], + 'V+Jets': measurement_config.VJets_category_templates_trees[category], + 'QCD': measurement_config.electron_QCD_MC_category_templates_trees[category], + 'SingleTop': measurement_config.SingleTop_category_templates_trees[category], + } + + # getting normalisations + normalisations_electron = { + # 'MET':get_fitted_normalisation( 'MET', 'electron', path_to_JSON, category, met_type ), + # 'HT':get_fitted_normalisation( 'HT', 'electron', path_to_JSON, category, met_type ), + # 'ST':get_fitted_normalisation( 'ST', 'electron', path_to_JSON, category, met_type ), + # 'MT':get_fitted_normalisation( 'MT', 'electron', path_to_JSON, category, met_type ), + # 'WPT':get_fitted_normalisation( 'WPT', 'electron', path_to_JSON, category, met_type ) + } + normalisations_muon = { + # 'MET':get_fitted_normalisation( 'MET', 'muon', path_to_JSON, category, met_type ), + # 'HT':get_fitted_normalisation( 'HT', 'muon', path_to_JSON, category, met_type ), + # 'ST':get_fitted_normalisation( 'ST', 'muon', path_to_JSON, category, met_type ), + # 'MT':get_fitted_normalisation( 'MT', 'muon', path_to_JSON, category, met_type ), + # 'WPT':get_fitted_normalisation( 'WPT', 'muon', path_to_JSON, category, met_type ) + } + preliminary = True + useQCDControl = True + b_tag_bin = '2orMoreBtags' + norm_variable = 'MET' + # comment out plots you don't want + include_plots = [ + 'HT', + # 'MET', + 'METNoHF', + 'ST', + 'WPT', + # 'NVertex', + # 'NVertexNoWeight', + 'LeptonPt', + 'AbsLeptonEta', + # # # 'Mjj', + # # # 'M3', + # # # 'angle_bl', + 'NJets', + # 'NBJets', + # 'NBJetsNoWeight', + # 'JetPt', + # 'AbsLeptonEta', + # 'RelIso', + # 'sigmaietaieta' + ] + + additional_qcd_plots = [ + # 'QCDHT', + # 'QCDMET', + # 'QCDST', + # 'QCDWPT', + # 'QCDAbsLeptonEta', + # 'QCDLeptonPt', + # 'QCDNJets', + # 'QCDsigmaietaieta', + # 'QCDRelIso', + # 'QCDHT_dataControl_mcSignal', + ] + + if make_additional_QCD_plots: + include_plots.extend( additional_qcd_plots ) + + + for channel, label in { + 'electron' : 'EPlusJets', + 'muon' : 'MuPlusJets' + }.iteritems() : + b_tag_bin = '2orMoreBtags' + + # Set folder for this batch of plots + output_folder = output_folder_base + "/Variables/" + make_folder_if_not_exists(output_folder) + + ################################################### + # HT + ################################################### + norm_variable = 'HT' + if 'HT' in include_plots: + print '---> HT' + make_plot( channel, + x_axis_title = '$%s$ [GeV]' % variables_latex['HT'], + y_axis_title = 'Events/(%i GeV)' % binWidth(control_plots_bins['HT']), + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + branchName = 'HT', + name_prefix = '%s_HT_' % label, + x_limits = control_plots_bins['HT'], + nBins = 20, + rebin = 1, + legend_location = ( 0.9, 0.83 ), + cms_logo_location = 'left', + use_qcd_data_region = useQCDControl, + ) + +################################################### + # MET + ################################################### + norm_variable = 'MET' + if 'METNoHF' in include_plots: + print '---> METNoHF' + make_plot( channel, + x_axis_title = '$%s$ [GeV]' % variables_latex['MET'], + y_axis_title = 'Events/(%i GeV)' % binWidth(control_plots_bins['MET']), + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + branchName = 'METNoHF', + name_prefix = '%s_METNoHF_' % label, + x_limits = control_plots_bins['MET'], + nBins = len(control_plots_bins['MET'])-1, + rebin = 1, + legend_location = ( 0.9, 0.83 ), + cms_logo_location = 'left', + use_qcd_data_region = useQCDControl, + ) + + ################################################### + # ST + ################################################### + norm_variable = 'ST' + if 'ST' in include_plots: + print '---> ST' + make_plot( channel, + x_axis_title = '$%s$ [GeV]' % variables_latex['ST'], + y_axis_title = 'Events/(%i GeV)' % binWidth(control_plots_bins['ST']), + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + branchName = 'STNoHF', + name_prefix = '%s_STNoHF_' % label, + x_limits = control_plots_bins['ST'], + nBins = 20, + rebin = 1, + legend_location = ( 0.9, 0.83 ), + cms_logo_location = 'left', + use_qcd_data_region = useQCDControl, + ) + + ################################################### + # WPT + ################################################### + norm_variable = 'WPT' + if 'WPT' in include_plots: + print '---> WPT' + make_plot( channel, + x_axis_title = '$%s$ [GeV]' % variables_latex['WPT'], + y_axis_title = 'Events/(%i GeV)' % binWidth(control_plots_bins['WPT']), + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + branchName = 'WPTNoHF', + name_prefix = '%s_WPTNoHF_' % label, + x_limits = control_plots_bins['WPT'], + nBins = 16, + rebin = 1, + legend_location = ( 0.9, 0.83 ), + cms_logo_location = 'left', + use_qcd_data_region = useQCDControl, + ) + + # Set folder for this batch of plots + output_folder = output_folder_base + "/FitVariables/" + make_folder_if_not_exists(output_folder) + + ################################################### + # NJets + ################################################### + if 'NJets' in include_plots: + print '---> NJets' + make_plot( channel, + x_axis_title = '$%s$' % control_plots_latex['NJets'], + y_axis_title = 'Events / 1', + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + branchName = 'NJets', + name_prefix = '%s_NJets_' % label, + x_limits = control_plots_bins['NJets'], + nBins = len(control_plots_bins['NJets'])-1, + rebin = 1, + legend_location = ( 0.9, 0.83 ), + cms_logo_location = 'left', + use_qcd_data_region = useQCDControl, + ) + + ################################################### + # Lepton Pt + ################################################### + if 'LeptonPt' in include_plots: + print '---> Lepton Pt' + binsLabel = 'ElectronPt' + if channel == 'muon': + binsLabel = 'MuonPt' + make_plot( channel, + x_axis_title = '$%s$' % control_plots_latex['pt'], + y_axis_title = 'Events/(%i GeV)' % binWidth(control_plots_bins[binsLabel]), + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % ( label ), + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % (label ), + branchName = 'lepton_pt', + name_prefix = '%s_LeptonPt_' % label, + x_limits = control_plots_bins[binsLabel], + nBins = len(control_plots_bins[binsLabel])-1, + rebin = 1, + legend_location = ( 0.9, 0.83 ), + cms_logo_location = 'left', + use_qcd_data_region = useQCDControl, + ) + ################################################### + # Lepton Eta + ################################################### + if 'LeptonEta' in include_plots: + print '---> Lepton Eta' + treeName = 'Electron/Electrons' + if channel == 'muon': + treeName = 'Muon/Muons' + + make_plot( channel, + x_axis_title = '$%s$' % control_plots_latex['eta'], + y_axis_title = 'Events/(%.1f)' % binWidth(control_plots_bins['LeptonEta']), + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/%s' % ( label, treeName), + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/%s' % (label, treeName), + branchName = 'lepton_eta', + name_prefix = '%s_LeptonEta_' % label, + x_limits = control_plots_bins['LeptonEta'], + nBins = len(control_plots_bins['LeptonEta'])-1, + rebin = 1, + legend_location = ( 0.95, 0.78 ), + cms_logo_location = 'right', + use_qcd_data_region = False, + ) + + ################################################### + # AbsLepton Eta + ################################################### + if 'AbsLeptonEta' in include_plots: + print '---> Abs Lepton Eta' + + make_plot( channel, + x_axis_title = '$%s$' % control_plots_latex['eta'], + y_axis_title = 'Events/(%.1f)' % binWidth(control_plots_bins['AbsLeptonEta']), + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % ( label ), + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % (label ), + branchName = 'abs(lepton_eta)', + name_prefix = '%s_AbsLeptonEta_' % label, + x_limits = control_plots_bins['AbsLeptonEta'], + nBins = len(control_plots_bins['AbsLeptonEta'])-1, + rebin = 1, + legend_location = ( 0.9, 0.88 ), + cms_logo_location = 'left', + use_qcd_data_region = 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 8a3ebf5e..bb4accf2 100644 --- a/src/cross_section_measurement/create_measurement.py +++ b/src/cross_section_measurement/create_measurement.py @@ -9,7 +9,7 @@ Example: python src/cross_section_measurement/create_measurement.py -c ''' - +from __future__ import print_function from optparse import OptionParser import tools.measurement from config import XSectionConfig, variable_binning @@ -41,7 +41,7 @@ def main(): categories.extend(measurement_config.categories_and_prefixes.keys()) categories.extend(measurement_config.rate_changing_systematics_names) categories.extend([measurement_config.vjets_theory_systematic_prefix + - systematic for systematic in measurement_config.generator_systematics if not ('mass' in systematic or 'hadronisation' in systematic)]) + systematic for systematic in measurement_config.generator_systematics if not ('mass' in systematic or 'hadronisation' in systematic or 'NLO' in systematic)]) for variable in variable_binning.bin_edges.keys(): for category in categories: @@ -406,7 +406,6 @@ def create_input(config, sample, variable, category, channel, template, if sample != 'data': if category in config.met_systematics_suffixes and not variable in config.variables_no_met: - # print variable, category branch = template.split('/')[-1] branch += '_METUncertainties[%s]' % config.met_systematics[ category] @@ -414,7 +413,19 @@ def create_input(config, sample, variable, category, channel, template, if 'JES_down' in category or 'JES_up' in category or 'JER_down' in category or 'JER_up' in category: tree += config.categories_and_prefixes[category] + if not sample == 'data': + if 'JES_down' in category: + input_file = input_file.replace('tree','minusJES_tree') + elif 'JES_up' in category: + input_file = input_file.replace('tree','plusJES_tree') + elif 'JER_up' in category: + input_file = input_file.replace('tree','plusJER_tree') + elif 'JER_down' in category: + input_file = input_file.replace('tree','minusJER_tree') + selection = '{0} >= 0'.format(branch) + if variable == 'abs_lepton_eta': + selection += ' && {0} <= 3'.format(branch) else: hist = template @@ -441,7 +452,8 @@ def create_input(config, sample, variable, category, channel, template, weight_branches.append('1') else: weight_branches.append('EventWeight') - weight_branches.append('PUWeight') + if category != 'PileUpSystematic': + weight_branches.append('PUWeight') if category == 'BJet_down': weight_branches.append('BJetDownWeight') elif category == 'BJet_up': diff --git a/src/cross_section_measurement/make_control_plots_fromTrees.py b/src/cross_section_measurement/make_control_plots_fromTrees.py index 2abfdf75..b0d5459f 100644 --- a/src/cross_section_measurement/make_control_plots_fromTrees.py +++ b/src/cross_section_measurement/make_control_plots_fromTrees.py @@ -39,9 +39,9 @@ def getPUWeights(histograms_to_draw, histogram_lables) : mcValues = list(mcHist.y()) weights = [ data / mc for data, mc in zip(dataValues, mcValues)] - print 'PU weights' - print 'Bin edges :',list(dataHist.xedges()) - print 'Weights : ',weights +# print 'PU weights' +# print 'Bin edges :',list(dataHist.xedges()) +# print 'Weights : ',weights def make_plot( channel, x_axis_title, y_axis_title, signal_region_tree, @@ -125,8 +125,9 @@ def make_plot( channel, x_axis_title, y_axis_title, selection = '1' if branchName == 'abs(lepton_eta)' : selection = 'lepton_eta > -10' + else: + selection = '%s >= 0' % branchName histograms = get_histograms_from_trees( trees = [signal_region_tree, control_region_tree], branch = branchName, weightBranch = weightBranchSignalRegion, files = histogram_files, nBins = nBins, xMin = x_limits[0], xMax = x_limits[-1], selection = selection ) - histograms_QCDControlRegion = None if use_qcd_data_region: qcd_control_region = signal_region_tree.replace( 'Ref selection', qcd_data_region ) @@ -216,6 +217,14 @@ def make_plot( channel, x_axis_title, y_axis_title, histogram_lables = ['data', 'QCD', 'V+Jets', 'Single-Top', samples_latex['TTJet']] histogram_colors = [colours['data'], colours['QCD'], colours['V+Jets'], colours['Single-Top'], colours['TTJet'] ] + + print 'Normalisation after selection' + print 'Data :',signal_region_hists['data'].integral(overflow=True) + print 'TTJet :',signal_region_hists['TTJet'].integral(overflow=True) + print 'Single Top :',signal_region_hists['SingleTop'].integral(overflow=True) + print 'V+Jets :',signal_region_hists['V+Jets'].integral(overflow=True) + print 'QCD :',qcd_from_data.integral(overflow=True) + histogram_properties = Histogram_properties() histogram_properties.name = name_prefix + b_tag_bin if category != 'central': @@ -344,34 +353,36 @@ def make_plot( channel, x_axis_title, y_axis_title, norm_variable = 'MET' # comment out plots you don't want include_plots = [ - 'HT', - 'MET', - 'METNoHF', - 'ST', - 'WPT', - 'NVertex', - 'NVertexNoWeight', - 'LeptonPt', - 'AbsLeptonEta', - # # 'Mjj', - # # 'M3', - # # 'angle_bl', + # 'HT', + # 'MET', + # 'METNoHF', + # 'ST', + # 'WPT', + # 'NVertex', + # 'NVertexNoWeight', + # 'LeptonPt', + # 'AbsLeptonEta', + # # # 'Mjj', + # # # 'M3', + # # # 'angle_bl', 'NJets', - 'NBJets', - 'NBJetsNoWeight' + # 'NBJets', + # 'NBJetsNoWeight', # 'JetPt', # 'AbsLeptonEta', # 'RelIso', + # 'sigmaietaieta' ] additional_qcd_plots = [ # 'QCDHT', - 'QCDMET', - 'QCDST', - 'QCDWPT', - 'QCDAbsLeptonEta', - 'QCDLeptonPt', - 'QCDNJets', + # 'QCDMET', + # 'QCDST', + # 'QCDWPT', + # 'QCDAbsLeptonEta', + # 'QCDLeptonPt', + # 'QCDNJets', + # 'QCDsigmaietaieta', # 'QCDRelIso', # 'QCDHT_dataControl_mcSignal', ] @@ -753,7 +764,7 @@ def make_plot( channel, x_axis_title, y_axis_title, rebin = 1, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', - use_qcd_data_region = useQCDControl, + use_qcd_data_region = False, ) ################################################### @@ -774,7 +785,7 @@ def make_plot( channel, x_axis_title, y_axis_title, rebin = 1, legend_location = ( 0.9, 0.88 ), cms_logo_location = 'left', - use_qcd_data_region = useQCDControl, + use_qcd_data_region = True, ) ################################################### @@ -801,13 +812,34 @@ def make_plot( channel, x_axis_title, y_axis_title, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', ) + ################################################### + # Sigma ieta ieta + ################################################### + + norm_variable = 'sigmaietaieta' + if 'sigmaietaieta' in include_plots and channel == 'electron': + print '---> sigmaietaieta' + make_plot( channel, + x_axis_title = '$%s$' % variables_latex['sigmaietaieta'], + y_axis_title = 'Events/(%i GeV)' % binWidth(control_plots_bins['sigmaietaieta']), + signal_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + control_region_tree = 'TTbar_plus_X_analysis/%s/Ref selection/FitVariables' % label, + branchName = 'sigmaIetaIeta', + name_prefix = '%s_sigmaIetaIeta_' % label, + x_limits = control_plots_bins['sigmaietaieta'], + nBins = len(control_plots_bins['sigmaietaieta'])-1, + rebin = 1, + legend_location = ( 0.9, 0.83 ), + cms_logo_location = 'left', + use_qcd_data_region = useQCDControl, + ) ################################################### # QCD Control Region ################################################### for channel, label in { - 'electronQCDNonIso' : 'EPlusJets/QCD non iso e+jets', - 'electronQCDConversions' : 'EPlusJets/QCDConversions', + # 'electronQCDNonIso' : 'EPlusJets/QCD non iso e+jets', + # 'electronQCDConversions' : 'EPlusJets/QCDConversions', 'muonQCDNonIso' : 'MuPlusJets/QCD iso > 0.3', 'muonQCDNonIso2' : 'MuPlusJets/QCD 0.12 < iso <= 0.3', }.iteritems() : @@ -815,6 +847,7 @@ def make_plot( channel, x_axis_title, y_axis_title, # Set folder for this batch of plots output_folder = output_folder_base + "QCDControl/Variables/%s/" % channel + # output_folder = output_folder_base + "QCDControl/Variables/%s/TightElectron/" % channel make_folder_if_not_exists(output_folder) print 'Control region :',label @@ -928,7 +961,8 @@ def make_plot( channel, x_axis_title, y_axis_title, ) # Set folder for this batch of plots - output_folder = output_folder_base + "QCDControl/Control/%s" % channel + output_folder = output_folder_base + "QCDControl/Control/%s/" % channel + # output_folder = output_folder_base + "QCDControl/Control/%s/TightElectron/" % channel make_folder_if_not_exists(output_folder) ################################################### # Abs Lepton Eta @@ -1019,4 +1053,26 @@ def make_plot( channel, x_axis_title, y_axis_title, rebin = 1, legend_location = ( 0.95, 0.78 ), cms_logo_location = 'right', + ) + + ################################################### + # Sigma ieta ieta + ################################################### + + norm_variable = 'sigmaietaieta' + if 'QCDsigmaietaieta' in include_plots and not 'MuPlusJets' in treeName: + print '---> sigmaietaieta' + make_plot( channel, + x_axis_title = '$%s$' % variables_latex['sigmaietaieta'], + y_axis_title = 'Events/(%i GeV)' % binWidth(control_plots_bins['sigmaietaieta']), + signal_region_tree = 'TTbar_plus_X_analysis/%s/FitVariables' % ( treeName ), + control_region_tree = 'TTbar_plus_X_analysis/%s/FitVariables' % ( treeName ), + branchName = 'sigmaIetaIeta', + name_prefix = '%s_sigmaIetaIeta_' % channel, + x_limits = control_plots_bins['sigmaietaieta'], + y_max_scale = 1.5, + nBins = len(control_plots_bins['sigmaietaieta'])-1, + rebin = 1, + legend_location = ( 0.95, 0.85 ), + cms_logo_location = 'left', ) \ No newline at end of file diff --git a/src/unfolding_tests/create_toy_mc.py b/src/unfolding_tests/create_toy_mc.py index f266da71..75a18d56 100644 --- a/src/unfolding_tests/create_toy_mc.py +++ b/src/unfolding_tests/create_toy_mc.py @@ -46,7 +46,7 @@ def main(): # variable = options.variable met_type = measurement_config.translate_options[options.metType] - create_toy_mc(input_file=measurement_config.unfolding_madgraph, + create_toy_mc(input_file=measurement_config.unfolding_central, output_folder=options.output_folder, # variable=variable, n_toy=options.n_toy_mc, diff --git a/src/unfolding_tests/create_unfolding_pull_data.py b/src/unfolding_tests/create_unfolding_pull_data.py index 80701a40..33d64f25 100644 --- a/src/unfolding_tests/create_unfolding_pull_data.py +++ b/src/unfolding_tests/create_unfolding_pull_data.py @@ -42,7 +42,7 @@ def main(): help="file with toy MC") parser.add_option("-v", "--variable", dest="variable", default='MET', help="set the variable to analyse (MET, HT, ST, MT, WPT)") - parser.add_option("-s", "--centre-of-mass-energy", dest="CoM", default=8, + parser.add_option("-s", "--centre-of-mass-energy", dest="CoM", default=13, help='''set the centre of mass energy for analysis. Default = 8 [TeV]''', type=int) parser.add_option("-c", "--channel", type='string',