ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.8
Committed: Thu Mar 8 12:41:29 2012 UTC (13 years, 2 months ago) by grchrist
Content type: text/x-python
Branch: MAIN
CVS Tags: V00-00-13
Changes since 1.7: +78 -14 lines
Log Message: 
now draws triggers from list, change of pkl in RateMonitor, new function (not yet called) to get additional lumi info in DataBaseParser

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