########################################################################
# -*- python -*-
# File Name .....: plotQCDMeasurementEffAndPur.py
# Original Author: Alexandros Attikis
# e-mail ........: attikis@cern.ch
# Institute .....: UCY
# Created .......: Friday, 16 Dec 2010 
# Last Edited ...: Friday, 7 Jan 2011
# Description ...: Python script that takes as input a dat or text file
#                  generated by the hplusPrintCounters.py script and
#                  uses as input the number of events passing specific
#                  cuts. It returns the efficiency and purity of the QCD
#                  samples, based on Bayesian methods.
#                  The efficiency is defined as: e = N/No
#                  where:
#                  N  = Events passing all the cuts
#                  No = Events passing the jet selection
#                  The purity is also calculated for each QCD pT-bin 
#                  by by the use of pseudo-experiments performed using
#                  the binomial distribution. Two parameters are required
#                  for this:
#                  i. the probability to pass the cut flow after jetSelection (fakeMET, Nbjets, e/mu veto, MET)
#                     i.e. the effiiciency of specific sample for the cut flow
#                  ii. the number of starting events (i.e. number of events after the Trigger, tauAntiIsolation and jetSelection)
########################################################################
from __future__ import division
#!/usr/bin env python
### to execute: |myLinux> python plotQCDMeasurementEffAndPur.py 
### OR
### to execute: |myLinux> python plotQCDMeasurementEffAndPur.py testFile.dat
### where testFile.dat is an example dat file provided for testing purposes

### System imports
import os, sys
from array import array
from math import fabs
### ROOT imports
from ROOT import *
### HiggsAnalysis imports
from HiggsAnalysis.HeavyChHiggsToTauNu.tools.dataset import *
from HiggsAnalysis.HeavyChHiggsToTauNu.tools.histograms import *
from HiggsAnalysis.HeavyChHiggsToTauNu.tools.tdrstyle import *
from HiggsAnalysis.HeavyChHiggsToTauNu.tools.histograms import *
import HiggsAnalysis.HeavyChHiggsToTauNu.tools.styles as styles
from FindFirstBinAbove import * 
from bayes import * 
from myArrays import *
#from HiggsAnalysis.HeavyChHiggsToTauNu.tools.styles import *
#from AtlasStyle import *


### Provide the name of the file to be opened 
if len(sys.argv) > 1:
    file = sys.argv[1]
### Otherwise use as default the name counters.dat
else:
    file = 'counters.dat'

### If the file already exists remove it since the script will automatically generate one.
if os.path.exists(file):
    os.system('rm %s'%(file))
os.system('hplusPrintCounters3.py --mode xsect --mainCounterOnly > %s ; echo "EOF" >> %s'%(file,file))

######################################################################################
### Options (script)
verbosity = False ### False to suppress most "print" commands

### Notes
# For some QCD samples the purity appears to be quite hight (almost 100%) which doesn't make sense! The reason is the very high Normalisation
# factors that these samples have which are much bigger than the rest. So, because the statistics are so low, if you try to set an upper limit
# on the purity say, that event when weighted by the Normalization factor could be either 1000 or 0, anything in between. This overwhelms the
# rest of the events considered (EWK = TTbar + WJets) and that's why you have so high purity. This is okay, BUT, you have to make sure that
# PURITY ERROR bar is also huge (e.g. spans from max to 0!)

######################################################################################
### Open file for processing
file = open('%s'%(file), 'r')

### Define lists to be used
DataSetName = []
CutNamesAndValues = []
DataSetAndXSectAndNormFact = []

######################################################################################
### Loop over lines and save the cut-flow info (Cross Section for MC, Events for data)
lim = True
while file:
    line = file.readline()
    # split the line and puts elements in list s as strings
    s = line.split()
    if 'EOF' in s:
        break
    n = len(s)
    # 10 is the number of columns
    ### Investigate lists which have three elements (i.e. 3 is the number of columns of each list)
    if n == 3:
            ### Save the names of the datasets, cross section and normalisation factor
            DataSetAndXSectAndNormFact.append([s[0], s[1], s[2]])
    if n == 10:
        if lim == True:
            for i in range(n):
                DataSetName.append(s[i])
            lim = False
        else:
            #print "CutNamesAndValues\tCut\tValue" 
            for i in range(1,n):
                CutNamesAndValues.append([DataSetName[i],s[0],float(s[i])])
                #print "%s\t%s\t%.0f"%(DataSetName[i],s[0],float(s[i]))

