ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.5
Committed: Wed Mar 7 12:21:26 2012 UTC (13 years, 1 month ago) by abrinke1
Content type: text/x-python
Branch: MAIN
CVS Tags: V00-00-09
Changes since 1.4: +428 -302 lines
Log Message: 
Many changes

File Contents

# User Rev Content
1 grchrist 1.3 #!/usr/bin/env python
2    
3 abrinke1 1.1 from DatabaseParser import *
4 abrinke1 1.5 from GetListOfRuns import *
5 abrinke1 1.1 import sys
6     import os
7     from numpy import *
8     import pickle
9    
10     from ROOT import gROOT, TCanvas, TF1, TGraph, TGraphErrors, TPaveStats, gPad, gStyle
11     from ROOT import TFile, TPaveText
12     from ROOT import gBenchmark
13     import array
14     import math
15    
16 abrinke1 1.5 from selectionParser import selectionParser
17    
18 abrinke1 1.1 def main():
19    
20 abrinke1 1.5 ## Can use any combination of LowestRunNumber, HighestRunNumber, and NumberOfRuns -
21     ## just modify "ExistingRuns.sort" and for "run in ExistingRuns" accordingly
22    
23     LowestRunNumber = 176000
24     HighestRunNumber = 180252
25     NumberOfRuns = 1
26    
27     run_list = []
28     ExistingRuns = GetLatestRunNumber(LowestRunNumber)
29     #ExistingRuns.sort(reverse = True) ##Allows you to count down from last run
30    
31     for run in ExistingRuns:
32     #if NumberOfRuns > 0:
33     if run <= HighestRunNumber:
34     run_list.append(run)
35     NumberOfRuns-=1
36    
37 abrinke1 1.1 ######## TO CREATE FITS #########
38 abrinke1 1.5 ## #run_list = [179497,179547,179558,179563,179889,179959,179977,180072,180076,180093,180241,180250,180252]
39    
40 abrinke1 1.1 ## trig_name = "HLT"
41 abrinke1 1.5 ## num_ls = 10
42     ## physics_active_psi = True ##Requires that physics and active be on, and that the prescale column is not 0
43     ## #JSON = [] ##To not use a JSON file, just leave the array empty
44     ## JSON = GetJSON("Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt") ##Returns array JSON[runs][ls_list]
45    
46 abrinke1 1.1 ## debug_print = False
47    
48 abrinke1 1.5 ## min_rate = 0.1
49 abrinke1 1.1 ## print_table = False
50 abrinke1 1.5 ## data_clean = True ##Gets rid of anomalous rate points, reqires physics_active_psi (PAP) and deadtime < 20%
51     ## ##plot_properties = [varX, varY, do_fit, save_root, save_png, fit_file]
52     ## plot_properties = [["delivered", "rate", True, True, False, ""]]
53    
54 abrinke1 1.1 ## masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
55 abrinke1 1.5 ## save_fits = False
56 abrinke1 1.1
57    
58 abrinke1 1.5 ## ######## TO SEE RATE VS PREDICTION ########
59     run_list = [180241]
60 abrinke1 1.1
61 abrinke1 1.5 trig_name = "IsoMu"
62 abrinke1 1.1 num_ls = 1
63 abrinke1 1.5 physics_active_psi = True
64     JSON = []
65 abrinke1 1.1 debug_print = False
66    
67 grchrist 1.3 min_rate = 10.0
68 abrinke1 1.1 print_table = False
69     data_clean = False
70 abrinke1 1.5 ##plot_properties = [varX, varY, do_fit, save_root, save_png, fit_file]
71     plot_properties = [["ls", "xsec", False, True, False, "Fits/2011/Fit_HLT_10LS_Run176023to180252.pkl"]]
72 abrinke1 1.1 masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
73     save_fits = False
74    
75 abrinke1 1.2
76 abrinke1 1.5 ######## END PARAMETERS - CALL FUNCTIONS ##########
77     Rates = GetDBRates(run_list, trig_name, num_ls, physics_active_psi, JSON, debug_print)
78 abrinke1 1.1 MakePlots(Rates, run_list, trig_name, num_ls, min_rate, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print)
79    
80    
81 abrinke1 1.5 def GetDBRates(run_list,trig_name,num_ls,physics_active_psi,JSON,debug_print):
82 abrinke1 1.1
83     Rates = {}
84 abrinke1 1.5 ## Save in RefRuns with name dependent on trig_name, num_ls, JSON, and physics_active_psi
85     if JSON:
86     if physics_active_psi:
87     RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_JPAP.pkl"
88     else:
89     RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_JSON.pkl"
90     else:
91     if physics_active_psi:
92     RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_PAP.pkl"
93     else:
94     RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS.pkl"
95    
96 abrinke1 1.1 RefRunFile = RefRunNameTemplate % (trig_name,num_ls)
97 abrinke1 1.5 RefRunFileHLT = RefRunNameTemplate % ("HLT",num_ls)
98     try: ##Open an existing RefRun file with the same parameters and trigger name
99 abrinke1 1.1 pkl_file = open(RefRunFile, 'rb')
100     Rates = pickle.load(pkl_file)
101     pkl_file.close()
102     os.remove(RefRunFile)
103     except:
104 abrinke1 1.5 try: ##Open an existing RefRun file with the same parameters and HLT for trigger name
105     pkl_file = open(RefRunFileHLT)
106     HLTRates = pickle.load(pkl_file)
107     for key in HLTRates:
108     if trig_name in str(key):
109     Rates[key] = HLTRates[key]
110     print str(RefRunFile)+" does not exist. Creating ..."
111     except:
112     print str(RefRunFile)+" does not exist. Creating ..."
113    
114     for RefRunNum in run_list:
115 abrinke1 1.1
116 abrinke1 1.5 if JSON:
117     if not RefRunNum in JSON:
118     continue
119 abrinke1 1.1
120     try:
121     ExistsAlready = False
122     for key in Rates:
123     if RefRunNum in Rates[key]["run"]:
124     ExistsAlready = True
125     break
126     if ExistsAlready:
127     continue
128     except:
129     print "Getting info for run "+str(RefRunNum)
130    
131     if RefRunNum < 1:
132     continue
133    
134 abrinke1 1.5 if True: ##Placeholder
135 abrinke1 1.1 if True: #May replace with "try" - for now it's good to know when problems happen
136     RefParser = DatabaseParser()
137     RefParser.RunNumber = RefRunNum
138     RefParser.ParseRunSetup()
139 abrinke1 1.5 RefLumiRangePhysicsActive = RefParser.GetLSRange(1,9999) ##Gets array of all LS with physics and active on
140     RefLumiArray = RefParser.GetLumiInfo() ##Gets array of all existing LS and their lumi info
141 abrinke1 1.2 RefLumiRange = []
142 abrinke1 1.1
143 abrinke1 1.5 for iterator in RefLumiArray[0]: ##Makes array of LS with proper PAP and JSON properties
144     if not physics_active_psi or (RefLumiArray[5][iterator] == 1 and RefLumiArray[6][iterator] == 1 and RefLumiArray[0][iterator] > 0):
145     if not JSON or RefRunNum in JSON:
146     if not JSON or iterator in JSON[RefRunNum]:
147     RefLumiRange.append(iterator)
148    
149     try:
150     nls = RefLumiRange[0]
151     LSRange = {}
152     except:
153     print "Run "+str(RefRunNum)+" has no good LS"
154     continue
155     if num_ls > len(RefLumiRange):
156 abrinke1 1.1 print "Run "+str(RefRunNum)+" is too short: from "+str(nls)+" to "+str(RefLumiRange[-1])+", while num_ls = "+str(num_ls)
157     continue
158     while nls < RefLumiRange[-1]-num_ls:
159 abrinke1 1.5 LSRange[nls] = []
160     counter = 0
161     for iterator in RefLumiRange:
162     if iterator >= nls and counter < num_ls:
163     LSRange[nls].append(iterator)
164     counter += 1
165 abrinke1 1.1 nls = LSRange[nls][-1]+1
166 abrinke1 1.5
167 abrinke1 1.1 print "Run "+str(RefRunNum)+" contains LS from "+str(min(LSRange))+" to "+str(max(LSRange))
168     for nls in LSRange:
169     TriggerRates = RefParser.GetHLTRates(LSRange[nls])
170 abrinke1 1.5
171     [inst, live, delivered, dead, pscols] = RefParser.GetAvLumiInfo(LSRange[nls])
172    
173 abrinke1 1.2 physics = 1
174     active = 1
175 abrinke1 1.5 psi = 99
176     for iterator in LSRange[nls]: ##Gets lowest value of physics, active, and psi in the set of lumisections
177 abrinke1 1.2 if RefLumiArray[5][iterator] == 0:
178     physics = 0
179     if RefLumiArray[6][iterator] == 0:
180     active = 0
181 abrinke1 1.5 if RefLumiArray[0][iterator] < psi:
182     psi = RefLumiArray[0][iterator]
183 abrinke1 1.2
184 abrinke1 1.5 if inst < 0 or live < 0 or delivered < 0:
185     print "Run "+str(RefRunNum)+" LS "+str(nls)+" inst lumi = "+str(inst)+" live lumi = "+str(live)+", delivered = "+str(delivered)+", physics = "+str(physics)+", active = "+str(active)
186 abrinke1 1.2
187 abrinke1 1.1 for key in TriggerRates:
188     if not trig_name in key:
189     continue
190     name = key
191 abrinke1 1.5 if re.match('.*_v[0-9]+',name): ##Removes _v#
192 abrinke1 1.1 name = name[:name.rfind('_')]
193    
194     if not Rates.has_key(name):
195     Rates[name] = {}
196     Rates[name]["run"] = []
197     Rates[name]["ls"] = []
198     Rates[name]["ps"] = []
199     Rates[name]["inst_lumi"] = []
200     Rates[name]["live_lumi"] = []
201     Rates[name]["delivered_lumi"] = []
202     Rates[name]["deadtime"] = []
203     Rates[name]["rawrate"] = []
204     Rates[name]["rate"] = []
205     Rates[name]["rawxsec"] = []
206     Rates[name]["xsec"] = []
207 abrinke1 1.2 Rates[name]["physics"] = []
208     Rates[name]["active"] = []
209 abrinke1 1.5 Rates[name]["psi"] = []
210 abrinke1 1.1 [avps, ps, rate, psrate] = TriggerRates[key]
211     Rates[name]["run"].append(RefRunNum)
212     Rates[name]["ls"].append(nls)
213     Rates[name]["ps"].append(ps)
214     Rates[name]["inst_lumi"].append(inst)
215     Rates[name]["live_lumi"].append(live)
216     Rates[name]["delivered_lumi"].append(delivered)
217     Rates[name]["deadtime"].append(dead)
218     Rates[name]["rawrate"].append(rate)
219 abrinke1 1.2 if live == 0:
220 abrinke1 1.5 Rates[name]["rate"].append(0)
221 abrinke1 1.2 Rates[name]["rawxsec"].append(0.0)
222     Rates[name]["xsec"].append(0.0)
223     else:
224 abrinke1 1.5 Rates[name]["rate"].append(psrate/(1.0-dead))
225 abrinke1 1.2 Rates[name]["rawxsec"].append(rate/live)
226     Rates[name]["xsec"].append(psrate/live)
227     Rates[name]["physics"].append(physics)
228     Rates[name]["active"].append(active)
229 abrinke1 1.5 Rates[name]["psi"].append(psi)
230     #except: ##If we replace "if True:" with "try:"
231 abrinke1 1.1 #print "Failed to parse run "+str(RefRunNum)
232    
233 abrinke1 1.5 RateOutput = open(RefRunFile, 'wb') ##Save new Rates[] to RefRuns
234 abrinke1 1.1 pickle.dump(Rates, RateOutput, 2)
235     RateOutput.close()
236     return Rates
237    
238     def MakePlots(Rates, run_list, trig_name, num_ls, min_rate, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print):
239     min_run = min(run_list)
240     max_run = max(run_list)
241    
242 abrinke1 1.5 InputFit = {}
243     OutputFit = {}
244 abrinke1 1.1
245     RootNameTemplate = "%s_%sLS_Run%sto%s.root"
246     RootFile = RootNameTemplate % (trig_name, num_ls, min_run, max_run)
247    
248 abrinke1 1.5 for varX, varY, do_fit, save_root, save_png, fit_file in plot_properties:
249     if fit_file:
250 abrinke1 1.1 pkl_file = open(fit_file, 'rb')
251 abrinke1 1.5 InputFit = pickle.load(pkl_file)
252 abrinke1 1.1 pkl_file.close()
253     if save_root:
254     try:
255     os.remove(RootFile)
256     except:
257     break
258    
259     for print_trigger in Rates:
260     meanrawrate = sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])
261     if not trig_name in print_trigger:
262     continue
263     if meanrawrate < min_rate:
264     continue
265     masked_trig = False
266     for mask in masked_triggers:
267     if str(mask) in print_trigger:
268     masked_trig = True
269     if masked_trig:
270     continue
271    
272 abrinke1 1.5 OutputFit[print_trigger] = {}
273 abrinke1 1.1
274     lowlumi = 0
275 abrinke1 1.2 numzeroes = 0
276     for live_lumi in Rates[print_trigger]["live_lumi"]:
277     if live_lumi < 1:
278     numzeroes+=1
279     meanlumi_init = sum(Rates[print_trigger]["live_lumi"])/(len(Rates[print_trigger]["live_lumi"])-numzeroes)
280 abrinke1 1.1 meanlumi = 0
281     highlumi = 0
282     lowxsec = 0
283     meanxsec = 0
284     highxsec = 0
285     nlow = 0
286     nhigh = 0
287     for iterator in range(len(Rates[print_trigger]["rate"])):
288     if not Rates[print_trigger]["run"][iterator] in run_list:
289     continue
290     if Rates[print_trigger]["live_lumi"][iterator] <= meanlumi_init:
291 abrinke1 1.5 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 and Rates[print_trigger]["deadtime"][iterator] < 0.20 and Rates[print_trigger]["psi"][iterator] > 0):
292 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
293     lowxsec+=Rates[print_trigger]["xsec"][iterator]
294     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
295     lowlumi+=Rates[print_trigger]["live_lumi"][iterator]
296     nlow+=1
297     if Rates[print_trigger]["live_lumi"][iterator] > meanlumi_init:
298 abrinke1 1.5 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 and Rates[print_trigger]["deadtime"][iterator] < 0.20 and Rates[print_trigger]["psi"][iterator] > 0):
299 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
300     highxsec+=Rates[print_trigger]["xsec"][iterator]
301     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
302     highlumi+=Rates[print_trigger]["live_lumi"][iterator]
303     nhigh+=1
304     meanxsec = meanxsec/(nlow+nhigh)
305     meanlumi = meanlumi/(nlow+nhigh)
306     slopexsec = ( (highxsec/nhigh) - (lowxsec/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
307    
308 abrinke1 1.5 [run_t,ls_t,ps_t,inst_t,live_t,delivered_t,deadtime_t,rawrate_t,rate_t,rawxsec_t,xsec_t,psi_t,e_run_t,e_ls_t,e_ps_t,e_inst_t,e_live_t,e_delivered_t,e_deadtime_t,e_rawrate_t,e_rate_t,e_rawxsec_t,e_xsec_t,e_psi_t,rawrate_fit_t,rate_fit_t,rawxsec_fit_t,xsec_fit_t,e_rawrate_fit_t,e_rate_fit_t,e_rawxsec_fit_t,e_xsec_fit_t] = MakePlotArrays()
309    
310     if fit_file:
311     FitType = InputFit[print_trigger][0]
312     X0 = InputFit[print_trigger][1]
313     X1 = InputFit[print_trigger][2]
314     X2 = InputFit[print_trigger][3]
315     X3 = InputFit[print_trigger][4]
316     Chi2 = InputFit[print_trigger][5]
317    
318 abrinke1 1.1 for iterator in range(len(Rates[print_trigger]["rate"])):
319     if not Rates[print_trigger]["run"][iterator] in run_list:
320     continue
321     prediction = meanxsec + slopexsec * (Rates[print_trigger]["live_lumi"][iterator] - meanlumi)
322     realvalue = Rates[print_trigger]["xsec"][iterator]
323 abrinke1 1.5 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 and Rates[print_trigger]["deadtime"][iterator] < 0.20 and Rates[print_trigger]["psi"][iterator] > 0):
324 abrinke1 1.1 run_t.append(Rates[print_trigger]["run"][iterator])
325     ls_t.append(Rates[print_trigger]["ls"][iterator])
326     ps_t.append(Rates[print_trigger]["ps"][iterator])
327     inst_t.append(Rates[print_trigger]["inst_lumi"][iterator])
328     live_t.append(Rates[print_trigger]["live_lumi"][iterator])
329     delivered_t.append(Rates[print_trigger]["delivered_lumi"][iterator])
330     deadtime_t.append(Rates[print_trigger]["deadtime"][iterator])
331     rawrate_t.append(Rates[print_trigger]["rawrate"][iterator])
332     rate_t.append(Rates[print_trigger]["rate"][iterator])
333     rawxsec_t.append(Rates[print_trigger]["rawxsec"][iterator])
334     xsec_t.append(Rates[print_trigger]["xsec"][iterator])
335 abrinke1 1.5 psi_t.append(Rates[print_trigger]["psi"][iterator])
336 abrinke1 1.1
337     e_run_t.append(0.0)
338     e_ls_t.append(0.0)
339     e_ps_t.append(0.0)
340     e_inst_t.append(14.14)
341     e_live_t.append(14.14)
342     e_delivered_t.append(14.14)
343 abrinke1 1.2 e_deadtime_t.append(0.01)
344 abrinke1 1.1 e_rawrate_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
345     e_rate_t.append(Rates[print_trigger]["ps"][iterator]*math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
346 abrinke1 1.5 e_psi_t.append(0.0)
347 abrinke1 1.2 if live_t[-1] == 0:
348     e_rawxsec_t.append(0)
349     e_xsec_t.append(0)
350     else:
351     e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["live_lumi"][iterator])
352     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])
353 abrinke1 1.5
354     if fit_file:
355     if FitType == "expo":
356     rate_prediction = X0 + X1*math.exp(X2*delivered_t[-1])
357     else:
358     rate_prediction = X0 + X1*delivered_t[-1] + X2*delivered_t[-1]*delivered_t[-1] + X3*delivered_t[-1]*delivered_t[-1]*delivered_t[-1]
359     ## if rate_t[-1] < 0.7 * rate_prediction or rate_t[-1] > 1.4 * rate_prediction:
360     ## 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])
361    
362 abrinke1 1.2 if live_t[-1] == 0:
363 abrinke1 1.5 rawrate_fit_t.append(0)
364     rate_fit_t.append(0)
365 abrinke1 1.2 rawxsec_fit_t.append(0)
366     xsec_fit_t.append(0)
367     else:
368 abrinke1 1.5 rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
369     rate_fit_t.append(rate_prediction)
370     rawxsec_fit_t.append(rawrate_fit_t[-1]/live_t[-1])
371     xsec_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/live_t[-1])
372     e_rawrate_fit_t.append(math.sqrt(Chi2)*rawrate_fit_t[-1]/rate_fit_t[-1])
373     e_rate_fit_t.append(math.sqrt(Chi2))
374     e_rawxsec_fit_t.append(math.sqrt(Chi2)*rawxsec_fit_t[-1]/rate_fit_t[-1])
375     e_xsec_fit_t.append(math.sqrt(Chi2)*xsec_fit_t[-1]/rate_fit_t[-1])
376    
377     else: ##If the data point does not pass the data_clean filter
378 abrinke1 1.1 if debug_print:
379     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)
380    
381 abrinke1 1.5 ## End "for iterator in range(len(Rates[print_trigger]["rate"])):" loop
382 abrinke1 1.1
383 abrinke1 1.5 AllPlotArrays = [run_t,ls_t,ps_t,inst_t,live_t,delivered_t,deadtime_t,rawrate_t,rate_t,rawxsec_t,xsec_t,psi_t,e_run_t,e_ls_t,e_ps_t,e_inst_t,e_live_t,e_delivered_t,e_deadtime_t,e_rawrate_t,e_rate_t,e_rawxsec_t,e_xsec_t,e_psi_t,rawrate_fit_t,rate_fit_t,rawxsec_fit_t,xsec_fit_t,e_rawrate_fit_t,e_rate_fit_t,e_rawxsec_fit_t,e_xsec_fit_t]
384     [VX, VXE, VY, VYE, VF, VFE] = GetVXVY(plot_properties, fit_file, AllPlotArrays)
385    
386 abrinke1 1.1
387 abrinke1 1.5 if save_root or save_png:
388     c1 = TCanvas(str(varX),str(varY))
389     c1.SetName(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
390    
391     gr1 = TGraphErrors(len(VX), VX, VY, VXE, VYE)
392     gr1.SetName("Graph_"+str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
393     gr1.GetXaxis().SetTitle(varX)
394     gr1.GetYaxis().SetTitle(varY)
395     gr1.SetTitle(str(print_trigger))
396     gr1.SetMinimum(0)
397     gr1.SetMaximum(1.2*max(VY))
398     #gr1.GetXaxis().SetLimits(min(VX)-0.2*max(VX),1.2*max(VX))
399     gr1.GetXaxis().SetLimits(0,1.2*max(VX))
400     gr1.SetMarkerStyle(8)
401     if fit_file:
402     gr1.SetMarkerSize(0.8)
403     else:
404     gr1.SetMarkerSize(0.5)
405     gr1.SetMarkerColor(2)
406    
407     if fit_file:
408     gr3 = TGraphErrors(len(VX), VX, VF, VXE, VFE)
409     gr3.SetMarkerStyle(8)
410     gr3.SetMarkerSize(0.4)
411     gr3.SetMarkerColor(4)
412     gr3.SetFillColor(4)
413     gr3.SetFillStyle(3003)
414 abrinke1 1.1
415 abrinke1 1.5 if do_fit:
416     if "rate" in varY:
417     f1a = TF1("f1a","pol2",0,8000)
418     f1a.SetLineColor(4)
419     f1a.SetLineWidth(2)
420     f1a.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
421     f1a.SetParLimits(1,0,2.0*max(VY)/(max(VX)*max(VX)))
422     #gr1.Fit("f1a","B","Q")
423     gr1.Fit("f1a","Q","rob=0.90")
424    
425     f1b = 0
426     f1c = 0
427     if True:
428     f1b = TF1("f1b","pol3",0,8000)
429     f1b.SetLineColor(2)
430     f1b.SetLineWidth(2)
431     f1b.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
432     f1b.SetParLimits(1,0,f1a.GetParameter(1)+0.0000001)
433     f1b.SetParLimits(2,0,f1a.GetParameter(2)+0.0000000001)
434     f1b.SetParLimits(3,0,2.0*max(VY)/(max(VX)*max(VX)*max(VX)))
435     gr1.Fit("f1b","Q","rob=0.90")
436     #if f1b.GetChisquare()/f1b.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
437     print "X0 = "+str(f1a.GetParameter(0))+" X1 = "+str(f1a.GetParameter(1))+" X2 = "+str(f1a.GetParameter(2))
438     print str(print_trigger)+" f1a Chi2 = "+str(10*f1a.GetChisquare()*math.sqrt(len(VY))/(math.sqrt(sum(VY))*num_ls*f1a.GetNDF()))+", f1b Chi2 = "+str(10*f1b.GetChisquare()*math.sqrt(len(VY))/(math.sqrt(sum(VY))*num_ls*f1b.GetNDF()))
439     print "X0 = "+str(f1b.GetParameter(0))+" X1 = "+str(f1b.GetParameter(1))+" X2 = "+str(f1b.GetParameter(2))+" X3 = "+str(f1b.GetParameter(3))
440 abrinke1 1.1
441 abrinke1 1.5 f1c = TF1("f1c","[0]+[1]*expo(2)",0,8000)
442     f1c.SetLineColor(3)
443     f1c.SetLineWidth(2)
444     f1c.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
445     f1c.SetParLimits(1,max(VY)/math.exp(10.0),max(VY)/math.exp(2.0))
446     f1c.SetParLimits(2,0.0,0.0000000001)
447     f1c.SetParLimits(3,2.0/max(VX),10.0/max(VX))
448     print str(max(VY)/math.exp(2.0))+" "+str(10.0/max(VX))
449     gr1.Fit("f1c","Q","rob=0.90")
450     #if f1c.GetChisquare()/f1c.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
451     print str(print_trigger)+" f1a Chi2 = "+str(10*f1a.GetChisquare()*math.sqrt(len(VY))/(math.sqrt(sum(VY))*num_ls*f1a.GetNDF()))+", f1c Chi2 = "+str(10*f1c.GetChisquare()*math.sqrt(len(VY))/(math.sqrt(sum(VY))*num_ls*f1c.GetNDF()))
452     print "X0 = "+str(f1c.GetParameter(0))+" X1 = "+str(f1c.GetParameter(1))+" X2 = "+str(f1c.GetParameter(2))+" X3 = "+str(f1c.GetParameter(3))
453 abrinke1 1.2
454 abrinke1 1.5 else: ##If this is not a rate plot
455     f1a = TF1("f1a","pol1",0,8000)
456     f1a.SetLineColor(4)
457     f1a.SetLineWidth(2)
458     if "xsec" in varY:
459     f1a.SetParLimits(0,0,meanxsec*1.5)
460     if slopexsec > 0:
461     f1a.SetParLimits(1,0,max(VY)/max(VX))
462     else:
463     f1a.SetParLimits(1,2*slopexsec,-2*slopexsec)
464 abrinke1 1.1 else:
465 abrinke1 1.5 f1a.SetParLimits(0,-1000,1000)
466     gr1.Fit("f1a","Q","rob=0.80")
467 abrinke1 1.1
468 abrinke1 1.5 if save_root or save_png:
469     gr1.Draw("APZ")
470     ## ##Option to draw stats box
471     ## p1 = TPaveStats()
472     ## p1 = gr1.GetListOfFunctions().FindObject("stats")
473     ## print p1
474     ## gr1.PaintStats(f1b).Draw("same")
475     if fit_file:
476     gr3.Draw("P3")
477     if do_fit:
478     f1a.Draw("same")
479     try:
480     f1b.Draw("same")
481     f1c.Draw("same")
482     except:
483     True
484     c1.Update()
485     if save_root:
486     myfile = TFile( RootFile, 'UPDATE' )
487     c1.Write()
488     myfile.Close()
489     if save_png:
490     c1.SaveAs(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX)+".png")
491 abrinke1 1.1
492    
493     if print_table or save_fits:
494 abrinke1 1.5 if not do_fit:
495     print "Can't have save_fits = True and do_fit = False"
496     continue
497     if f1c.GetChisquare()/f1c.GetNDF() < 0.95*f1a.GetChisquare()/f1a.GetNDF() and f1c.GetChisquare()/f1c.GetNDF() < 0.95*f1b.GetChisquare()/f1b.GetNDF():
498     OutputFit[print_trigger] = ["expo", f1c.GetParameter(0), f1c.GetParameter(1), f1c.GetParameter(3), 0.0, f1c.GetChisquare()/f1c.GetNDF(), meanrawrate]
499     elif f1b.GetChisquare()/f1b.GetNDF() < 0.95*f1a.GetChisquare()/f1a.GetNDF():
500     OutputFit[print_trigger] = ["poly", f1b.GetParameter(0), f1b.GetParameter(1), f1b.GetParameter(2), f1b.GetParameter(3), f1b.GetChisquare()/f1b.GetNDF(), meanrawrate]
501     else:
502     OutputFit[print_trigger] = ["poly", f1a.GetParameter(0), f1a.GetParameter(1), f1a.GetParameter(2), 0.0, f1a.GetChisquare()/f1a.GetNDF(), meanrawrate]
503 abrinke1 1.1
504     if save_root:
505     print "Output root file is "+str(RootFile)
506    
507     if save_fits:
508     FitNameTemplate = "Fits/2011/Fit_%s_%sLS_Run%sto%s.pkl"
509     FitFile = FitNameTemplate % (trig_name, num_ls, min_run, max_run)
510     if os.path.exists(FitFile):
511     os.remove(FitFile)
512 abrinke1 1.5 FitOutputFile = open(FitFile, 'wb')
513     pickle.dump(OutputFit, FitOutputFile, 2)
514     FitOutputFile.close()
515 abrinke1 1.1
516     if print_table:
517 abrinke1 1.5 print "The expo fit is of the form p0+p1*e^(p2*x), poly is p0+(p1/10^3)*x+(p2/10^6)*x^2+(p3/10^9)*x^3, where x is Deliv. Lumi."
518     print '%60s%10s%10s%10s%10s%10s%10s%10s' % ("Trig", "fit", "p0", "p1", "p2", "p3", "Chi2", "Av raw")
519     for print_trigger in OutputFit:
520     _trigger = (print_trigger[:56] + '...') if len(print_trigger) > 59 else print_trigger
521 abrinke1 1.1 try:
522 abrinke1 1.5 if OutputFit[print_trigger][0] == "poly":
523     print '%60s%10s%10s%10s%10s%10s%10s' % (_trigger, OutputFit[print_trigger][0], round(OutputFit[print_trigger][1],3), round(OutputFit[print_trigger][2],6)*1000, round(OutputFit[print_trigger][3],9)*1000000, round(OutputFit[print_trigger][4],12)*1000000000, round(OutputFit[print_trigger][5],2), round(OutputFit[print_trigger][6],3))
524     else:
525     print '%60s%10s%10s%10s%10s%10s%10s' % (_trigger, OutputFit[print_trigger][0], OutputFit[print_trigger][1], OutputFit[print_trigger][2], OutputFit[print_trigger][3], OutputFit[print_trigger][4], round(OutputFit[print_trigger][5],2), round(OutputFit[print_trigger][6],3))
526 abrinke1 1.1 except:
527     print str(print_trigger)+" is somehow broken"
528    
529 abrinke1 1.5
530     ############# SUPPORTING FUNCTIONS ################
531    
532    
533     def GetJSON(json_file):
534    
535     input_file = open(json_file)
536     file_content = input_file.read()
537     inputRange = selectionParser(file_content)
538     JSON = inputRange.runsandls()
539     return JSON
540     ##JSON is an array: JSON[run_number] = [1st ls, 2nd ls, 3rd ls ... nth ls]
541    
542     def MakePlotArrays():
543     run_t = array.array('f')
544     ls_t = array.array('f')
545     ps_t = array.array('f')
546     inst_t = array.array('f')
547     live_t = array.array('f')
548     delivered_t = array.array('f')
549     deadtime_t = array.array('f')
550     rawrate_t = array.array('f')
551     rate_t = array.array('f')
552     rawxsec_t = array.array('f')
553     xsec_t = array.array('f')
554     psi_t = array.array('f')
555    
556     e_run_t = array.array('f')
557     e_ls_t = array.array('f')
558     e_ps_t = array.array('f')
559     e_inst_t = array.array('f')
560     e_live_t = array.array('f')
561     e_delivered_t = array.array('f')
562     e_deadtime_t = array.array('f')
563     e_rawrate_t = array.array('f')
564     e_rate_t = array.array('f')
565     e_rawxsec_t = array.array('f')
566     e_xsec_t = array.array('f')
567     e_psi_t = array.array('f')
568    
569     rawrate_fit_t = array.array('f')
570     rate_fit_t = array.array('f')
571     rawxsec_fit_t = array.array('f')
572     xsec_fit_t = array.array('f')
573     e_rawrate_fit_t = array.array('f')
574     e_rate_fit_t = array.array('f')
575     e_rawxsec_fit_t = array.array('f')
576     e_xsec_fit_t = array.array('f')
577    
578     return [run_t,ls_t,ps_t,inst_t,live_t,delivered_t,deadtime_t,rawrate_t,rate_t,rawxsec_t,xsec_t,psi_t,e_run_t,e_ls_t,e_ps_t,e_inst_t,e_live_t,e_delivered_t,e_deadtime_t,e_rawrate_t,e_rate_t,e_rawxsec_t,e_xsec_t,e_psi_t,rawrate_fit_t,rate_fit_t,rawxsec_fit_t,xsec_fit_t,e_rawrate_fit_t,e_rate_fit_t,e_rawxsec_fit_t,e_xsec_fit_t]
579    
580    
581     def GetVXVY(plot_properties, fit_file, AllPlotArrays):
582    
583     VF = "0"
584     VFE = "0"
585    
586     [run_t,ls_t,ps_t,inst_t,live_t,delivered_t,deadtime_t,rawrate_t,rate_t,rawxsec_t,xsec_t,psi_t,e_run_t,e_ls_t,e_ps_t,e_inst_t,e_live_t,e_delivered_t,e_deadtime_t,e_rawrate_t,e_rate_t,e_rawxsec_t,e_xsec_t,e_psi_t,rawrate_fit_t,rate_fit_t,rawxsec_fit_t,xsec_fit_t,e_rawrate_fit_t,e_rate_fit_t,e_rawxsec_fit_t,e_xsec_fit_t] = AllPlotArrays
587     for varX, varY, do_fit, save_root, save_png, fit_file in plot_properties:
588     if varX == "run":
589     VX = run_t
590     VXE = run_t_e
591     elif varX == "ls":
592     VX = ls_t
593     VXE = e_ls_t
594     elif varX == "ps":
595     VX = ps_t
596     VXE = e_ps_t
597     elif varX == "inst":
598     VX = inst_t
599     VXE = e_inst_t
600     elif varX == "live":
601     VX = live_t
602     VXE = e_live_t
603     elif varX == "delivered":
604     VX = delivered_t
605     VXE = e_delivered_t
606     elif varX == "deadtime":
607     VX = deadtime_t
608     VXE = e_deadtime_t
609     elif varX == "rawrate":
610     VX = rawrate_t
611     VXE = e_rawrate_t
612     elif varX == "rate":
613     VX = rate_t
614     VXE = e_rate_t
615     elif varX == "rawxsec":
616     VX = rawxsec_t
617     VXE = e_rawxsec_t
618     elif varX == "xsec":
619     VX = xsec_t
620     VXE = e_xsec_t
621     elif varX == "psi":
622     VX = psi_t
623     VXE = e_psi_t
624     else:
625     print "No valid variable entered for X"
626     continue
627     if varY == "run":
628     VY = run_t
629     VYE = run_t_e
630     elif varY == "ls":
631     VY = ls_t
632     VYE = e_ls_t
633     elif varY == "ps":
634     VY = ps_t
635     VYE = e_ps_t
636     elif varY == "inst":
637     VY = inst_t
638     VYE = e_inst_t
639     elif varY == "live":
640     VY = live_t
641     VYE = e_live_t
642     elif varY == "delivered":
643     VY = delivered_t
644     VYE = e_delivered_t
645     elif varY == "deadtime":
646     VY = deadtime_t
647     VYE = e_deadtime_t
648     elif varY == "rawrate":
649     VY = rawrate_t
650     VYE = e_rawrate_t
651     if fit_file:
652     VF = rawrate_fit_t
653     VFE = e_rawrate_fit_t
654     elif varY == "rate":
655     VY = rate_t
656     VYE = e_rate_t
657     if fit_file:
658     VF = rate_fit_t
659     VFE = e_rate_fit_t
660     elif varY == "rawxsec":
661     VY = rawxsec_t
662     VYE = e_rawxsec_t
663     if fit_file:
664     VF = rawxsec_fit_t
665     VFE = e_rawxsec_fit_t
666     elif varY == "xsec":
667     VY = xsec_t
668     VYE = e_xsec_t
669     if fit_file:
670     VF = xsec_fit_t
671     VFE = e_xsec_fit_t
672     elif varY == "psi":
673     VY = psi_t
674     VYE = e_psi_t
675     else:
676     print "No valid variable entered for Y"
677     continue
678    
679     return [VX, VXE, VY, VYE, VF, VFE]
680    
681    
682 abrinke1 1.1 if __name__=='__main__':
683     main()