ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.7
Committed: Wed Mar 7 21:56:46 2012 UTC (13 years, 1 month ago) by abrinke1
Content type: text/x-python
Branch: MAIN
CVS Tags: V00-00-12, V00-00-11
Changes since 1.6: +42 -27 lines
Log Message: 
Fixes for divide by zero, other tweaks

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.7 ## ###### TO CREATE FITS #########
38 abrinke1 1.5 ## #run_list = [179497,179547,179558,179563,179889,179959,179977,180072,180076,180093,180241,180250,180252]
39 abrinke1 1.7 ## #run_list = [180241]
40     ## trig_name = "IsoMu"
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 abrinke1 1.7 ## plot_properties = [["delivered", "rate", True, False, False, ""]]
53 abrinke1 1.5
54 abrinke1 1.1 ## masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
55 abrinke1 1.7 ## save_fits = True
56 abrinke1 1.1
57    
58 abrinke1 1.7 ###### TO SEE RATE VS PREDICTION ########
59     run_list = [180252]
60 abrinke1 1.1
61 abrinke1 1.7 trig_name = "Mu"
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 abrinke1 1.7 min_rate = 1.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 abrinke1 1.7 plot_properties = [["ls", "rawrate", 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 abrinke1 1.7 ##Limits Rates[] to runs in run_list
261     NewTrigger = {}
262     for key in Rates[print_trigger]:
263     NewTrigger[key] = []
264     for iterator in range (len(Rates[print_trigger]["run"])):
265     if Rates[print_trigger]["run"][iterator] in run_list:
266     for key in Rates[print_trigger]:
267     NewTrigger[key].append(Rates[print_trigger][key][iterator])
268     Rates[print_trigger] = NewTrigger
269    
270 abrinke1 1.1 meanrawrate = sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])
271     if not trig_name in print_trigger:
272     continue
273     if meanrawrate < min_rate:
274     continue
275     masked_trig = False
276     for mask in masked_triggers:
277     if str(mask) in print_trigger:
278     masked_trig = True
279     if masked_trig:
280     continue
281    
282 abrinke1 1.5 OutputFit[print_trigger] = {}
283 abrinke1 1.1
284     lowlumi = 0
285 abrinke1 1.7 meanlumi_init = median(Rates[print_trigger]["live_lumi"])
286 abrinke1 1.1 meanlumi = 0
287     highlumi = 0
288     lowxsec = 0
289     meanxsec = 0
290     highxsec = 0
291     nlow = 0
292     nhigh = 0
293     for iterator in range(len(Rates[print_trigger]["rate"])):
294     if Rates[print_trigger]["live_lumi"][iterator] <= meanlumi_init:
295 abrinke1 1.7 if ( 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 and Rates[print_trigger]["live_lumi"] > 500):
296 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
297     lowxsec+=Rates[print_trigger]["xsec"][iterator]
298     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
299     lowlumi+=Rates[print_trigger]["live_lumi"][iterator]
300     nlow+=1
301     if Rates[print_trigger]["live_lumi"][iterator] > meanlumi_init:
302 abrinke1 1.7 if ( 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 and Rates[print_trigger]["live_lumi"] > 500):
303 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
304     highxsec+=Rates[print_trigger]["xsec"][iterator]
305     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
306     highlumi+=Rates[print_trigger]["live_lumi"][iterator]
307     nhigh+=1
308 abrinke1 1.7 try:
309     meanxsec = meanxsec/(nlow+nhigh)
310     meanlumi = meanlumi/(nlow+nhigh)
311     slopexsec = ( (highxsec/nhigh) - (lowxsec/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
312     except:
313     print str(print_trigger)+" has no good datapoints - setting initial xsec slope estimate to 0"
314     meanxsec = median(Rates[print_trigger]["xsec"])
315     meanlumi = median(Rates[print_trigger]["live_lumi"])
316     slopexsec = 0
317 abrinke1 1.1
318 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()
319    
320     if fit_file:
321     FitType = InputFit[print_trigger][0]
322     X0 = InputFit[print_trigger][1]
323     X1 = InputFit[print_trigger][2]
324     X2 = InputFit[print_trigger][3]
325     X3 = InputFit[print_trigger][4]
326     Chi2 = InputFit[print_trigger][5]
327 abrinke1 1.7 ##print str(print_trigger)+" "+str(FitType)+" "+str(X0)+" "+str(X1)+" "+str(X2)+" "+str(X3)
328 abrinke1 1.5
329 abrinke1 1.1 for iterator in range(len(Rates[print_trigger]["rate"])):
330     if not Rates[print_trigger]["run"][iterator] in run_list:
331     continue
332     prediction = meanxsec + slopexsec * (Rates[print_trigger]["live_lumi"][iterator] - meanlumi)
333     realvalue = Rates[print_trigger]["xsec"][iterator]
334 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):
335 abrinke1 1.1 run_t.append(Rates[print_trigger]["run"][iterator])
336     ls_t.append(Rates[print_trigger]["ls"][iterator])
337     ps_t.append(Rates[print_trigger]["ps"][iterator])
338     inst_t.append(Rates[print_trigger]["inst_lumi"][iterator])
339     live_t.append(Rates[print_trigger]["live_lumi"][iterator])
340     delivered_t.append(Rates[print_trigger]["delivered_lumi"][iterator])
341     deadtime_t.append(Rates[print_trigger]["deadtime"][iterator])
342     rawrate_t.append(Rates[print_trigger]["rawrate"][iterator])
343     rate_t.append(Rates[print_trigger]["rate"][iterator])
344     rawxsec_t.append(Rates[print_trigger]["rawxsec"][iterator])
345     xsec_t.append(Rates[print_trigger]["xsec"][iterator])
346 abrinke1 1.5 psi_t.append(Rates[print_trigger]["psi"][iterator])
347 abrinke1 1.1
348     e_run_t.append(0.0)
349     e_ls_t.append(0.0)
350     e_ps_t.append(0.0)
351     e_inst_t.append(14.14)
352     e_live_t.append(14.14)
353     e_delivered_t.append(14.14)
354 abrinke1 1.2 e_deadtime_t.append(0.01)
355 abrinke1 1.1 e_rawrate_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
356     e_rate_t.append(Rates[print_trigger]["ps"][iterator]*math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
357 abrinke1 1.5 e_psi_t.append(0.0)
358 abrinke1 1.2 if live_t[-1] == 0:
359     e_rawxsec_t.append(0)
360     e_xsec_t.append(0)
361     else:
362     e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["live_lumi"][iterator])
363     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])
364 abrinke1 1.5
365     if fit_file:
366     if FitType == "expo":
367     rate_prediction = X0 + X1*math.exp(X2*delivered_t[-1])
368     else:
369     rate_prediction = X0 + X1*delivered_t[-1] + X2*delivered_t[-1]*delivered_t[-1] + X3*delivered_t[-1]*delivered_t[-1]*delivered_t[-1]
370     ## if rate_t[-1] < 0.7 * rate_prediction or rate_t[-1] > 1.4 * rate_prediction:
371     ## 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])
372    
373 abrinke1 1.2 if live_t[-1] == 0:
374 abrinke1 1.5 rawrate_fit_t.append(0)
375     rate_fit_t.append(0)
376 abrinke1 1.2 rawxsec_fit_t.append(0)
377     xsec_fit_t.append(0)
378 abrinke1 1.7 e_rawrate_fit_t.append(0)
379     e_rate_fit_t.append(math.sqrt(Chi2))
380     e_rawxsec_fit_t.append(0)
381     e_xsec_fit_t.append(0)
382 abrinke1 1.2 else:
383 abrinke1 1.5 rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
384     rate_fit_t.append(rate_prediction)
385     rawxsec_fit_t.append(rawrate_fit_t[-1]/live_t[-1])
386     xsec_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/live_t[-1])
387 abrinke1 1.7 e_rawrate_fit_t.append(math.sqrt(Chi2)*rawrate_fit_t[-1]/rate_fit_t[-1])
388     e_rate_fit_t.append(math.sqrt(Chi2))
389     e_rawxsec_fit_t.append(math.sqrt(Chi2)*rawxsec_fit_t[-1]/rate_fit_t[-1])
390     e_xsec_fit_t.append(math.sqrt(Chi2)*xsec_fit_t[-1]/rate_fit_t[-1])
391 abrinke1 1.5
392     else: ##If the data point does not pass the data_clean filter
393 abrinke1 1.1 if debug_print:
394     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)
395    
396 abrinke1 1.5 ## End "for iterator in range(len(Rates[print_trigger]["rate"])):" loop
397 abrinke1 1.1
398 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]
399     [VX, VXE, VY, VYE, VF, VFE] = GetVXVY(plot_properties, fit_file, AllPlotArrays)
400 abrinke1 1.1
401 abrinke1 1.5 if save_root or save_png:
402     c1 = TCanvas(str(varX),str(varY))
403     c1.SetName(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
404    
405     gr1 = TGraphErrors(len(VX), VX, VY, VXE, VYE)
406     gr1.SetName("Graph_"+str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
407     gr1.GetXaxis().SetTitle(varX)
408     gr1.GetYaxis().SetTitle(varY)
409     gr1.SetTitle(str(print_trigger))
410     gr1.SetMinimum(0)
411     gr1.SetMaximum(1.2*max(VY))
412     #gr1.GetXaxis().SetLimits(min(VX)-0.2*max(VX),1.2*max(VX))
413     gr1.GetXaxis().SetLimits(0,1.2*max(VX))
414     gr1.SetMarkerStyle(8)
415     if fit_file:
416     gr1.SetMarkerSize(0.8)
417     else:
418     gr1.SetMarkerSize(0.5)
419     gr1.SetMarkerColor(2)
420    
421     if fit_file:
422     gr3 = TGraphErrors(len(VX), VX, VF, VXE, VFE)
423     gr3.SetMarkerStyle(8)
424     gr3.SetMarkerSize(0.4)
425     gr3.SetMarkerColor(4)
426     gr3.SetFillColor(4)
427     gr3.SetFillStyle(3003)
428 abrinke1 1.1
429 abrinke1 1.5 if do_fit:
430     if "rate" in varY:
431     f1a = TF1("f1a","pol2",0,8000)
432     f1a.SetLineColor(4)
433     f1a.SetLineWidth(2)
434     f1a.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
435     f1a.SetParLimits(1,0,2.0*max(VY)/(max(VX)*max(VX)))
436     #gr1.Fit("f1a","B","Q")
437     gr1.Fit("f1a","Q","rob=0.90")
438    
439     f1b = 0
440     f1c = 0
441     if True:
442     f1b = TF1("f1b","pol3",0,8000)
443     f1b.SetLineColor(2)
444     f1b.SetLineWidth(2)
445     f1b.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
446     f1b.SetParLimits(1,0,f1a.GetParameter(1)+0.0000001)
447     f1b.SetParLimits(2,0,f1a.GetParameter(2)+0.0000000001)
448     f1b.SetParLimits(3,0,2.0*max(VY)/(max(VX)*max(VX)*max(VX)))
449     gr1.Fit("f1b","Q","rob=0.90")
450     #if f1b.GetChisquare()/f1b.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
451     print "X0 = "+str(f1a.GetParameter(0))+" X1 = "+str(f1a.GetParameter(1))+" X2 = "+str(f1a.GetParameter(2))
452     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()))
453     print "X0 = "+str(f1b.GetParameter(0))+" X1 = "+str(f1b.GetParameter(1))+" X2 = "+str(f1b.GetParameter(2))+" X3 = "+str(f1b.GetParameter(3))
454 abrinke1 1.1
455 abrinke1 1.5 f1c = TF1("f1c","[0]+[1]*expo(2)",0,8000)
456     f1c.SetLineColor(3)
457     f1c.SetLineWidth(2)
458     f1c.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
459     f1c.SetParLimits(1,max(VY)/math.exp(10.0),max(VY)/math.exp(2.0))
460     f1c.SetParLimits(2,0.0,0.0000000001)
461     f1c.SetParLimits(3,2.0/max(VX),10.0/max(VX))
462     print str(max(VY)/math.exp(2.0))+" "+str(10.0/max(VX))
463     gr1.Fit("f1c","Q","rob=0.90")
464     #if f1c.GetChisquare()/f1c.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
465     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()))
466     print "X0 = "+str(f1c.GetParameter(0))+" X1 = "+str(f1c.GetParameter(1))+" X2 = "+str(f1c.GetParameter(2))+" X3 = "+str(f1c.GetParameter(3))
467 abrinke1 1.2
468 abrinke1 1.5 else: ##If this is not a rate plot
469     f1a = TF1("f1a","pol1",0,8000)
470     f1a.SetLineColor(4)
471     f1a.SetLineWidth(2)
472     if "xsec" in varY:
473     f1a.SetParLimits(0,0,meanxsec*1.5)
474     if slopexsec > 0:
475     f1a.SetParLimits(1,0,max(VY)/max(VX))
476     else:
477     f1a.SetParLimits(1,2*slopexsec,-2*slopexsec)
478 abrinke1 1.1 else:
479 abrinke1 1.5 f1a.SetParLimits(0,-1000,1000)
480     gr1.Fit("f1a","Q","rob=0.80")
481 abrinke1 1.1
482 abrinke1 1.5 if save_root or save_png:
483     gr1.Draw("APZ")
484     ## ##Option to draw stats box
485     ## p1 = TPaveStats()
486     ## p1 = gr1.GetListOfFunctions().FindObject("stats")
487     ## print p1
488     ## gr1.PaintStats(f1b).Draw("same")
489     if fit_file:
490     gr3.Draw("P3")
491     if do_fit:
492     f1a.Draw("same")
493     try:
494     f1b.Draw("same")
495     f1c.Draw("same")
496     except:
497     True
498     c1.Update()
499     if save_root:
500     myfile = TFile( RootFile, 'UPDATE' )
501     c1.Write()
502     myfile.Close()
503     if save_png:
504     c1.SaveAs(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX)+".png")
505 abrinke1 1.1
506    
507     if print_table or save_fits:
508 abrinke1 1.5 if not do_fit:
509     print "Can't have save_fits = True and do_fit = False"
510     continue
511     if f1c.GetChisquare()/f1c.GetNDF() < 0.95*f1a.GetChisquare()/f1a.GetNDF() and f1c.GetChisquare()/f1c.GetNDF() < 0.95*f1b.GetChisquare()/f1b.GetNDF():
512     OutputFit[print_trigger] = ["expo", f1c.GetParameter(0), f1c.GetParameter(1), f1c.GetParameter(3), 0.0, f1c.GetChisquare()/f1c.GetNDF(), meanrawrate]
513     elif f1b.GetChisquare()/f1b.GetNDF() < 0.95*f1a.GetChisquare()/f1a.GetNDF():
514     OutputFit[print_trigger] = ["poly", f1b.GetParameter(0), f1b.GetParameter(1), f1b.GetParameter(2), f1b.GetParameter(3), f1b.GetChisquare()/f1b.GetNDF(), meanrawrate]
515     else:
516     OutputFit[print_trigger] = ["poly", f1a.GetParameter(0), f1a.GetParameter(1), f1a.GetParameter(2), 0.0, f1a.GetChisquare()/f1a.GetNDF(), meanrawrate]
517 abrinke1 1.1
518     if save_root:
519     print "Output root file is "+str(RootFile)
520    
521     if save_fits:
522     FitNameTemplate = "Fits/2011/Fit_%s_%sLS_Run%sto%s.pkl"
523     FitFile = FitNameTemplate % (trig_name, num_ls, min_run, max_run)
524     if os.path.exists(FitFile):
525     os.remove(FitFile)
526 abrinke1 1.5 FitOutputFile = open(FitFile, 'wb')
527     pickle.dump(OutputFit, FitOutputFile, 2)
528     FitOutputFile.close()
529 abrinke1 1.7 print "Output fit file is "+str(FitFile)
530 abrinke1 1.1
531     if print_table:
532 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."
533     print '%60s%10s%10s%10s%10s%10s%10s%10s' % ("Trig", "fit", "p0", "p1", "p2", "p3", "Chi2", "Av raw")
534     for print_trigger in OutputFit:
535     _trigger = (print_trigger[:56] + '...') if len(print_trigger) > 59 else print_trigger
536 abrinke1 1.1 try:
537 abrinke1 1.5 if OutputFit[print_trigger][0] == "poly":
538     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))
539     else:
540     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))
541 abrinke1 1.1 except:
542     print str(print_trigger)+" is somehow broken"
543    
544 abrinke1 1.5
545     ############# SUPPORTING FUNCTIONS ################
546    
547    
548     def GetJSON(json_file):
549    
550     input_file = open(json_file)
551     file_content = input_file.read()
552     inputRange = selectionParser(file_content)
553     JSON = inputRange.runsandls()
554     return JSON
555     ##JSON is an array: JSON[run_number] = [1st ls, 2nd ls, 3rd ls ... nth ls]
556    
557     def MakePlotArrays():
558     run_t = array.array('f')
559     ls_t = array.array('f')
560     ps_t = array.array('f')
561     inst_t = array.array('f')
562     live_t = array.array('f')
563     delivered_t = array.array('f')
564     deadtime_t = array.array('f')
565     rawrate_t = array.array('f')
566     rate_t = array.array('f')
567     rawxsec_t = array.array('f')
568     xsec_t = array.array('f')
569     psi_t = array.array('f')
570    
571     e_run_t = array.array('f')
572     e_ls_t = array.array('f')
573     e_ps_t = array.array('f')
574     e_inst_t = array.array('f')
575     e_live_t = array.array('f')
576     e_delivered_t = array.array('f')
577     e_deadtime_t = array.array('f')
578     e_rawrate_t = array.array('f')
579     e_rate_t = array.array('f')
580     e_rawxsec_t = array.array('f')
581     e_xsec_t = array.array('f')
582     e_psi_t = array.array('f')
583    
584     rawrate_fit_t = array.array('f')
585     rate_fit_t = array.array('f')
586     rawxsec_fit_t = array.array('f')
587     xsec_fit_t = array.array('f')
588     e_rawrate_fit_t = array.array('f')
589     e_rate_fit_t = array.array('f')
590     e_rawxsec_fit_t = array.array('f')
591     e_xsec_fit_t = array.array('f')
592    
593     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]
594    
595    
596     def GetVXVY(plot_properties, fit_file, AllPlotArrays):
597    
598     VF = "0"
599     VFE = "0"
600    
601     [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
602     for varX, varY, do_fit, save_root, save_png, fit_file in plot_properties:
603     if varX == "run":
604     VX = run_t
605     VXE = run_t_e
606     elif varX == "ls":
607     VX = ls_t
608     VXE = e_ls_t
609     elif varX == "ps":
610     VX = ps_t
611     VXE = e_ps_t
612     elif varX == "inst":
613     VX = inst_t
614     VXE = e_inst_t
615     elif varX == "live":
616     VX = live_t
617     VXE = e_live_t
618     elif varX == "delivered":
619     VX = delivered_t
620     VXE = e_delivered_t
621     elif varX == "deadtime":
622     VX = deadtime_t
623     VXE = e_deadtime_t
624     elif varX == "rawrate":
625     VX = rawrate_t
626     VXE = e_rawrate_t
627     elif varX == "rate":
628     VX = rate_t
629     VXE = e_rate_t
630     elif varX == "rawxsec":
631     VX = rawxsec_t
632     VXE = e_rawxsec_t
633     elif varX == "xsec":
634     VX = xsec_t
635     VXE = e_xsec_t
636     elif varX == "psi":
637     VX = psi_t
638     VXE = e_psi_t
639     else:
640     print "No valid variable entered for X"
641     continue
642     if varY == "run":
643     VY = run_t
644     VYE = run_t_e
645     elif varY == "ls":
646     VY = ls_t
647     VYE = e_ls_t
648     elif varY == "ps":
649     VY = ps_t
650     VYE = e_ps_t
651     elif varY == "inst":
652     VY = inst_t
653     VYE = e_inst_t
654     elif varY == "live":
655     VY = live_t
656     VYE = e_live_t
657     elif varY == "delivered":
658     VY = delivered_t
659     VYE = e_delivered_t
660     elif varY == "deadtime":
661     VY = deadtime_t
662     VYE = e_deadtime_t
663     elif varY == "rawrate":
664     VY = rawrate_t
665     VYE = e_rawrate_t
666     if fit_file:
667     VF = rawrate_fit_t
668     VFE = e_rawrate_fit_t
669     elif varY == "rate":
670     VY = rate_t
671     VYE = e_rate_t
672     if fit_file:
673     VF = rate_fit_t
674     VFE = e_rate_fit_t
675     elif varY == "rawxsec":
676     VY = rawxsec_t
677     VYE = e_rawxsec_t
678     if fit_file:
679     VF = rawxsec_fit_t
680     VFE = e_rawxsec_fit_t
681     elif varY == "xsec":
682     VY = xsec_t
683     VYE = e_xsec_t
684     if fit_file:
685     VF = xsec_fit_t
686     VFE = e_xsec_fit_t
687     elif varY == "psi":
688     VY = psi_t
689     VYE = e_psi_t
690     else:
691     print "No valid variable entered for Y"
692     continue
693    
694     return [VX, VXE, VY, VYE, VF, VFE]
695    
696    
697 abrinke1 1.1 if __name__=='__main__':
698     main()