if verbosity == True:
    ### Print the DataSet Name, Cross Section and MC Normalisation Factor (in this order)
    print
    print "********************************************************************"
    print "* Printing DataSet Name, Cross Section and MC Normalisation Factor *" 
    print "********************************************************************"
    print "DataSetAndXSectAndNormFact = ", DataSetAndXSectAndNormFact
    print "For Wjets  we have ..."
    print "DataSetAndXSectAndNormFact[0][0] = ", DataSetAndXSectAndNormFact[0][0]
    print "DataSetAndXSectAndNormFact[0][1] = ", DataSetAndXSectAndNormFact[0][1]
    print "DataSetAndXSectAndNormFact[0][2] = ", DataSetAndXSectAndNormFact[0][2]
    print "For QCD_Pt30to50 we have ..."
    print "DataSetAndXSectAndNormFact[1][0] = ", DataSetAndXSectAndNormFact[1][0]
    print "DataSetAndXSectAndNormFact[1][1] = ", DataSetAndXSectAndNormFact[1][1]
    print "DataSetAndXSectAndNormFact[1][2] = ", DataSetAndXSectAndNormFact[1][2]

### Declare the Normalisation Factors for each DataSet
WJets_NormFactor        = float(DataSetAndXSectAndNormFact[0][2])
qcd30to50_NormFactor    = float(DataSetAndXSectAndNormFact[1][2])
qcd50to80_NormFactor    = float(DataSetAndXSectAndNormFact[2][2])
qcd80to120_NormFactor   = float(DataSetAndXSectAndNormFact[3][2])
qcd120to170_NormFactor  = float(DataSetAndXSectAndNormFact[4][2])
qcd170to230_NormFactor  = float(DataSetAndXSectAndNormFact[5][2])
qcd230to300_NormFactor  = float(DataSetAndXSectAndNormFact[6][2])
TTbarJets_NormFactor    = float(DataSetAndXSectAndNormFact[7][2])
TTbar_NormFactor        = float(DataSetAndXSectAndNormFact[8][2])
### Define lists that will hold the information (Cut Name and Xsection/Events surviving Cut)
WJets        = []
qcd30to50    = []
qcd50to80    = []
qcd80to120   = []
qcd120to170  = []
qcd170to230  = []
qcd230to300  = []
TTbarJets    = []
TTbar        = []
######################################################################################
### Loop over all cut names and create new lists according to the DataSet Name
for item in CutNamesAndValues:
    name = item[0]
    tmp = [item[1],item[2]]
    if name == 'WJets':
        WJets.append(tmp)
    if name == 'QCD_Pt30to50':
        qcd30to50.append(tmp)
    if name == 'QCD_Pt50to80':
        qcd50to80.append(tmp)
    if name == 'QCD_Pt80to120':
        qcd80to120.append(tmp)
    if name == 'QCD_Pt120to170':
        qcd120to170.append(tmp)
    if name == 'QCD_Pt170to230':
        qcd170to230.append(tmp)
    if name == 'QCD_Pt230to300':
        qcd230to300.append(tmp)
    if name == 'TTbarJets':
        TTbarJets.append(tmp)
    if name == 'TTbar':
        TTbar.append(tmp)
######################################################################################
### Print lists 
if verbosity == True:
    print
    print "*******************************************"
    print "* Printing cut-flow for each DataSet Name *"
    print "*******************************************"
    print 'qcd30to50', qcd30to50
    print
    print 'qcd50to80', qcd50to80
    print
    print 'qcd80to120', qcd80to120
    print
    print 'qcd120to170', qcd120to170
    print
    print 'qcd170to230', qcd170to230
    print
    print 'qcd230to300', qcd230to300
    print
    print 'TTbarJets', TTbarJets
    print
    print 'TTbar', TTbar
    print
    print 'W+Jets', WJets
    print "***************************************************************************************************************"
######################################################################################
### Calculate cross sections before and after cuts
# ElectroWeak (TTbar + WJets)
XSectTotalTTbar  = TTbar[4][1]
XSectTotalWJets  = WJets[4][1]
XSectPassedTTbar = TTbar[-1][1]
XSectPassedWJets = WJets[-1][1]
XSectTotalEWK    = XSectTotalTTbar + XSectTotalWJets
XSectPassedEWK   = XSectPassedTTbar + XSectPassedWJets
### QCD 
XSectTotalQCD  = qcd30to50[4][1]+qcd50to80[4][1]+qcd80to120[4][1]+qcd120to170[4][1]+qcd170to230[4][1]+qcd230to300[4][1]
XSectPassedQCD = qcd30to50[-1][1]+qcd50to80[-1][1]+qcd80to120[-1][1]+qcd120to170[-1][1]+qcd170to230[-1][1]+qcd230to300[-1][1]
XSectPassedQCD_UpLimit = qcd30to50_NormFactor+qcd50to80_NormFactor+qcd80to120[-1][1]+qcd120to170[-1][1]+qcd170to230[-1][1]+qcd230to300[-1][1]
### Save Events before and after cuts
# ElectroWeak (TTbar + WJets)
EvtsTotalTTbar  = XSectTotalTTbar/TTbar_NormFactor
EvtsTotalWJets  = XSectTotalWJets/WJets_NormFactor
EvtsTotalEWK    = EvtsTotalTTbar + EvtsTotalWJets
#
EvtsPassedTTbar = XSectPassedTTbar/TTbar_NormFactor
EvtsPassedWJets = XSectPassedWJets/WJets_NormFactor
EvtsPassedEWK   = EvtsPassedTTbar + EvtsPassedWJets
### QCD 
EvtsTotalQCD  = qcd30to50[4][1]/qcd30to50_NormFactor + qcd50to80[4][1]/qcd50to80_NormFactor + qcd80to120[4][1]/qcd80to120_NormFactor + qcd120to170[4][1]/qcd120to170_NormFactor + qcd170to230[4][1]/qcd170to230_NormFactor + qcd230to300[4][1]/qcd230to300_NormFactor
EvtsPassedQCD = qcd30to50[-1][1]/qcd30to50_NormFactor + qcd50to80[-1][1]/qcd50to80_NormFactor + qcd80to120[-1][1]/qcd80to120_NormFactor + qcd120to170[-1][1]/qcd120to170_NormFactor + qcd170to230[-1][1]/qcd170to230_NormFactor + qcd230to300[-1][1]/qcd230to300_NormFactor

