UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py

from DatabaseParser import *
import sys
import os
from numpy import *
import pickle

from ROOT import gROOT, TCanvas, TF1, TGraph, TGraphErrors, TPaveStats, gPad, gStyle
from ROOT import TFile, TPaveText
from ROOT import gBenchmark
import array
import math

def main():

    ######## TO CREATE FITS #########
##     run_list = [179497,179547,179558,179563,179889,179959,179977,180072,180076,180093,180241,180250,180252]
##     trig_name = "HLT"
##     num_ls = 20
##     debug_print = False

##     min_rate = 1.0
##     print_table = False
##     data_clean = True
##     ##plot_properties = [varX, varY, do_fit, save_root, save_png, overlay_fit, fit_file]
##     plot_properties = [["live", "rate", True, True, False, False, ""]]
##     masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
##     save_fits = True
    

    ######## TO SEE RATE VS PREDICTION ########
    run_list = [180241]

    trig_name = "Mu"
    num_ls = 1
    debug_print = False

    min_rate = 3.0
    print_table = False
    data_clean = False
    ##plot_properties = [varX, varY, do_fit, save_root, save_png, overlay_fit, fit_file]
    plot_properties = [["ls", "rawrate", False, True, False, True,"Fits/2011/Fit_HLT_20LS_Run179497to180252.pkl"]]
    masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
    save_fits = False
    
    ########  END PARAMETERS ##########
    
    Rates = GetDBRates(run_list, trig_name, num_ls, debug_print)
    MakePlots(Rates, run_list, trig_name, num_ls, min_rate, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print)
    

def GetDBRates(run_list,trig_name,num_ls, debug_print):
    
    Rates = {}
    RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS.pkl"
    RefRunFile = RefRunNameTemplate % (trig_name,num_ls)
    try:
        pkl_file = open(RefRunFile, 'rb')
        Rates = pickle.load(pkl_file)
        pkl_file.close()
        os.remove(RefRunFile)
    except:
        print str(RefRunFile)+" does not exist. Creating ..."

    for RefRunNum in run_list: #Will change to a "count back runs" scheme, or something like that

        try:
            ExistsAlready = False
            for key in Rates:
                if RefRunNum in Rates[key]["run"]:
                    ExistsAlready = True
                    break
            if ExistsAlready:
                continue
        except:
            print "Getting info for run "+str(RefRunNum)
        
        if RefRunNum < 1:
            continue

        if True:
            if True: #May replace with "try" - for now it's good to know when problems happen
                RefParser = DatabaseParser()
                RefParser.RunNumber = RefRunNum
                RefParser.ParseRunSetup()
                RefLumiRange = RefParser.GetLSRange(1,9999)
                #RefLumiArray = RefParser.GetLumiInfo()

                nls = RefLumiRange[0]
                LSRange = {}
                if nls >= RefLumiRange[-1]-num_ls:
                    print "Run "+str(RefRunNum)+" is too short: from "+str(nls)+" to "+str(RefLumiRange[-1])+", while num_ls = "+str(num_ls)
                    continue
                while nls < RefLumiRange[-1]-num_ls:
                    if num_ls > 1:
                        LSRange[nls] = RefParser.GetLSRange(nls,num_ls)
                    else:
                        LSRange[nls] = [nls]
                    nls = LSRange[nls][-1]+1
                print "Run "+str(RefRunNum)+" contains LS from "+str(min(LSRange))+" to "+str(max(LSRange))
                for nls in LSRange:
                    TriggerRates = RefParser.GetHLTRates(LSRange[nls])
                    [inst, live, delivered, dead, pscols] = RefParser.GetAvLumiInfo(LSRange[nls])
                    for key in TriggerRates:
                        if not trig_name in key:
                            continue
                        name = key
                        if re.match('.*_v[0-9]+',name):
                            name = name[:name.rfind('_')]

                        if not Rates.has_key(name):
                            Rates[name] = {}
                            Rates[name]["run"] = []
                            Rates[name]["ls"] = []
                            Rates[name]["ps"] = []
                            Rates[name]["inst_lumi"] = []
                            Rates[name]["live_lumi"] = []
                            Rates[name]["delivered_lumi"] = []
                            Rates[name]["deadtime"] = []
                            Rates[name]["rawrate"] = []
                            Rates[name]["rate"] = []
                            Rates[name]["rawxsec"] = []
                            Rates[name]["xsec"] = []

                        if delivered == 0:
                            continue
                        [avps, ps, rate, psrate] = TriggerRates[key]
                        Rates[name]["run"].append(RefRunNum)
                        Rates[name]["ls"].append(nls)
                        Rates[name]["ps"].append(ps)
                        Rates[name]["inst_lumi"].append(inst)
                        Rates[name]["live_lumi"].append(live)
                        Rates[name]["delivered_lumi"].append(delivered)
                        Rates[name]["deadtime"].append(dead)
                        Rates[name]["rawrate"].append(rate)
                        Rates[name]["rate"].append(psrate/(1.0-dead))
                        Rates[name]["rawxsec"].append(rate/delivered)
                        Rates[name]["xsec"].append(psrate/((1.0-dead)*delivered))
            #except:
                #print "Failed to parse run "+str(RefRunNum)

    RateOutput = open(RefRunFile, 'wb')
    pickle.dump(Rates, RateOutput, 2)
    RateOutput.close()
    #print Rates
    return Rates

