--- UserCode/OSUT3Analysis/Configuration/scripts/fitMCToData.py 2013/04/07 00:57:26 1.1 +++ UserCode/OSUT3Analysis/Configuration/scripts/fitMCToData.py 2013/08/23 16:22:50 1.5 @@ -2,46 +2,64 @@ import sys import os import re -from optparse import OptionParser +import time +from math import * from array import * from decimal import * - +from optparse import OptionParser from OSUT3Analysis.Configuration.configurationOptions import * from OSUT3Analysis.Configuration.processingUtilities import * +from OSUT3Analysis.Configuration.formattingUtilities import * + -##set default plotting options -line_width = 2 -plotting_options = "" + +### parse the command-line options parser = OptionParser() parser = set_commandline_arguments(parser) -parser.add_option("--hist", action="store_true", dest="plot_hist", default=False, - help="plot as hollow histograms instead of error bar crosses") -parser.add_option("--line-width", dest="line_width", - help="set line width (default is 2)") -parser.add_option("--pdf", action="store_true", dest="plot_savePdf", default=False, - help="save plot as pdf in addition") +parser.add_option("-f", "--fancy", action="store_true", dest="makeFancy", default=False, + help="removes the title and replaces it with the official CMS plot heading") +parser.add_option("--ylog", action="store_true", dest="setLogY", default=False, + help="Set logarithmic scale on vertical axis on all plots") +parser.add_option("--ymin", dest="setYMin", + help="Minimum of y axis") +parser.add_option("--ymax", dest="setYMax", + help="Maximum of y axis") +parser.add_option("-E", "--ratioRelErrMax", dest="ratioRelErrMax", + help="maximum error used in rebinning the ratio histogram") (arguments, args) = parser.parse_args() + if arguments.localConfig: - sys.path.insert(0,os.getcwd()) + sys.path.append(os.getcwd()) exec("from " + arguments.localConfig.rstrip('.py') + " import *") -outputFileName = "mc_fit_to_data.root" -if arguments.outputFileName: - outputFileName = arguments.outputFileName -pdfFileName = outputFileName[:-5] + ".pdf" - -condor_dir = set_condor_output_dir(arguments) - +#### deal with conflicting arguments +if arguments.normalizeToData and arguments.normalizeToUnitArea: + print "Conflicting normalizations requsted, will normalize to unit area" + arguments.normalizeToData = False +if arguments.normalizeToData and arguments.noStack: + print "You have asked to scale non-stacked backgrounds to data. This is a very strange request. Will normalize to unit area instead" + arguments.normalizeToData = False + arguments.normalizeToUnitArea = True if arguments.makeRatioPlots and arguments.makeDiffPlots: print "You have requested both ratio and difference plots. Will make just ratio plots instead" arguments.makeRatioPlots = False +if arguments.makeRatioPlots and arguments.noStack: + 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." + arguments.makeRatioPlots = False +if arguments.makeDiffPlots and arguments.noStack: + 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." + arguments.makeDiffPlots = False + + +from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, THStack, TH1F, TCanvas, TString, TLegend, TLegendEntry, THStack, TIter, TKey, TPaveLabel, TPaveText, TF1, gPad + + +### setting ROOT options so our plots will look awesome and everyone will love us -from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, THStack, TH1F, TCanvas, TString, TLegend, TArrow, THStack, TIter, TKey, TPaveLabel, TPaveText, TF1, gPad -sys.argv = [] gROOT.SetBatch() gStyle.SetOptStat(0) gStyle.SetCanvasBorderMode(0) @@ -49,281 +67,577 @@ gStyle.SetPadBorderMode(0) gStyle.SetPadColor(0) gStyle.SetCanvasColor(0) gStyle.SetTextFont(42) +gStyle.SetCanvasDefH(600) +gStyle.SetCanvasDefW(600) +gStyle.SetCanvasDefX(0) +gStyle.SetCanvasDefY(0) +gStyle.SetPadTopMargin(0.07) +gStyle.SetPadBottomMargin(0.13) +gStyle.SetPadLeftMargin(0.15) +gStyle.SetPadRightMargin(0.05) +gStyle.SetTitleColor(1, "XYZ") +gStyle.SetTitleFont(42, "XYZ") +gStyle.SetTitleSize(0.04, "XYZ") +gStyle.SetTitleXOffset(1.1) +gStyle.SetTitleYOffset(2) +gStyle.SetTextAlign(12) +gStyle.SetLabelColor(1, "XYZ") +gStyle.SetLabelFont(42, "XYZ") +gStyle.SetLabelOffset(0.007, "XYZ") +gStyle.SetLabelSize(0.04, "XYZ") +gStyle.SetAxisColor(1, "XYZ") +gStyle.SetStripDecimals(True) +gStyle.SetTickLength(0.03, "XYZ") +gStyle.SetNdivisions(510, "XYZ") +gStyle.SetPadTickX(1) +gStyle.SetPadTickY(1) gROOT.ForceStyle() -outputFile = TFile(outputFileName, "RECREATE") -datasets_needed = [] -for histogram in input_histograms: - for dataset in histogram['datasets']: - if dataset not in datasets_needed: - datasets_needed.append(dataset) - if histogram['target_dataset'] not in datasets_needed: - datasets_needed.append(histogram['target_dataset']) - -weight = intLumi / 10000.0 -for dataset in datasets_needed: - dataset_file = "%s/%s.root" % (condor_dir,dataset) - fin = TFile (dataset_file) - flags = fin.Get ("flags") - noWeights = flags and flags.GetBinContent (1) - fin.Close () +#set the text for the luminosity label +if(intLumi < 1000.): + LumiInPb = intLumi + LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}" + LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInPb) + " pb^{-1}" +else: + LumiInFb = intLumi/1000. + LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInFb) + " fb^{-1}" + +#bestest place for lumi. label, in top left corner +topLeft_x_left = 0.1375839 +topLeft_x_right = 0.4580537 +topLeft_y_bottom = 0.8479021 +topLeft_y_top = 0.9475524 +topLeft_y_offset = 0.035 + +#set the text for the fancy heading +HeaderText = "CMS Preliminary: " + LumiText + " at #sqrt{s} = 8 TeV" + +#position for header +header_x_left = 0.2181208 +header_x_right = 0.9562937 +header_y_bottom = 0.9479866 +header_y_top = 0.9947552 - if types[dataset] != "data" and not noWeights: - os.system("mergeTFileServiceHistograms -i %s -o %s -w %g" % (dataset_file, dataset_file + "_tmp", weight)) - else: - os.system("mergeTFileServiceHistograms -i %s -o %s -w %g" % (dataset_file, dataset_file + "_tmp", 1.0)) -for histogram in input_histograms: - HistogramsToFit = [] - HistogramDatasets = [] - TargetDataset = histogram['target_dataset'] - Stack = [] - Stack.append (THStack("stack_before",histogram['name'])) - Stack.append (THStack("stack_after",histogram['name'])) +########################################################################################################################################## +########################################################################################################################################## +########################################################################################################################################## + +# some fancy-ass code from Andrzej Zuranski to merge bins in the ratio plot until the error goes below some threshold +def ratioHistogram( dataHist, mcHist, relErrMax=0.10): + + if not dataHist: + print "Error: trying to run ratioHistogram but dataHist is invalid" + return + + if not mcHist: + print "Error: trying to run ratioHistogram but mcHist is invalid" + return + + def groupR(group): + Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]] + return (Data-MC)/MC if MC else 0 + + def groupErr(group): + Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]] + dataErr2,mcErr2 = [sum(hist.GetBinError(i)**2 for i in group) for hist in [dataHist,mcHist]] + if Data > 0 and MC > 0 and Data != MC: + return abs(math.sqrt( (dataErr2+mcErr2)/(Data-MC)**2 + mcErr2/MC**2 ) * (Data-MC)/MC) + else: + return 0 - if(intLumi < 1000.): - LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}" + def regroup(groups): + err,iG = max( (groupErr(g),groups.index(g)) for g in groups ) + if err < relErrMax or len(groups)<3 : return groups + iH = max( [iG-1,iG+1], key = lambda i: groupErr(groups[i]) if 0<=i 0): + Target.Scale(1./Target.Integral()) + + ### formatting target histogram and adding to legend + legendIndex = 0 + Legend.AddEntry(Target,labels[TargetDataset],"LEP") + legendIndex = legendIndex+1 + if not outputFile.Get ("OSUAnalysis"): outputFile.mkdir ("OSUAnalysis") - if not outputFile.Get ("OSUAnalysis/" + histogram['channel']): - outputFile.Get ("OSUAnalysis").mkdir (histogram['channel']) - - for dataset in histogram['datasets']: - fileName = condor_dir + "/" + dataset + ".root_tmp" - if not os.path.exists(fileName): - continue - inputFile = TFile(fileName) - if inputFile.IsZombie() or not inputFile.GetNkeys(): - continue - - Histogram = inputFile.Get("OSUAnalysis/"+histogram['channel']+"/"+histogram['name']).Clone() + if not outputFile.Get ("OSUAnalysis/" + distribution['channel']): + outputFile.Get ("OSUAnalysis").mkdir (distribution['channel']) + + for sample in distribution['datasets']: # loop over different samples requested to be fit + + dataset_file = "%s/%s.root" % (condor_dir,sample) + inputFile = TFile(dataset_file) + HistogramObj = inputFile.Get(pathToDir+"/"+distribution['channel']+"/"+distribution['name']) + if not HistogramObj: + print "WARNING: Could not find histogram " + pathToDir + "/" + distribution['name'] + " in file " + dataset_file + ". Will skip it and continue." + continue + Histogram = HistogramObj.Clone() Histogram.SetDirectory(0) inputFile.Close() - if arguments.rebinFactor: RebinFactor = int(arguments.rebinFactor) - if Histogram.GetNbinsX() >= RebinFactor*10: + #don't rebin histograms which will have less than 5 bins or any gen-matching histograms + if Histogram.GetNbinsX() >= RebinFactor*5 and Histogram.GetName().find("GenMatch") is -1: Histogram.Rebin(RebinFactor) - numBgMCSamples += 1 + xAxisLabel = Histogram.GetXaxis().GetTitle() + unitBeginIndex = xAxisLabel.find("[") + unitEndIndex = xAxisLabel.find("]") + + if unitBeginIndex is not -1 and unitEndIndex is not -1: #x axis has a unit + yAxisLabel = "Entries / " + str(Histogram.GetXaxis().GetBinWidth(1)) + " " + xAxisLabel[unitBeginIndex+1:unitEndIndex] + else: + yAxisLabel = "Entries per bin (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " width)" + + if not arguments.makeFancy: + histoTitle = Histogram.GetTitle() + else: + histoTitle = "" + - if(arguments.noStack): + legLabel = labels[sample] + if (arguments.printYields): + yieldHist = Histogram.Integral() + legLabel = legLabel + " (%.1f)" % yieldHist + FittingLegendEntries.append(legLabel) + + if( types[sample] == "bgMC"): + + numFittingSamples += 1 + fittingIntegral += Histogram.Integral() + + Histogram.SetLineStyle(1) + if(arguments.noStack): + Histogram.SetFillStyle(0) + Histogram.SetLineColor(colors[sample]) + Histogram.SetLineWidth(2) + else: + Histogram.SetFillStyle(1001) + Histogram.SetFillColor(colors[sample]) + Histogram.SetLineColor(1) + Histogram.SetLineWidth(1) + + elif( types[sample] == "signalMC"): + + numFittingSamples += 1 + Histogram.SetFillStyle(0) - Histogram.SetLineColor(colors[dataset]) + Histogram.SetLineColor(colors[sample]) + Histogram.SetLineStyle(1) Histogram.SetLineWidth(2) - BgMCLegend.AddEntry(Histogram,labels[dataset],"L") - else: - Histogram.SetFillStyle(1001) - Histogram.SetFillColor(colors[dataset]) - Histogram.SetLineColor(1) - Histogram.SetLineWidth(1) - BgMCLegend.AddEntry(Histogram,labels[dataset],"F") - Histogram.SetLineStyle(1) - + if(arguments.normalizeToUnitArea and Histogram.Integral() > 0): + Histogram.Scale(1./Histogram.Integral()) + HistogramsToFit.append(Histogram) - HistogramDatasets.append(dataset) + FittingHistogramDatasets.append(sample) + + #scaling histograms as per user's specifications + if targetIntegral > 0 and fittingIntegral > 0: + scaleFactor = targetIntegral/fittingIntegral + for fittingHist in HistogramsToFit: + if arguments.normalizeToData: + fittingHist.Scale(scaleFactor) + + if arguments.normalizeToUnitArea and not arguments.noStack and fittingIntegral > 0: + fittingHist.Scale(1./fittingIntegral) + elif arguments.normalizeToUnitArea and arguments.noStack and fittingHist.Integral() > 0: + fittingHist.Scale(1./fittingHist.Integral()) + def fitf (x, par): xBin = HistogramsToFit[0].FindBin (x[0]) value = 0.0 - + sumOfWeights = 0.0 + for i in range (0, len (HistogramsToFit)): - value += par[i] * HistogramsToFit[i].GetBinContent (xBin) - + weight = 1.0 / (HistogramsToFit[i].GetBinError (xBin) * HistogramsToFit[i].GetBinError (xBin)) + sumOfWeights += weight + value += weight * par[i] * HistogramsToFit[i].GetBinContent (xBin) + value /= sumOfWeights + return value lowerLimit = Target.GetBinLowEdge (1) upperLimit = Target.GetBinLowEdge (Target.GetNbinsX ()) + Target.GetBinWidth (Target.GetNbinsX ()) - if 'lowerLimit' in histogram: - lowerLimit = histogram['lowerLimit'] - if 'upperLimit' in histogram: - upperLimit = histogram['upperLimit'] + if 'lowerLimit' in distribution: + lowerLimit = distribution['lowerLimit'] + if 'upperLimit' in distribution: + upperLimit = distribution['upperLimit'] func = TF1 ("fit", fitf, lowerLimit, upperLimit, len (HistogramsToFit)) for i in range (0, len (HistogramsToFit)): func.SetParameter (i, 1.0) - func.SetParLimits (i, 0.0, 10.0) - func.SetParName (i, labels[HistogramDatasets[i]]) + func.SetParName (i, labels[FittingHistogramDatasets[i]]) - for i in range (0, histogram['iterations'] - 1): + for i in range (0, distribution['iterations'] - 1): print "Iteration " + str (i + 1) + "..." Target.Fit ("fit", "QEMR0") Target.Fit ("fit", "VEMR0") - func.SetLineWidth (line_width) - func.SetLineColor (632) - - Target.SetMarkerStyle (20) - Target.SetMarkerColor (1) - Target.SetLineColor (1) - Canvas = TCanvas(histogram['name']) + finalMax = 0 + if not arguments.noStack: + for fittingHist in HistogramsToFit: + finalMax += fittingHist.GetMaximum() + else: + for fittingHist in HistogramsToFit: + if(fittingHist.GetMaximum() > finalMax): + finalMax = fittingHist.GetMaximum() + if(Target.GetMaximum() > finalMax): + finalMax = Target.GetMaximum() + + Target.SetMaximum(1.1*finalMax) + Target.SetMinimum(0.0001) + + Canvas = TCanvas(distribution['name'] + "_FitFunction") Canvas.cd (1) Target.Draw () - func.Draw (plotting_options + "same") - - outputFile.cd ("OSUAnalysis/" + histogram['channel']) + func.Draw ("same") + + outputFile.cd ("OSUAnalysis/" + distribution['channel']) Canvas.Write () - if arguments.plot_savePdf: - if histogram == input_histograms[0]: - Canvas.Print (pdfFileName + "(", "pdf") - else: - Canvas.Print (pdfFileName, "pdf") + if arguments.savePDFs: + if histogram == input_histograms[0]: + Canvas.Print (pdfFileName + "(", "pdf") + else: + Canvas.Print (pdfFileName, "pdf") Target.SetStats (0) - for i in range (0, 2): + + + + ### formatting bgMC histograms and adding to legend + legendIndex = numFittingSamples-1 + for Histogram in reversed(HistogramsToFit): + if(arguments.noStack): + Legend.AddEntry(Histogram,FittingLegendEntries[legendIndex],"L") + else: + Legend.AddEntry(Histogram,FittingLegendEntries[legendIndex],"F") + legendIndex = legendIndex-1 + + + ### Drawing histograms to canvas + + makeRatioPlots = arguments.makeRatioPlots + makeDiffPlots = arguments.makeDiffPlots + + yAxisMin = 0.0001 + if arguments.setYMin: + yAxisMin = float(arguments.setYMin) + + + ### Draw everything to the canvases !!!! + + for i in range (0, 2): # 0 => before, 1 => after if i == 1: + ratios = [] for j in range (0, len (HistogramsToFit)): HistogramsToFit[j].Scale (func.GetParameter (j)) + ratios.append(func.GetParameter (j)) - for bgMCHist in HistogramsToFit: + for fittingHist in HistogramsToFit: if not arguments.noStack: - Stack[i].Add(bgMCHist) + Stack_list[i].Add(fittingHist) - finalMax = 0 - if not arguments.noStack: - finalMax = Stack[i].GetMaximum() - else: - for bgMCHist in HistogramsToFit: - if(bgMCHist.GetMaximum() > finalMax): - finalMax = bgMCHist.GetMaximum() - if(Target.GetMaximum() > finalMax): - finalMax = Target.GetMaximum() - makeRatioPlots = arguments.makeRatioPlots - makeDiffPlots = arguments.makeDiffPlots + #creating the histogram to represent the statistical errors on the stack + if not arguments.noStack: + ErrorHisto = HistogramsToFit[0].Clone("errors") + ErrorHisto.SetFillStyle(3001) + ErrorHisto.SetFillColor(13) + ErrorHisto.SetLineWidth(0) + if i == 1: + Legend.AddEntry(ErrorHisto,"Stat. Errors","F") + for Histogram in HistogramsToFit: + if Histogram is not HistogramsToFit[0]: + ErrorHisto.Add(Histogram) + if i == 0: - Canvas = TCanvas(histogram['name'] + "_Before") + Canvas = TCanvas(distribution['name'] + "_Before") if i == 1: - Canvas = TCanvas(histogram['name'] + "_After") + Canvas = TCanvas(distribution['name'] + "_After") + + if makeRatioPlots or makeDiffPlots: Canvas.SetFillStyle(0) Canvas.Divide(1,2) Canvas.cd(1) - gPad.SetPad(0.01,0.25,0.99,0.99) - gPad.SetMargin(0.1,0.05,0.02,0.07) + gPad.SetPad(0,0.25,1,1) + gPad.SetMargin(0.15,0.05,0.01,0.07) gPad.SetFillStyle(0) gPad.Update() gPad.Draw() + if arguments.setLogY: + gPad.SetLogy() Canvas.cd(2) - gPad.SetPad(0.01,0.01,0.99,0.25) - #format: gPad.SetMargin(l,r,b,t) - gPad.SetMargin(0.1,0.05,0.4,0.02) + gPad.SetPad(0,0,1,0.25) + # format: gPad.SetMargin(l,r,b,t) + gPad.SetMargin(0.15,0.05,0.4,0.01) gPad.SetFillStyle(0) gPad.SetGridy(1) gPad.Update() gPad.Draw() - + Canvas.cd(1) - if not arguments.noStack: - Stack[i].SetTitle(histoTitle) - Stack[i].Draw("HIST") - Stack[i].GetXaxis().SetTitle(xAxisLabel) - Stack[i].SetMaximum(1.1*finalMax) - Stack[i].SetMinimum(0.0001) - if makeRatioPlots or makeDiffPlots: - Stack[i].GetHistogram().GetXaxis().SetLabelSize(0) + ### finding the maximum value of anything going on the canvas, so we know how to set the y-axis + finalMax = 0 + if numFittingSamples is not 0 and not arguments.noStack: + finalMax = ErrorHisto.GetMaximum() + ErrorHisto.GetBinError(ErrorHisto.GetMaximumBin()) else: + for bgMCHist in HistogramsToFit: + if(bgMCHist.GetMaximum() > finalMax): + finalMax = bgMCHist.GetMaximum() + if(Target.GetMaximum() > finalMax): + finalMax = Target.GetMaximum() + Target.GetBinError(Target.GetMaximumBin()) + finalMax = 1.15*finalMax + if arguments.setYMax: + finalMax = float(arguments.setYMax) + + + if not arguments.noStack: # draw stacked background samples + Stack_list[i].SetTitle(histoTitle) + Stack_list[i].Draw("HIST") + Stack_list[i].GetXaxis().SetTitle(xAxisLabel) + Stack_list[i].GetYaxis().SetTitle(yAxisLabel) + Stack_list[i].SetMaximum(finalMax) + Stack_list[i].SetMinimum(yAxisMin) + if makeRatioPlots or makeDiffPlots: + Stack_list[i].GetHistogram().GetXaxis().SetLabelSize(0) + #draw shaded error bands + ErrorHisto.Draw("A E2 SAME") + + else: #draw the unstacked backgrounds HistogramsToFit[0].SetTitle(histoTitle) HistogramsToFit[0].Draw("HIST") HistogramsToFit[0].GetXaxis().SetTitle(xAxisLabel) - HistogramsToFit[0].SetMaximum(1.1*finalMax) - HistogramsToFit[0].SetMinimum(0.0001) + HistogramsToFit[0].GetYaxis().SetTitle(yAxisLabel) + HistogramsToFit[0].SetMaximum(finalMax) + HistogramsToFit[0].SetMinimum(yAxisMin) for bgMCHist in HistogramsToFit: - bgMCHist.Draw("HIST SAME") + bgMCHist.Draw("A HIST SAME") - dataYield = Target.Integral (1, Target.GetNbinsX ()) - mcYield = 0.0 - for bgMCHist in HistogramsToFit: - mcYield += bgMCHist.Integral (1, bgMCHist.GetNbinsX ()) - Label.Clear () - if i == 0: - Label.AddText ("Before Fit to Data") - if i == 1: - Label.AddText ("After Fit to Data") - Label.AddText ("data yield: " + str (dataYield)) - Label.AddText ("MC yield: " + str (mcYield)) - - Target.Draw("E SAME") - BgMCLegend.Draw() - LumiLabel.Draw() - Label.Draw() + Target.Draw("A E X0 SAME") - if makeRatioPlots or makeDiffPlots: + + + #legend coordinates, empirically determined :-) + x_left = 0.6761745 + x_right = 0.9328859 + x_width = x_right - x_left + y_max = 0.9 + entry_height = 0.05 + + if(numFittingSamples is not 0): #then draw the data & bgMC legend + + numExtraEntries = 2 # count the target and (lack of) title + Legend.SetX1NDC(x_left) + numExtraEntries = numExtraEntries + 1 # count the stat. errors entry + + Legend.SetY1NDC(y_max-entry_height*(numExtraEntries+numFittingSamples)) + Legend.SetX2NDC(x_right) + Legend.SetY2NDC(y_max) + Legend.Draw() + + RatiosLabel.SetX1NDC(x_left - 0.1) + RatiosLabel.SetX2NDC(x_right) + RatiosLabel.SetY2NDC(Legend.GetY1NDC() - 0.1) + RatiosLabel.SetY1NDC(RatiosLabel.GetY2NDC() - entry_height*(numFittingSamples)) + + # Deciding which text labels to draw and drawing them + drawLumiLabel = False + drawNormLabel = False + offsetNormLabel = False + drawHeaderLabel = False + + if not arguments.normalizeToUnitArea: #don't draw the lumi label if there's no data and it's scaled to unit area + drawLumiLabel = True + # move the normalization label down before drawing if we drew the lumi. label + offsetNormLabel = True + if arguments.normalizeToUnitArea or arguments.normalizeToData: + drawNormLabel = True + if arguments.makeFancy: + drawHeaderLabel = True + drawLumiLabel = False + + # now that flags are set, draw the appropriate labels + + if drawLumiLabel: + LumiLabel.Draw() + + if drawNormLabel: + if offsetNormLabel: + NormLabel.SetY1NDC(topLeft_y_bottom-topLeft_y_offset) + NormLabel.SetY2NDC(topLeft_y_top-topLeft_y_offset) + else: + NormLabel.SetY1NDC(topLeft_y_bottom) + NormLabel.SetY2NDC(topLeft_y_top) + NormLabel.Draw() + + if drawHeaderLabel: + HeaderLabel.Draw() + + YieldsLabel.Clear() + mcYield = Stack_list[i].GetStack().Last().Integral() + dataYield = Target.Integral() + if i == 0: + YieldsLabel.AddText ("Before Fit to Data") + if i == 1: + YieldsLabel.AddText ("After Fit to Data") + YieldsLabel.AddText ("data yield: " + '%.1f' % dataYield) + YieldsLabel.AddText ("MC yield: " + '%.1f' % mcYield) + if i == 1: + for j in range(0,len(FittingLegendEntries)): + RatiosLabel.AddText (FittingLegendEntries[j]+" ratio: " + '%.2f' % ratios[j]) + YieldsLabel.Draw() + RatiosLabel.Draw() + + # drawing the ratio or difference plot if requested + if (makeRatioPlots or makeDiffPlots): Canvas.cd(2) - BgSum = Stack[i].GetStack().Last() - Comparison = Target.Clone() - Comparison.Add(BgSum,-1) - if not makeDiffPlots: - Comparison.Divide(BgSum) - Comparison.SetTitle("") - Comparison.GetXaxis().SetTitle(xAxisLabel) + BgSum = Stack_list[i].GetStack().Last() if makeRatioPlots: - Comparison.GetYaxis().SetTitle("#frac{Data-MC}{MC}") + if arguments.ratioRelErrMax: + Comparison = ratioHistogram(Target,BgSum,arguments.ratioRelErrMax) + else: + Comparison = ratioHistogram(Target,BgSum) elif makeDiffPlots: + Comparison = Target.Clone("diff") + Comparison.Add(BgSum,-1) + Comparison.SetTitle("") Comparison.GetYaxis().SetTitle("Data-MC") + Comparison.GetXaxis().SetTitle(xAxisLabel) Comparison.GetYaxis().CenterTitle() Comparison.GetYaxis().SetTitleSize(0.1) - Comparison.GetYaxis().SetTitleOffset(0.35) + Comparison.GetYaxis().SetTitleOffset(0.5) Comparison.GetXaxis().SetTitleSize(0.15) Comparison.GetYaxis().SetLabelSize(0.1) Comparison.GetXaxis().SetLabelSize(0.15) if makeRatioPlots: - Comparison.GetYaxis().SetRangeUser(-1,1) + RatioYRange = 1.15 + if arguments.ratioYRange: + RatioYRange = float(arguments.ratioYRange) + Comparison.GetYaxis().SetRangeUser(-1*RatioYRange, RatioYRange) elif makeDiffPlots: YMax = Comparison.GetMaximum() YMin = Comparison.GetMinimum() @@ -336,25 +650,55 @@ for histogram in input_histograms: Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax) else: Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin) - + Comparison.GetYaxis().SetNdivisions(205) Comparison.Draw() - outputFile.cd ("OSUAnalysis/" + histogram['channel']) + + if i == 0: - Canvas.Write (histogram['name'] + "_Before") - if arguments.plot_savePdf: - Canvas.Print (pdfFileName, "pdf") + Canvas.Write (distribution['name'] + "_Before") + if arguments.savePDFs: + pathToDirString = plainTextString(pathToDir) + Canvas.SaveAs(condor_dir+"/fitting_histogram_pdfs/"+pathToDirString+"/"+distribution['name']+"_Before.pdf") + if i == 1: - Canvas.Write (histogram['name'] + "_After") - if arguments.plot_savePdf: - if histogram == input_histograms[-1]: - Canvas.Print (pdfFileName + ")", "pdf") - else: - Canvas.Print (pdfFileName, "pdf") + Canvas.Write (distribution['name'] + "_After") + if arguments.savePDFs: + pathToDirString = plainTextString(pathToDir) + Canvas.SaveAs(condor_dir+"/fitting_histogram_pdfs/"+pathToDirString+"/"+distribution['name']+"_After.pdf") + + + + +########################################################################################################################################## +########################################################################################################################################## +########################################################################################################################################## + + +########################################################################################################################################## +########################################################################################################################################## +########################################################################################################################################## + + +condor_dir = set_condor_output_dir(arguments) + + +# make output file +outputFileName = "mc_fit_to_data.root" +if arguments.outputFileName: + outputFileName = arguments.outputFileName + +outputFile = TFile(condor_dir + "/" + outputFileName, "RECREATE") + + +if arguments.savePDFs: + os.system("rm -rf %s/fitting_histograms_pdfs" % (condor_dir)) + os.system("mkdir %s/fitting_histograms_pdfs" % (condor_dir)) + -outputFile.Close () +#get root directory in the first layer, generally "OSUAnalysis" +for distribution in input_distributions: + MakeOneDHist("OSUAnalysis",distribution) -for dataset in datasets_needed: - dataset_file = "%s/%s.root_tmp" % (condor_dir,dataset) - os.remove(dataset_file) +outputFile.Close()