### Calculate Efficiencies
# QCD combined: Using Bayes
qcdObj_eff = bayesDivide(EvtsPassedQCD, EvtsTotalQCD)
QCD_efficiency_bayes = qcdObj_eff.GetY()[0]
QCD_eff_errUp = qcdObj_eff.GetEYhigh()[0]
QCD_eff_errDown = qcdObj_eff.GetEYlow()[0]
# QCD combined: Using Standard Efficiency Definition
QCD_efficiency_estimate = XSectPassedQCD/XSectTotalQCD

# TTbar: Using Bayes
ttbarObj_eff = bayesDivide(EvtsPassedTTbar, EvtsTotalTTbar)
TTbar_efficiency_bayes = ttbarObj_eff.GetY()[0]
TTbar_eff_errUp = ttbarObj_eff.GetEYhigh()[0]
TTbar_eff_errDown = ttbarObj_eff.GetEYlow()[0]
# TTbar: Using Standard Efficiency Definition
TTbar_efficiency_estimate = XSectPassedTTbar/XSectTotalTTbar

# WJets: Using Bayes
wjetsObj_eff = bayesDivide(EvtsPassedWJets, EvtsTotalWJets)
WJets_efficiency_bayes = wjetsObj_eff.GetY()[0]
WJets_eff_errUp = wjetsObj_eff.GetEYhigh()[0]
WJets_eff_errDown = wjetsObj_eff.GetEYlow()[0]
# WJets: Using Standard Efficiency Definition
WJets_efficiency_estimate = XSectPassedWJets/XSectTotalWJets

### Calculate Purity of QCD
# Using Bayes
qcdObj_pur = bayesDivide(EvtsPassedQCD, EvtsPassedEWK+EvtsPassedQCD)
QCD_purity_bayes = qcdObj_pur.GetY()[0]
QCD_pur_errUp = qcdObj_pur.GetEYhigh()[0]
QCD_pur_errDown = qcdObj_pur.GetEYlow()[0]
# Using Standard Purity Definition
QCD_purity_estimate = XSectPassedQCD/(XSectPassedEWK+XSectPassedQCD)
# Try to roughly place an Upper Limit
QCD_purity_UpLimit = XSectPassedQCD_UpLimit/(XSectPassedEWK+XSectPassedQCD_UpLimit)
######################################################################################
### Calculate Efficiency using the ROOT class "bayesDivide", for each QCD DataSet and write to a file
# Declare (double) arrays to be filled with efficiencies and purities
effArray         = array("d",[ ])
effArray_errDown =  array("d",[ ])
effArray_errUp   =  array("d",[ ])

#filename = "../efficiencies.dat"
filename = "efficiencies.dat"
EfficienciesFile = open(filename,"a")

effObj = bayesDivide(qcd30to50[-1][1], qcd30to50[4][1])
eff = effObj.GetY()[0]
eff_errUp = effObj.GetEYhigh()[0]
eff_errDown = effObj.GetEYlow()[0]
effArray.append(eff)
effArray_errDown.append(eff_errDown)
effArray_errUp.append(eff_errUp)
EfficienciesFile.write(str(eff) + " ")
EfficienciesFile.write(str(eff_errDown) + " ")
EfficienciesFile.write(str(eff_errUp) + " ")
EfficienciesFile.write("\n")

effObj = bayesDivide(qcd50to80[-1][1], qcd50to80[4][1])
eff = effObj.GetY()[0]
eff_errUp = effObj.GetEYhigh()[0]
eff_errDown = effObj.GetEYlow()[0]
effArray.append(eff)
effArray_errDown.append(eff_errDown)
effArray_errUp.append(eff_errUp)
EfficienciesFile.write(str(eff) + " ")
EfficienciesFile.write(str(eff_errDown) + " ")
EfficienciesFile.write(str(eff_errUp) + " ")
EfficienciesFile.write("\n")