def MakePlots(Rates, run_list, trig_name, num_ls, min_rate, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print):
    min_run = min(run_list)
    max_run = max(run_list)

    Input = {}
    Output = {}

    RootNameTemplate = "%s_%sLS_Run%sto%s.root"
    RootFile = RootNameTemplate % (trig_name, num_ls, min_run, max_run)

    for varX, varY, do_fit, save_root, save_png, overlay_fit, fit_file in plot_properties:
        if overlay_fit:
            pkl_file = open(fit_file, 'rb')
            Input = pickle.load(pkl_file)
            pkl_file.close()
        if save_root:
            try:
                os.remove(RootFile)
            except:
                break

    for print_trigger in Rates:
        meanrawrate = sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])
        if not trig_name in print_trigger:
            continue
        if meanrawrate < min_rate:
            continue
        masked_trig = False
        for mask in masked_triggers:
            if str(mask) in print_trigger:
                masked_trig = True
        if masked_trig:
            continue
                
        Output[print_trigger] = {}

        lowlumi = 0
        meanlumi_init = sum(Rates[print_trigger]["live_lumi"])/len(Rates[print_trigger]["live_lumi"])
        meanlumi = 0
        highlumi = 0
        lowxsec = 0
        meanxsec = 0
        highxsec = 0
        nlow = 0
        nhigh = 0
        for iterator in range(len(Rates[print_trigger]["rate"])):
            if not Rates[print_trigger]["run"][iterator] in run_list:
                continue
            if Rates[print_trigger]["live_lumi"][iterator] <= meanlumi_init:
                if not data_clean or Rates[print_trigger]["rawrate"][iterator] > 0.04:
                    meanxsec+=Rates[print_trigger]["xsec"][iterator]
                    lowxsec+=Rates[print_trigger]["xsec"][iterator]
                    meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
                    lowlumi+=Rates[print_trigger]["live_lumi"][iterator]
                    nlow+=1
            if Rates[print_trigger]["live_lumi"][iterator] > meanlumi_init:
                if not data_clean or Rates[print_trigger]["rawrate"][iterator] > 0.04:
                    meanxsec+=Rates[print_trigger]["xsec"][iterator]
                    highxsec+=Rates[print_trigger]["xsec"][iterator]
                    meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
                    highlumi+=Rates[print_trigger]["live_lumi"][iterator]
                    nhigh+=1
        meanxsec = meanxsec/(nlow+nhigh)
        meanlumi = meanlumi/(nlow+nhigh)
        slopexsec = ( (highxsec/nhigh) - (lowxsec/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )

        run_t = array.array('f')
        ls_t = array.array('f')
        ps_t = array.array('f')
        inst_t = array.array('f')
        live_t = array.array('f')
        delivered_t = array.array('f')
        deadtime_t = array.array('f')
        rawrate_t = array.array('f')
        rate_t = array.array('f')
        rawxsec_t = array.array('f')
        xsec_t = array.array('f')

        e_run_t = array.array('f')
        e_ls_t = array.array('f')
        e_ps_t = array.array('f')
        e_inst_t = array.array('f')
        e_live_t = array.array('f')
        e_delivered_t = array.array('f')
        e_deadtime_t = array.array('f')
        e_rawrate_t = array.array('f')
        e_rate_t = array.array('f')
        e_rawxsec_t = array.array('f')
        e_xsec_t = array.array('f')

        rawrate_fit_t = array.array('f')
        rate_fit_t = array.array('f')
        rawxsec_fit_t = array.array('f')
        xsec_fit_t = array.array('f')
        e_rawrate_fit_t = array.array('f')
        e_rate_fit_t = array.array('f')
        e_rawxsec_fit_t = array.array('f')
        e_xsec_fit_t = array.array('f')

        if overlay_fit:
            X0 = Input[print_trigger][0]
            X1 = Input[print_trigger][1]
            X2 = Input[print_trigger][2]
            Chi2 = Input[print_trigger][3]
            
        for iterator in range(len(Rates[print_trigger]["rate"])):
            if not Rates[print_trigger]["run"][iterator] in run_list:
                continue
            prediction = meanxsec + slopexsec * (Rates[print_trigger]["live_lumi"][iterator] - meanlumi)
            realvalue = Rates[print_trigger]["xsec"][iterator]
            if not data_clean or (realvalue > 0.4*prediction and realvalue < 2.5*prediction) or (realvalue > 0.4*meanxsec and realvalue < 2.5*meanxsec) or prediction < 0:
                run_t.append(Rates[print_trigger]["run"][iterator])
                ls_t.append(Rates[print_trigger]["ls"][iterator])
                ps_t.append(Rates[print_trigger]["ps"][iterator])
                inst_t.append(Rates[print_trigger]["inst_lumi"][iterator])
                live_t.append(Rates[print_trigger]["live_lumi"][iterator])
                delivered_t.append(Rates[print_trigger]["delivered_lumi"][iterator])
                deadtime_t.append(Rates[print_trigger]["deadtime"][iterator])
                rawrate_t.append(Rates[print_trigger]["rawrate"][iterator])
                rate_t.append(Rates[print_trigger]["rate"][iterator])
                rawxsec_t.append(Rates[print_trigger]["rawxsec"][iterator])
                xsec_t.append(Rates[print_trigger]["xsec"][iterator])

                e_run_t.append(0.0)
                e_ls_t.append(0.0)
                e_ps_t.append(0.0)
                e_inst_t.append(14.14)
                e_live_t.append(14.14)
                e_delivered_t.append(14.14)
                e_deadtime_t.append(1.0)
                e_rawrate_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
                e_rate_t.append(Rates[print_trigger]["ps"][iterator]*math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
                e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["delivered_lumi"][iterator])
                e_xsec_t.append(Rates[print_trigger]["ps"][iterator]*math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["live_lumi"][iterator])

                if overlay_fit:
                    rate_prediction = X0 + X1*live_t[-1] + X2*live_t[-1]*live_t[-1] 
                    rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
                    rate_fit_t.append(rate_prediction)
                    rawxsec_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]*delivered_t[-1]))
                    xsec_fit_t.append(rate_prediction/delivered_t[-1])
                    e_rawrate_fit_t.append(e_rawrate_t[-1]*sqrt(Chi2))
                    e_rate_fit_t.append(e_rate_t[-1]*sqrt(Chi2))
                    e_rawxsec_fit_t.append(e_rawxsec_t[-1]*sqrt(Chi2))
                    e_xsec_fit_t.append(e_xsec_t[-1]*sqrt(Chi2))
 
            else:
                if debug_print:
                    print str(print_trigger)+" has xsec "+str(round(Rates[print_trigger]["xsec"][iterator],6))+" at lumi "+str(round(Rates[print_trigger]["live_lumi"][iterator],2))+" where the expected value is "+str(prediction)

        for varX, varY, do_fit, save_root, save_png, overlay_fit, fit_file in plot_properties:
            if varX == "run":
                VX = run_t
                VXE = run_t_e
            elif varX == "ls":
                VX = ls_t
                VXE = e_ls_t
            elif varX == "ps":
                VX = ps_t
                VXE = e_ps_t
            elif varX == "inst":
                VX = inst_t
                VXE = e_inst_t
            elif varX == "live":
                VX = live_t
                VXE = e_live_t
            elif varX == "delivered":
                VX = delivered_t
                VXE = e_delivered_t
            elif varX == "rawrate":
                VX = rawrate_t
                VXE = e_rawrate_t
            elif varX == "rate":
                VX = rate_t
                VXE = e_rate_t
            elif varX == "rawxsec":
                VX = rawxsec_t
                VXE = e_rawxsec_t
            elif varX == "xsec":
                VX = xsec_t
                VXE = e_xsec_t
            else:
                print "No valid variable entered for X"
                continue
            if varY == "run":
                VY = run_t
                VYE = run_t_e
            elif varY == "ls":
                VY = ls_t
                VYE = e_ls_t
            elif varY == "ps":
                VY = ps_t
                VYE = e_ps_t
            elif varY == "inst":
                VY = inst_t
                VYE = e_inst_t
            elif varY == "live":
                VY = live_t
                VYE = e_live_t
            elif varY == "delivered":
                VY = delivered_t
                VYE = e_delivered_t
            elif varY == "rawrate":
                VY = rawrate_t
                VYE = e_rawrate_t
                if overlay_fit:
                    VF = rawrate_fit_t
                    VFE = e_rawrate_fit_t
            elif varY == "rate":
                VY = rate_t
                VYE = e_rate_t
                if overlay_fit:
                    VF = rate_fit_t
                    VFE = e_rate_fit_t
            elif varY == "rawxsec":
                VY = rawxsec_t
                VYE = e_rawxsec_t
                if overlay_fit:
                    VF = rawxsec_fit_t
                    VFE = e_rawxsec_fit_t
            elif varY == "xsec":
                VY = xsec_t
                VYE = e_xsec_t
                if overlay_fit:
                    VF = xsec_fit_t
                    VFE = e_xsec_fit_t
            else:
                print "No valid variable entered for Y"
                continue


            if save_root or save_png:
                c1 = TCanvas(str(varX),str(varY))
                c1.SetName(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))

            gr1 = TGraphErrors(len(VX), VX, VY, VXE, VYE)
            gr1.SetName("Graph_"+str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
            gr1.GetXaxis().SetTitle(varX)
            gr1.GetYaxis().SetTitle(varY)
            gr1.SetTitle(str(print_trigger))
            gr1.SetMinimum(0)
            gr1.SetMaximum(1.2*max(VY))
            gr1.GetXaxis().SetLimits(min(VX)-0.2*max(VX),1.2*max(VX))
            gr1.SetMarkerStyle(8)
            if overlay_fit:
                gr1.SetMarkerSize(0.8)
            else:
                gr1.SetMarkerSize(0.5)
            gr1.SetMarkerColor(2)

##             gr2 = TGraphErrors(len(VX), live_t, xsec_t, e_live_t, e_xsec_t)

            if overlay_fit:
                gr3 = TGraphErrors(len(VX), VX, VF, VXE, VFE)
                gr3.SetMarkerStyle(8)
                gr3.SetMarkerSize(0.4)
                gr3.SetMarkerColor(4)
                gr3.SetFillColor(4)
                gr3.SetFillStyle(3003)
            
            if do_fit:
                if "rate" in varY:
##                     f2a = TF1("f2a","pol1",0,8000)
##                     f2a.SetParLimits(0,0,meanxsec*1.5)
##                     if slopexsec > 0:
##                         f2a.SetParLimits(1,0,max(xsec_t)/max(live_t))
##                     else:
##                         f2a.SetParLimits(1,2*slopexsec,-2*slopexsec)
##                     gr2.Fit("f2a","Q","rob=0.80")
                    
                    f1a = TF1("f1a","pol2",0,8000)
                    f1a.SetLineColor(4)
                    f1a.SetLineWidth(2)
                    f1a.SetParLimits(0,0,1000)
                    f1a.SetParLimits(1,0,1000)
                    #gr1.Fit("f1a","B","Q")
                    gr1.Fit("f1a","Q","rob=0.80")

##                     f1b = TF1("f1b","pol2",0,8000)
##                     f1b.SetParameters(0.0,f2a.GetParameter(0),f2a.GetParameter(1))
##                     f1b.SetLineColor(3)
##                     f1b.SetLineWidth(2)
                    
                else:
                    f1a = TF1("f1a","pol1",0,8000)
                    f1a.SetLineColor(4)
                    f1a.SetLineWidth(2)
                    if "xsec" in varY:
                        f1a.SetParLimits(0,0,meanxsec*1.5)
                        if slopexsec > 0:
                            f1a.SetParLimits(1,0,max(VY)/max(VX))
                        else:
                            f1a.SetParLimits(1,2*slopexsec,-2*slopexsec)
                    else:
                        f1a.SetParLimits(0,0,1000)
                    gr1.Fit("f1a","Q","rob=0.80")

            if save_root or save_png:
                gr1.Draw("APZ")
                if overlay_fit:
                    gr3.Draw("P3")
                if do_fit:
                    f1a.Draw("same")
##                     if "rate" in varY:
##                         f1b.Draw("same")
                c1.Update()
                if save_root:
                    myfile = TFile( RootFile, 'UPDATE' )
                    c1.Write()
                    myfile.Close()
                if save_png:
                    c1.SaveAs(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX)+".png")
                
        
        ##p1 = TPaveStats()                                                                                                                              
        ##p1 = gr1.GetListOfFunctions().FindObject("stats")                                                                                              
        ##print p1                                                                                                                                       
        ##gr1.PaintStats(f1b).Draw("same")                                                                                                               

        if print_table or save_fits:
            Output[print_trigger] = [f1a.GetParameter(0), f1a.GetParameter(1), f1a.GetParameter(2), f1a.GetChisquare()/f1a.GetNDF(), f1a.GetParameter(0)+5000*(f1a.GetParameter(1)+f1a.GetParameter(2)*5000), meanrawrate]

    if save_root:
        print "Output root file is "+str(RootFile)

    if save_fits:
        FitNameTemplate = "Fits/2011/Fit_%s_%sLS_Run%sto%s.pkl"
        FitFile = FitNameTemplate % (trig_name, num_ls, min_run, max_run)
        if os.path.exists(FitFile):
            os.remove(FitFile)
        FitOutput = open(FitFile, 'wb')
        pickle.dump(Output, FitOutput, 2)
        FitOutput.close()

    if print_table:
        print '%60s%10s%10s%10s%10s%10s%10s' % ("Trig", "p0", "p1", "p2", "Chi2", "5e33 pred", "Av raw")
        for print_trigger in Output:
            _trigger = (print_trigger[:56] + '...') if len(print_trigger) > 59 else print_trigger
            try:
                print '%60s%10s%10s%10s%10s%10s%10s' % (_trigger, round(Output[print_trigger][0],3), round(Output[print_trigger][1],6)*1000, round(Output[print_trigger][2],9)*1000000, round(Output[print_trigger][3],2), round(Output[print_trigger][4],2) , round(Output[print_trigger][5],3))
            except:
                print str(print_trigger)+" is somehow broken"

