kiesel/plotTree/qcdClosure.py

#! /usr/bin/env python2
# -*- coding: utf-8 -*-

import ROOT
import argparse
import ConfigParser
import os.path

# to use user defined help message, sys.arv has to be sent to python and not
# to TApplication
ROOT.PyConfig.IgnoreCommandLineOptions = True

# use tdr style defined in own macro
import Styles
Styles.tdrStyle()

# this is the configuration file for plots including label, etc.
axisConf = ConfigParser.SafeConfigParser()
axisConf.read("axis.cfg")

# global variables:
integratedLuminosity = 19.3 #fb

def readTree( filename, treename = "susyTree" ):
        """
        filename: name of file containing the tree
        treename: name of the tree
        returns: TTree Object
        """
        if not os.path.isfile(filename):
                print( "File %s does not exist"%filename)
        tree = ROOT.TChain( treename )
        tree.AddFile( filename )
        return tree

def createHistoFromTree2D(tree, variable, weight="", nBin=50):
        """
        tree: tree to create histo from
        variable: variable to plot (must be a branch of the tree)
        weight: weights to apply (e.g. "var1*(var2 > 15)" will use weights from var1 and cut on var2 > 15
        nBins, firstBin, lastBin: number of bins, first bin and last bin (same as in TH1F constructor)
        nEvents: number of events to process (-1 = all)
        returns: histogram
        """
        from random import randint
        from sys import maxint
        #make a random name you could give something meaningfull here,
        #but that would make this less readable
        name = "%x"%(randint(0, maxint))

        # automatic binning if lastbin < firstbin
        result = ROOT.TH2F(name, variable, nBin, 0, -1, nBin, 0, -1)
        result.Sumw2()
        tree.Draw("%s>>%s"%(variable, name), weight, "goff,colz")
        return result

def createHistoFromTree(tree, variable, weight="", nBins=100, firstBin=None, lastBin=None, nEvents=-1):
        """
        tree: tree to create histo from
        variable: variable to plot (must be a branch of the tree)
        weight: weights to apply (e.g. "var1*(var2 > 15)" will use weights from var1 and cut on var2 > 15
        nBins, firstBin, lastBin: number of bins, first bin and last bin (same as in TH1F constructor)
        nEvents: number of events to process (-1 = all)
        returns: histogram
        """
        if ":" in variable:
                return createHistoFromTree2D(tree, variable, weight="")
        from ROOT import TH1F
        from random import randint
        from sys import maxint
        if nEvents < 0:
                nEvents = maxint
        if firstBin == None:
                firstBin = tree.GetMinimum( variable )
        if lastBin == None:
                lastBin = tree.GetMaximum( variable )
        #make a random name you could give something meaningfull here,
        #but that would make this less readable
        name = "%x"%(randint(0, maxint))
        result = TH1F(name, variable, nBins, firstBin, lastBin)
        result.Sumw2()
        tree.Draw("%s>>%s"%(variable, name), weight, "goff", nEvents)
        return result

def plot2D( tree, plot, cut, save=False ):
        # read axis config file
        var = plot.split(":")
        try:
                xlabel = axisConf.get( var[0], "label" )
        except:
                xlabel = ""
        try:
                ylabel = axisConf.get( var[1], "label" )
        except:
                ylabel = ""
        try:
                xunit = axisConf.get( var[0], "unit" )
        except:
                xunit = ""
        try:
                yunit = axisConf.get( var[1], "unit" )
        except:
                yunit = ""
        if xunit != "":
                xunit = " ["+xunit+"]"
        if yunit != "":
                yunit = " ["+yunit+"]"

        # modify histo
        histo = createHistoFromTree2D( tree, plot, cut )
        histo.SetTitle(";%s%s;%s%s"%(ylabel,yunit,xlabel,xunit))

        # draw histo
        canvas = ROOT.TCanvas()
        canvas.cd()
        histo.Draw("colz")

        # draw information
        cutText = ROOT.TPaveText(.5,.8,.9,.9,"ndc")
        cutText.SetBorderSize(0)
        cutText.SetFillColor(0)
        for line in cut.split("&&"):
                cutText.AddText( line )
        cutText.Draw()

        topInfo = ROOT.TPaveText(.1,.95,1,1,"ndc")
        topInfo.SetBorderSize(0)
        topInfo.SetFillColor(0)
        topInfo.AddText("Work in progress, %sfb^{-1}, #sqrt{s}=8TeV"%"?")
        topInfo.Draw()

        if save:
                canvas.SaveAs("%.pdf"%variable)
        else:
                raw_input()