effObj = bayesDivide(qcd80to120[-1][1], qcd80to120[4][1])
eff = effObj.GetY()[0]
eff_errUp = effObj.GetEYhigh()[0]
eff_errDown = effObj.GetEYlow()[0]
effArray.append(eff)
effArray_errDown.append(eff_errDown)
effArray_errUp.append(eff_errUp)
EfficienciesFile.write(str(eff) + " ")
EfficienciesFile.write(str(eff_errDown) + " ")
EfficienciesFile.write(str(eff_errUp) + " ")
EfficienciesFile.write("\n")

effObj = bayesDivide(qcd120to170[-1][1], qcd120to170[4][1])
eff = effObj.GetY()[0]
eff_errUp = effObj.GetEYhigh()[0]
eff_errDown = effObj.GetEYlow()[0]
effArray.append(eff)
effArray_errDown.append(eff_errDown)
effArray_errUp.append(eff_errUp)
EfficienciesFile.write(str(eff) + " ")
EfficienciesFile.write(str(eff_errDown) + " ")
EfficienciesFile.write(str(eff_errUp) + " ")
EfficienciesFile.write("\n")

effObj = bayesDivide(qcd170to230[-1][1], qcd170to230[4][1])
eff = effObj.GetY()[0]
eff_errUp = effObj.GetEYhigh()[0]
eff_errDown = effObj.GetEYlow()[0]
effArray.append(eff)
effArray_errDown.append(eff_errDown)
effArray_errUp.append(eff_errUp)
EfficienciesFile.write(str(eff) + " ")
EfficienciesFile.write(str(eff_errDown) + " ")
EfficienciesFile.write(str(eff_errUp) + " ")
EfficienciesFile.write("\n")

effObj = bayesDivide(qcd230to300[-1][1], qcd230to300[4][1])
eff = effObj.GetY()[0]
eff_errUp = effObj.GetEYhigh()[0]
eff_errDown = effObj.GetEYlow()[0]
effArray.append(eff)
effArray_errDown.append(eff_errDown)
effArray_errUp.append(eff_errUp)
EfficienciesFile.write(str(eff) + " ")
EfficienciesFile.write(str(eff_errDown) + " ")
EfficienciesFile.write(str(eff_errUp) + " ")
EfficienciesFile.write("\n")

######################################################################################
if verbosity == True:
    print
    print "**************"
    print "* Efficiency *"
    print "**************"
    print "effArray = ", effArray
    print
    print "effArray_errDown = ", effArray_errDown
    print
    print "effArray_errUp = ", effArray_errUp

    print "*******************"
    print "* All QCD Samples *"
    print "*******************"
    print "--> Purity (Bayes) = ", QCD_purity_bayes
    print "--> Purity (Estimate) = ", QCD_purity_estimate
    print "--> Purity (Upper Limit) = ", QCD_purity_UpLimit
    print
    print "--> Efficiency (Bayes) = ", QCD_efficiency_bayes, " + ", QCD_eff_errUp, " - ", QCD_eff_errDown
    print "--> Efficiency (Estimate) = ", QCD_efficiency_estimate
    print "****************"
    print "* TTbar Sample *"
    print "****************"
    print "--> Efficiency (Bayes) = ", TTbar_efficiency_bayes, " + ", TTbar_eff_errUp, " - ", TTbar_eff_errDown
    print "--> Efficiency (Estimate) = ", TTbar_efficiency_estimate
    print "****************"
    print "* WJets Sample *"
    print "****************"
    print "--> Efficiency (Bayes) = ", WJets_efficiency_bayes, " + ", WJets_eff_errUp, " - ", WJets_eff_errDown
    print "--> Efficiency (Estimate) = ", WJets_efficiency_estimate
    print "***************************************************************************************************************"
######################################################################################

### Customise ROOT
gROOT.Reset()
ROOT.gROOT.SetBatch(False)
style = TDRStyle()
style.setPalettePretty()
style.setWide(True)
gStyle.SetOptStat(0000)
#gROOT.SetStyle("Plain");

### Create TCanvas and TGraph with Asymmetric Error Bars using Bayesian Statistical Tools
cEff = TCanvas ("Efficiencies", "Efficiencies", 20, 20, 800, 800)
cEff.cd()

qcdBin  = array("d",[50, 80, 120, 170, 230, 300 ])
qcdBin_down = array("d",[0, 0, 0, 0, 0, 0 ])
qcdBin_up  = array("d",[0, 0, 0, 0, 0, 0 ])
graph = TGraphAsymmErrors(6, qcdBin, effArray, qcdBin_down,  qcdBin_up, effArray_errDown, effArray_errUp) ### works

graph.Draw("AP")
graph.SetMarkerStyle(kFullCircle)
graph.SetMarkerColor(kRed)
graph.SetMarkerSize(1)
gPad.SetLogy(0)
gPad.SetLogx(0)
graph.GetYaxis().SetTitle("Efficiency (#epsilon)")
graph.GetXaxis().SetTitle("#hat{p}_{T}")
### Finally, re-draw graph and update canvas and gPad
graph.Draw("AP")
cEff.Update()
#gPad.Update()
##########################################################################################################################################
### Purity Measurement