if __name__=='__main__':
    main()
Revision:	1.1
Committed:	Thu Feb 23 16:44:36 2012 UTC (13 years, 2 months ago) by abrinke1
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	V00-00-07
Branch point for:	V00-00-06
Log Message:	Creates and stores rate histories, predictions, and overlays
#	User	Rev	Content
1	abrinke1	1.1	from DatabaseParser import *
2			import sys
3			import os
4			from numpy import *
5			import pickle
6
7			from ROOT import gROOT, TCanvas, TF1, TGraph, TGraphErrors, TPaveStats, gPad, gStyle
8			from ROOT import TFile, TPaveText
9			from ROOT import gBenchmark
10			import array
11			import math
12
13			def main():
14
15			######## TO CREATE FITS #########
16			## run_list = [179497,179547,179558,179563,179889,179959,179977,180072,180076,180093,180241,180250,180252]
17			## trig_name = "HLT"
18			## num_ls = 20
19			## debug_print = False
20
21			## min_rate = 1.0
22			## print_table = False
23			## data_clean = True
24			## ##plot_properties = [varX, varY, do_fit, save_root, save_png, overlay_fit, fit_file]
25			## plot_properties = [["live", "rate", True, True, False, False, ""]]
26			## masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
27			## save_fits = True
28
29
30			######## TO SEE RATE VS PREDICTION ########
31			run_list = [180241]
32
33			trig_name = "Mu"
34			num_ls = 1
35			debug_print = False
36
37			min_rate = 3.0
38			print_table = False
39			data_clean = False
40			##plot_properties = [varX, varY, do_fit, save_root, save_png, overlay_fit, fit_file]
41			plot_properties = [["ls", "rawrate", False, True, False, True,"Fits/2011/Fit_HLT_20LS_Run179497to180252.pkl"]]
42			masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
43			save_fits = False
44
45			######## END PARAMETERS ##########
46
47			Rates = GetDBRates(run_list, trig_name, num_ls, debug_print)
48			MakePlots(Rates, run_list, trig_name, num_ls, min_rate, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print)
49
50
51			def GetDBRates(run_list,trig_name,num_ls, debug_print):
52
53			Rates = {}
54			RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS.pkl"
55			RefRunFile = RefRunNameTemplate % (trig_name,num_ls)
56			try:
57			pkl_file = open(RefRunFile, 'rb')
58			Rates = pickle.load(pkl_file)
59			pkl_file.close()
60			os.remove(RefRunFile)
61			except:
62			print str(RefRunFile)+" does not exist. Creating ..."
63
64			for RefRunNum in run_list: #Will change to a "count back runs" scheme, or something like that
65
66			try:
67			ExistsAlready = False
68			for key in Rates:
69			if RefRunNum in Rates[key]["run"]:
70			ExistsAlready = True
71			break
72			if ExistsAlready:
73			continue
74			except:
75			print "Getting info for run "+str(RefRunNum)
76
77			if RefRunNum < 1:
78			continue
79
80			if True:
81			if True: #May replace with "try" - for now it's good to know when problems happen
82			RefParser = DatabaseParser()
83			RefParser.RunNumber = RefRunNum
84			RefParser.ParseRunSetup()
85			RefLumiRange = RefParser.GetLSRange(1,9999)
86			#RefLumiArray = RefParser.GetLumiInfo()
87
88			nls = RefLumiRange[0]
89			LSRange = {}
90			if nls >= RefLumiRange[-1]-num_ls:
91			print "Run "+str(RefRunNum)+" is too short: from "+str(nls)+" to "+str(RefLumiRange[-1])+", while num_ls = "+str(num_ls)
92			continue
93			while nls < RefLumiRange[-1]-num_ls:
94			if num_ls > 1:
95			LSRange[nls] = RefParser.GetLSRange(nls,num_ls)
96			else:
97			LSRange[nls] = [nls]
98			nls = LSRange[nls][-1]+1
99			print "Run "+str(RefRunNum)+" contains LS from "+str(min(LSRange))+" to "+str(max(LSRange))
100			for nls in LSRange:
101			TriggerRates = RefParser.GetHLTRates(LSRange[nls])
102			[inst, live, delivered, dead, pscols] = RefParser.GetAvLumiInfo(LSRange[nls])
103			for key in TriggerRates:
104			if not trig_name in key:
105			continue
106			name = key
107			if re.match('.*_v[0-9]+',name):
108			name = name[:name.rfind('_')]
109
110			if not Rates.has_key(name):
111			Rates[name] = {}
112			Rates[name]["run"] = []
113			Rates[name]["ls"] = []
114			Rates[name]["ps"] = []
115			Rates[name]["inst_lumi"] = []
116			Rates[name]["live_lumi"] = []
117			Rates[name]["delivered_lumi"] = []
118			Rates[name]["deadtime"] = []
119			Rates[name]["rawrate"] = []
120			Rates[name]["rate"] = []
121			Rates[name]["rawxsec"] = []
122			Rates[name]["xsec"] = []
123
124			if delivered == 0:
125			continue
126			[avps, ps, rate, psrate] = TriggerRates[key]
127			Rates[name]["run"].append(RefRunNum)
128			Rates[name]["ls"].append(nls)
129			Rates[name]["ps"].append(ps)
130			Rates[name]["inst_lumi"].append(inst)
131			Rates[name]["live_lumi"].append(live)
132			Rates[name]["delivered_lumi"].append(delivered)
133			Rates[name]["deadtime"].append(dead)
134			Rates[name]["rawrate"].append(rate)
135			Rates[name]["rate"].append(psrate/(1.0-dead))
136			Rates[name]["rawxsec"].append(rate/delivered)
137			Rates[name]["xsec"].append(psrate/((1.0-dead)*delivered))
138			#except:
139			#print "Failed to parse run "+str(RefRunNum)
140
141			RateOutput = open(RefRunFile, 'wb')
142			pickle.dump(Rates, RateOutput, 2)
143			RateOutput.close()
144			#print Rates
145			return Rates
146
147			def MakePlots(Rates, run_list, trig_name, num_ls, min_rate, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print):
148			min_run = min(run_list)
149			max_run = max(run_list)
150
151			Input = {}
152			Output = {}
153
154			RootNameTemplate = "%s_%sLS_Run%sto%s.root"
155			RootFile = RootNameTemplate % (trig_name, num_ls, min_run, max_run)
156
157			for varX, varY, do_fit, save_root, save_png, overlay_fit, fit_file in plot_properties:
158			if overlay_fit:
159			pkl_file = open(fit_file, 'rb')
160			Input = pickle.load(pkl_file)
161			pkl_file.close()
162			if save_root:
163			try:
164			os.remove(RootFile)
165			except:
166			break
167
168			for print_trigger in Rates:
169			meanrawrate = sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])
170			if not trig_name in print_trigger:
171			continue
172			if meanrawrate < min_rate:
173			continue
174			masked_trig = False
175			for mask in masked_triggers:
176			if str(mask) in print_trigger:
177			masked_trig = True
178			if masked_trig:
179			continue
180
181			Output[print_trigger] = {}
182
183			lowlumi = 0
184			meanlumi_init = sum(Rates[print_trigger]["live_lumi"])/len(Rates[print_trigger]["live_lumi"])
185			meanlumi = 0
186			highlumi = 0
187			lowxsec = 0
188			meanxsec = 0
189			highxsec = 0
190			nlow = 0
191			nhigh = 0
192			for iterator in range(len(Rates[print_trigger]["rate"])):
193			if not Rates[print_trigger]["run"][iterator] in run_list:
194			continue
195			if Rates[print_trigger]["live_lumi"][iterator] <= meanlumi_init:
196			if not data_clean or Rates[print_trigger]["rawrate"][iterator] > 0.04:
197			meanxsec+=Rates[print_trigger]["xsec"][iterator]
198			lowxsec+=Rates[print_trigger]["xsec"][iterator]
199			meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
200			lowlumi+=Rates[print_trigger]["live_lumi"][iterator]
201			nlow+=1
202			if Rates[print_trigger]["live_lumi"][iterator] > meanlumi_init:
203			if not data_clean or Rates[print_trigger]["rawrate"][iterator] > 0.04:
204			meanxsec+=Rates[print_trigger]["xsec"][iterator]
205			highxsec+=Rates[print_trigger]["xsec"][iterator]
206			meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
207			highlumi+=Rates[print_trigger]["live_lumi"][iterator]
208			nhigh+=1
209			meanxsec = meanxsec/(nlow+nhigh)
210			meanlumi = meanlumi/(nlow+nhigh)
211			slopexsec = ( (highxsec/nhigh) - (lowxsec/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
212
213			run_t = array.array('f')
214			ls_t = array.array('f')
215			ps_t = array.array('f')
216			inst_t = array.array('f')
217			live_t = array.array('f')
218			delivered_t = array.array('f')
219			deadtime_t = array.array('f')
220			rawrate_t = array.array('f')
221			rate_t = array.array('f')
222			rawxsec_t = array.array('f')
223			xsec_t = array.array('f')
224
225			e_run_t = array.array('f')
226			e_ls_t = array.array('f')
227			e_ps_t = array.array('f')
228			e_inst_t = array.array('f')
229			e_live_t = array.array('f')
230			e_delivered_t = array.array('f')
231			e_deadtime_t = array.array('f')
232			e_rawrate_t = array.array('f')
233			e_rate_t = array.array('f')
234			e_rawxsec_t = array.array('f')
235			e_xsec_t = array.array('f')
236
237			rawrate_fit_t = array.array('f')
238			rate_fit_t = array.array('f')
239			rawxsec_fit_t = array.array('f')
240			xsec_fit_t = array.array('f')
241			e_rawrate_fit_t = array.array('f')
242			e_rate_fit_t = array.array('f')
243			e_rawxsec_fit_t = array.array('f')
244			e_xsec_fit_t = array.array('f')
245
246			if overlay_fit:
247			X0 = Input[print_trigger][0]
248			X1 = Input[print_trigger][1]
249			X2 = Input[print_trigger][2]
250			Chi2 = Input[print_trigger][3]
251
252			for iterator in range(len(Rates[print_trigger]["rate"])):
253			if not Rates[print_trigger]["run"][iterator] in run_list:
254			continue
255			prediction = meanxsec + slopexsec * (Rates[print_trigger]["live_lumi"][iterator] - meanlumi)
256			realvalue = Rates[print_trigger]["xsec"][iterator]
257			if not data_clean or (realvalue > 0.4prediction and realvalue < 2.5prediction) or (realvalue > 0.4meanxsec and realvalue < 2.5meanxsec) or prediction < 0:
258			run_t.append(Rates[print_trigger]["run"][iterator])
259			ls_t.append(Rates[print_trigger]["ls"][iterator])
260			ps_t.append(Rates[print_trigger]["ps"][iterator])
261			inst_t.append(Rates[print_trigger]["inst_lumi"][iterator])
262			live_t.append(Rates[print_trigger]["live_lumi"][iterator])
263			delivered_t.append(Rates[print_trigger]["delivered_lumi"][iterator])
264			deadtime_t.append(Rates[print_trigger]["deadtime"][iterator])
265			rawrate_t.append(Rates[print_trigger]["rawrate"][iterator])
266			rate_t.append(Rates[print_trigger]["rate"][iterator])
267			rawxsec_t.append(Rates[print_trigger]["rawxsec"][iterator])
268			xsec_t.append(Rates[print_trigger]["xsec"][iterator])
269
270			e_run_t.append(0.0)
271			e_ls_t.append(0.0)
272			e_ps_t.append(0.0)
273			e_inst_t.append(14.14)
274			e_live_t.append(14.14)
275			e_delivered_t.append(14.14)
276			e_deadtime_t.append(1.0)
277			e_rawrate_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
278			e_rate_t.append(Rates[print_trigger]["ps"][iterator]math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls23.3)))
279			e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["delivered_lumi"][iterator])
280			e_xsec_t.append(Rates[print_trigger]["ps"][iterator]math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls23.3))/Rates[print_trigger]["live_lumi"][iterator])
281
282			if overlay_fit:
283			rate_prediction = X0 + X1live_t[-1] + X2live_t[-1]*live_t[-1]
284			rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
285			rate_fit_t.append(rate_prediction)
286			rawxsec_fit_t.append(rate_prediction(1.0-deadtime_t[-1])/(ps_t[-1]delivered_t[-1]))
287			xsec_fit_t.append(rate_prediction/delivered_t[-1])
288			e_rawrate_fit_t.append(e_rawrate_t[-1]*sqrt(Chi2))
289			e_rate_fit_t.append(e_rate_t[-1]*sqrt(Chi2))
290			e_rawxsec_fit_t.append(e_rawxsec_t[-1]*sqrt(Chi2))
291			e_xsec_fit_t.append(e_xsec_t[-1]*sqrt(Chi2))
292
293			else:
294			if debug_print:
295			print str(print_trigger)+" has xsec "+str(round(Rates[print_trigger]["xsec"][iterator],6))+" at lumi "+str(round(Rates[print_trigger]["live_lumi"][iterator],2))+" where the expected value is "+str(prediction)
296
297			for varX, varY, do_fit, save_root, save_png, overlay_fit, fit_file in plot_properties:
298			if varX == "run":
299			VX = run_t
300			VXE = run_t_e
301			elif varX == "ls":
302			VX = ls_t
303			VXE = e_ls_t
304			elif varX == "ps":
305			VX = ps_t
306			VXE = e_ps_t
307			elif varX == "inst":
308			VX = inst_t
309			VXE = e_inst_t
310			elif varX == "live":
311			VX = live_t
312			VXE = e_live_t
313			elif varX == "delivered":
314			VX = delivered_t
315			VXE = e_delivered_t
316			elif varX == "rawrate":
317			VX = rawrate_t
318			VXE = e_rawrate_t
319			elif varX == "rate":
320			VX = rate_t
321			VXE = e_rate_t
322			elif varX == "rawxsec":
323			VX = rawxsec_t
324			VXE = e_rawxsec_t
325			elif varX == "xsec":
326			VX = xsec_t
327			VXE = e_xsec_t
328			else:
329			print "No valid variable entered for X"
330			continue
331			if varY == "run":
332			VY = run_t
333			VYE = run_t_e
334			elif varY == "ls":
335			VY = ls_t
336			VYE = e_ls_t
337			elif varY == "ps":
338			VY = ps_t
339			VYE = e_ps_t
340			elif varY == "inst":
341			VY = inst_t
342			VYE = e_inst_t
343			elif varY == "live":
344			VY = live_t
345			VYE = e_live_t
346			elif varY == "delivered":
347			VY = delivered_t
348			VYE = e_delivered_t
349			elif varY == "rawrate":
350			VY = rawrate_t
351			VYE = e_rawrate_t
352			if overlay_fit:
353			VF = rawrate_fit_t
354			VFE = e_rawrate_fit_t
355			elif varY == "rate":
356			VY = rate_t
357			VYE = e_rate_t
358			if overlay_fit:
359			VF = rate_fit_t
360			VFE = e_rate_fit_t
361			elif varY == "rawxsec":
362			VY = rawxsec_t
363			VYE = e_rawxsec_t
364			if overlay_fit:
365			VF = rawxsec_fit_t
366			VFE = e_rawxsec_fit_t
367			elif varY == "xsec":
368			VY = xsec_t
369			VYE = e_xsec_t
370			if overlay_fit:
371			VF = xsec_fit_t
372			VFE = e_xsec_fit_t
373			else:
374			print "No valid variable entered for Y"
375			continue
376
377
378			if save_root or save_png:
379			c1 = TCanvas(str(varX),str(varY))
380			c1.SetName(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
381
382			gr1 = TGraphErrors(len(VX), VX, VY, VXE, VYE)
383			gr1.SetName("Graph_"+str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
384			gr1.GetXaxis().SetTitle(varX)
385			gr1.GetYaxis().SetTitle(varY)
386			gr1.SetTitle(str(print_trigger))
387			gr1.SetMinimum(0)
388			gr1.SetMaximum(1.2*max(VY))
389			gr1.GetXaxis().SetLimits(min(VX)-0.2max(VX),1.2max(VX))
390			gr1.SetMarkerStyle(8)
391			if overlay_fit:
392			gr1.SetMarkerSize(0.8)
393			else:
394			gr1.SetMarkerSize(0.5)
395			gr1.SetMarkerColor(2)
396
397			## gr2 = TGraphErrors(len(VX), live_t, xsec_t, e_live_t, e_xsec_t)
398
399			if overlay_fit:
400			gr3 = TGraphErrors(len(VX), VX, VF, VXE, VFE)
401			gr3.SetMarkerStyle(8)
402			gr3.SetMarkerSize(0.4)
403			gr3.SetMarkerColor(4)
404			gr3.SetFillColor(4)
405			gr3.SetFillStyle(3003)
406
407			if do_fit:
408			if "rate" in varY:
409			## f2a = TF1("f2a","pol1",0,8000)
410			## f2a.SetParLimits(0,0,meanxsec*1.5)
411			## if slopexsec > 0:
412			## f2a.SetParLimits(1,0,max(xsec_t)/max(live_t))
413			## else:
414			## f2a.SetParLimits(1,2slopexsec,-2slopexsec)
415			## gr2.Fit("f2a","Q","rob=0.80")
416
417			f1a = TF1("f1a","pol2",0,8000)
418			f1a.SetLineColor(4)
419			f1a.SetLineWidth(2)
420			f1a.SetParLimits(0,0,1000)
421			f1a.SetParLimits(1,0,1000)
422			#gr1.Fit("f1a","B","Q")
423			gr1.Fit("f1a","Q","rob=0.80")
424
425			## f1b = TF1("f1b","pol2",0,8000)
426			## f1b.SetParameters(0.0,f2a.GetParameter(0),f2a.GetParameter(1))
427			## f1b.SetLineColor(3)
428			## f1b.SetLineWidth(2)
429
430			else:
431			f1a = TF1("f1a","pol1",0,8000)
432			f1a.SetLineColor(4)
433			f1a.SetLineWidth(2)
434			if "xsec" in varY:
435			f1a.SetParLimits(0,0,meanxsec*1.5)
436			if slopexsec > 0:
437			f1a.SetParLimits(1,0,max(VY)/max(VX))
438			else:
439			f1a.SetParLimits(1,2slopexsec,-2slopexsec)
440			else:
441			f1a.SetParLimits(0,0,1000)
442			gr1.Fit("f1a","Q","rob=0.80")
443
444			if save_root or save_png:
445			gr1.Draw("APZ")
446			if overlay_fit:
447			gr3.Draw("P3")
448			if do_fit:
449			f1a.Draw("same")
450			## if "rate" in varY:
451			## f1b.Draw("same")
452			c1.Update()
453			if save_root:
454			myfile = TFile( RootFile, 'UPDATE' )
455			c1.Write()
456			myfile.Close()
457			if save_png:
458			c1.SaveAs(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX)+".png")
459
460
461			##p1 = TPaveStats()
462			##p1 = gr1.GetListOfFunctions().FindObject("stats")
463			##print p1
464			##gr1.PaintStats(f1b).Draw("same")
465
466			if print_table or save_fits:
467			Output[print_trigger] = [f1a.GetParameter(0), f1a.GetParameter(1), f1a.GetParameter(2), f1a.GetChisquare()/f1a.GetNDF(), f1a.GetParameter(0)+5000(f1a.GetParameter(1)+f1a.GetParameter(2)5000), meanrawrate]
468
469			if save_root:
470			print "Output root file is "+str(RootFile)
471
472			if save_fits:
473			FitNameTemplate = "Fits/2011/Fit_%s_%sLS_Run%sto%s.pkl"
474			FitFile = FitNameTemplate % (trig_name, num_ls, min_run, max_run)
475			if os.path.exists(FitFile):
476			os.remove(FitFile)
477			FitOutput = open(FitFile, 'wb')
478			pickle.dump(Output, FitOutput, 2)
479			FitOutput.close()
480
481			if print_table:
482			print '%60s%10s%10s%10s%10s%10s%10s' % ("Trig", "p0", "p1", "p2", "Chi2", "5e33 pred", "Av raw")
483			for print_trigger in Output:
484			_trigger = (print_trigger[:56] + '...') if len(print_trigger) > 59 else print_trigger
485			try:
486			print '%60s%10s%10s%10s%10s%10s%10s' % (_trigger, round(Output[print_trigger][0],3), round(Output[print_trigger][1],6)1000, round(Output[print_trigger][2],9)1000000, round(Output[print_trigger][3],2), round(Output[print_trigger][4],2) , round(Output[print_trigger][5],3))
487			except:
488			print str(print_trigger)+" is somehow broken"
489
490			if __name__=='__main__':
491			main()