ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.24
Committed: Mon Mar 26 16:19:34 2012 UTC (13 years, 1 month ago) by grchrist
Content type: text/x-python
Branch: MAIN
Changes since 1.23: +4 -1 lines
Log Message:
fixed bug where we couldn't find fit file

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