### Set Options
xMin  = 0.0
xMax  = 1.0
nBins = 200
nPseudoExperiments = 100;
confidenceLevel = 0.683;

print "*** Creating Histograms"
# Create Histograms. Number of bins = 200, Max = 100, Min = -100
hBinomial_qcd30to50   =  TH1F("qcd30to50", "qcd30to50", 100000, 0, 1000000)
hBinomial_qcd50to80   =  TH1F("qcd50to80", "qcd50to80", 100000, 0, 1000000)
hBinomial_qcd80to120  =  TH1F("qcd80to120", "qcd80to120", 100000, 0, 1000000)
hBinomial_qcd120to170 =  TH1F("qcd120to170", "qcd120to170", 100000, 0, 1000000)
hBinomial_qcd170to230 =  TH1F("qcd170to230", "qcd170to230", 100000, 0, 1000000)
hBinomial_qcd230to300 =  TH1F("qcd230to300", "qcd230to300", 100000, 0, 1000000)
hBinomial_TTbar       =  TH1F("TTbar", "TTbar", 100000, 0, 1000000)
hBinomial_WJets       =  TH1F("WJets", "Wjets", 100000, 0, 1000000 )
hQCDPurity            =  TH1F("QCDPurity", "QCDPurity", nBins, xMin, xMax)

#print "effArray = ", effArray
### Calculate Events at given Integrated Luminosity
IntLumi = 40
qcd30to50_EvtsPassed   = qcd30to50[4][1]*IntLumi
qcd50to80_EvtsPassed   = qcd50to80[4][1]*IntLumi
qcd80to120_EvtsPassed  = qcd80to120[4][1]*IntLumi
qcd120to170_EvtsPassed = qcd120to170[4][1]*IntLumi
qcd170to230_EvtsPassed = qcd170to230[4][1]*IntLumi
qcd230to300_EvtsPassed = qcd230to300[4][1]*IntLumi
### Calculate Events at given Integrated Luminosity, before the Cut-Flow (but after Trigger, Anti-Isolation and Njets cuts)
qcd30to50_Evts   = qcd30to50[4][1]*IntLumi*qcd30to50_NormFactor
qcd50to80_Evts   = qcd50to80[4][1]*IntLumi*qcd50to80_NormFactor
qcd80to120_Evts  = qcd80to120[4][1]*IntLumi*qcd80to120_NormFactor
qcd120to170_Evts = qcd120to170[4][1]*IntLumi*qcd120to170_NormFactor
qcd170to230_Evts = qcd170to230[4][1]*IntLumi*qcd170to230_NormFactor
qcd230to300_Evts = qcd230to300[4][1]*IntLumi*qcd230to300_NormFactor
ttbar_Evts       = TTbar[4][1]*IntLumi*TTbarJets_NormFactor
wjets_Evts       = WJets[4][1]*IntLumi*WJets_NormFactor
qcd_Evts         = qcd30to50_Evts+qcd50to80_Evts+qcd80to120_Evts+qcd120to170_Evts+qcd170to230_Evts+qcd230to300_Evts
### Calculate the QCD Purity before the Cut-Flow (but after Trigger, Anti-Isolation and Njets cuts)
QCD_Purity_beforeCutFlow = (qcd_Evts)/(qcd_Evts+ttbar_Evts+wjets_Evts)
### Calculate xSections of samples
qcd30to50_xSectionPassedJetSelection   = qcd30to50[4][1]
qcd50to80_xSectionPassedJetSelection   = qcd50to80[4][1]
qcd80to120_xSectionPassedJetSelection  = qcd80to120[4][1]
qcd120to170_xSectionPassedJetSelection = qcd120to170[4][1]
qcd170to230_xSectionPassedJetSelection = qcd170to230[4][1]
qcd230to300_xSectionPassedJetSelection = qcd230to300[4][1]
ttbar_xSectionPassedJetSelection       = 38.03 #TTbar[4][1]
wjets_xSectionPassedJetSelection       = 19.14 #WJets[4][1]
    
### Define efficiencies and get values from the efficiency array
qcd30to50_eff   = effArray[0]
qcd50to80_eff   = effArray[1]
qcd80to120_eff  = effArray[2]
qcd120to170_eff = effArray[3]
qcd170to230_eff = effArray[4]
qcd230to300_eff = effArray[5]
ttbar_eff       = TTbar[-1][1]/TTbar[4][1]
wjets_eff       = WJets[-1][1]/WJets[4][1]

print "*** Filling Histograms"

### Loop to fill the histograms using the random number generator "gRandom" 
for i in range(nPseudoExperiments):
### Note: gRandom->Binomial(Total Number of generated Events for sample, probability that a given event will pass the given cut flow)
### Estimate the number of events that will pass all cut-flow according to the binomial distribution with parameters:
### par1 = Total number of events before cuts, par2 = probability that they will pass selection (i.e. efficiency of cut)
### Using number of events at specific integrated luminosity (IntLumi)
    percent = float(i)/float(float(nPseudoExperiments)/100.0);
    #percent = float(i) / float(nPseudoExperiments/100.0)
    if( percent%10 == 0):
        print " Progress: ", str(percent), "%"

