UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py

#!/usr/bin/env python

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 = [180250]

    trig_name = "Mu"
    num_ls = 1
    debug_print = False

    min_rate = 10.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()
                RefLumiRangePhysicsActive = RefParser.GetLSRange(1,9999)
                RefLumiArray = RefParser.GetLumiInfo()
                RefLumiRange = []
                for iterator in RefLumiArray[0]:
                    RefLumiRange.append(iterator)

                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)
                        LSRange[nls] = []
                        for iterator in RefLumiRange:
                            if iterator >= nls and iterator < nls+num_ls:
                                LSRange[nls].append(iterator)
                    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,lastpscol] = RefParser.GetAvLumiInfo(LSRange[nls])
                    physics = 1
                    active = 1
                    for iterator in LSRange[nls]:
                        if RefLumiArray[5][iterator] == 0:
                            physics = 0
                        if RefLumiArray[6][iterator] == 0:
                            active = 0

                    if live < 0:
                        print "Run "+str(RefRunNum)+" LS "+str(nls)+" live lumi = "+str(live)+", delivered = "+str(delivered)+", physics = "+str(physics)+", active = "+str(active)

                    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"] = []
                            Rates[name]["physics"] = []
                            Rates[name]["active"] = []
                        [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))
                        if live == 0:
                            Rates[name]["rawxsec"].append(0.0)
                            Rates[name]["xsec"].append(0.0)
                        else:
                            Rates[name]["rawxsec"].append(rate/live)
                            Rates[name]["xsec"].append(psrate/live)
                        Rates[name]["physics"].append(physics)
                        Rates[name]["active"].append(active)
            #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
        numzeroes = 0
        for live_lumi in Rates[print_trigger]["live_lumi"]:
            if live_lumi < 1:
                numzeroes+=1
        meanlumi_init = sum(Rates[print_trigger]["live_lumi"])/(len(Rates[print_trigger]["live_lumi"])-numzeroes)
        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 and Rates[print_trigger]["physics"][iterator] == 1 and Rates[print_trigger]["active"][iterator] == 1):
                    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 and Rates[print_trigger]["physics"][iterator] == 1 and Rates[print_trigger]["active"][iterator] == 1):
                    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 ) and Rates[print_trigger]["physics"][iterator] == 1 and Rates[print_trigger]["active"][iterator] == 1 ):
                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(0.01)
                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)))
                
                if live_t[-1] == 0:
                    e_rawxsec_t.append(0)
                    e_xsec_t.append(0)
                else:
                    e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["live_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*delivered_t[-1] + X2*delivered_t[-1]*delivered_t[-1]
                    if rate_t[-1] < 0.7 * rate_prediction or rate_t[-1] > 1.4 * rate_prediction:
                        pass
                        #print str(run_t[-1])+"  "+str(ls_t[-1])+"  "+str(print_trigger)+"  "+str(ps_t[-1])+"  "+str(deadtime_t[-1])+"  "+str(rate_prediction)+"  "+str(rate_t[-1])+"  "+str(rawrate_t[-1])
                    rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
                    rate_fit_t.append(rate_prediction)
                    if live_t[-1] == 0:
                        rawxsec_fit_t.append(0)
                        xsec_fit_t.append(0)
                    else:
                        rawxsec_fit_t.append(rate_prediction/(ps_t[-1]*live_t[-1]))
                        xsec_fit_t.append(rate_prediction/live_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 == "deadtime":
                VX = deadtime_t
                VXE = e_deadtime_t
                print deadtime_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 == "deadtime":
                VY = deadtime_t
                VYE = e_deadtime_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")

                    if f1a.GetChisquare()/f1a.GetNDF() > 20:
                        f1b = TF1("f1b","pol3",0,8000)
                        f1b.SetLineColor(4)
                        f1b.SetLineWidth(2)
                        f1b.SetParLimits(0,0,1000)
                        f1b.SetParLimits(1,0,1000)
                        f1b.SetParLimits(2,0,1000)
                        gr1.Fit("f1b","Q","rob=0.80")
                        if f1b.GetChisquare()/f1b.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
                            print str(print_trigger)+" f1a Chi2 = "+str(f1a.GetChisquare()/f1a.GetNDF())+", f1b Chi2 = "+str(f1b.GetChisquare()/f1b.GetNDF())
                        

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