ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.25
Committed: Mon Mar 26 20:11:55 2012 UTC (13 years, 1 month ago) by grchrist
Content type: text/x-python
Branch: MAIN
Changes since 1.24: +39 -6 lines
Log Message:
a few more changes to make things clearer

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 amott 1.18 import getopt
10 abrinke1 1.1
11     from ROOT import gROOT, TCanvas, TF1, TGraph, TGraphErrors, TPaveStats, gPad, gStyle
12 amott 1.18 from ROOT import TFile, TPaveText, TBrowser
13 abrinke1 1.1 from ROOT import gBenchmark
14     import array
15     import math
16 grchrist 1.8 from ReadConfig import RateMonConfig
17 abrinke1 1.1
18 abrinke1 1.5 from selectionParser import selectionParser
19    
20 amott 1.18 def usage():
21     print sys.argv[0]+" [options] <list of runs>"
22     print "This script is used to generate fits and do secondary shifter validation"
23     print "<list of runs> this is a list of the form: a b c-d e f-g, specifying individual runs and/or run ranges"
24     print " be careful with using ranges (c-d), it is highly recommended to use a JSON in this case"
25     print "options: "
26     print "--makeFits run in fit making mode"
27     print "--secondary run in secondary shifter mode"
28     print "--fitFile=<path> path to the fit file"
29     print "--json=<path> path to the JSON file"
30     print "--TriggerList=<path> path to the trigger list (without versions!)"
31    
32     class Modes:
33     none,fits,secondary = range(3)
34    
35 abrinke1 1.1 def main():
36 amott 1.18 try:
37 grchrist 1.22 try:
38     opt, args = getopt.getopt(sys.argv[1:],"",["makeFits","secondary","fitFile=","json=","TriggerList="])
39    
40     except getopt.GetoptError, err:
41     print str(err)
42     usage()
43     sys.exit(2)
44    
45 grchrist 1.21 ## if len(args)<1:
46     ## print "\nPlease specify at least 1 run to look at\n"
47     ## usage()
48     ## sys.exit(0)
49 amott 1.18
50 grchrist 1.25
51     ##### RUN LIST ########
52 grchrist 1.22 run_list=[]
53 grchrist 1.23
54 grchrist 1.22 if len(args)<1:
55     inputrunlist=[]
56     print "No runs specified"
57     runinput=raw_input("Enter run range in form <run1> <run2> <run3> or <run1>-<run2>:")
58     inputrunlist.append(runinput)
59    
60    
61     if runinput.find(' ')!=-1:
62     args=runinput.split(' ')
63     else:
64     args.append(runinput)
65    
66 grchrist 1.23
67    
68 grchrist 1.22 for r in args:
69     if r.find('-')!=-1: # r is a run range
70     rrange = r.split('-')
71     if len(rrange)!=2:
72     print "Invalid run range %s" % (r,)
73     sys.exit(0)
74     try:
75     for rr in range(int(rrange[0]),int(rrange[1])+1):
76     run_list.append(rr)
77     except:
78     print "Invalid run range %s" % (r,)
79     sys.exit(0)
80     else: # r is not a run range
81     try:
82     run_list.append(int(r))
83     except:
84     print "Invalid run %s" % (r,)
85 grchrist 1.23 print "run list=",run_list
86 grchrist 1.21
87 grchrist 1.25 ##### READ CMD LINE ARGS #########
88 grchrist 1.22 mode = Modes.none
89     fitFile = ""
90     jsonfile = ""
91     trig_list = []
92 grchrist 1.25 max_dt=-1.0
93     subsys=-1.0
94     SubSystemOff={'All':False,'Mu':False,'HCal':False,'ECal':False,'Tracker':False,'EndCap':False,'Beam':False}
95 grchrist 1.22 for o,a in opt:
96     if o == "--makeFits":
97     mode = Modes.fits
98     elif o == "--secondary":
99     mode = Modes.secondary
100     elif o == "--fitFile":
101     fitFile = str(a)
102     elif o == "--json":
103     jsonfile = a
104 grchrist 1.25 elif o== "--maxdt":
105     max_dt = a
106    
107     elif o=="--All":
108     subsys=1
109     SubSystemOff["All"]=True
110     elif o=="--Mu":
111     subsys=1
112     SubSystemOff["Mu"]=True
113     elif o=="--HCal":
114     SubSystemOff["HCal"]=True
115     subsys=1
116     elif o=="--Tracker":
117     SubSystemOff["Tracker"]=True
118     subsys=1
119     elif o=="--EndCap":
120     SubSystemOff["EndCap"]=True
121     subsys=1
122     elif o=="--Beam":
123     SubSystemOff["Beam"]=True
124     subsys=1
125 grchrist 1.22 elif o == "--TriggerList":
126     try:
127     f = open(a)
128     for entry in f:
129     if entry.startswith('#'):
130     continue
131     if entry.find(':')!=-1:
132     entry = entry[:entry.find(':')] ## We can point this to the existing monitor list, just remove everything after ':'!
133     if entry.find('#')!=-1:
134     entry = entry[:entry.find('#')] ## We can point this to the existing monitor list, just remove everything after ':'!
135     trig_list.append( entry.rstrip('\n'))
136     except:
137     print "\nInvalid Trigger List\n"
138     sys.exit(0)
139     else:
140     print "\nInvalid Option %s\n" % (str(o),)
141     usage()
142     sys.exit(2)
143    
144     print "\n\n"
145 grchrist 1.25 ###### MODES #########
146    
147 grchrist 1.22 if mode == Modes.none: ## no mode specified
148     print "\nNo operation mode specified!\n"
149     modeinput=raw_input("Enter mode, --makeFits or --secondary:")
150     print "modeinput=",modeinput
151     if not (modeinput=="--makeFits" or modeinput=="--secondary"):
152     print "not either"
153     usage()
154 amott 1.18 sys.exit(0)
155 grchrist 1.22 elif modeinput == "--makeFits":
156     mode=Modes.fits
157     elif modeinput =="--secondary":
158     mode=Modes.secondary
159     else:
160     print "FATAL ERROR: No Mode specified"
161 amott 1.18 sys.exit(0)
162 grchrist 1.22
163     if mode == Modes.fits:
164     print "Running in Fit Making mode\n\n"
165     elif mode == Modes.secondary:
166     print "Running in Secondary Shifter mode\n\n"
167     else: ## should never get here, but exit if we do
168 grchrist 1.21 print "FATAL ERROR: No Mode specified"
169     sys.exit(0)
170 grchrist 1.22
171     if fitFile=="" and not mode==Modes.fits:
172     print "\nPlease specify fit file. These are available:\n"
173     path="Fits/2011/" # insert the path to the directory of interest
174     dirList=os.listdir(path)
175     for fname in dirList:
176     print fname
177 grchrist 1.24 fitFile = path+raw_input("Enter fit file in format Fit_HLT_10LS_Run176023to180252.pkl: ")
178    
179 grchrist 1.21 ##usage()
180     ##sys.exit(0)
181 grchrist 1.22 elif fitFile=="":
182     fitFile="Fits/2011/Fit_HLT_10LS_Run%sto%s.pkl" % (min(run_list),max(run_list))
183     print "fit file=",fitFile
184 amott 1.18
185 grchrist 1.24 print "fitFile=",fitFile
186 grchrist 1.25
187     ###### TRIGGER LIST #######
188 grchrist 1.24
189 grchrist 1.22 if trig_list == []:
190    
191     print "\nPlease specify list of triggers\n"
192     print "Available lists are:"
193     dirList=os.listdir(".")
194     for fname in dirList:
195     entry=fname
196     if entry.find('.')!=-1:
197     extension = entry[entry.find('.'):] ## We can point this to the existing monitor list, just remove everything after ':'!
198     if extension==".list":
199     print fname
200     trig_input=raw_input("\nEnter triggers in format HLT_IsoMu24_eta2p1 or a .list file: ")
201 grchrist 1.21
202 grchrist 1.22 if trig_input.find('.')!=-1:
203     extension = trig_input[trig_input.find('.'):]
204 grchrist 1.21 if extension==".list":
205 grchrist 1.22 try:
206     fl=open(trig_input)
207     except:
208     print "Cannot open file"
209     usage()
210     sys.exit(0)
211 grchrist 1.21
212 grchrist 1.22 for line in fl:
213     if line.startswith('#'):
214     continue
215     if len(line)<1:
216     continue
217 grchrist 1.21
218 grchrist 1.22 if len(line)>=2:
219     arg=line.rstrip('\n').rstrip(' ').lstrip(' ')
220     trig_list.append(arg)
221     else:
222     arg=''
223     else:
224     trig_list.append(trig_input)
225 grchrist 1.21
226    
227    
228    
229     ##usage()
230     ##sys.exit(0)
231    
232    
233 amott 1.18
234 abrinke1 1.5 ## Can use any combination of LowestRunNumber, HighestRunNumber, and NumberOfRuns -
235     ## just modify "ExistingRuns.sort" and for "run in ExistingRuns" accordingly
236    
237 grchrist 1.22 if jsonfile=="":
238     JSON=[]
239     else:
240     print "Using JSON: %s" % (jsonfile,)
241     JSON = GetJSON(jsonfile) ##Returns array JSON[runs][ls_list]
242 abrinke1 1.5
243 grchrist 1.21
244 abrinke1 1.5
245 grchrist 1.22 ###### TO CREATE FITS #########
246     if mode == Modes.fits:
247     trig_name = "HLT"
248     num_ls = 10
249     physics_active_psi = True ##Requires that physics and active be on, and that the prescale column is not 0
250     #JSON = [] ##To not use a JSON file, just leave the array empty
251     debug_print = False
252     no_versions=False
253     min_rate = 0.1
254     print_table = False
255     data_clean = True ##Gets rid of anomalous rate points, reqires physics_active_psi (PAP) and deadtime < 20%
256     ##plot_properties = [varX, varY, do_fit, save_root, save_png, fit_file]
257     plot_properties = [["delivered", "rate", True, True, False, fitFile]]
258 amott 1.18
259 grchrist 1.22 masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
260     save_fits = True
261 grchrist 1.25 if max_dt==-1.0:
262     max_dt=0.08 ## no deadtime cutuse 2.0
263 grchrist 1.22 force_new=True
264     print_info=True
265 grchrist 1.25 if subsys==-1.0:
266     SubSystemOff={'All':False,'Mu':False,'HCal':False,'ECal':False,'Tracker':False,'EndCap':False,'Beam':True}
267 grchrist 1.22
268    
269     ###### TO SEE RATE VS PREDICTION ########
270     if mode == Modes.secondary:
271     trig_name = "HLT"
272     num_ls = 1
273     physics_active_psi = True
274     debug_print = False
275     no_versions=False
276     min_rate = 1.0
277     print_table = False
278     data_clean = True
279     ##plot_properties = [varX, varY, do_fit, save_root, save_png, fit_file]
280     ##plot_properties = [["ls", "rawrate", False, True, False, "Fits/2011/Fit_HLT_10LS_Run176023to180252.pkl"]]
281     ##plot_properties = [["ls", "rawrate", False, True, False, "Fits/2011/Fit_HLT_10LS_Run179497to180252.pkl"]]
282     plot_properties = [["ls", "rawrate", False, True, False,fitFile]]
283    
284     masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
285     save_fits = False
286 grchrist 1.25 if max_dt==-1.0:
287     max_dt=2.0 ## no deadtime cut=2.0
288 grchrist 1.22 force_new=True
289     print_info=True
290 grchrist 1.25 if subsys==-1.0:
291     SubSystemOff={'All':True,'Mu':False,'HCal':False,'ECal':False,'Tracker':False,'EndCap':False,'Beam':True}
292 grchrist 1.22
293    
294 grchrist 1.13
295    
296 grchrist 1.10
297 grchrist 1.22 ######## END PARAMETERS - CALL FUNCTIONS ##########
298     [Rates,LumiPageInfo]= GetDBRates(run_list, trig_name, trig_list, num_ls, max_dt, physics_active_psi, JSON, debug_print, force_new, SubSystemOff)
299     ##if not checkLS(Rates,LumiPageInfo,trig_list):
300     ## print "Missing LS!"
301    
302    
303     ##for iterator in range(len(Rates["HLT_IsoMu30_eta2p1"]["rawrate"])):
304     ## print iterator, "ls=",Rates["HLT_IsoMu30_eta2p1"]["ls"][iterator],"rate=",round(Rates["HLT_IsoMu30_eta2p1"]["rawrate"][iterator],2)
305    
306     rootFileName = MakePlots(Rates, LumiPageInfo, run_list, trig_name, trig_list, num_ls, min_rate, max_dt, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print,SubSystemOff, print_info)
307     except KeyboardInterrupt:
308     print "Wait... come back..."
309 abrinke1 1.1
310 amott 1.18 def GetDBRates(run_list,trig_name,trig_list_noV, num_ls, max_dt, physics_active_psi,JSON,debug_print, force_new, SubSystemOff):
311 abrinke1 1.1
312     Rates = {}
313 grchrist 1.10 LumiPageInfo={}
314 abrinke1 1.5 ## Save in RefRuns with name dependent on trig_name, num_ls, JSON, and physics_active_psi
315     if JSON:
316 amott 1.18 #print "Using JSON file"
317 abrinke1 1.5 if physics_active_psi:
318     RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_JPAP.pkl"
319     else:
320     RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_JSON.pkl"
321     else:
322 grchrist 1.14 print "Using Physics and Active ==1"
323 abrinke1 1.5 if physics_active_psi:
324     RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_PAP.pkl"
325     else:
326     RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS.pkl"
327 grchrist 1.10
328    
329 abrinke1 1.1 RefRunFile = RefRunNameTemplate % (trig_name,num_ls)
330 abrinke1 1.5 RefRunFileHLT = RefRunNameTemplate % ("HLT",num_ls)
331 grchrist 1.10
332     print "RefRun=",RefRunFile
333     print "RefRunFileHLT",RefRunFileHLT
334 grchrist 1.11 if not force_new:
335     try: ##Open an existing RefRun file with the same parameters and trigger name
336     pkl_file = open(RefRunFile, 'rb')
337     Rates = pickle.load(pkl_file)
338     pkl_file.close()
339     os.remove(RefRunFile)
340     print "using",RefRunFile
341    
342    
343 abrinke1 1.5 except:
344 grchrist 1.11 try: ##Open an existing RefRun file with the same parameters and HLT for trigger name
345     pkl_file = open(RefRunFileHLT)
346     HLTRates = pickle.load(pkl_file)
347     for key in HLTRates:
348     if trig_name in str(key):
349     Rates[key] = HLTRates[key]
350     #print str(RefRunFile)+" does not exist. Creating ..."
351     except:
352     print str(RefRunFile)+" does not exist. Creating ..."
353 grchrist 1.10
354     ## try the lumis file
355     RefLumiNameTemplate = "RefRuns/2011/Lumis_%s_%sLS.pkl"
356 grchrist 1.11 RefLumiFile= RefLumiNameTemplate % ("HLT",num_ls)
357     if not force_new:
358     try:
359     pkl_lumi_file = open(RefLumiFile, 'rb')
360     LumiPageInfo = pickle.load(pkl_lumi_file)
361     pkl_lumi_file.close()
362     os.remove(RefLumiFile)
363     print "using",RefLumiFile
364     except:
365     print str(RefLumiFile)+" doesn't exist. Make it..."
366    
367 abrinke1 1.5 for RefRunNum in run_list:
368 grchrist 1.11
369 abrinke1 1.5 if JSON:
370     if not RefRunNum in JSON:
371     continue
372 abrinke1 1.1 try:
373 grchrist 1.10
374 abrinke1 1.1 ExistsAlready = False
375     for key in Rates:
376     if RefRunNum in Rates[key]["run"]:
377     ExistsAlready = True
378     break
379 grchrist 1.10
380    
381     LumiExistsLAready=False
382     for v in LumiPageInfo.itervalues():
383     #print RefRunNum, v["Run"]
384    
385     if RefRunNum == v["Run"]:
386     LumiExistsAlready=True
387     break
388     if ExistsAlready and LumiExistsAlready:
389 abrinke1 1.1 continue
390 grchrist 1.10
391    
392 abrinke1 1.1 except:
393     print "Getting info for run "+str(RefRunNum)
394    
395     if RefRunNum < 1:
396     continue
397 grchrist 1.17 print "calculating rates and green lumis for run ",RefRunNum
398 grchrist 1.10
399 abrinke1 1.5 if True: ##Placeholder
400 abrinke1 1.1 if True: #May replace with "try" - for now it's good to know when problems happen
401     RefParser = DatabaseParser()
402     RefParser.RunNumber = RefRunNum
403     RefParser.ParseRunSetup()
404 abrinke1 1.5 RefLumiRangePhysicsActive = RefParser.GetLSRange(1,9999) ##Gets array of all LS with physics and active on
405     RefLumiArray = RefParser.GetLumiInfo() ##Gets array of all existing LS and their lumi info
406 abrinke1 1.2 RefLumiRange = []
407 grchrist 1.17 RefMoreLumiArray = RefParser.GetMoreLumiInfo()#dict with keys as bits from lumisections WBM page and values are dicts with key=LS:value=bit
408 amott 1.18
409    
410     ## We have specified the trig list without version numbers, we add them specific to this run
411 grchrist 1.22 ##print "Processing Triggers: "
412 amott 1.18 trig_list=[]
413     for entry in trig_list_noV:
414     trig_list.append(RefParser.GetTriggerVersion(entry))
415     if trig_list[-1]=="":
416     print ">> WARNING: could not find version for trigger %s, SKIPPING" % (entry,)
417     else:
418 grchrist 1.22 ##print ">> %s " % (trig_list[-1],)
419     pass
420 grchrist 1.17 #DeadTimeBeamActive=RefParser.GetDeadTimeBeamActive()
421     #print "deadtime ls run 180250=",DeadTimeBeamActive
422 abrinke1 1.5 for iterator in RefLumiArray[0]: ##Makes array of LS with proper PAP and JSON properties
423 grchrist 1.11 ##cheap way of getting PSCol None-->0
424 amott 1.18 if RefLumiArray[0][iterator] not in range(1,9):
425 grchrist 1.11 RefLumiArray[0][iterator]=0
426 grchrist 1.15
427 grchrist 1.11
428 grchrist 1.17 if not physics_active_psi or (RefLumiArray[5][iterator] == 1 and RefLumiArray[6][iterator] == 1 and RefMoreLumiArray["b1pres"][iterator]==1 and RefMoreLumiArray["b2pres"][iterator]==1 and RefMoreLumiArray["b1stab"][iterator] and RefMoreLumiArray["b2stab"][iterator]==1):
429 abrinke1 1.5 if not JSON or RefRunNum in JSON:
430     if not JSON or iterator in JSON[RefRunNum]:
431     RefLumiRange.append(iterator)
432 grchrist 1.15 #print iterator, RefLumiArray[0][iterator], "active=",RefLumiArray[5][iterator],"physics=",RefLumiArray[6][iterator]
433 grchrist 1.9 #print "hbhea=",RefMoreLumiArray['hbhea'][iterator]
434    
435 abrinke1 1.5 try:
436     nls = RefLumiRange[0]
437     LSRange = {}
438     except:
439     print "Run "+str(RefRunNum)+" has no good LS"
440     continue
441     if num_ls > len(RefLumiRange):
442 abrinke1 1.1 print "Run "+str(RefRunNum)+" is too short: from "+str(nls)+" to "+str(RefLumiRange[-1])+", while num_ls = "+str(num_ls)
443     continue
444     while nls < RefLumiRange[-1]-num_ls:
445 abrinke1 1.5 LSRange[nls] = []
446     counter = 0
447     for iterator in RefLumiRange:
448     if iterator >= nls and counter < num_ls:
449     LSRange[nls].append(iterator)
450     counter += 1
451 abrinke1 1.1 nls = LSRange[nls][-1]+1
452 abrinke1 1.5
453 grchrist 1.17 #print "Run "+str(RefRunNum)+" contains LS from "+str(min(LSRange))+" to "+str(max(LSRange))
454 grchrist 1.8 for nls in sorted(LSRange.iterkeys()):
455 grchrist 1.9
456 abrinke1 1.1 TriggerRates = RefParser.GetHLTRates(LSRange[nls])
457 grchrist 1.16 #print nls, RefLumiArray[1][nls], RefLumiArray[2][nls]
458 abrinke1 1.5
459     [inst, live, delivered, dead, pscols] = RefParser.GetAvLumiInfo(LSRange[nls])
460 grchrist 1.17 deadtimebeamactive=RefParser.GetDeadTimeBeamActive(LSRange[nls])
461     #print LSRange,deadtimebeamactive
462 abrinke1 1.2 physics = 1
463     active = 1
464 abrinke1 1.5 psi = 99
465     for iterator in LSRange[nls]: ##Gets lowest value of physics, active, and psi in the set of lumisections
466 abrinke1 1.2 if RefLumiArray[5][iterator] == 0:
467     physics = 0
468     if RefLumiArray[6][iterator] == 0:
469     active = 0
470 abrinke1 1.5 if RefLumiArray[0][iterator] < psi:
471     psi = RefLumiArray[0][iterator]
472 abrinke1 1.2
473 grchrist 1.17 MoreLumiMulti=LumiRangeGreens(RefMoreLumiArray,LSRange,nls,RefRunNum,deadtimebeamactive)
474 grchrist 1.10
475     #print MoreLumiMulti.keys()
476     # print "\n\n\n"
477    
478    
479 grchrist 1.9
480 abrinke1 1.5 if inst < 0 or live < 0 or delivered < 0:
481     print "Run "+str(RefRunNum)+" LS "+str(nls)+" inst lumi = "+str(inst)+" live lumi = "+str(live)+", delivered = "+str(delivered)+", physics = "+str(physics)+", active = "+str(active)
482 abrinke1 1.2
483 grchrist 1.10
484    
485     LumiPageInfo[nls]=MoreLumiMulti
486     ##print LumiPageInfo[nls]
487     ## try:
488     ## LumiPageInfo[keys].append(values)
489     ## except:
490     ## print "Failed",RefRunNum, nls, keys
491    
492 abrinke1 1.1 for key in TriggerRates:
493 grchrist 1.12 ##if not key in trig_list:
494     ## continue
495    
496     ##if not trig_name in key:
497     ## continue
498 amott 1.18 name = StripVersion(key)
499 grchrist 1.12
500 grchrist 1.14 ##if re.match('.*_v[0-9]+',name): ##Removes _v#
501     ## name = name[:name.rfind('_')]
502 amott 1.18 if not name in trig_list_noV:
503 grchrist 1.12 continue
504     #print "trigger=",name, trig_list
505    
506 abrinke1 1.1 if not Rates.has_key(name):
507     Rates[name] = {}
508     Rates[name]["run"] = []
509     Rates[name]["ls"] = []
510     Rates[name]["ps"] = []
511     Rates[name]["inst_lumi"] = []
512     Rates[name]["live_lumi"] = []
513     Rates[name]["delivered_lumi"] = []
514     Rates[name]["deadtime"] = []
515     Rates[name]["rawrate"] = []
516     Rates[name]["rate"] = []
517     Rates[name]["rawxsec"] = []
518     Rates[name]["xsec"] = []
519 abrinke1 1.2 Rates[name]["physics"] = []
520     Rates[name]["active"] = []
521 abrinke1 1.5 Rates[name]["psi"] = []
522 grchrist 1.9
523 grchrist 1.10 #for keys, values in MoreLumiMulti.iteritems():
524     # Rates[name][keys] = []
525    
526 grchrist 1.9
527 abrinke1 1.1 [avps, ps, rate, psrate] = TriggerRates[key]
528 grchrist 1.11 #print "TriggerRates=",TriggerRates[key], "key=",key
529 abrinke1 1.1 Rates[name]["run"].append(RefRunNum)
530     Rates[name]["ls"].append(nls)
531     Rates[name]["ps"].append(ps)
532     Rates[name]["inst_lumi"].append(inst)
533     Rates[name]["live_lumi"].append(live)
534     Rates[name]["delivered_lumi"].append(delivered)
535 grchrist 1.17 Rates[name]["deadtime"].append(deadtimebeamactive)
536 abrinke1 1.1 Rates[name]["rawrate"].append(rate)
537 abrinke1 1.2 if live == 0:
538 abrinke1 1.5 Rates[name]["rate"].append(0)
539 abrinke1 1.2 Rates[name]["rawxsec"].append(0.0)
540     Rates[name]["xsec"].append(0.0)
541     else:
542 grchrist 1.17 Rates[name]["rate"].append(psrate/(1.0-deadtimebeamactive))
543 abrinke1 1.2 Rates[name]["rawxsec"].append(rate/live)
544     Rates[name]["xsec"].append(psrate/live)
545     Rates[name]["physics"].append(physics)
546     Rates[name]["active"].append(active)
547 abrinke1 1.5 Rates[name]["psi"].append(psi)
548 grchrist 1.17 #print iterator, "LS=", nls, "dt=",round(deadtimebeamactive,2), "deliv=",delivered, "live=",live
549 grchrist 1.9
550 grchrist 1.10 #for keys, values in MoreLumiMulti.iteritems():
551     # Rates[name][keys].append(values)
552 grchrist 1.9 #print nls, name, keys, values
553     #print " "
554 abrinke1 1.5 #except: ##If we replace "if True:" with "try:"
555 abrinke1 1.1 #print "Failed to parse run "+str(RefRunNum)
556    
557 grchrist 1.10 RateOutput = open(RefRunFile, 'wb') ##Save new Rates[] to RefRuns
558     pickle.dump(Rates, RateOutput, 2)
559     RateOutput.close()
560     LumiOutput = open(RefLumiFile,'wb')
561     pickle.dump(LumiPageInfo,LumiOutput, 2)
562     LumiOutput.close()
563    
564    
565     return [Rates,LumiPageInfo]
566 abrinke1 1.1
567 grchrist 1.11 def MakePlots(Rates, LumiPageInfo, run_list, trig_name, trig_list, num_ls, min_rate, max_dt, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print, SubSystemOff, print_info):
568 abrinke1 1.1 min_run = min(run_list)
569     max_run = max(run_list)
570 grchrist 1.11
571 abrinke1 1.5 InputFit = {}
572     OutputFit = {}
573 grchrist 1.14 first_trigger=True
574 abrinke1 1.1
575     RootNameTemplate = "%s_%sLS_Run%sto%s.root"
576     RootFile = RootNameTemplate % (trig_name, num_ls, min_run, max_run)
577    
578 amott 1.20 [[varX, varY, do_fit, save_root, save_png, fit_file]] = plot_properties
579     if not do_fit:
580     try:
581 abrinke1 1.1 pkl_file = open(fit_file, 'rb')
582 abrinke1 1.5 InputFit = pickle.load(pkl_file)
583 amott 1.20 except:
584     print "ERROR: could not open fit file: %s" % (fit_file,)
585     if save_root:
586     try:
587     os.remove(RootFile)
588     except:
589     pass
590    
591     ## check that all the triggers we ask to plot are in the input fit
592     if not save_fits:
593     goodtrig_list = []
594     for trig in trig_list:
595     if not InputFit.has_key(trig):
596     print "WARNING: No Fit Prediction for Trigger %s, SKIPPING" % (trig,)
597     else:
598     goodtrig_list.append(trig)
599     trig_list = goodtrig_list
600 amott 1.18
601 abrinke1 1.1 for print_trigger in Rates:
602 abrinke1 1.7 ##Limits Rates[] to runs in run_list
603     NewTrigger = {}
604 grchrist 1.8 if not print_trigger in trig_list:
605     continue
606 abrinke1 1.7 for key in Rates[print_trigger]:
607     NewTrigger[key] = []
608     for iterator in range (len(Rates[print_trigger]["run"])):
609     if Rates[print_trigger]["run"][iterator] in run_list:
610     for key in Rates[print_trigger]:
611     NewTrigger[key].append(Rates[print_trigger][key][iterator])
612     Rates[print_trigger] = NewTrigger
613    
614 abrinke1 1.1 meanrawrate = sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])
615     if not trig_name in print_trigger:
616     continue
617     if meanrawrate < min_rate:
618     continue
619     masked_trig = False
620     for mask in masked_triggers:
621     if str(mask) in print_trigger:
622     masked_trig = True
623     if masked_trig:
624     continue
625    
626 abrinke1 1.5 OutputFit[print_trigger] = {}
627 abrinke1 1.1
628     lowlumi = 0
629 abrinke1 1.7 meanlumi_init = median(Rates[print_trigger]["live_lumi"])
630 abrinke1 1.1 meanlumi = 0
631     highlumi = 0
632     lowxsec = 0
633     meanxsec = 0
634     highxsec = 0
635     nlow = 0
636     nhigh = 0
637 grchrist 1.15
638 abrinke1 1.1 for iterator in range(len(Rates[print_trigger]["rate"])):
639     if Rates[print_trigger]["live_lumi"][iterator] <= meanlumi_init:
640 grchrist 1.11 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] < max_dt and Rates[print_trigger]["psi"][iterator] > 0 and Rates[print_trigger]["live_lumi"] > 500):
641 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
642     lowxsec+=Rates[print_trigger]["xsec"][iterator]
643     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
644     lowlumi+=Rates[print_trigger]["live_lumi"][iterator]
645     nlow+=1
646     if Rates[print_trigger]["live_lumi"][iterator] > meanlumi_init:
647 grchrist 1.11 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] < max_dt and Rates[print_trigger]["psi"][iterator] > 0 and Rates[print_trigger]["live_lumi"] > 500):
648 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
649     highxsec+=Rates[print_trigger]["xsec"][iterator]
650     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
651     highlumi+=Rates[print_trigger]["live_lumi"][iterator]
652     nhigh+=1
653 abrinke1 1.7 try:
654     meanxsec = meanxsec/(nlow+nhigh)
655     meanlumi = meanlumi/(nlow+nhigh)
656     slopexsec = ( (highxsec/nhigh) - (lowxsec/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
657     except:
658     print str(print_trigger)+" has no good datapoints - setting initial xsec slope estimate to 0"
659     meanxsec = median(Rates[print_trigger]["xsec"])
660     meanlumi = median(Rates[print_trigger]["live_lumi"])
661     slopexsec = 0
662 abrinke1 1.1
663 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()
664    
665 amott 1.20 if not do_fit:
666 abrinke1 1.5 FitType = InputFit[print_trigger][0]
667     X0 = InputFit[print_trigger][1]
668     X1 = InputFit[print_trigger][2]
669     X2 = InputFit[print_trigger][3]
670     X3 = InputFit[print_trigger][4]
671     Chi2 = InputFit[print_trigger][5]
672 grchrist 1.14 #print str(print_trigger)+" "+str(FitType)+" "+str(X0)+" "+str(X1)+" "+str(X2)+" "+str(X3)
673 amott 1.20 #if (first_trigger):
674     # print '%20s % 10s % 6s % 5s % 5s % 3s % 4s' % ('trigger', 'fit type ', 'cubic', 'quad', ' linear', ' c ', 'Chi2')
675     # first_trigger=False
676     #print '%20s % 10s % 2.2g % 2.2g % 2.2g % 2.2g % 2.2g' % (print_trigger, FitType, X3, X2, X1, X0, Chi2)
677 grchrist 1.14 #print '{}, {}, {:02.2g}, {:02.2g}, {:02.2g}, {:02.2g} '.format(print_trigger, FitType, X0, X1, X2, X3)
678 grchrist 1.8 ## we are 2 lumis off when we start! -gets worse when we skip lumis
679 grchrist 1.17 it_offset=0
680 abrinke1 1.1 for iterator in range(len(Rates[print_trigger]["rate"])):
681     if not Rates[print_trigger]["run"][iterator] in run_list:
682     continue
683     prediction = meanxsec + slopexsec * (Rates[print_trigger]["live_lumi"][iterator] - meanlumi)
684     realvalue = Rates[print_trigger]["xsec"][iterator]
685 grchrist 1.8
686 grchrist 1.11
687     if pass_cuts(data_clean, realvalue, prediction, meanxsec, Rates, print_trigger, iterator, num_ls,LumiPageInfo,SubSystemOff,max_dt,print_info, trig_list):
688 grchrist 1.8
689 abrinke1 1.1 run_t.append(Rates[print_trigger]["run"][iterator])
690     ls_t.append(Rates[print_trigger]["ls"][iterator])
691     ps_t.append(Rates[print_trigger]["ps"][iterator])
692     inst_t.append(Rates[print_trigger]["inst_lumi"][iterator])
693 grchrist 1.17 live_t.append(Rates[print_trigger]["live_lumi"][iterator])
694     delivered_t.append(Rates[print_trigger]["delivered_lumi"][iterator])
695     deadtime_t.append(Rates[print_trigger]["deadtime"][iterator])
696 abrinke1 1.1 rawrate_t.append(Rates[print_trigger]["rawrate"][iterator])
697     rate_t.append(Rates[print_trigger]["rate"][iterator])
698     rawxsec_t.append(Rates[print_trigger]["rawxsec"][iterator])
699     xsec_t.append(Rates[print_trigger]["xsec"][iterator])
700 abrinke1 1.5 psi_t.append(Rates[print_trigger]["psi"][iterator])
701 abrinke1 1.1
702     e_run_t.append(0.0)
703     e_ls_t.append(0.0)
704     e_ps_t.append(0.0)
705     e_inst_t.append(14.14)
706     e_live_t.append(14.14)
707     e_delivered_t.append(14.14)
708 abrinke1 1.2 e_deadtime_t.append(0.01)
709 abrinke1 1.1 e_rawrate_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
710     e_rate_t.append(Rates[print_trigger]["ps"][iterator]*math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
711 abrinke1 1.5 e_psi_t.append(0.0)
712 abrinke1 1.2 if live_t[-1] == 0:
713     e_rawxsec_t.append(0)
714     e_xsec_t.append(0)
715     else:
716 grchrist 1.16 try:
717     e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["live_lumi"][iterator])
718     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])
719     except:
720     e_rawxsec_t.append(0.)
721     e_xsec_t.append(0.)
722 amott 1.20 if not do_fit:
723 abrinke1 1.5 if FitType == "expo":
724     rate_prediction = X0 + X1*math.exp(X2*delivered_t[-1])
725     else:
726     rate_prediction = X0 + X1*delivered_t[-1] + X2*delivered_t[-1]*delivered_t[-1] + X3*delivered_t[-1]*delivered_t[-1]*delivered_t[-1]
727     ## if rate_t[-1] < 0.7 * rate_prediction or rate_t[-1] > 1.4 * rate_prediction:
728     ## 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])
729    
730 abrinke1 1.2 if live_t[-1] == 0:
731 abrinke1 1.5 rawrate_fit_t.append(0)
732     rate_fit_t.append(0)
733 abrinke1 1.2 rawxsec_fit_t.append(0)
734     xsec_fit_t.append(0)
735 abrinke1 1.7 e_rawrate_fit_t.append(0)
736     e_rate_fit_t.append(math.sqrt(Chi2))
737     e_rawxsec_fit_t.append(0)
738     e_xsec_fit_t.append(0)
739 grchrist 1.11 #print "live_t=0", ls_t[-1], rawrate_fit_t[-1]
740 abrinke1 1.2 else:
741 grchrist 1.11 if ps_t[-1]>0.0:
742     rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
743     else:
744     rawrate_fit_t.append(0.0)
745    
746 abrinke1 1.5 rate_fit_t.append(rate_prediction)
747     rawxsec_fit_t.append(rawrate_fit_t[-1]/live_t[-1])
748     xsec_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/live_t[-1])
749 abrinke1 1.7 e_rawrate_fit_t.append(math.sqrt(Chi2)*rawrate_fit_t[-1]/rate_fit_t[-1])
750     e_rate_fit_t.append(math.sqrt(Chi2))
751     e_rawxsec_fit_t.append(math.sqrt(Chi2)*rawxsec_fit_t[-1]/rate_fit_t[-1])
752     e_xsec_fit_t.append(math.sqrt(Chi2)*xsec_fit_t[-1]/rate_fit_t[-1])
753 grchrist 1.11 #print "live_t>0", ls_t[-1], rawrate_fit_t[-1]
754 abrinke1 1.5
755 grchrist 1.17 ##print iterator, 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"])
756 grchrist 1.16
757 grchrist 1.8
758    
759 abrinke1 1.5 else: ##If the data point does not pass the data_clean filter
760 grchrist 1.11 #print "not passed", iterator, ls_t[-1], rawrate_fit_t[-1]
761 abrinke1 1.1 if debug_print:
762     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)
763    
764 abrinke1 1.5 ## End "for iterator in range(len(Rates[print_trigger]["rate"])):" loop
765 abrinke1 1.1
766 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]
767     [VX, VXE, VY, VYE, VF, VFE] = GetVXVY(plot_properties, fit_file, AllPlotArrays)
768 abrinke1 1.1
769 abrinke1 1.5 if save_root or save_png:
770     c1 = TCanvas(str(varX),str(varY))
771     c1.SetName(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
772    
773     gr1 = TGraphErrors(len(VX), VX, VY, VXE, VYE)
774     gr1.SetName("Graph_"+str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
775     gr1.GetXaxis().SetTitle(varX)
776     gr1.GetYaxis().SetTitle(varY)
777     gr1.SetTitle(str(print_trigger))
778     gr1.SetMinimum(0)
779     gr1.SetMaximum(1.2*max(VY))
780     #gr1.GetXaxis().SetLimits(min(VX)-0.2*max(VX),1.2*max(VX))
781     gr1.GetXaxis().SetLimits(0,1.2*max(VX))
782     gr1.SetMarkerStyle(8)
783     if fit_file:
784     gr1.SetMarkerSize(0.8)
785     else:
786     gr1.SetMarkerSize(0.5)
787     gr1.SetMarkerColor(2)
788    
789 amott 1.20 if not do_fit:
790 abrinke1 1.5 gr3 = TGraphErrors(len(VX), VX, VF, VXE, VFE)
791     gr3.SetMarkerStyle(8)
792     gr3.SetMarkerSize(0.4)
793     gr3.SetMarkerColor(4)
794     gr3.SetFillColor(4)
795     gr3.SetFillStyle(3003)
796 abrinke1 1.1
797 abrinke1 1.5 if do_fit:
798     if "rate" in varY:
799     f1a = TF1("f1a","pol2",0,8000)
800     f1a.SetLineColor(4)
801     f1a.SetLineWidth(2)
802     f1a.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
803     f1a.SetParLimits(1,0,2.0*max(VY)/(max(VX)*max(VX)))
804     #gr1.Fit("f1a","B","Q")
805     gr1.Fit("f1a","Q","rob=0.90")
806    
807     f1b = 0
808     f1c = 0
809     if True:
810     f1b = TF1("f1b","pol3",0,8000)
811     f1b.SetLineColor(2)
812     f1b.SetLineWidth(2)
813     f1b.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
814     f1b.SetParLimits(1,0,f1a.GetParameter(1)+0.0000001)
815     f1b.SetParLimits(2,0,f1a.GetParameter(2)+0.0000000001)
816     f1b.SetParLimits(3,0,2.0*max(VY)/(max(VX)*max(VX)*max(VX)))
817     gr1.Fit("f1b","Q","rob=0.90")
818     #if f1b.GetChisquare()/f1b.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
819 grchrist 1.15 #print "X0 = "+str(f1a.GetParameter(0))+" X1 = "+str(f1a.GetParameter(1))+" X2 = "+str(f1a.GetParameter(2))
820     #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()))
821     #print "X0 = "+str(f1b.GetParameter(0))+" X1 = "+str(f1b.GetParameter(1))+" X2 = "+str(f1b.GetParameter(2))+" X3 = "+str(f1b.GetParameter(3))
822     if (first_trigger):
823 grchrist 1.16 print '%-60s %4s x0 x1 x2 x3 chi2 ndf chi2/ndf' % ('trigger', 'type')
824 grchrist 1.15
825     first_trigger=False
826    
827    
828 abrinke1 1.1
829 abrinke1 1.5 f1c = TF1("f1c","[0]+[1]*expo(2)",0,8000)
830     f1c.SetLineColor(3)
831     f1c.SetLineWidth(2)
832     f1c.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
833     f1c.SetParLimits(1,max(VY)/math.exp(10.0),max(VY)/math.exp(2.0))
834     f1c.SetParLimits(2,0.0,0.0000000001)
835     f1c.SetParLimits(3,2.0/max(VX),10.0/max(VX))
836 grchrist 1.15 #print str(max(VY)/math.exp(2.0))+" "+str(10.0/max(VX))
837 abrinke1 1.5 gr1.Fit("f1c","Q","rob=0.90")
838     #if f1c.GetChisquare()/f1c.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
839 grchrist 1.15 #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()))
840     #print "X0 = "+str(f1c.GetParameter(0))+" X1 = "+str(f1c.GetParameter(1))+" X2 = "+str(f1c.GetParameter(2))+" X3 = "+str(f1c.GetParameter(3))
841    
842    
843    
844    
845    
846     if (f1c.GetChisquare()/f1c.GetNDF() < f1b.GetChisquare()/f1b.GetNDF() and f1c.GetChisquare()/f1c.GetNDF() < f1a.GetChisquare()/f1a.GetNDF()):
847 grchrist 1.16 print '%-60s expo % .2f+/-%.2f % .2e+/-%.1e % .2e+/-%.1e % .2e+/-%.1e %7.2f %4.0f %5.3f ' % (print_trigger, f1c.GetParameter(0), f1c.GetParError(0), f1c.GetParameter(1), f1c.GetParError(1), 0 , 0 , 0 , 0 , f1c.GetChisquare(), f1c.GetNDF(), f1c.GetChisquare()/f1c.GetNDF())
848 grchrist 1.15 elif (f1b.GetChisquare()/f1b.GetNDF() < f1a.GetChisquare()/f1a.GetNDF()):
849 grchrist 1.16 print '%-60s cube % .2f+/-%.2f % .2e+/-%.1e % .2e+/-%.1e % .2e+/-%.1e %7.2f %4.0f %5.3f ' % (print_trigger, f1b.GetParameter(0), f1b.GetParError(0), f1b.GetParameter(1), f1b.GetParError(1), f1b.GetParameter(2), f1b.GetParError(2), f1b.GetParameter(3), f1b.GetParError(3), f1b.GetChisquare(), f1b.GetNDF(), f1b.GetChisquare()/f1b.GetNDF())
850 grchrist 1.15 else:
851 grchrist 1.16 print '%-60s quad % .2f+/-%.2f % .2e+/-%.1e % .2e+/-%.1e % .2e+/-%.1e %7.2f %4.0f %5.3f ' % (print_trigger, f1a.GetParameter(0), f1a.GetParError(0), f1a.GetParameter(1), f1a.GetParError(1), f1a.GetParameter(2), f1a.GetParError(2), 0 , 0 , f1a.GetChisquare(), f1a.GetNDF(), f1a.GetChisquare()/f1a.GetNDF())
852 abrinke1 1.2
853 grchrist 1.15
854    
855 abrinke1 1.5 else: ##If this is not a rate plot
856     f1a = TF1("f1a","pol1",0,8000)
857     f1a.SetLineColor(4)
858     f1a.SetLineWidth(2)
859     if "xsec" in varY:
860     f1a.SetParLimits(0,0,meanxsec*1.5)
861     if slopexsec > 0:
862     f1a.SetParLimits(1,0,max(VY)/max(VX))
863     else:
864     f1a.SetParLimits(1,2*slopexsec,-2*slopexsec)
865 abrinke1 1.1 else:
866 abrinke1 1.5 f1a.SetParLimits(0,-1000,1000)
867     gr1.Fit("f1a","Q","rob=0.80")
868 abrinke1 1.1
869 abrinke1 1.5 if save_root or save_png:
870     gr1.Draw("APZ")
871     ## ##Option to draw stats box
872     ## p1 = TPaveStats()
873     ## p1 = gr1.GetListOfFunctions().FindObject("stats")
874     ## print p1
875     ## gr1.PaintStats(f1b).Draw("same")
876 amott 1.20 if not do_fit:
877 abrinke1 1.5 gr3.Draw("P3")
878     if do_fit:
879     f1a.Draw("same")
880     try:
881     f1b.Draw("same")
882     f1c.Draw("same")
883     except:
884     True
885     c1.Update()
886     if save_root:
887     myfile = TFile( RootFile, 'UPDATE' )
888     c1.Write()
889     myfile.Close()
890     if save_png:
891     c1.SaveAs(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX)+".png")
892 abrinke1 1.1
893    
894     if print_table or save_fits:
895 abrinke1 1.5 if not do_fit:
896     print "Can't have save_fits = True and do_fit = False"
897     continue
898     if f1c.GetChisquare()/f1c.GetNDF() < 0.95*f1a.GetChisquare()/f1a.GetNDF() and f1c.GetChisquare()/f1c.GetNDF() < 0.95*f1b.GetChisquare()/f1b.GetNDF():
899 grchrist 1.16 OutputFit[print_trigger] = ["expo", f1c.GetParameter(0), f1c.GetParameter(1), f1c.GetParameter(3), 0.0, f1c.GetChisquare()/f1c.GetNDF(), meanrawrate, f1c.GetParError(0), f1c.GetParError(1), f1c.GetParError(2), f1c.GetParError(3)]
900 abrinke1 1.5 elif f1b.GetChisquare()/f1b.GetNDF() < 0.95*f1a.GetChisquare()/f1a.GetNDF():
901 grchrist 1.16 OutputFit[print_trigger] = ["poly", f1b.GetParameter(0), f1b.GetParameter(1), f1b.GetParameter(2), f1b.GetParameter(3), f1b.GetChisquare()/f1b.GetNDF(), meanrawrate,f1b.GetParError(0), f1b.GetParError(1), f1b.GetParError(2), f1b.GetParError(3)]
902 abrinke1 1.5 else:
903 grchrist 1.16 OutputFit[print_trigger] = ["poly", f1a.GetParameter(0), f1a.GetParameter(1), f1a.GetParameter(2), 0.0, f1a.GetChisquare()/f1a.GetNDF(), meanrawrate, f1a.GetParError(0), f1a.GetParError(1), f1a.GetParError(2), 0.0]
904 abrinke1 1.1
905     if save_root:
906     print "Output root file is "+str(RootFile)
907    
908     if save_fits:
909 amott 1.20 #FitNameTemplate = "Fits/2011/Fit_%s_%sLS_Run%sto%s.pkl"
910     #FitFile = FitNameTemplate % (trig_name, num_ls, min_run, max_run)
911     if os.path.exists(fit_file):
912     os.remove(fit_file)
913     FitOutputFile = open(fit_file, 'wb')
914 abrinke1 1.5 pickle.dump(OutputFit, FitOutputFile, 2)
915     FitOutputFile.close()
916 amott 1.20 print "Output fit file is "+str(fit_file)
917 abrinke1 1.1
918     if print_table:
919 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."
920     print '%60s%10s%10s%10s%10s%10s%10s%10s' % ("Trig", "fit", "p0", "p1", "p2", "p3", "Chi2", "Av raw")
921     for print_trigger in OutputFit:
922     _trigger = (print_trigger[:56] + '...') if len(print_trigger) > 59 else print_trigger
923 abrinke1 1.1 try:
924 abrinke1 1.5 if OutputFit[print_trigger][0] == "poly":
925     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))
926     else:
927     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))
928 abrinke1 1.1 except:
929     print str(print_trigger)+" is somehow broken"
930 amott 1.18 return RootFile
931 abrinke1 1.5
932     ############# SUPPORTING FUNCTIONS ################
933    
934    
935     def GetJSON(json_file):
936    
937     input_file = open(json_file)
938     file_content = input_file.read()
939     inputRange = selectionParser(file_content)
940     JSON = inputRange.runsandls()
941     return JSON
942     ##JSON is an array: JSON[run_number] = [1st ls, 2nd ls, 3rd ls ... nth ls]
943    
944     def MakePlotArrays():
945     run_t = array.array('f')
946     ls_t = array.array('f')
947     ps_t = array.array('f')
948     inst_t = array.array('f')
949     live_t = array.array('f')
950     delivered_t = array.array('f')
951     deadtime_t = array.array('f')
952     rawrate_t = array.array('f')
953     rate_t = array.array('f')
954     rawxsec_t = array.array('f')
955     xsec_t = array.array('f')
956     psi_t = array.array('f')
957    
958     e_run_t = array.array('f')
959     e_ls_t = array.array('f')
960     e_ps_t = array.array('f')
961     e_inst_t = array.array('f')
962     e_live_t = array.array('f')
963     e_delivered_t = array.array('f')
964     e_deadtime_t = array.array('f')
965     e_rawrate_t = array.array('f')
966     e_rate_t = array.array('f')
967     e_rawxsec_t = array.array('f')
968     e_xsec_t = array.array('f')
969     e_psi_t = array.array('f')
970    
971     rawrate_fit_t = array.array('f')
972     rate_fit_t = array.array('f')
973     rawxsec_fit_t = array.array('f')
974     xsec_fit_t = array.array('f')
975     e_rawrate_fit_t = array.array('f')
976     e_rate_fit_t = array.array('f')
977     e_rawxsec_fit_t = array.array('f')
978     e_xsec_fit_t = array.array('f')
979    
980     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]
981    
982    
983     def GetVXVY(plot_properties, fit_file, AllPlotArrays):
984    
985     VF = "0"
986     VFE = "0"
987    
988     [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
989     for varX, varY, do_fit, save_root, save_png, fit_file in plot_properties:
990     if varX == "run":
991     VX = run_t
992     VXE = run_t_e
993     elif varX == "ls":
994     VX = ls_t
995     VXE = e_ls_t
996     elif varX == "ps":
997     VX = ps_t
998     VXE = e_ps_t
999     elif varX == "inst":
1000     VX = inst_t
1001     VXE = e_inst_t
1002     elif varX == "live":
1003     VX = live_t
1004     VXE = e_live_t
1005     elif varX == "delivered":
1006     VX = delivered_t
1007     VXE = e_delivered_t
1008     elif varX == "deadtime":
1009     VX = deadtime_t
1010     VXE = e_deadtime_t
1011     elif varX == "rawrate":
1012     VX = rawrate_t
1013     VXE = e_rawrate_t
1014     elif varX == "rate":
1015     VX = rate_t
1016     VXE = e_rate_t
1017     elif varX == "rawxsec":
1018     VX = rawxsec_t
1019     VXE = e_rawxsec_t
1020     elif varX == "xsec":
1021     VX = xsec_t
1022     VXE = e_xsec_t
1023     elif varX == "psi":
1024     VX = psi_t
1025     VXE = e_psi_t
1026     else:
1027     print "No valid variable entered for X"
1028     continue
1029     if varY == "run":
1030     VY = run_t
1031     VYE = run_t_e
1032     elif varY == "ls":
1033     VY = ls_t
1034     VYE = e_ls_t
1035     elif varY == "ps":
1036     VY = ps_t
1037     VYE = e_ps_t
1038     elif varY == "inst":
1039     VY = inst_t
1040     VYE = e_inst_t
1041     elif varY == "live":
1042     VY = live_t
1043     VYE = e_live_t
1044     elif varY == "delivered":
1045     VY = delivered_t
1046     VYE = e_delivered_t
1047     elif varY == "deadtime":
1048     VY = deadtime_t
1049     VYE = e_deadtime_t
1050     elif varY == "rawrate":
1051     VY = rawrate_t
1052     VYE = e_rawrate_t
1053     if fit_file:
1054     VF = rawrate_fit_t
1055     VFE = e_rawrate_fit_t
1056     elif varY == "rate":
1057     VY = rate_t
1058     VYE = e_rate_t
1059     if fit_file:
1060     VF = rate_fit_t
1061     VFE = e_rate_fit_t
1062     elif varY == "rawxsec":
1063     VY = rawxsec_t
1064     VYE = e_rawxsec_t
1065     if fit_file:
1066     VF = rawxsec_fit_t
1067     VFE = e_rawxsec_fit_t
1068     elif varY == "xsec":
1069     VY = xsec_t
1070     VYE = e_xsec_t
1071     if fit_file:
1072     VF = xsec_fit_t
1073     VFE = e_xsec_fit_t
1074     elif varY == "psi":
1075     VY = psi_t
1076     VYE = e_psi_t
1077     else:
1078     print "No valid variable entered for Y"
1079     continue
1080    
1081     return [VX, VXE, VY, VYE, VF, VFE]
1082    
1083 grchrist 1.11 def pass_cuts(data_clean, realvalue, prediction, meanxsec, Rates, print_trigger, iterator, num_ls,LumiPageInfo,SubSystemOff, max_dt, print_info, trig_list):
1084 grchrist 1.17 it_offset=0
1085 grchrist 1.15 Passed=True
1086     subsystemfailed=[]
1087 grchrist 1.11
1088 grchrist 1.8 if num_ls==1:
1089     ##fit is 2 ls ahead of real rate
1090 grchrist 1.17
1091 grchrist 1.8
1092 grchrist 1.10 LS=Rates[print_trigger]["ls"][iterator]
1093     #print "ls=",LS,
1094     LSRange=LumiPageInfo[LS]["LSRange"]
1095     #print LSRange,
1096     LS2=LSRange[-1]
1097     #LS2=LSRange.pop()
1098     #print "LS2=",LS2
1099    
1100     #print LumiPageInfo[LS]
1101     lumidict={}
1102     lumidict=LumiPageInfo[LS]
1103 grchrist 1.15
1104 grchrist 1.11
1105    
1106    
1107    
1108     if print_info:
1109     if (iterator==0 and print_trigger==trig_list[0]):
1110     print '%10s%10s%10s%10s%10s%10s%10s%15s%20s' % ("Status", "Run", "LS", "Physics", "Active", "Deadtime", " MaxDeadTime", " Passed all subsystems?", " List of Subsystems failed")
1111    
1112     ## if SubSystemOff["All"]:
1113     ## for keys in LumiPageInfo[LS]:
1114     ## #print LS, keys, LumiPageInfo[LS][keys]
1115     ## if not LumiPageInfo[LS][keys]:
1116     ## Passed=False
1117     ## subsystemfailed.append(keys)
1118     ## break
1119     ## else:
1120     if SubSystemOff["Mu"] or SubSystemOff["All"]:
1121     if not (LumiPageInfo[LS]["rpc"] or LumiPageInfo[LS]["dt0"] or LumiPageInfo[LS]["dtp"] or LumiPageInfo[LS]["dtm"] or LumiPageInfo["cscp"] or LumiPageInfo["cscm"]):
1122     Passed=False
1123     subsystemfailed.append("Mu")
1124     if SubSystemOff["HCal"] or SubSystemOff["All"]:
1125     if not (LumiPageInfo[LS]["hbhea"] and LumiPageInfo[LS]["hbheb"] and LumiPageInfo[LS]["hbhec"]):
1126     Passed=False
1127     subsystemfailed.append("HCal")
1128     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["hf"]):
1129     Passed=False
1130     subsystemfailed.append("HCal-EndCap")
1131     if SubSystemOff["ECal"] or SubSystemOff["All"]:
1132     if not (LumiPageInfo[LS]["ebp"] and LumiPageInfo[LS]["ebm"]):
1133     Passed=False
1134     subsystemfailed.append("ECal")
1135     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["eep"] and LumiPageInfo[LS]["eem"] and LumiPageInfo[LS]["esp"] or LumiPageInfo[LS]["esm"]):
1136     Passed=False
1137     subsystemfailed.append("ECal-EndCap")
1138     if SubSystemOff["Tracker"] or SubSystemOff["All"]:
1139     if not (LumiPageInfo[LS]["tob"] and LumiPageInfo[LS]["tibtid"] and LumiPageInfo[LS]["bpix"] and LumiPageInfo[LS]["fpix"]):
1140     Passed=False
1141     subsystemfailed.append("Tracker")
1142     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["tecp"] and LumiPageInfo[LS]["tecm"]):
1143     Passed=False
1144     subsystemfailed.append("Tracker-EndCap")
1145     if SubSystemOff["Beam"] or SubSystemOff["All"]:
1146     if not(LumiPageInfo[LS]["b1pres"] and LumiPageInfo[LS]["b2pres"] and LumiPageInfo[LS]["b1stab"] and LumiPageInfo[LS]["b2stab"]):
1147     Passed=False
1148     subsystemfailed.append("Beam")
1149 grchrist 1.12 else:
1150 grchrist 1.17
1151 grchrist 1.12 Passed=True
1152 grchrist 1.10
1153 grchrist 1.17 #print "LS",LS, Passed, round(Rates[print_trigger]["deadtime"][iterator],2), max_dt
1154 grchrist 1.10
1155 grchrist 1.8 if not data_clean or (
1156 grchrist 1.11 ## (
1157 grchrist 1.8
1158 grchrist 1.11 ## (
1159     ## realvalue > 0.4*prediction
1160     ## and realvalue < 2.5*prediction
1161     ## )
1162 grchrist 1.8
1163 grchrist 1.11 ## or
1164 grchrist 1.8
1165 grchrist 1.11 ## (
1166     ## realvalue > 0.4*meanxsec
1167     ## and realvalue < 2.5*meanxsec
1168     ## )
1169 grchrist 1.8
1170 grchrist 1.11 ## or prediction < 0
1171     ## )
1172 grchrist 1.8
1173 grchrist 1.11 Rates[print_trigger]["physics"][iterator] == 1
1174 grchrist 1.8 and Rates[print_trigger]["active"][iterator] == 1
1175 grchrist 1.17 and Rates[print_trigger]["deadtime"][iterator] < max_dt
1176 grchrist 1.11 #and Rates[print_trigger]["psi"][iterator] > 0
1177 grchrist 1.10 and Passed
1178 grchrist 1.8 ):
1179 grchrist 1.11 #print LS, "True"
1180 grchrist 1.8 return True
1181     else:
1182 grchrist 1.11
1183 grchrist 1.13 if (print_info and print_trigger==trig_list[0] and num_ls==1):
1184 grchrist 1.11
1185 grchrist 1.17 print '%10s%10s%10s%10s%10s%10s%10s%15s%20s' % ("Failed", Rates[print_trigger]["run"][iterator], LS, Rates[print_trigger]["physics"][iterator], Rates[print_trigger]["active"][iterator], round(Rates[print_trigger]["deadtime"][iterator],2), max_dt, Passed, subsystemfailed)
1186 grchrist 1.15 ##elif(print_info and print_trigger==trig_list[0]):
1187 grchrist 1.17 ## print '%10s%10s%10s%10s%10s%10s%10s%15s%20s' % ("Failed", Rates[print_trigger]["run"][iterator], LumiPageInfo[Rates[print_trigger]["ls"][iterator]]["LSRange"], Rates[print_trigger]["physics"][iterator], Rates[print_trigger]["active"][iterator], round(Rates[print_trigger]["deadtime"][iterator],2), max_dt, Passed, subsystemfailed)
1188 grchrist 1.8 return False
1189 abrinke1 1.5
1190 grchrist 1.10
1191     #### LumiRangeGreens ####
1192     ####inputs: RefMoreLumiArray --dict with lumi page info in LS by LS blocks,
1193     #### LRange --list range over lumis,
1194     #### nls --number of lumisections
1195     #### RefRunNum --run number
1196     ####
1197     ####outputs RangeMoreLumi --lumi page info in dict LSRange blocks with lumi, added items Run and LSRange
1198 grchrist 1.17 def LumiRangeGreens(RefMoreLumiArray,LSRange,nls,RefRunNum,deadtimebeamactive):
1199 grchrist 1.9
1200     RangeMoreLumi={}
1201     for keys,values in RefMoreLumiArray.iteritems():
1202     RangeMoreLumi[keys]=1
1203 grchrist 1.10
1204 grchrist 1.9 for iterator in LSRange[nls]:
1205     for keys, values in RefMoreLumiArray.iteritems():
1206     if RefMoreLumiArray[keys][iterator]==0:
1207     RangeMoreLumi[keys]=0
1208 grchrist 1.10 RangeMoreLumi['LSRange']=LSRange[nls]
1209     RangeMoreLumi['Run']=RefRunNum
1210 grchrist 1.17 RangeMoreLumi['DeadTimeBeamActive']=deadtimebeamactive
1211 grchrist 1.9 return RangeMoreLumi
1212    
1213 grchrist 1.10 #### CheckLumis ####
1214     ####inputs:
1215     #### PageLumiInfo --dict of LS with dict of some lumipage info
1216     #### Rates --dict of triggernames with dict of info
1217     def checkLS(Rates, PageLumiInfo,trig_list):
1218     rateslumis=Rates[trig_list[-1]]["ls"]
1219     keys=PageLumiInfo.keys()
1220     print "lumi run=",PageLumiInfo[keys[-1]]["Run"]
1221     ll=0
1222     for ls in keys:
1223     print ls,rateslumis[ll]
1224     ll=ll+1
1225     return False
1226    
1227    
1228 grchrist 1.9
1229 abrinke1 1.1 if __name__=='__main__':
1230     main()