### Generate Binomial distributions for all samples, given the efficiencies of the cut flow.
    qcd30to50_EvtsAfterJetSelBinomial   = gRandom.Binomial( int(qcd30to50_xSectionPassedJetSelection*IntLumi), qcd30to50_eff)
    qcd50to80_EvtsAfterJetSelBinomial   = gRandom.Binomial( int(qcd50to80_xSectionPassedJetSelection*IntLumi), qcd50to80_eff)
    qcd80to120_EvtsAfterJetSelBinomial  = gRandom.Binomial( int(qcd80to120_xSectionPassedJetSelection*IntLumi), qcd80to120_eff)
    qcd120to170_EvtsAfterJetSelBinomial = gRandom.Binomial( int(qcd120to170_xSectionPassedJetSelection*IntLumi), qcd120to170_eff)
    qcd170to230_EvtsAfterJetSelBinomial = gRandom.Binomial( int(qcd170to230_xSectionPassedJetSelection*IntLumi), qcd170to230_eff)
    qcd230to300_EvtsAfterJetSelBinomial = gRandom.Binomial( int(qcd230to300_xSectionPassedJetSelection*IntLumi), qcd230to300_eff)
    ttbar_EvtsAfterJetSelBinomial       = gRandom.Binomial( int(ttbar_xSectionPassedJetSelection*IntLumi), ttbar_eff)
    wjets_EvtsAfterJetSelBinomial       = gRandom.Binomial( int(wjets_xSectionPassedJetSelection*IntLumi), wjets_eff)

### Fill Histos with binomial numbers
    hBinomial_qcd30to50.Fill(qcd30to50_EvtsAfterJetSelBinomial)
    hBinomial_qcd50to80.Fill(qcd50to80_EvtsAfterJetSelBinomial)
    hBinomial_qcd80to120.Fill(qcd80to120_EvtsAfterJetSelBinomial)
    hBinomial_qcd120to170.Fill(qcd120to170_EvtsAfterJetSelBinomial)
    hBinomial_qcd170to230.Fill(qcd170to230_EvtsAfterJetSelBinomial)
    hBinomial_qcd230to300.Fill(qcd230to300_EvtsAfterJetSelBinomial)
    hBinomial_TTbar.Fill(ttbar_EvtsAfterJetSelBinomial)
    hBinomial_WJets.Fill(wjets_EvtsAfterJetSelBinomial)
    
### Calculate combined events passed for QCD (all Pt bins)
    qcd_EvtsAfterJetSelBinomial = qcd30to50_EvtsAfterJetSelBinomial + qcd50to80_EvtsAfterJetSelBinomial + qcd80to120_EvtsAfterJetSelBinomial + qcd120to170_EvtsAfterJetSelBinomial + qcd170to230_EvtsAfterJetSelBinomial + qcd230to300_EvtsAfterJetSelBinomial  

    QCD_Purity = (qcd_EvtsAfterJetSelBinomial)/(qcd_EvtsAfterJetSelBinomial+ttbar_EvtsAfterJetSelBinomial+wjets_EvtsAfterJetSelBinomial)

### Fill purity histos
    hQCDPurity.Fill(QCD_Purity)
    #binMax = hQCDPurity.GetMaximumBin()
    #xMaxNonEmpty = hQCDPurity.GetXaxis().GetBinCenter(binMax)
    #xMinNonEmpty = hQCDPurity.GetXaxis().GetBinCenter(FindFirstBinAbove(hQCDPurity, 0, 1))
    #hQCDPurity.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.5) , xMaxNonEmpty*(1+0.5))

