From c8e3a48ea2786cb586857567ea051cfe0a5dd8e4 Mon Sep 17 00:00:00 2001 From: Salvatore Rappoccio Date: Tue, 19 Nov 2013 10:14:01 -0500 Subject: [PATCH 1/2] Adding first dummy script to unfold top xs --- IHeartNY/test/unfoldTopPt.py | 142 +++++++++++++++++++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 IHeartNY/test/unfoldTopPt.py diff --git a/IHeartNY/test/unfoldTopPt.py b/IHeartNY/test/unfoldTopPt.py new file mode 100644 index 00000000..4d3b6320 --- /dev/null +++ b/IHeartNY/test/unfoldTopPt.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python +# ============================================================================== +# File and Version Information: +# $Id: RooUnfoldExample.py 302 2011-09-30 20:39:20Z T.J.Adye $ +# +# Description: +# Simple example usage of the RooUnfold package using toy MC. +# +# Author: Tim Adye +# +# ============================================================================== + +from ROOT import gRandom, TH1, TH1D, cout, TFile, gSystem, TCanvas, TPad, gROOT, gStyle, THStack, TLegend, TLatex + +gROOT.Macro("rootlogon.C") + + +gStyle.SetOptStat(000000) + +gStyle.SetTitleFont(43) +#gStyle.SetTitleFontSize(0.05) +gStyle.SetTitleFont(43, "XYZ") +gStyle.SetTitleSize(30, "XYZ") +gStyle.SetTitleOffset(3.5, "X") +gStyle.SetLabelFont(43, "XYZ") +gStyle.SetLabelSize(24, "XYZ") + + +gSystem.Load("RooUnfold-1.1.1/libRooUnfold.so") + +from ROOT import RooUnfoldResponse +from ROOT import RooUnfold +from ROOT import RooUnfoldBayes +from ROOT import RooUnfoldSvd +# from ROOT import RooUnfoldTUnfold + +# Constants +lum = 19.7 # fb-1 +sigma_ttbar_NNLO = 245.8 * 1000. # fb, from http://arxiv.org/pdf/1303.6254.pdf +Nmc = 21675970 +SF_b = 0.97 +SF_t = 0.94 + +# ============================================================================== +# Example Unfolding +# ============================================================================== + +fdata = TFile("TTSemilepAnalyzer_unfolding_data_type1.root") +fmc = TFile("TTSemilepAnalyzer_antibtag_w_mucut_type1.root") + +response= fmc.Get("response_pt") + +print "==================================== Get Hists =====================================" +hTrue= fmc.Get("ptGenTop") +hMeas= fdata.Get("ptRecoTop") + +hTrue.Scale( sigma_ttbar_NNLO * lum / float(Nmc) * SF_b * SF_t ) + +print "==================================== UNFOLD ===================================" +unfold= RooUnfoldBayes (response, hMeas, 10); # OR +#unfold= RooUnfoldSvd (response, hMeas, 20); # OR +# unfold= RooUnfoldTUnfold (response, hMeas); + + + +c1 = TCanvas("c", "c", 700, 700) +pad1 = TPad("pad1","pad1",0,0.3,1,1) +pad1.SetBottomMargin(0.05); +pad1.SetLeftMargin(0.2) +pad1.Draw(); +pad1.cd(); + + +hReco= unfold.Hreco(); + + + +hFrac = hReco.Clone() +hFrac.SetName("hFrac") +hFrac.SetTitle(";p_{T} (GeV);Measured/Truth") +hFrac.Divide(hTrue) +hMeas.Sumw2() + + +#unfold.PrintTable (cout, hTrue); +hReco.SetMarkerStyle(21) +hMeas.SetMarkerStyle( 25); + +hTrue.SetTitle(";;d #sigma / d p_{T} (fb / GeV)") +hTrue.Draw('hist') +hReco.Draw('same') +hMeas.Draw('same') +hTrue.UseCurrentStyle() +hTrue.SetLineColor(4); +hTrue.GetYaxis().SetTitleSize(25) +hTrue.GetXaxis().SetLabelSize(0) +hTrue.GetXaxis().SetRangeUser(400., 10000.) + + +# Correct for bin width +for ibin in range(1, hTrue.GetXaxis().GetNbins() ) : + width = hTrue.GetBinWidth(ibin) + hTrue.SetBinContent(ibin, hTrue.GetBinContent(ibin) / width ) + hMeas.SetBinContent(ibin, hMeas.GetBinContent(ibin) / width ) + hReco.SetBinContent(ibin, hReco.GetBinContent(ibin) / width ) + hTrue.SetBinError(ibin, hTrue.GetBinError(ibin) / width ) + hMeas.SetBinError(ibin, hMeas.GetBinError(ibin) / width ) + hReco.SetBinError(ibin, hReco.GetBinError(ibin) / width ) + +leg = TLegend(0.6, 0.6, 0.84, 0.84) +leg.SetFillColor(0) +leg.SetBorderSize(0) +leg.AddEntry( hReco, 'Unfolded data', 'p') +leg.AddEntry( hTrue, 'Generated', 'l') +leg.AddEntry( hMeas, 'Raw data', 'p') +leg.Draw() + +text1 = TLatex() +text1.SetNDC() +text1.SetTextFont(42) +text1.DrawLatex(0.20,0.87, "#scale[1.0]{CMS Preliminary, L = 19.7 fb^{-1} at #sqrt{s} = 8 TeV}") + + +c1.cd(); +pad2 = TPad("pad2","pad2",0,0.0,1,0.28) +pad2.SetTopMargin(0.05); +pad2.SetBottomMargin(0.4); +pad2.SetLeftMargin(0.2) +pad2.Draw(); +pad2.cd(); +hFrac.SetMaximum(2.0) +hFrac.SetMinimum(0.0) +hFrac.UseCurrentStyle() +hFrac.GetYaxis().SetTitleSize(25) +hFrac.GetYaxis().SetTitleOffset(1.7) +hFrac.GetXaxis().SetLabelSize(20) +hFrac.GetYaxis().SetNdivisions(2,4,0,False) + +hFrac.Draw("e") +hFrac.GetXaxis().SetRangeUser(400., 10000.) + +c1.Update() From b9568f282389fba429ecbcb420d6c221d9ebe7cb Mon Sep 17 00:00:00 2001 From: Salvatore Rappoccio Date: Tue, 19 Nov 2013 10:28:27 -0500 Subject: [PATCH 2/2] Updating unfolding script to run closure test --- IHeartNY/test/unfoldTopPt.py | 86 +++++++++++++++++++++++++++++------- 1 file changed, 71 insertions(+), 15 deletions(-) diff --git a/IHeartNY/test/unfoldTopPt.py b/IHeartNY/test/unfoldTopPt.py index 4d3b6320..206b0176 100644 --- a/IHeartNY/test/unfoldTopPt.py +++ b/IHeartNY/test/unfoldTopPt.py @@ -10,6 +10,29 @@ # # ============================================================================== +from optparse import OptionParser + + +parser = OptionParser() + + + +parser.add_option('--closureTest', metavar='F', action='store_true', + default=False, + dest='closureTest', + help='Run Closure Test') + + +parser.add_option('--plotFullRange', metavar='F', action='store_true', + default=False, + dest='plotFullRange', + help='Plot the full pt range (pt > 300 GeV)') + + +(options, args) = parser.parse_args() + +argv = [] + from ROOT import gRandom, TH1, TH1D, cout, TFile, gSystem, TCanvas, TPad, gROOT, gStyle, THStack, TLegend, TLatex gROOT.Macro("rootlogon.C") @@ -35,11 +58,18 @@ # from ROOT import RooUnfoldTUnfold # Constants -lum = 19.7 # fb-1 -sigma_ttbar_NNLO = 245.8 * 1000. # fb, from http://arxiv.org/pdf/1303.6254.pdf -Nmc = 21675970 -SF_b = 0.97 -SF_t = 0.94 +if options.closureTest == False : + lum = 19.7 # fb-1 + sigma_ttbar_NNLO = 245.8 * 1000. # fb, from http://arxiv.org/pdf/1303.6254.pdf + Nmc = 21675970 + SF_b = 0.97 + SF_t = 0.94 +else : + lum = 1.0 + sigma_ttbar_NNLO = 1.0 + Nmc = 1.0 + SF_b = 1.0 + SF_t = 1.0 # ============================================================================== # Example Unfolding @@ -52,7 +82,10 @@ print "==================================== Get Hists =====================================" hTrue= fmc.Get("ptGenTop") -hMeas= fdata.Get("ptRecoTop") +if options.closureTest == False : + hMeas= fdata.Get("ptRecoTop") +else : + hMeas= fmc.Get("ptRecoTop") hTrue.Scale( sigma_ttbar_NNLO * lum / float(Nmc) * SF_b * SF_t ) @@ -86,15 +119,19 @@ hReco.SetMarkerStyle(21) hMeas.SetMarkerStyle( 25); -hTrue.SetTitle(";;d #sigma / d p_{T} (fb / GeV)") -hTrue.Draw('hist') -hReco.Draw('same') +hReco.SetTitle(";;d #sigma / d p_{T} (fb / GeV)") +hReco.SetMinimum(0.0) +hReco.Draw() +hTrue.Draw('hist same') hMeas.Draw('same') hTrue.UseCurrentStyle() hTrue.SetLineColor(4); hTrue.GetYaxis().SetTitleSize(25) hTrue.GetXaxis().SetLabelSize(0) -hTrue.GetXaxis().SetRangeUser(400., 10000.) + +if options.plotFullRange == False : + hReco.GetXaxis().SetRangeUser(400., 10000.) + # Correct for bin width @@ -107,12 +144,20 @@ hMeas.SetBinError(ibin, hMeas.GetBinError(ibin) / width ) hReco.SetBinError(ibin, hReco.GetBinError(ibin) / width ) -leg = TLegend(0.6, 0.6, 0.84, 0.84) +leg = TLegend(0.4, 0.6, 0.84, 0.84) leg.SetFillColor(0) leg.SetBorderSize(0) -leg.AddEntry( hReco, 'Unfolded data', 'p') -leg.AddEntry( hTrue, 'Generated', 'l') -leg.AddEntry( hMeas, 'Raw data', 'p') + +if options.closureTest == False : + leg.AddEntry( hReco, 'Unfolded data', 'p') + leg.AddEntry( hTrue, 'Generated', 'l') + leg.AddEntry( hMeas, 'Raw data', 'p') +else : + leg.AddEntry( hReco, 'Unfolded MC : Closure test', 'p') + leg.AddEntry( hTrue, 'Generated', 'l') + leg.AddEntry( hMeas, 'Raw MC', 'p') + + leg.Draw() text1 = TLatex() @@ -137,6 +182,17 @@ hFrac.GetYaxis().SetNdivisions(2,4,0,False) hFrac.Draw("e") -hFrac.GetXaxis().SetRangeUser(400., 10000.) +if options.plotFullRange == False : + hFrac.GetXaxis().SetRangeUser(400., 10000.) + c1.Update() + +append = "" +if options.closureTest : + append += "_closure" +if options.plotFullRange : + append += "_fullrange" + +c1.Print("unfolded_ttbar_xs" + append + ".pdf", "pdf") +c1.Print("unfolded_ttbar_xs" + append + ".png", "png")