def fillWeights( tree, jetCut, photonCut ):

        # this is the binning with witch the fo reweighting will be done and has to be studied
        nBins = 50
        xMin = 80
        xMax = 300
        h_photonPt = createHistoFromTree( tree, "photon[0].pt", photonCut, nBins, xMin, xMax )
        h_jetPt = createHistoFromTree( tree, "photon[0].pt", jetCut, nBins, xMin, xMax )

        # h_jetPt is now histogram with w^{-1} = h_jetPt / h_photonPt
        h_jetPt.Divide( h_photonPt )
        import array
        weight = array.array( "f", [0] )
        tree.Branch( "weight2", weight, "weight2/F" )

        for event in tree:
                pt = tree.photon.at(0).pt
                #print "p_T = %i GeV"%pt
                weight = h_jetPt.GetBinContent( h_jetPt.FindBin( pt ) )
                print weight
                #print tree.weight2
                #tree.Fill()


def plot( tree, plot, cut, save=False ):
        # read axis config file
        try:
                label = axisConf.get( plot, "label" )
                unit = axisConf.get( plot, "unit" )
                if unit != "":
                        unit = " ["+unit+"]"
        except:
                print "please specify label and unit in axis configuration file"
                label = plot
                unit = ""

        # modify histo
        histo = createHistoFromTree( tree, plot, cut )
        histo.SetTitle(";%s%s;Entries"%(label,unit))
        histo.SetLineColor(1)

        # draw histo
        canvas = ROOT.TCanvas()
        canvas.cd()
        canvas.SetLogy()
        histo.Draw("hist")

        # draw information
        cutText = ROOT.TPaveText(.5,.8,.9,.9,"ndc")
        cutText.SetBorderSize(0)
        cutText.SetFillColor(0)
        for line in cut.split("&&"):
                cutText.AddText( line )
        cutText.Draw()

        topInfo = ROOT.TPaveText(.1,.95,1,1,"ndc")
        topInfo.SetBorderSize(0)
        topInfo.SetFillColor(0)
        topInfo.AddText("Work in progress, %sfb^{-1}, #sqrt{s}=8TeV"%"?")
        topInfo.Draw()

        if save:
                canvas.SaveAs("%s.pdf"%plot)
        else:
                raw_input()


if __name__ == "__main__":
        # include knowledge about objects saved in the tree
        ROOT.gSystem.Load("../treeWriter/libTreeObjects.so")

        arguments = argparse.ArgumentParser( description="Calculate weighting "
                        +"factors for QCD background estimation." )
        arguments.add_argument("-f", "--file", default="../treeWriter/myTree.root",
                        help="ROOT file containing a TTree produced by 'treeWriter'.")
        arguments.add_argument("-c", "--cut", default="", help="TCut string.")
        arguments.add_argument("-d", "--distribution", default =['met'], nargs="+",
                        help="Distributions to plot.")
        arguments.add_argument("--save", action="store_true", help="Save canvas as pdf.")
        opts = arguments.parse_args()

        tree = readTree( opts.file )

        # see https://twiki.cern.ch/twiki/bin/viewauth/CMS/CutBasedPhotonID2012
        # for more information on 2012 Photon ID
        # TODO: Single tower H/E
        # TODO: correct isolation with ρ
        photonCut2012ID = "photon.sigmaIetaIeta<0.012 \
&& photon.chargedIso < 2.6 \
&& photon.neutralIso < 3.5 + 0.04*photon.pt \
&& photon.photonIso < 1.3 + 0.005*photon.pt"

        jetPhotonCut = "photon.sigmaIetaIeta > 0.012 && photon.sigmaIetaIeta<0.02 \
&& photon.chargedIso < 2.6 \
&& photon.neutralIso < 3.5 + 0.04*photon.pt \
&& photon.photonIso < 1.3 + 0.005*photon.pt"

        fillWeights( tree, jetPhotonCut, photonCut2012ID )

        if opts.cut == "photon":
                opts.cut = photonCut2012ID
        elif opts.cut == "jetphoton":
                opts.cut = jetPhotonCut

        if "all" in opts.distribution:
                opts.distribution = axisConf.sections()

        for plots in opts.distribution:
                if ":" in plots:
                        plot2D( tree, plots, opts.cut, opts.save )
                else:
                        plot( tree, plots, opts.cut, opts.save )

