ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/fitMCToData.py

Comparing UserCode/OSUT3Analysis/Configuration/scripts/fitMCToData.py (file contents):
Revision 1.1 by ahart, Sun Apr 7 00:57:26 2013 UTC vs.
Revision 1.5 by ahart, Fri Aug 23 16:22:50 2013 UTC

#	Line 2 | Line 2
2		import sys
3		import os
4		import re
5	<	from optparse import OptionParser
5	>	import time
6	>	from math import *
7		from array import *
8		from decimal import *
9	<
9	>	from optparse import OptionParser
10		from OSUT3Analysis.Configuration.configurationOptions import *
11		from OSUT3Analysis.Configuration.processingUtilities import *
12	+	from OSUT3Analysis.Configuration.formattingUtilities import *
13	+
14
15	<	##set default plotting options
16	<	line_width = 2
14	<	plotting_options = ""
15	>
16	>	### parse the command-line options
17
18		parser = OptionParser()
19		parser = set_commandline_arguments(parser)
20
21	<	parser.add_option("--hist", action="store_true", dest="plot_hist", default=False,
22	<	help="plot as hollow histograms instead of error bar crosses")
23	<	parser.add_option("--line-width", dest="line_width",
24	<	help="set line width (default is 2)")
25	<	parser.add_option("--pdf", action="store_true", dest="plot_savePdf", default=False,
26	<	help="save plot as pdf in addition")
21	>	parser.add_option("-f", "--fancy", action="store_true", dest="makeFancy", default=False,
22	>	help="removes the title and replaces it with the official CMS plot heading")
23	>	parser.add_option("--ylog", action="store_true", dest="setLogY", default=False,
24	>	help="Set logarithmic scale on vertical axis on all plots")
25	>	parser.add_option("--ymin", dest="setYMin",
26	>	help="Minimum of y axis")
27	>	parser.add_option("--ymax", dest="setYMax",
28	>	help="Maximum of y axis")
29	>	parser.add_option("-E", "--ratioRelErrMax", dest="ratioRelErrMax",
30	>	help="maximum error used in rebinning the ratio histogram")
31
32		(arguments, args) = parser.parse_args()
33
34	+
35		if arguments.localConfig:
36	<	sys.path.insert(0,os.getcwd())
36	>	sys.path.append(os.getcwd())
37		exec("from " + arguments.localConfig.rstrip('.py') + " import *")
38
39	<	outputFileName = "mc_fit_to_data.root"
40	<	if arguments.outputFileName:
41	<	outputFileName = arguments.outputFileName
42	<	pdfFileName = outputFileName[:-5] + ".pdf"
43	<
44	<	condor_dir = set_condor_output_dir(arguments)
45	<
39	>	#### deal with conflicting arguments
40	>	if arguments.normalizeToData and arguments.normalizeToUnitArea:
41	>	print "Conflicting normalizations requsted, will normalize to unit area"
42	>	arguments.normalizeToData = False
43	>	if arguments.normalizeToData and arguments.noStack:
44	>	print "You have asked to scale non-stacked backgrounds to data.  This is a very strange request.  Will normalize to unit area instead"
45	>	arguments.normalizeToData = False
46	>	arguments.normalizeToUnitArea = True
47		if arguments.makeRatioPlots and arguments.makeDiffPlots:
48		print "You have requested both ratio and difference plots.  Will make just ratio plots instead"
49		arguments.makeRatioPlots = False
50	+	if arguments.makeRatioPlots and arguments.noStack:
51	+	print "You have asked to make a ratio plot and to not stack the backgrounds.  This is a very strange request.  Will skip making the ratio plot."
52	+	arguments.makeRatioPlots = False
53	+	if arguments.makeDiffPlots and arguments.noStack:
54	+	print "You have asked to make a difference plot and to not stack the backgrounds.  This is a very strange request.  Will skip making the difference plot."
55	+	arguments.makeDiffPlots = False
56	+
57	+
58	+	from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, THStack, TH1F, TCanvas, TString, TLegend, TLegendEntry, THStack, TIter, TKey, TPaveLabel, TPaveText, TF1, gPad
59	+
60	+
61	+	### setting ROOT options so our plots will look awesome and everyone will love us
62
43	–	from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, THStack, TH1F, TCanvas, TString, TLegend, TArrow, THStack, TIter, TKey, TPaveLabel, TPaveText, TF1, gPad
44	–	sys.argv = []
63		gROOT.SetBatch()
64		gStyle.SetOptStat(0)
65		gStyle.SetCanvasBorderMode(0)
#	Line 49 | Line 67 | gStyle.SetPadBorderMode(0)
67		gStyle.SetPadColor(0)
68		gStyle.SetCanvasColor(0)
69		gStyle.SetTextFont(42)
70	+	gStyle.SetCanvasDefH(600)
71	+	gStyle.SetCanvasDefW(600)
72	+	gStyle.SetCanvasDefX(0)
73	+	gStyle.SetCanvasDefY(0)
74	+	gStyle.SetPadTopMargin(0.07)
75	+	gStyle.SetPadBottomMargin(0.13)
76	+	gStyle.SetPadLeftMargin(0.15)
77	+	gStyle.SetPadRightMargin(0.05)
78	+	gStyle.SetTitleColor(1, "XYZ")
79	+	gStyle.SetTitleFont(42, "XYZ")
80	+	gStyle.SetTitleSize(0.04, "XYZ")
81	+	gStyle.SetTitleXOffset(1.1)
82	+	gStyle.SetTitleYOffset(2)
83	+	gStyle.SetTextAlign(12)
84	+	gStyle.SetLabelColor(1, "XYZ")
85	+	gStyle.SetLabelFont(42, "XYZ")
86	+	gStyle.SetLabelOffset(0.007, "XYZ")
87	+	gStyle.SetLabelSize(0.04, "XYZ")
88	+	gStyle.SetAxisColor(1, "XYZ")
89	+	gStyle.SetStripDecimals(True)
90	+	gStyle.SetTickLength(0.03, "XYZ")
91	+	gStyle.SetNdivisions(510, "XYZ")
92	+	gStyle.SetPadTickX(1)
93	+	gStyle.SetPadTickY(1)
94		gROOT.ForceStyle()
95
54	–	outputFile = TFile(outputFileName, "RECREATE")
96
97	<	datasets_needed = []
98	<	for histogram in input_histograms:
99	<	for dataset in histogram['datasets']:
100	<	if dataset not in datasets_needed:
101	<	datasets_needed.append(dataset)
102	<	if histogram['target_dataset'] not in datasets_needed:
103	<	datasets_needed.append(histogram['target_dataset'])
104	<
105	<	weight = intLumi / 10000.0
106	<	for dataset in datasets_needed:
107	<	dataset_file = "%s/%s.root" % (condor_dir,dataset)
108	<	fin = TFile (dataset_file)
109	<	flags = fin.Get ("flags")
110	<	noWeights = flags and flags.GetBinContent (1)
111	<	fin.Close ()
97	>	#set the text for the luminosity label
98	>	if(intLumi < 1000.):
99	>	LumiInPb = intLumi
100	>	LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
101	>	LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInPb) + " pb^{-1}"
102	>	else:
103	>	LumiInFb = intLumi/1000.
104	>	LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInFb) + " fb^{-1}"
105	>
106	>	#bestest place for lumi. label, in top left corner
107	>	topLeft_x_left    = 0.1375839
108	>	topLeft_x_right   = 0.4580537
109	>	topLeft_y_bottom  = 0.8479021
110	>	topLeft_y_top     = 0.9475524
111	>	topLeft_y_offset  = 0.035
112	>
113	>	#set the text for the fancy heading
114	>	HeaderText = "CMS Preliminary: " + LumiText + " at #sqrt{s} = 8 TeV"
115	>
116	>	#position for header
117	>	header_x_left    = 0.2181208
118	>	header_x_right   = 0.9562937
119	>	header_y_bottom  = 0.9479866
120	>	header_y_top     = 0.9947552
121
72	–	if types[dataset] != "data" and not noWeights:
73	–	os.system("mergeTFileServiceHistograms -i %s -o %s -w %g" % (dataset_file, dataset_file + "_tmp", weight))
74	–	else:
75	–	os.system("mergeTFileServiceHistograms -i %s -o %s -w %g" % (dataset_file, dataset_file + "_tmp", 1.0))
122
77	–	for histogram in input_histograms:
78	–	HistogramsToFit = []
79	–	HistogramDatasets = []
80	–	TargetDataset = histogram['target_dataset']
123
124	<	Stack = []
125	<	Stack.append (THStack("stack_before",histogram['name']))
126	<	Stack.append (THStack("stack_after",histogram['name']))
124	>	##########################################################################################################################################
125	>	##########################################################################################################################################
126	>	##########################################################################################################################################
127	>
128	>	# some fancy-ass code from Andrzej Zuranski to merge bins in the ratio plot until the error goes below some threshold
129	>	def ratioHistogram( dataHist, mcHist, relErrMax=0.10):
130	>
131	>	if not dataHist:
132	>	print "Error:  trying to run ratioHistogram but dataHist is invalid"
133	>	return
134	>
135	>	if not mcHist:
136	>	print "Error:  trying to run ratioHistogram but mcHist is invalid"
137	>	return
138	>
139	>	def groupR(group):
140	>	Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]]
141	>	return (Data-MC)/MC if MC else 0
142	>
143	>	def groupErr(group):
144	>	Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]]
145	>	dataErr2,mcErr2 = [sum(hist.GetBinError(i)**2 for i in group) for hist in [dataHist,mcHist]]
146	>	if Data > 0 and MC > 0 and Data != MC:
147	>	return abs(math.sqrt( (dataErr2+mcErr2)/(Data-MC)**2 + mcErr2/MC**2 ) * (Data-MC)/MC)
148	>	else:
149	>	return 0
150
151	<	if(intLumi < 1000.):
152	<	LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
151	>	def regroup(groups):
152	>	err,iG = max( (groupErr(g),groups.index(g)) for g in groups )
153	>	if err < relErrMax or len(groups)<3 : return groups
154	>	iH = max( [iG-1,iG+1], key = lambda i: groupErr(groups[i]) if 0<=i<len(groups) else -1 )
155	>	iLo,iHi = sorted([iG,iH])
156	>	return regroup(groups[:iLo] + [groups[iLo]+groups[iHi]] + groups[iHi+1:])
157	>
158	>	#don't rebin the histograms of the number of a given object (except for the pileup ones)
159	>	if ((dataHist.GetName().find("num") is not -1 and dataHist.GetName().find("Primaryvertexs") is -1) or
160	>	dataHist.GetName().find("CutFlow")  is not -1 or
161	>	dataHist.GetName().find("GenMatch") is not -1):
162	>	ratio = dataHist.Clone()
163	>	ratio.Add(mcHist,-1)
164	>	ratio.Divide(mcHist)
165	>	ratio.SetTitle("")
166		else:
167	<	getcontext().prec = 2
168	<	LumiInFb = intLumi/1000.
169	<	LumiText = "L_{int} = " + str(LumiInFb) + " fb^{-1}"
167	>	groups = regroup( [(i,) for i in range(1,1+dataHist.GetNbinsX())] )
168	>	ratio = TH1F("ratio","",len(groups), array('d', [dataHist.GetBinLowEdge(min(g)) for g in groups ] + [dataHist.GetXaxis().GetBinUpEdge(dataHist.GetNbinsX())]) )
169	>	for i,g in enumerate(groups) :
170	>	ratio.SetBinContent(i+1,groupR(g))
171	>	ratio.SetBinError(i+1,groupErr(g))
172	>
173	>	ratio.GetYaxis().SetTitle("#frac{Data-MC}{MC}")
174	>	ratio.SetLineColor(1)
175	>	ratio.SetLineWidth(2)
176	>	return ratio
177	>
178	>	##########################################################################################################################################
179	>	##########################################################################################################################################
180	>	##########################################################################################################################################
181	>
182	>
183	>	def MakeOneDHist(pathToDir,distribution):
184	>
185	>	numFittingSamples = 0
186	>
187	>	HeaderLabel = TPaveLabel(header_x_left,header_y_bottom,header_x_right,header_y_top,HeaderText,"NDC")
188	>	HeaderLabel.SetTextAlign(32)
189	>	HeaderLabel.SetBorderSize(0)
190	>	HeaderLabel.SetFillColor(0)
191	>	HeaderLabel.SetFillStyle(0)
192
193	<	LumiLabel = TPaveLabel(0.1,0.8,0.34,0.9,LumiText,"NDC")
193	>	LumiLabel = TPaveLabel(topLeft_x_left,topLeft_y_bottom,topLeft_x_right,topLeft_y_top,LumiText,"NDC")
194		LumiLabel.SetBorderSize(0)
195		LumiLabel.SetFillColor(0)
196		LumiLabel.SetFillStyle(0)
197	+
198	+	NormLabel = TPaveLabel()
199	+	NormLabel.SetDrawOption("NDC")
200	+	NormLabel.SetX1NDC(topLeft_x_left)
201	+	NormLabel.SetX2NDC(topLeft_x_right)
202	+
203	+	NormLabel.SetBorderSize(0)
204	+	NormLabel.SetFillColor(0)
205	+	NormLabel.SetFillStyle(0)
206	+
207	+	NormText = ""
208	+	if arguments.normalizeToUnitArea:
209	+	NormText = "Scaled to unit area"
210	+	elif arguments.normalizeToData:
211	+	NormText = "MC scaled to data"
212	+	NormLabel.SetLabel(NormText)
213	+
214	+	YieldsLabel = TPaveText(0.39, 0.7, 0.59, 0.9,"NDC")
215	+	YieldsLabel.SetBorderSize(0)
216	+	YieldsLabel.SetFillColor(0)
217	+	YieldsLabel.SetFillStyle(0)
218	+	YieldsLabel.SetTextAlign(12)
219	+
220	+	RatiosLabel = TPaveText()
221	+	RatiosLabel.SetDrawOption("NDC")
222	+	RatiosLabel.SetBorderSize(0)
223	+	RatiosLabel.SetFillColor(0)
224	+	RatiosLabel.SetFillStyle(0)
225	+	RatiosLabel.SetTextAlign(32)
226	+
227	+
228	+	Legend = TLegend()
229	+	Legend.SetBorderSize(0)
230	+	Legend.SetFillColor(0)
231	+	Legend.SetFillStyle(0)
232	+
233	+
234	+	#    outputFile.cd(pathToDir)
235	+
236	+	fittingIntegral = 0
237	+	scaleFactor = 1
238
239	<	Label = TPaveText(0.34, 0.7, 0.64, 0.9,"NDC")
240	<	Label.SetBorderSize(0)
100	<	Label.SetFillColor(0)
101	<	Label.SetFillStyle(0)
102	<	Label.SetTextAlign(12)
103	<
104	<	BgMCLegend = TLegend(0.70,0.65,0.94,0.89, "Data & Bkgd. MC")
105	<	BgMCLegend.SetBorderSize(0)
106	<	BgMCLegend.SetFillColor(0)
107	<	BgMCLegend.SetFillStyle(0)
239	>	HistogramsToFit = []
240	>	TargetDataset = distribution['target_dataset']
241
242	<	scaleFactor = 1
242	>	FittingLegendEntries = []
243	>	DataLegendEntries = []
244
245	<	numBgMCSamples = 0
112	<	numDataSamples = 0
113	<	numSignalSamples = 0
245	>	FittingHistogramDatasets = []
246
247	<	fileName = condor_dir + "/" + histogram['target_dataset'] + ".root_tmp"
247	>
248	>	Stack_list = []
249	>	Stack_list.append (THStack("stack_before",distribution['name']))
250	>	Stack_list.append (THStack("stack_after",distribution['name']))
251	>
252	>	fileName = condor_dir + "/" + distribution['target_dataset'] + ".root"
253		if not os.path.exists(fileName):
254	<	continue
254	>	return
255		inputFile = TFile(fileName)
256		if inputFile.IsZombie() or not inputFile.GetNkeys():
257	<	continue
257	>	return
258
122	–	Target = inputFile.Get("OSUAnalysis/"+histogram['channel']+"/"+histogram['name']).Clone()
123	–	Target.SetDirectory(0)
124	–	inputFile.Close()
259
126	–	numDataSamples += 1
260
261	+	Target = inputFile.Get("OSUAnalysis/"+distribution['channel']+"/"+distribution['name']).Clone()
262	+	Target.SetDirectory(0)
263	+	inputFile.Close()
264	+
265	+	Target.SetMarkerStyle(20)
266	+	Target.SetMarkerSize(0.8)
267		Target.SetFillStyle(0)
268		Target.SetLineColor(colors[TargetDataset])
269		Target.SetLineStyle(1)
270		Target.SetLineWidth(2)
271	<
272	<	BgMCLegend.AddEntry(Target,labels[TargetDataset],"LEP")
273	<
274	<	xAxisLabel = Target.GetXaxis().GetTitle()
275	<	histoTitle = Target.GetTitle()
276	<
271	>	targetIntegral = Target.Integral()
272	>	if(arguments.normalizeToUnitArea and Target.Integral() > 0):
273	>	Target.Scale(1./Target.Integral())
274	>
275	>	### formatting target histogram and adding to legend
276	>	legendIndex = 0
277	>	Legend.AddEntry(Target,labels[TargetDataset],"LEP")
278	>	legendIndex = legendIndex+1
279	>
280		if not outputFile.Get ("OSUAnalysis"):
281		outputFile.mkdir ("OSUAnalysis")
282	<	if not outputFile.Get ("OSUAnalysis/" + histogram['channel']):
283	<	outputFile.Get ("OSUAnalysis").mkdir (histogram['channel'])
284	<
285	<	for dataset in histogram['datasets']:
286	<	fileName = condor_dir + "/" + dataset + ".root_tmp"
287	<	if not os.path.exists(fileName):
288	<	continue
289	<	inputFile = TFile(fileName)
290	<	if inputFile.IsZombie() or not inputFile.GetNkeys():
291	<	continue
292	<
293	<	Histogram = inputFile.Get("OSUAnalysis/"+histogram['channel']+"/"+histogram['name']).Clone()
282	>	if not outputFile.Get ("OSUAnalysis/" + distribution['channel']):
283	>	outputFile.Get ("OSUAnalysis").mkdir (distribution['channel'])
284	>
285	>	for sample in distribution['datasets']: # loop over different samples requested to be fit
286	>
287	>	dataset_file = "%s/%s.root" % (condor_dir,sample)
288	>	inputFile = TFile(dataset_file)
289	>	HistogramObj = inputFile.Get(pathToDir+"/"+distribution['channel']+"/"+distribution['name'])
290	>	if not HistogramObj:
291	>	print "WARNING:  Could not find histogram " + pathToDir + "/" + distribution['name'] + " in file " + dataset_file + ".  Will skip it and continue."
292	>	continue
293	>	Histogram = HistogramObj.Clone()
294		Histogram.SetDirectory(0)
295		inputFile.Close()
154	–
296		if arguments.rebinFactor:
297		RebinFactor = int(arguments.rebinFactor)
298	<	if Histogram.GetNbinsX() >= RebinFactor*10:
298	>	#don't rebin histograms which will have less than 5 bins or any gen-matching histograms
299	>	if Histogram.GetNbinsX() >= RebinFactor*5 and Histogram.GetName().find("GenMatch") is -1:
300		Histogram.Rebin(RebinFactor)
301
302	<	numBgMCSamples += 1
302	>	xAxisLabel = Histogram.GetXaxis().GetTitle()
303	>	unitBeginIndex = xAxisLabel.find("[")
304	>	unitEndIndex = xAxisLabel.find("]")
305	>
306	>	if unitBeginIndex is not -1 and unitEndIndex is not -1: #x axis has a unit
307	>	yAxisLabel = "Entries / " + str(Histogram.GetXaxis().GetBinWidth(1)) + " " + xAxisLabel[unitBeginIndex+1:unitEndIndex]
308	>	else:
309	>	yAxisLabel = "Entries per bin (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " width)"
310	>
311	>	if not arguments.makeFancy:
312	>	histoTitle = Histogram.GetTitle()
313	>	else:
314	>	histoTitle = ""
315	>
316
317	<	if(arguments.noStack):
317	>	legLabel = labels[sample]
318	>	if (arguments.printYields):
319	>	yieldHist = Histogram.Integral()
320	>	legLabel = legLabel + " (%.1f)" % yieldHist
321	>	FittingLegendEntries.append(legLabel)
322	>
323	>	if( types[sample] == "bgMC"):
324	>
325	>	numFittingSamples += 1
326	>	fittingIntegral += Histogram.Integral()
327	>
328	>	Histogram.SetLineStyle(1)
329	>	if(arguments.noStack):
330	>	Histogram.SetFillStyle(0)
331	>	Histogram.SetLineColor(colors[sample])
332	>	Histogram.SetLineWidth(2)
333	>	else:
334	>	Histogram.SetFillStyle(1001)
335	>	Histogram.SetFillColor(colors[sample])
336	>	Histogram.SetLineColor(1)
337	>	Histogram.SetLineWidth(1)
338	>
339	>	elif( types[sample] == "signalMC"):
340	>
341	>	numFittingSamples += 1
342	>
343		Histogram.SetFillStyle(0)
344	<	Histogram.SetLineColor(colors[dataset])
344	>	Histogram.SetLineColor(colors[sample])
345	>	Histogram.SetLineStyle(1)
346		Histogram.SetLineWidth(2)
347	<	BgMCLegend.AddEntry(Histogram,labels[dataset],"L")
348	<	else:
349	<	Histogram.SetFillStyle(1001)
169	<	Histogram.SetFillColor(colors[dataset])
170	<	Histogram.SetLineColor(1)
171	<	Histogram.SetLineWidth(1)
172	<	BgMCLegend.AddEntry(Histogram,labels[dataset],"F")
173	<	Histogram.SetLineStyle(1)
174	<
347	>	if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
348	>	Histogram.Scale(1./Histogram.Integral())
349	>
350		HistogramsToFit.append(Histogram)
351	<	HistogramDatasets.append(dataset)
351	>	FittingHistogramDatasets.append(sample)
352	>
353	>	#scaling histograms as per user's specifications
354	>	if targetIntegral > 0 and fittingIntegral > 0:
355	>	scaleFactor = targetIntegral/fittingIntegral
356	>	for fittingHist in HistogramsToFit:
357	>	if arguments.normalizeToData:
358	>	fittingHist.Scale(scaleFactor)
359	>
360	>	if arguments.normalizeToUnitArea and not arguments.noStack and fittingIntegral > 0:
361	>	fittingHist.Scale(1./fittingIntegral)
362	>	elif arguments.normalizeToUnitArea and arguments.noStack and fittingHist.Integral() > 0:
363	>	fittingHist.Scale(1./fittingHist.Integral())
364	>
365
366		def fitf (x, par):
367		xBin = HistogramsToFit[0].FindBin (x[0])
368		value = 0.0
369	<
369	>	sumOfWeights = 0.0
370	>
371		for i in range (0, len (HistogramsToFit)):
372	<	value += par[i] * HistogramsToFit[i].GetBinContent (xBin)
373	<
372	>	weight = 1.0 / (HistogramsToFit[i].GetBinError (xBin) * HistogramsToFit[i].GetBinError (xBin))
373	>	sumOfWeights += weight
374	>	value += weight * par[i] * HistogramsToFit[i].GetBinContent (xBin)
375	>	value /= sumOfWeights
376	>
377		return value
378
379		lowerLimit = Target.GetBinLowEdge (1)
380		upperLimit = Target.GetBinLowEdge (Target.GetNbinsX ()) + Target.GetBinWidth (Target.GetNbinsX ())
381	<	if 'lowerLimit' in histogram:
382	<	lowerLimit = histogram['lowerLimit']
383	<	if 'upperLimit' in histogram:
384	<	upperLimit = histogram['upperLimit']
381	>	if 'lowerLimit' in distribution:
382	>	lowerLimit = distribution['lowerLimit']
383	>	if 'upperLimit' in distribution:
384	>	upperLimit = distribution['upperLimit']
385		func = TF1 ("fit", fitf, lowerLimit, upperLimit, len (HistogramsToFit))
386
387		for i in range (0, len (HistogramsToFit)):
388		func.SetParameter (i, 1.0)
389	<	func.SetParLimits (i, 0.0, 10.0)
198	<	func.SetParName (i, labels[HistogramDatasets[i]])
389	>	func.SetParName (i, labels[FittingHistogramDatasets[i]])
390
391	<	for i in range (0, histogram['iterations'] - 1):
391	>	for i in range (0, distribution['iterations'] - 1):
392		print "Iteration " + str (i + 1) + "..."
393		Target.Fit ("fit", "QEMR0")
394		Target.Fit ("fit", "VEMR0")
395
205	–	func.SetLineWidth (line_width)
206	–	func.SetLineColor (632)
207	–
208	–	Target.SetMarkerStyle (20)
209	–	Target.SetMarkerColor (1)
210	–	Target.SetLineColor (1)
396
397	<	Canvas = TCanvas(histogram['name'])
397	>	finalMax = 0
398	>	if not arguments.noStack:
399	>	for fittingHist in HistogramsToFit:
400	>	finalMax += fittingHist.GetMaximum()
401	>	else:
402	>	for fittingHist in HistogramsToFit:
403	>	if(fittingHist.GetMaximum() > finalMax):
404	>	finalMax = fittingHist.GetMaximum()
405	>	if(Target.GetMaximum() > finalMax):
406	>	finalMax = Target.GetMaximum()
407	>
408	>	Target.SetMaximum(1.1*finalMax)
409	>	Target.SetMinimum(0.0001)
410	>
411	>	Canvas = TCanvas(distribution['name'] + "_FitFunction")
412		Canvas.cd (1)
413		Target.Draw ()
414	<	func.Draw (plotting_options + "same")
415	<
416	<	outputFile.cd ("OSUAnalysis/" + histogram['channel'])
414	>	func.Draw ("same")
415	>
416	>	outputFile.cd ("OSUAnalysis/" + distribution['channel'])
417		Canvas.Write ()
418	<	if arguments.plot_savePdf:
419	<	if histogram == input_histograms[0]:
420	<	Canvas.Print (pdfFileName + "(", "pdf")
421	<	else:
422	<	Canvas.Print (pdfFileName, "pdf")
418	>	if arguments.savePDFs:
419	>	if histogram == input_histograms[0]:
420	>	Canvas.Print (pdfFileName + "(", "pdf")
421	>	else:
422	>	Canvas.Print (pdfFileName, "pdf")
423		Target.SetStats (0)
424
425	<	for i in range (0, 2):
425	>
426	>
427	>
428	>	### formatting bgMC histograms and adding to legend
429	>	legendIndex = numFittingSamples-1
430	>	for Histogram in reversed(HistogramsToFit):
431	>	if(arguments.noStack):
432	>	Legend.AddEntry(Histogram,FittingLegendEntries[legendIndex],"L")
433	>	else:
434	>	Legend.AddEntry(Histogram,FittingLegendEntries[legendIndex],"F")
435	>	legendIndex = legendIndex-1
436	>
437	>
438	>	### Drawing histograms to canvas
439	>
440	>	makeRatioPlots = arguments.makeRatioPlots
441	>	makeDiffPlots = arguments.makeDiffPlots
442	>
443	>	yAxisMin = 0.0001
444	>	if arguments.setYMin:
445	>	yAxisMin = float(arguments.setYMin)
446	>
447	>
448	>	### Draw everything to the canvases !!!!
449	>
450	>	for i in range (0, 2): # 0 => before, 1 => after
451
452		if i == 1:
453	+	ratios = []
454		for j in range (0, len (HistogramsToFit)):
455		HistogramsToFit[j].Scale (func.GetParameter (j))
456	+	ratios.append(func.GetParameter (j))
457
458	<	for bgMCHist in HistogramsToFit:
458	>	for fittingHist in HistogramsToFit:
459		if not arguments.noStack:
460	<	Stack[i].Add(bgMCHist)
460	>	Stack_list[i].Add(fittingHist)
461
236	–	finalMax = 0
237	–	if not arguments.noStack:
238	–	finalMax = Stack[i].GetMaximum()
239	–	else:
240	–	for bgMCHist in HistogramsToFit:
241	–	if(bgMCHist.GetMaximum() > finalMax):
242	–	finalMax = bgMCHist.GetMaximum()
243	–	if(Target.GetMaximum() > finalMax):
244	–	finalMax = Target.GetMaximum()
462
463	<	makeRatioPlots = arguments.makeRatioPlots
464	<	makeDiffPlots = arguments.makeDiffPlots
463	>	#creating the histogram to represent the statistical errors on the stack
464	>	if not arguments.noStack:
465	>	ErrorHisto = HistogramsToFit[0].Clone("errors")
466	>	ErrorHisto.SetFillStyle(3001)
467	>	ErrorHisto.SetFillColor(13)
468	>	ErrorHisto.SetLineWidth(0)
469	>	if i == 1:
470	>	Legend.AddEntry(ErrorHisto,"Stat. Errors","F")
471	>	for Histogram in HistogramsToFit:
472	>	if Histogram is not HistogramsToFit[0]:
473	>	ErrorHisto.Add(Histogram)
474	>
475		if i == 0:
476	<	Canvas = TCanvas(histogram['name'] + "_Before")
476	>	Canvas = TCanvas(distribution['name'] + "_Before")
477		if i == 1:
478	<	Canvas = TCanvas(histogram['name'] + "_After")
478	>	Canvas = TCanvas(distribution['name'] + "_After")
479	>
480	>
481		if makeRatioPlots or makeDiffPlots:
482		Canvas.SetFillStyle(0)
483		Canvas.Divide(1,2)
484		Canvas.cd(1)
485	<	gPad.SetPad(0.01,0.25,0.99,0.99)
486	<	gPad.SetMargin(0.1,0.05,0.02,0.07)
485	>	gPad.SetPad(0,0.25,1,1)
486	>	gPad.SetMargin(0.15,0.05,0.01,0.07)
487		gPad.SetFillStyle(0)
488		gPad.Update()
489		gPad.Draw()
490	+	if arguments.setLogY:
491	+	gPad.SetLogy()
492		Canvas.cd(2)
493	<	gPad.SetPad(0.01,0.01,0.99,0.25)
494	<	#format: gPad.SetMargin(l,r,b,t)
495	<	gPad.SetMargin(0.1,0.05,0.4,0.02)
493	>	gPad.SetPad(0,0,1,0.25)
494	>	# format: gPad.SetMargin(l,r,b,t)
495	>	gPad.SetMargin(0.15,0.05,0.4,0.01)
496		gPad.SetFillStyle(0)
497		gPad.SetGridy(1)
498		gPad.Update()
499		gPad.Draw()
500	<
500	>
501		Canvas.cd(1)
502
503	<	if not arguments.noStack:
504	<	Stack[i].SetTitle(histoTitle)
505	<	Stack[i].Draw("HIST")
506	<	Stack[i].GetXaxis().SetTitle(xAxisLabel)
276	<	Stack[i].SetMaximum(1.1*finalMax)
277	<	Stack[i].SetMinimum(0.0001)
278	<	if makeRatioPlots or makeDiffPlots:
279	<	Stack[i].GetHistogram().GetXaxis().SetLabelSize(0)
503	>	### finding the maximum value of anything going on the canvas, so we know how to set the y-axis
504	>	finalMax = 0
505	>	if numFittingSamples is not 0 and not arguments.noStack:
506	>	finalMax = ErrorHisto.GetMaximum() + ErrorHisto.GetBinError(ErrorHisto.GetMaximumBin())
507		else:
508	+	for bgMCHist in HistogramsToFit:
509	+	if(bgMCHist.GetMaximum() > finalMax):
510	+	finalMax = bgMCHist.GetMaximum()
511	+	if(Target.GetMaximum() > finalMax):
512	+	finalMax = Target.GetMaximum() + Target.GetBinError(Target.GetMaximumBin())
513	+	finalMax = 1.15*finalMax
514	+	if arguments.setYMax:
515	+	finalMax = float(arguments.setYMax)
516	+
517	+
518	+	if not arguments.noStack: # draw stacked background samples
519	+	Stack_list[i].SetTitle(histoTitle)
520	+	Stack_list[i].Draw("HIST")
521	+	Stack_list[i].GetXaxis().SetTitle(xAxisLabel)
522	+	Stack_list[i].GetYaxis().SetTitle(yAxisLabel)
523	+	Stack_list[i].SetMaximum(finalMax)
524	+	Stack_list[i].SetMinimum(yAxisMin)
525	+	if makeRatioPlots or makeDiffPlots:
526	+	Stack_list[i].GetHistogram().GetXaxis().SetLabelSize(0)
527	+	#draw shaded error bands
528	+	ErrorHisto.Draw("A E2 SAME")
529	+
530	+	else: #draw the unstacked backgrounds
531		HistogramsToFit[0].SetTitle(histoTitle)
532		HistogramsToFit[0].Draw("HIST")
533		HistogramsToFit[0].GetXaxis().SetTitle(xAxisLabel)
534	<	HistogramsToFit[0].SetMaximum(1.1*finalMax)
535	<	HistogramsToFit[0].SetMinimum(0.0001)
534	>	HistogramsToFit[0].GetYaxis().SetTitle(yAxisLabel)
535	>	HistogramsToFit[0].SetMaximum(finalMax)
536	>	HistogramsToFit[0].SetMinimum(yAxisMin)
537		for bgMCHist in HistogramsToFit:
538	<	bgMCHist.Draw("HIST SAME")
538	>	bgMCHist.Draw("A HIST SAME")
539
540	<	dataYield = Target.Integral (1, Target.GetNbinsX ())
290	<	mcYield = 0.0
291	<	for bgMCHist in HistogramsToFit:
292	<	mcYield += bgMCHist.Integral (1, bgMCHist.GetNbinsX ())
293	<	Label.Clear ()
294	<	if i == 0:
295	<	Label.AddText ("Before Fit to Data")
296	<	if i == 1:
297	<	Label.AddText ("After Fit to Data")
298	<	Label.AddText ("data yield: " + str (dataYield))
299	<	Label.AddText ("MC yield: " + str (mcYield))
300	<
301	<	Target.Draw("E SAME")
302	<	BgMCLegend.Draw()
303	<	LumiLabel.Draw()
304	<	Label.Draw()
540	>	Target.Draw("A E X0 SAME")
541
542	<	if makeRatioPlots or makeDiffPlots:
542	>
543	>
544	>	#legend coordinates, empirically determined :-)
545	>	x_left = 0.6761745
546	>	x_right = 0.9328859
547	>	x_width = x_right - x_left
548	>	y_max = 0.9
549	>	entry_height = 0.05
550	>
551	>	if(numFittingSamples is not 0): #then draw the data & bgMC legend
552	>
553	>	numExtraEntries = 2 # count the target and (lack of) title
554	>	Legend.SetX1NDC(x_left)
555	>	numExtraEntries = numExtraEntries + 1 # count the stat. errors entry
556	>
557	>	Legend.SetY1NDC(y_max-entry_height*(numExtraEntries+numFittingSamples))
558	>	Legend.SetX2NDC(x_right)
559	>	Legend.SetY2NDC(y_max)
560	>	Legend.Draw()
561	>
562	>	RatiosLabel.SetX1NDC(x_left - 0.1)
563	>	RatiosLabel.SetX2NDC(x_right)
564	>	RatiosLabel.SetY2NDC(Legend.GetY1NDC() - 0.1)
565	>	RatiosLabel.SetY1NDC(RatiosLabel.GetY2NDC() - entry_height*(numFittingSamples))
566	>
567	>	# Deciding which text labels to draw and drawing them
568	>	drawLumiLabel = False
569	>	drawNormLabel = False
570	>	offsetNormLabel = False
571	>	drawHeaderLabel = False
572	>
573	>	if not arguments.normalizeToUnitArea: #don't draw the lumi label if there's no data and it's scaled to unit area
574	>	drawLumiLabel = True
575	>	# move the normalization label down before drawing if we drew the lumi. label
576	>	offsetNormLabel = True
577	>	if arguments.normalizeToUnitArea or arguments.normalizeToData:
578	>	drawNormLabel = True
579	>	if arguments.makeFancy:
580	>	drawHeaderLabel = True
581	>	drawLumiLabel = False
582	>
583	>	# now that flags are set, draw the appropriate labels
584	>
585	>	if drawLumiLabel:
586	>	LumiLabel.Draw()
587	>
588	>	if drawNormLabel:
589	>	if offsetNormLabel:
590	>	NormLabel.SetY1NDC(topLeft_y_bottom-topLeft_y_offset)
591	>	NormLabel.SetY2NDC(topLeft_y_top-topLeft_y_offset)
592	>	else:
593	>	NormLabel.SetY1NDC(topLeft_y_bottom)
594	>	NormLabel.SetY2NDC(topLeft_y_top)
595	>	NormLabel.Draw()
596	>
597	>	if drawHeaderLabel:
598	>	HeaderLabel.Draw()
599	>
600	>	YieldsLabel.Clear()
601	>	mcYield = Stack_list[i].GetStack().Last().Integral()
602	>	dataYield = Target.Integral()
603	>	if i == 0:
604	>	YieldsLabel.AddText ("Before Fit to Data")
605	>	if i == 1:
606	>	YieldsLabel.AddText ("After Fit to Data")
607	>	YieldsLabel.AddText ("data yield: " + '%.1f' % dataYield)
608	>	YieldsLabel.AddText ("MC yield: " + '%.1f' % mcYield)
609	>	if i == 1:
610	>	for j in range(0,len(FittingLegendEntries)):
611	>	RatiosLabel.AddText (FittingLegendEntries[j]+" ratio: " + '%.2f' % ratios[j])
612	>	YieldsLabel.Draw()
613	>	RatiosLabel.Draw()
614	>
615	>	# drawing the ratio or difference plot if requested
616	>	if (makeRatioPlots or makeDiffPlots):
617		Canvas.cd(2)
618	<	BgSum = Stack[i].GetStack().Last()
309	<	Comparison = Target.Clone()
310	<	Comparison.Add(BgSum,-1)
311	<	if not makeDiffPlots:
312	<	Comparison.Divide(BgSum)
313	<	Comparison.SetTitle("")
314	<	Comparison.GetXaxis().SetTitle(xAxisLabel)
618	>	BgSum = Stack_list[i].GetStack().Last()
619		if makeRatioPlots:
620	<	Comparison.GetYaxis().SetTitle("#frac{Data-MC}{MC}")
620	>	if arguments.ratioRelErrMax:
621	>	Comparison = ratioHistogram(Target,BgSum,arguments.ratioRelErrMax)
622	>	else:
623	>	Comparison = ratioHistogram(Target,BgSum)
624		elif makeDiffPlots:
625	+	Comparison = Target.Clone("diff")
626	+	Comparison.Add(BgSum,-1)
627	+	Comparison.SetTitle("")
628		Comparison.GetYaxis().SetTitle("Data-MC")
629	+	Comparison.GetXaxis().SetTitle(xAxisLabel)
630		Comparison.GetYaxis().CenterTitle()
631		Comparison.GetYaxis().SetTitleSize(0.1)
632	<	Comparison.GetYaxis().SetTitleOffset(0.35)
632	>	Comparison.GetYaxis().SetTitleOffset(0.5)
633		Comparison.GetXaxis().SetTitleSize(0.15)
634		Comparison.GetYaxis().SetLabelSize(0.1)
635		Comparison.GetXaxis().SetLabelSize(0.15)
636		if makeRatioPlots:
637	<	Comparison.GetYaxis().SetRangeUser(-1,1)
637	>	RatioYRange = 1.15
638	>	if arguments.ratioYRange:
639	>	RatioYRange = float(arguments.ratioYRange)
640	>	Comparison.GetYaxis().SetRangeUser(-1*RatioYRange, RatioYRange)
641		elif makeDiffPlots:
642		YMax = Comparison.GetMaximum()
643		YMin = Comparison.GetMinimum()
#	Line 336 | Line 650 | for histogram in input_histograms:
650		Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax)
651		else:
652		Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin)
653	<
653	>
654		Comparison.GetYaxis().SetNdivisions(205)
655		Comparison.Draw()
656
657	<	outputFile.cd ("OSUAnalysis/" + histogram['channel'])
657	>
658	>
659		if i == 0:
660	<	Canvas.Write (histogram['name'] + "_Before")
661	<	if arguments.plot_savePdf:
662	<	Canvas.Print (pdfFileName, "pdf")
660	>	Canvas.Write (distribution['name'] + "_Before")
661	>	if arguments.savePDFs:
662	>	pathToDirString = plainTextString(pathToDir)
663	>	Canvas.SaveAs(condor_dir+"/fitting_histogram_pdfs/"+pathToDirString+"/"+distribution['name']+"_Before.pdf")
664	>
665		if i == 1:
666	<	Canvas.Write (histogram['name'] + "_After")
667	<	if arguments.plot_savePdf:
668	<	if histogram == input_histograms[-1]:
669	<	Canvas.Print (pdfFileName + ")", "pdf")
670	<	else:
671	<	Canvas.Print (pdfFileName, "pdf")
666	>	Canvas.Write (distribution['name'] + "_After")
667	>	if arguments.savePDFs:
668	>	pathToDirString = plainTextString(pathToDir)
669	>	Canvas.SaveAs(condor_dir+"/fitting_histogram_pdfs/"+pathToDirString+"/"+distribution['name']+"_After.pdf")
670	>
671	>
672	>
673	>
674	>	##########################################################################################################################################
675	>	##########################################################################################################################################
676	>	##########################################################################################################################################
677	>
678	>
679	>	##########################################################################################################################################
680	>	##########################################################################################################################################
681	>	##########################################################################################################################################
682	>
683	>
684	>	condor_dir = set_condor_output_dir(arguments)
685	>
686	>
687	>	# make output file
688	>	outputFileName = "mc_fit_to_data.root"
689	>	if arguments.outputFileName:
690	>	outputFileName = arguments.outputFileName
691	>
692	>	outputFile = TFile(condor_dir + "/" + outputFileName, "RECREATE")
693	>
694	>
695	>	if arguments.savePDFs:
696	>	os.system("rm -rf %s/fitting_histograms_pdfs" % (condor_dir))
697	>	os.system("mkdir %s/fitting_histograms_pdfs" % (condor_dir))
698	>
699
700	<	outputFile.Close ()
700	>	#get root directory in the first layer, generally "OSUAnalysis"
701	>	for distribution in input_distributions:
702	>	MakeOneDHist("OSUAnalysis",distribution)
703
704	<	for dataset in datasets_needed:
359	<	dataset_file = "%s/%s.root_tmp" % (condor_dir,dataset)
360	<	os.remove(dataset_file)
704	>	outputFile.Close()

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines