ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.21
Committed: Mon Mar 26 12:40:43 2012 UTC (13 years, 1 month ago) by grchrist
Content type: text/x-python
Branch: MAIN
Changes since 1.20: +79 -17 lines
Log Message:
a bit more user friendly

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