Revision:	1.2
Committed:	Mon Mar 25 12:02:29 2013 UTC (12 years, 1 month ago) by kiesel
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.1:	+32 -14 lines
Log Message:	help message fixed, library included
#	User	Rev	Content
1	kiesel	1.1	#! /usr/bin/env python2
2			# -- coding: utf-8 --
3
4			import ROOT
5			import argparse
6			import ConfigParser
7			import os.path
8	kiesel	1.2
9			# to use user defined help message, sys.arv has to be sent to python and not
10			# to TApplication
11			ROOT.PyConfig.IgnoreCommandLineOptions = True
12
13			# use tdr style defined in own macro
14	kiesel	1.1	import Styles
15			Styles.tdrStyle()
16
17	kiesel	1.2	# this is the configuration file for plots including label, etc.
18	kiesel	1.1	axisConf = ConfigParser.SafeConfigParser()
19			axisConf.read("axis.cfg")
20
21	kiesel	1.2	# global variables:
22	kiesel	1.1	integratedLuminosity = 19.3 #fb
23
24	kiesel	1.2	def readTree( filename, treename = "susyTree" ):
25	kiesel	1.1	"""
26			filename: name of file containing the tree
27			treename: name of the tree
28			returns: TTree Object
29			"""
30			if not os.path.isfile(filename):
31			print( "File %s does not exist"%filename)
32			tree = ROOT.TChain( treename )
33			tree.AddFile( filename )
34			return tree
35
36			def createHistoFromTree2D(tree, variable, weight="", nBin=50):
37			"""
38			tree: tree to create histo from
39			variable: variable to plot (must be a branch of the tree)
40			weight: weights to apply (e.g. "var1*(var2 > 15)" will use weights from var1 and cut on var2 > 15
41			nBins, firstBin, lastBin: number of bins, first bin and last bin (same as in TH1F constructor)
42			nEvents: number of events to process (-1 = all)
43			returns: histogram
44			"""
45			from random import randint
46			from sys import maxint
47			#make a random name you could give something meaningfull here,
48			#but that would make this less readable
49			name = "%x"%(randint(0, maxint))
50
51			# automatic binning if lastbin < firstbin
52			result = ROOT.TH2F(name, variable, nBin, 0, -1, nBin, 0, -1)
53			result.Sumw2()
54			tree.Draw("%s>>%s"%(variable, name), weight, "goff,colz")
55			return result
56
57			def createHistoFromTree(tree, variable, weight="", nBins=100, firstBin=None, lastBin=None, nEvents=-1):
58			"""
59			tree: tree to create histo from
60			variable: variable to plot (must be a branch of the tree)
61			weight: weights to apply (e.g. "var1*(var2 > 15)" will use weights from var1 and cut on var2 > 15
62			nBins, firstBin, lastBin: number of bins, first bin and last bin (same as in TH1F constructor)
63			nEvents: number of events to process (-1 = all)
64			returns: histogram
65			"""
66			if ":" in variable:
67			return createHistoFromTree2D(tree, variable, weight="")
68			from ROOT import TH1F
69			from random import randint
70			from sys import maxint
71			if nEvents < 0:
72			nEvents = maxint
73			if firstBin == None:
74			firstBin = tree.GetMinimum( variable )
75			if lastBin == None:
76			lastBin = tree.GetMaximum( variable )
77			#make a random name you could give something meaningfull here,
78			#but that would make this less readable
79			name = "%x"%(randint(0, maxint))
80			result = TH1F(name, variable, nBins, firstBin, lastBin)
81			result.Sumw2()
82			tree.Draw("%s>>%s"%(variable, name), weight, "goff", nEvents)
83			return result
84
85			def plot2D( tree, plot, cut, save=False ):
86			# read axis config file
87			var = plot.split(":")
88			try:
89			xlabel = axisConf.get( var[0], "label" )
90			except:
91			xlabel = ""
92			try:
93			ylabel = axisConf.get( var[1], "label" )
94			except:
95			ylabel = ""
96			try:
97			xunit = axisConf.get( var[0], "unit" )
98			except:
99			xunit = ""
100			try:
101			yunit = axisConf.get( var[1], "unit" )
102			except:
103			yunit = ""
104			if xunit != "":
105			xunit = " ["+xunit+"]"
106			if yunit != "":
107			yunit = " ["+yunit+"]"
108
109			# modify histo
110			histo = createHistoFromTree2D( tree, plot, cut )
111			histo.SetTitle(";%s%s;%s%s"%(ylabel,yunit,xlabel,xunit))
112
113			# draw histo
114			canvas = ROOT.TCanvas()
115			canvas.cd()
116			histo.Draw("colz")
117
118			# draw information
119			cutText = ROOT.TPaveText(.5,.8,.9,.9,"ndc")
120			cutText.SetBorderSize(0)
121			cutText.SetFillColor(0)
122			for line in cut.split("&&"):
123			cutText.AddText( line )
124			cutText.Draw()
125
126			topInfo = ROOT.TPaveText(.1,.95,1,1,"ndc")
127			topInfo.SetBorderSize(0)
128			topInfo.SetFillColor(0)
129			topInfo.AddText("Work in progress, %sfb^{-1}, #sqrt{s}=8TeV"%"?")
130			topInfo.Draw()
131
132			if save:
133			canvas.SaveAs("%.pdf"%variable)
134			else:
135			raw_input()
136
137			def fillWeights( tree, jetCut, photonCut ):
138
139			# this is the binning with witch the fo reweighting will be done and has to be studied
140			nBins = 50
141			xMin = 80
142			xMax = 300
143	kiesel	1.2	h_photonPt = createHistoFromTree( tree, "photon[0].pt", photonCut, nBins, xMin, xMax )
144			h_jetPt = createHistoFromTree( tree, "photon[0].pt", jetCut, nBins, xMin, xMax )
145	kiesel	1.1
146			# h_jetPt is now histogram with w^{-1} = h_jetPt / h_photonPt
147			h_jetPt.Divide( h_photonPt )
148	kiesel	1.2	import array
149			weight = array.array( "f", [0] )
150			tree.Branch( "weight2", weight, "weight2/F" )
151
152			for event in tree:
153			pt = tree.photon.at(0).pt
154			#print "p_T = %i GeV"%pt
155			weight = h_jetPt.GetBinContent( h_jetPt.FindBin( pt ) )
156			print weight
157			#print tree.weight2
158			#tree.Fill()
159	kiesel	1.1
160
161
162			def plot( tree, plot, cut, save=False ):
163			# read axis config file
164			try:
165			label = axisConf.get( plot, "label" )
166			unit = axisConf.get( plot, "unit" )
167			if unit != "":
168			unit = " ["+unit+"]"
169			except:
170			print "please specify label and unit in axis configuration file"
171			label = plot
172			unit = ""
173
174			# modify histo
175			histo = createHistoFromTree( tree, plot, cut )
176			histo.SetTitle(";%s%s;Entries"%(label,unit))
177			histo.SetLineColor(1)
178
179			# draw histo
180			canvas = ROOT.TCanvas()
181			canvas.cd()
182			canvas.SetLogy()
183			histo.Draw("hist")
184
185			# draw information
186			cutText = ROOT.TPaveText(.5,.8,.9,.9,"ndc")
187			cutText.SetBorderSize(0)
188			cutText.SetFillColor(0)
189			for line in cut.split("&&"):
190			cutText.AddText( line )
191			cutText.Draw()
192
193			topInfo = ROOT.TPaveText(.1,.95,1,1,"ndc")
194			topInfo.SetBorderSize(0)
195			topInfo.SetFillColor(0)
196			topInfo.AddText("Work in progress, %sfb^{-1}, #sqrt{s}=8TeV"%"?")
197			topInfo.Draw()
198
199			if save:
200	kiesel	1.2	canvas.SaveAs("%s.pdf"%plot)
201	kiesel	1.1	else:
202			raw_input()
203
204
205			if __name__ == "__main__":
206	kiesel	1.2	# include knowledge about objects saved in the tree
207			ROOT.gSystem.Load("../treeWriter/libTreeObjects.so")
208
209			arguments = argparse.ArgumentParser( description="Calculate weighting "
210			+"factors for QCD background estimation." )
211			arguments.add_argument("-f", "--file", default="../treeWriter/myTree.root",
212			help="ROOT file containing a TTree produced by 'treeWriter'.")
213			arguments.add_argument("-c", "--cut", default="", help="TCut string.")
214	kiesel	1.1	arguments.add_argument("-d", "--distribution", default =['met'], nargs="+",
215	kiesel	1.2	help="Distributions to plot.")
216	kiesel	1.1	arguments.add_argument("--save", action="store_true", help="Save canvas as pdf.")
217			opts = arguments.parse_args()
218
219			tree = readTree( opts.file )
220
221			# see https://twiki.cern.ch/twiki/bin/viewauth/CMS/CutBasedPhotonID2012
222			# for more information on 2012 Photon ID
223			# TODO: Single tower H/E
224			# TODO: correct isolation with ρ
225			photonCut2012ID = "photon.sigmaIetaIeta<0.012 \
226			&& photon.chargedIso < 2.6 \
227			&& photon.neutralIso < 3.5 + 0.04*photon.pt \
228			&& photon.photonIso < 1.3 + 0.005*photon.pt"
229
230			jetPhotonCut = "photon.sigmaIetaIeta > 0.012 && photon.sigmaIetaIeta<0.02 \
231			&& photon.chargedIso < 2.6 \
232			&& photon.neutralIso < 3.5 + 0.04*photon.pt \
233			&& photon.photonIso < 1.3 + 0.005*photon.pt"
234
235			fillWeights( tree, jetPhotonCut, photonCut2012ID )
236
237			if opts.cut == "photon":
238			opts.cut = photonCut2012ID
239			elif opts.cut == "jetphoton":
240			opts.cut = jetPhotonCut
241
242			if "all" in opts.distribution:
243			opts.distribution = axisConf.sections()
244
245			for plots in opts.distribution:
246			if ":" in plots:
247			plot2D( tree, plots, opts.cut, opts.save )
248			else:
249			plot( tree, plots, opts.cut, opts.save )
250