addCmsPreliminaryText() 
print "*** Creating Canvas"
### Create a canvas to put the histograms on.
cBinom =  TCanvas("Binomials", "Binomials", 1)
cPur =  TCanvas("Purity", "Purity", 1)
cBinom.Divide(2,4)
cBinom.cd(1)
hBinomial_qcd30to50.Draw()
hBinomial_qcd30to50.SetFillColor(kRed+1)
hBinomial_qcd30to50.GetXaxis().SetTitle( "Events" )
binMax = hBinomial_qcd30to50.GetMaximumBin()
xMaxNonEmpty = hBinomial_qcd30to50.GetXaxis().GetBinCenter(binMax)
xMinNonEmpty = hBinomial_qcd30to50.GetXaxis().GetBinCenter(FindFirstBinAbove(hBinomial_qcd30to50, 0, 1))
hBinomial_qcd30to50.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.2) , xMaxNonEmpty*(1+0.2))
addCmsPreliminaryText()
###
cBinom.cd(2)
hBinomial_qcd50to80.Draw()
hBinomial_qcd50to80.SetFillColor(kRed+2)
hBinomial_qcd50to80.GetXaxis().SetTitle( "Events" )
binMax = hBinomial_qcd50to80.GetMaximumBin()
xMaxNonEmpty = hBinomial_qcd50to80.GetXaxis().GetBinCenter(binMax)
xMinNonEmpty = hBinomial_qcd50to80.GetXaxis().GetBinCenter(FindFirstBinAbove(hBinomial_qcd50to80, 0, 1))
hBinomial_qcd50to80.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.2) , xMaxNonEmpty*(1+0.2))
addCmsPreliminaryText()
###
cBinom.cd(3)
hBinomial_qcd80to120.Draw()
hBinomial_qcd80to120.SetFillColor(kRed+3)
hBinomial_qcd80to120.GetXaxis().SetTitle( "Events" )
binMax = hBinomial_qcd80to120.GetMaximumBin()
xMaxNonEmpty = hBinomial_qcd80to120.GetXaxis().GetBinCenter(binMax)
xMinNonEmpty = hBinomial_qcd80to120.GetXaxis().GetBinCenter(FindFirstBinAbove(hBinomial_qcd80to120, 0, 1))
hBinomial_qcd80to120.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.2) , xMaxNonEmpty*(1+0.2))
addCmsPreliminaryText()
###
cBinom.cd(4)
hBinomial_qcd120to170.Draw()
hBinomial_qcd120to170.SetFillColor(kRed+4)
hBinomial_qcd120to170.GetXaxis().SetTitle( "Events" )
binMax = hBinomial_qcd120to170.GetMaximumBin()
xMaxNonEmpty = hBinomial_qcd120to170.GetXaxis().GetBinCenter(binMax)
xMinNonEmpty = hBinomial_qcd120to170.GetXaxis().GetBinCenter(FindFirstBinAbove(hBinomial_qcd120to170, 0, 1))
hBinomial_qcd120to170.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.2) , xMaxNonEmpty*(1+0.2))
addCmsPreliminaryText()
###
cBinom.cd(5)
hBinomial_qcd170to230.Draw()
hBinomial_qcd170to230.SetFillColor(kRed-1)
hBinomial_qcd170to230.GetXaxis().SetTitle( "Events" )
binMax = hBinomial_qcd170to230.GetMaximumBin()
xMaxNonEmpty = hBinomial_qcd170to230.GetXaxis().GetBinCenter(binMax)
xMinNonEmpty = hBinomial_qcd170to230.GetXaxis().GetBinCenter(FindFirstBinAbove(hBinomial_qcd170to230, 0, 1))
hBinomial_qcd170to230.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.2) , xMaxNonEmpty*(1+0.2))
addCmsPreliminaryText()
###
cBinom.cd(6)
hBinomial_qcd230to300.Draw()
hBinomial_qcd230to300.SetFillColor(kRed-2)
hBinomial_qcd230to300.GetXaxis().SetTitle( "Events" )
binMax = hBinomial_qcd230to300.GetMaximumBin()
xMaxNonEmpty = hBinomial_qcd230to300.GetXaxis().GetBinCenter(binMax)
xMinNonEmpty = hBinomial_qcd230to300.GetXaxis().GetBinCenter(FindFirstBinAbove(hBinomial_qcd230to300, 0, 1))
hBinomial_qcd230to300.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.2) , xMaxNonEmpty*(1+0.2))
addCmsPreliminaryText()
###
cBinom.cd(7)
hBinomial_TTbar.Draw()
hBinomial_TTbar.SetFillColor(kBlack)
hBinomial_TTbar.GetXaxis().SetTitle( "Events" )
binMax = hBinomial_TTbar.GetMaximumBin()
xMaxNonEmpty = hBinomial_TTbar.GetXaxis().GetBinCenter(binMax)
xMinNonEmpty = hBinomial_TTbar.GetXaxis().GetBinCenter(FindFirstBinAbove(hBinomial_TTbar, 0, 1))
hBinomial_TTbar.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.5) , xMaxNonEmpty*(1+0.5))
###
cBinom.cd(8)
hBinomial_WJets.Draw()
hBinomial_WJets.SetFillColor(kBlue)
hBinomial_WJets.GetXaxis().SetTitle( "Events" )
binMax = hBinomial_WJets.GetMaximumBin()
xMaxNonEmpty = hBinomial_WJets.GetXaxis().GetBinCenter(binMax)
xMinNonEmpty = hBinomial_WJets.GetXaxis().GetBinCenter(FindFirstBinAbove(hBinomial_WJets, 0, 1))
hBinomial_WJets.GetXaxis().SetRangeUser( xMinNonEmpty*(1-0.5) , xMaxNonEmpty*(1+0.5))
addCmsPreliminaryText()

### Purity plot
cPur.cd()
hQCDPurity.Draw()
hQCDPurity.SetFillColor(kOrange)
hQCDPurity.GetXaxis().SetTitle( "QCD Purity" )
hQCDPurity.GetYaxis().SetTitle( "Pseudo-Experiments" )
### Definitions
binWidth = hQCDPurity.GetBinWidth(1);  
TotalIntegral = hQCDPurity.Integral( hQCDPurity.FindBin(0.0), hQCDPurity.FindBin(1.0) );
### Define lists to be used
intervals = []
lowBin   = []
upBin    = []

### Determine lowest bin where an entry is found (to minimise the loop iterations)
lowLimitBin = -1
upLimitBin  = -1
for i in range(0,hQCDPurity.GetNbinsX()+1):
    if (hQCDPurity.GetBinContent(i) > 0.0 ):
        lowLimitBin = i
        break

### Determine max bin where an entry is found (to minimise the loop iterations). Reverse (decrement) loop
for j in range(hQCDPurity.GetNbinsX()+1, 0, -1):
    if (hQCDPurity.GetBinContent(j) > 0.0 ):
        upLimitBin = j
        break
### Perform a double loop to determine all the intervals that contain an area which is at least equal to the confidenceLevel
for i in range ( lowLimitBin, hQCDPurity.GetNbinsX()+1, +1 ):
### Second loop index set to first loop index to avoid double (or reverse) calculations
    for j in range ( i, upLimitBin+1, +1 ):
        tmpIntegral = hQCDPurity.Integral(i, j, "")
        #print "tmpIntegral = ", tmpIntegral
        tmpIntegralFraction = tmpIntegral/TotalIntegral
        if(tmpIntegralFraction >= confidenceLevel):
            a = i
            b = j
            interval = b-a
      ### Save all the intervals that contain the area >= confidence level. We will find the minimum later on.
            intervals.append(interval)
            lowBin.append(a)
            upBin.append(b)
            break
### Definition of variables
myMin = 9999999;
k = 0;
index = 0;
QCD_Purity_mean = hQCDPurity.GetMean()

### Loop over stored intervals to find the minimum interval containing an area >= to the confidenceLevel chosen.
for s in range( len(intervals) ):
### Determine the minimum by comparing all values. Store the vector index.
    if( intervals[s] < myMin):
        myMin = intervals[s]
        index = s

### Determine lowError (Less than mean) and highError (Greater than mean)
result   = QCD_Purity_mean*100
lowError = fabs(QCD_Purity_mean-(lowBin[index]*binWidth))*100
upError  = fabs(QCD_Purity_mean-(upBin[index]*binWidth))*100

### Then create the TString that contains them and draw it on the current activated canvas
myText  = "(" + str(format(result, ('.1f'))) + "^{+" + str(format(upError, ('.1f'))) + "}_{-" + str(format(lowError, ('.1f'))) + "})%"
myYCoordinate = hQCDPurity.GetMaximum()
myXCoordinate = (xMax-xMin)*(0.10)
#text = TText(myXCoordinate, myYCoordinate*(1-0.1), myText)
text = TLatex(myXCoordinate, myYCoordinate*(1-0.1), myText)
text.Draw()
gPad.Update()

### Let's create a file and write it to disk.
filename = "../purities.dat"
### Open file in "append" mode
PuritiesFile = open(filename,"a")
### Append to the file
PuritiesFile.write(str(result) + " ")
PuritiesFile.write(str(lowError) + " ")
PuritiesFile.write(str(upError) + " ")
PuritiesFile.write("\n")
PuritiesFile.close()
EfficienciesFile.close()

### Options
ROOT.gPad.SetLogy(False)
addCmsPreliminaryText()
SavePlots = True
SavePath  = "/afs/cern.ch/user/a/attikis/qcdMeasurementByIsolationVeto/eff/HltJet30u_MetCut_36X_MC_pattuple_v6_1/figures/"
SaveAs    = "pseudoExperiments_MetCut70"
### Save canvases to file
if(SavePlots):
    print "*** Saving canvases to path: ", SavePath
    cEff.SaveAs( SavePath + SaveAs + "_cEff_" + ".png");
    cEff.SaveAs( SavePath + SaveAs + "_cEff_" + ".C");
    cEff.SaveAs( SavePath + SaveAs + "_cEff_" + ".eps");
    cEff.SaveAs( SavePath + SaveAs + "_cEff_" + ".pdf");

    cPur.SaveAs( SavePath + SaveAs + "_cPur_" + ".png");
    cPur.SaveAs( SavePath + SaveAs + "_cPur_" + ".C");
    cPur.SaveAs( SavePath + SaveAs + "_cPur_" + ".eps");
    cPur.SaveAs( SavePath + SaveAs + "_cPur_" + ".pdf");

    cBinom.SaveAs( SavePath + SaveAs + "_cBinom_" + ".png");
    cBinom.SaveAs( SavePath + SaveAs + "_cBinom_" + ".C");
    cBinom.SaveAs( SavePath + SaveAs + "_cBinom_" + ".eps");
    cBinom.SaveAs( SavePath + SaveAs + "_cBinom_" + ".pdf");

### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
raw_input('*** Done! Press \"Enter\" to exit to terminal ...')
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
