ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.73
Committed: Fri Nov 16 10:48:53 2012 UTC (12 years, 5 months ago) by grchrist
Content type: text/x-python
Branch: MAIN
Changes since 1.72: +7 -5 lines
Log Message:
added info to parser to get L1 seed changes in correct format

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 awoodard 1.47 from StreamMonitor import StreamMonitor
11 awoodard 1.72 from itertools import groupby
12     from operator import itemgetter
13 grchrist 1.29
14 abrinke1 1.1 from ROOT import gROOT, TCanvas, TF1, TGraph, TGraphErrors, TPaveStats, gPad, gStyle
15 amott 1.18 from ROOT import TFile, TPaveText, TBrowser
16 abrinke1 1.1 from ROOT import gBenchmark
17     import array
18     import math
19 grchrist 1.8 from ReadConfig import RateMonConfig
20 muell149 1.46 from TablePrint import *
21 abrinke1 1.5 from selectionParser import selectionParser
22    
23 amott 1.18 def usage():
24     print sys.argv[0]+" [options] <list of runs>"
25     print "This script is used to generate fits and do secondary shifter validation"
26     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"
27     print " be careful with using ranges (c-d), it is highly recommended to use a JSON in this case"
28     print "options: "
29     print "--makeFits run in fit making mode"
30     print "--secondary run in secondary shifter mode"
31 grchrist 1.28 print "--TMD put in TMD predictions"
32 amott 1.18 print "--fitFile=<path> path to the fit file"
33     print "--json=<path> path to the JSON file"
34     print "--TriggerList=<path> path to the trigger list (without versions!)"
35 grchrist 1.26 print "--maxdt=<max deadtime> Mask LS above max deadtime threshold"
36     print "--All Mask LS with any red LS on WBM LS page (not inc castor zdc etc)"
37     print "--Mu Mask LS with Mu off"
38     print "--HCal Mask LS with HCal barrel off"
39     print "--Tracker Mask LS with Tracker barrel off"
40     print "--ECal Mask LS with ECal barrel off"
41     print "--EndCap Mask LS with EndCap sys off, used in combination with other subsys"
42     print "--Beam Mask LS with Beam off"
43 grchrist 1.34 print "--NoVersion Ignore version numbers"
44 grchrist 1.35 print "--linear Force Linear fits"
45 grchrist 1.38 print "--inst Fits using inst not delivered"
46     print "--TMDerr Use errors from TMD predictions"
47 muell149 1.46 print "--write Writes fit info into csv, for ranking nonlinear triggers"
48 awoodard 1.50 print "--AllTriggers Run for all triggers instead of specifying a trigger list"
49 muell149 1.46
50 amott 1.18 class Modes:
51 grchrist 1.29 none,fits,secondary = range(3)
52    
53     def pickYear():
54     global thisyear
55 grchrist 1.32 thisyear="2012"
56 grchrist 1.29 print "Year set to ",thisyear
57    
58 abrinke1 1.1 def main():
59 amott 1.18 try:
60 grchrist 1.32 ##set year to 2012
61 grchrist 1.29 pickYear()
62 awoodard 1.47
63 grchrist 1.22 try:
64 grchrist 1.61 opt, args = getopt.getopt(sys.argv[1:],"",["makeFits","secondary","fitFile=","json=","TriggerList=","maxdt=","All","Mu","HCal","Tracker","ECal","EndCap","Beam","NoVersion","linear","inst","TMDerr","write","AllTriggers","UsePSCol="])
65 grchrist 1.22
66     except getopt.GetoptError, err:
67     print str(err)
68     usage()
69     sys.exit(2)
70    
71 grchrist 1.25 ##### RUN LIST ########
72 grchrist 1.22 run_list=[]
73 grchrist 1.23
74 grchrist 1.22 if len(args)<1:
75     inputrunlist=[]
76     print "No runs specified"
77     runinput=raw_input("Enter run range in form <run1> <run2> <run3> or <run1>-<run2>:")
78     inputrunlist.append(runinput)
79    
80     if runinput.find(' ')!=-1:
81     args=runinput.split(' ')
82     else:
83     args.append(runinput)
84    
85     for r in args:
86     if r.find('-')!=-1: # r is a run range
87     rrange = r.split('-')
88     if len(rrange)!=2:
89     print "Invalid run range %s" % (r,)
90     sys.exit(0)
91     try:
92     for rr in range(int(rrange[0]),int(rrange[1])+1):
93     run_list.append(rr)
94     except:
95     print "Invalid run range %s" % (r,)
96     sys.exit(0)
97     else: # r is not a run range
98     try:
99     run_list.append(int(r))
100     except:
101     print "Invalid run %s" % (r,)
102 grchrist 1.21
103 grchrist 1.25 ##### READ CMD LINE ARGS #########
104 grchrist 1.22 mode = Modes.none
105     fitFile = ""
106     jsonfile = ""
107     trig_list = []
108 grchrist 1.25 max_dt=-1.0
109     subsys=-1.0
110 grchrist 1.61 NoVersion=True
111 grchrist 1.35 linear=False
112 grchrist 1.38 do_inst=False
113     TMDerr=False
114 muell149 1.46 wp_bool=False
115 grchrist 1.59 all_triggers=False
116 grchrist 1.73 DoL1=False
117 grchrist 1.61 UsePSCol=-1
118 grchrist 1.25 SubSystemOff={'All':False,'Mu':False,'HCal':False,'ECal':False,'Tracker':False,'EndCap':False,'Beam':False}
119 grchrist 1.22 for o,a in opt:
120     if o == "--makeFits":
121     mode = Modes.fits
122     elif o == "--secondary":
123     mode = Modes.secondary
124     elif o == "--fitFile":
125     fitFile = str(a)
126     elif o == "--json":
127     jsonfile = a
128 grchrist 1.26 elif o=="--maxdt":
129     max_dt = float(a)
130 grchrist 1.25 elif o=="--All":
131     subsys=1
132     SubSystemOff["All"]=True
133     elif o=="--Mu":
134     subsys=1
135     SubSystemOff["Mu"]=True
136     elif o=="--HCal":
137     SubSystemOff["HCal"]=True
138     subsys=1
139     elif o=="--Tracker":
140     SubSystemOff["Tracker"]=True
141     subsys=1
142 grchrist 1.26 elif o=="--ECal":
143     SubSystemOff["ECal"]=True
144     subsys=1
145 grchrist 1.25 elif o=="--EndCap":
146     SubSystemOff["EndCap"]=True
147     subsys=1
148     elif o=="--Beam":
149     SubSystemOff["Beam"]=True
150     subsys=1
151 grchrist 1.34 elif o=="--NoVersion":
152     NoVersion=True
153 grchrist 1.35 elif o=="--linear":
154     linear=True
155 grchrist 1.38 elif o=="--inst":
156     do_inst=True
157     elif o=="--TMDerr":
158     TMDerr=True
159 muell149 1.46 elif o=="--write":
160     wp_bool=True
161 awoodard 1.50 elif o=="--AllTriggers":
162 grchrist 1.61 all_triggers=True
163     elif o=="--UsePSCol":
164     UsePSCol=int(a)
165 grchrist 1.22 elif o == "--TriggerList":
166     try:
167     f = open(a)
168     for entry in f:
169     if entry.startswith('#'):
170     continue
171     if entry.find(':')!=-1:
172     entry = entry[:entry.find(':')] ## We can point this to the existing monitor list, just remove everything after ':'!
173     if entry.find('#')!=-1:
174     entry = entry[:entry.find('#')] ## We can point this to the existing monitor list, just remove everything after ':'!
175     trig_list.append( entry.rstrip('\n'))
176     except:
177     print "\nInvalid Trigger List\n"
178     sys.exit(0)
179     else:
180     print "\nInvalid Option %s\n" % (str(o),)
181     usage()
182     sys.exit(2)
183    
184     print "\n\n"
185 grchrist 1.25 ###### MODES #########
186    
187 grchrist 1.22 if mode == Modes.none: ## no mode specified
188     print "\nNo operation mode specified!\n"
189     modeinput=raw_input("Enter mode, --makeFits or --secondary:")
190     print "modeinput=",modeinput
191     if not (modeinput=="--makeFits" or modeinput=="--secondary"):
192     print "not either"
193     usage()
194 amott 1.18 sys.exit(0)
195 grchrist 1.22 elif modeinput == "--makeFits":
196     mode=Modes.fits
197     elif modeinput =="--secondary":
198     mode=Modes.secondary
199     else:
200     print "FATAL ERROR: No Mode specified"
201 amott 1.18 sys.exit(0)
202 grchrist 1.22
203     if mode == Modes.fits:
204     print "Running in Fit Making mode\n\n"
205     elif mode == Modes.secondary:
206     print "Running in Secondary Shifter mode\n\n"
207     else: ## should never get here, but exit if we do
208 grchrist 1.21 print "FATAL ERROR: No Mode specified"
209     sys.exit(0)
210 grchrist 1.22
211     if fitFile=="" and not mode==Modes.fits:
212     print "\nPlease specify fit file. These are available:\n"
213 grchrist 1.29 path="Fits/%s/" % (thisyear) # insert the path to the directory of interest
214 grchrist 1.22 dirList=os.listdir(path)
215     for fname in dirList:
216     print fname
217 grchrist 1.24 fitFile = path+raw_input("Enter fit file in format Fit_HLT_10LS_Run176023to180252.pkl: ")
218    
219 grchrist 1.22 elif fitFile=="":
220 grchrist 1.39 NoVstr=""
221     if NoVersion:
222     NoVstr="NoV_"
223 grchrist 1.38 if not do_inst:
224 grchrist 1.39 fitFile="Fits/%s/Fit_HLT_%s10LS_Run%sto%s.pkl" % (thisyear,NoVstr,min(run_list),max(run_list))
225 grchrist 1.38 else:
226 grchrist 1.39 fitFile="Fits/%s/Fit_inst_HLT_%s10LS_Run%sto%s.pkl" % (thisyear,NoVstr,min(run_list),max(run_list))
227 grchrist 1.26
228 grchrist 1.39 if "NoV" in fitFile:
229     NoVersion=True
230 grchrist 1.25
231     ###### TRIGGER LIST #######
232 grchrist 1.24
233 awoodard 1.50 if trig_list == [] and not all_triggers:
234 grchrist 1.22
235     print "\nPlease specify list of triggers\n"
236     print "Available lists are:"
237     dirList=os.listdir(".")
238     for fname in dirList:
239     entry=fname
240     if entry.find('.')!=-1:
241     extension = entry[entry.find('.'):] ## We can point this to the existing monitor list, just remove everything after ':'!
242     if extension==".list":
243     print fname
244 grchrist 1.26 trig_input=raw_input("\nEnter triggers in format HLT_IsoMu30_eta2p1 or a .list file: ")
245 grchrist 1.21
246 grchrist 1.22 if trig_input.find('.')!=-1:
247     extension = trig_input[trig_input.find('.'):]
248 grchrist 1.21 if extension==".list":
249 grchrist 1.22 try:
250     fl=open(trig_input)
251     except:
252     print "Cannot open file"
253     usage()
254     sys.exit(0)
255 grchrist 1.21
256 grchrist 1.22 for line in fl:
257     if line.startswith('#'):
258     continue
259     if len(line)<1:
260     continue
261 grchrist 1.21
262 grchrist 1.22 if len(line)>=2:
263     arg=line.rstrip('\n').rstrip(' ').lstrip(' ')
264     trig_list.append(arg)
265     else:
266     arg=''
267     else:
268     trig_list.append(trig_input)
269 grchrist 1.21
270 abrinke1 1.5 ## Can use any combination of LowestRunNumber, HighestRunNumber, and NumberOfRuns -
271     ## just modify "ExistingRuns.sort" and for "run in ExistingRuns" accordingly
272    
273 grchrist 1.22 if jsonfile=="":
274     JSON=[]
275     else:
276     print "Using JSON: %s" % (jsonfile,)
277     JSON = GetJSON(jsonfile) ##Returns array JSON[runs][ls_list]
278 abrinke1 1.5
279 grchrist 1.22 ###### TO CREATE FITS #########
280     if mode == Modes.fits:
281     trig_name = "HLT"
282 grchrist 1.34 num_ls = 10
283 grchrist 1.22 physics_active_psi = True ##Requires that physics and active be on, and that the prescale column is not 0
284     #JSON = [] ##To not use a JSON file, just leave the array empty
285     debug_print = False
286     no_versions=False
287 grchrist 1.34 min_rate = 0.0
288 grchrist 1.22 print_table = False
289     data_clean = True ##Gets rid of anomalous rate points, reqires physics_active_psi (PAP) and deadtime < 20%
290     ##plot_properties = [varX, varY, do_fit, save_root, save_png, fit_file]
291 grchrist 1.38 if not do_inst:
292     plot_properties = [["delivered", "rate", True, True, False, fitFile]]
293 awoodard 1.54 # plot_properties = [["delivered", "rawrate", True, True, False, fitFile]]
294 grchrist 1.38 else:
295     plot_properties = [["inst", "rate", True, True, False, fitFile]]
296 amott 1.18
297 grchrist 1.34 masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_Zero"]
298 grchrist 1.22 save_fits = True
299 grchrist 1.25 if max_dt==-1.0:
300     max_dt=0.08 ## no deadtime cutuse 2.0
301 grchrist 1.22 force_new=True
302     print_info=True
303 grchrist 1.25 if subsys==-1.0:
304     SubSystemOff={'All':False,'Mu':False,'HCal':False,'ECal':False,'Tracker':False,'EndCap':False,'Beam':True}
305 grchrist 1.22
306     ###### TO SEE RATE VS PREDICTION ########
307     if mode == Modes.secondary:
308     trig_name = "HLT"
309     num_ls = 1
310     physics_active_psi = True
311     debug_print = False
312     no_versions=False
313 grchrist 1.34 min_rate = 0.0
314 grchrist 1.22 print_table = False
315     data_clean = True
316     ##plot_properties = [varX, varY, do_fit, save_root, save_png, fit_file]
317 grchrist 1.70 plot_properties = [["ls", "rate", False, True, False,fitFile]]
318 awoodard 1.52 ## rate is calculated as: (measured rate, deadtime corrected) * prescale [prediction not dt corrected]
319     ## rawrate is calculated as: measured rate [prediction is dt corrected]
320 grchrist 1.22
321 grchrist 1.34 masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_Zero"]
322 grchrist 1.22 save_fits = False
323 grchrist 1.25 if max_dt==-1.0:
324     max_dt=2.0 ## no deadtime cut=2.0
325 grchrist 1.22 force_new=True
326     print_info=True
327 grchrist 1.25 if subsys==-1.0:
328     SubSystemOff={'All':True,'Mu':False,'HCal':False,'ECal':False,'Tracker':False,'EndCap':False,'Beam':True}
329 grchrist 1.13
330 grchrist 1.26 for k in SubSystemOff.iterkeys():
331     print k,"=",SubSystemOff[k]," ",
332     print " "
333 grchrist 1.59
334 grchrist 1.22 ######## END PARAMETERS - CALL FUNCTIONS ##########
335 grchrist 1.61 [Rates,LumiPageInfo, L1_trig_list]= GetDBRates(run_list, trig_name, trig_list, num_ls, max_dt, physics_active_psi, JSON, debug_print, force_new, SubSystemOff,NoVersion,all_triggers, DoL1,UsePSCol)
336 grchrist 1.59 if DoL1:
337     trig_list=L1_trig_list
338 grchrist 1.73 ##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,NoVersion, linear, do_inst, TMDerr,wp_bool,all_triggers)
339 awoodard 1.47
340 grchrist 1.22 except KeyboardInterrupt:
341     print "Wait... come back..."
342 abrinke1 1.1
343 awoodard 1.47
344 grchrist 1.61 def GetDBRates(run_list,trig_name,trig_list, num_ls, max_dt, physics_active_psi,JSON,debug_print, force_new, SubSystemOff,NoVersion,all_triggers, DoL1,UsePSCol):
345 abrinke1 1.1
346     Rates = {}
347 grchrist 1.10 LumiPageInfo={}
348 abrinke1 1.5 ## Save in RefRuns with name dependent on trig_name, num_ls, JSON, and physics_active_psi
349     if JSON:
350 amott 1.18 #print "Using JSON file"
351 abrinke1 1.5 if physics_active_psi:
352 grchrist 1.29 RefRunNameTemplate = "RefRuns/%s/Rates_%s_%sLS_JPAP.pkl"
353 abrinke1 1.5 else:
354 grchrist 1.29 RefRunNameTemplate = "RefRuns/%s/Rates_%s_%sLS_JSON.pkl"
355 abrinke1 1.5 else:
356 grchrist 1.14 print "Using Physics and Active ==1"
357 abrinke1 1.5 if physics_active_psi:
358 grchrist 1.29 RefRunNameTemplate = "RefRuns/%s/Rates_%s_%sLS_PAP.pkl"
359 abrinke1 1.5 else:
360 grchrist 1.29 RefRunNameTemplate = "RefRuns/%s/Rates_%s_%sLS.pkl"
361 grchrist 1.10
362 grchrist 1.29 RefRunFile = RefRunNameTemplate % (thisyear,trig_name,num_ls)
363     RefRunFileHLT = RefRunNameTemplate % (thisyear,"HLT",num_ls)
364 grchrist 1.10
365     print "RefRun=",RefRunFile
366     print "RefRunFileHLT",RefRunFileHLT
367 grchrist 1.11 if not force_new:
368     try: ##Open an existing RefRun file with the same parameters and trigger name
369     pkl_file = open(RefRunFile, 'rb')
370     Rates = pickle.load(pkl_file)
371     pkl_file.close()
372     os.remove(RefRunFile)
373     print "using",RefRunFile
374    
375 abrinke1 1.5 except:
376 grchrist 1.11 try: ##Open an existing RefRun file with the same parameters and HLT for trigger name
377     pkl_file = open(RefRunFileHLT)
378     HLTRates = pickle.load(pkl_file)
379     for key in HLTRates:
380     if trig_name in str(key):
381     Rates[key] = HLTRates[key]
382     #print str(RefRunFile)+" does not exist. Creating ..."
383     except:
384     print str(RefRunFile)+" does not exist. Creating ..."
385 grchrist 1.10
386     ## try the lumis file
387 grchrist 1.29 RefLumiNameTemplate = "RefRuns/%s/Lumis_%s_%sLS.pkl"
388     RefLumiFile= RefLumiNameTemplate % (thisyear,"HLT",num_ls)
389 grchrist 1.11 if not force_new:
390     try:
391     pkl_lumi_file = open(RefLumiFile, 'rb')
392     LumiPageInfo = pickle.load(pkl_lumi_file)
393     pkl_lumi_file.close()
394     os.remove(RefLumiFile)
395     print "using",RefLumiFile
396     except:
397     print str(RefLumiFile)+" doesn't exist. Make it..."
398 grchrist 1.34
399     trig_list_noV=[]
400     for trigs in trig_list:
401     trig_list_noV.append(StripVersion(trigs))
402 grchrist 1.40
403 grchrist 1.34 if NoVersion:
404     trig_list=trig_list_noV
405    
406 abrinke1 1.5 for RefRunNum in run_list:
407     if JSON:
408     if not RefRunNum in JSON:
409     continue
410 awoodard 1.47 try:
411 abrinke1 1.1 ExistsAlready = False
412     for key in Rates:
413     if RefRunNum in Rates[key]["run"]:
414     ExistsAlready = True
415     break
416 grchrist 1.10
417     LumiExistsLAready=False
418     for v in LumiPageInfo.itervalues():
419     if RefRunNum == v["Run"]:
420     LumiExistsAlready=True
421     break
422     if ExistsAlready and LumiExistsAlready:
423 abrinke1 1.1 continue
424 grchrist 1.10
425 abrinke1 1.1 except:
426     print "Getting info for run "+str(RefRunNum)
427    
428     if RefRunNum < 1:
429     continue
430 grchrist 1.30 ColRunNum,isCol,isGood = GetLatestRunNumber(RefRunNum)
431     if not isGood:
432     print "Run ",RefRunNum, " is not Collisions"
433    
434     continue
435 awoodard 1.50
436 grchrist 1.26 if not isCol:
437     print "Run ",RefRunNum, " is not Collisions"
438    
439     continue
440    
441 grchrist 1.17 print "calculating rates and green lumis for run ",RefRunNum
442 grchrist 1.73
443 abrinke1 1.5 if True: ##Placeholder
444 abrinke1 1.1 if True: #May replace with "try" - for now it's good to know when problems happen
445     RefParser = DatabaseParser()
446     RefParser.RunNumber = RefRunNum
447     RefParser.ParseRunSetup()
448 abrinke1 1.5 RefLumiRangePhysicsActive = RefParser.GetLSRange(1,9999) ##Gets array of all LS with physics and active on
449     RefLumiArray = RefParser.GetLumiInfo() ##Gets array of all existing LS and their lumi info
450 abrinke1 1.2 RefLumiRange = []
451 grchrist 1.17 RefMoreLumiArray = RefParser.GetMoreLumiInfo()#dict with keys as bits from lumisections WBM page and values are dicts with key=LS:value=bit
452 grchrist 1.73 L1HLTseeds=RefParser.GetL1HLTseeds()
453     HLTL1PS=RefParser.GetL1PSbyseed()
454 grchrist 1.59 if DoL1:
455 grchrist 1.73
456 grchrist 1.59 for HLTkey in trig_list:
457     #print name
458     if "L1" in HLTkey:
459     continue
460     else:
461     try:
462     for L1seed in L1HLTseeds[HLTkey]:
463     if L1seed not in trig_list:
464     trig_list.append(L1seed)
465     except:
466     pass
467 abrinke1 1.5 for iterator in RefLumiArray[0]: ##Makes array of LS with proper PAP and JSON properties
468 grchrist 1.11 ##cheap way of getting PSCol None-->0
469 amott 1.18 if RefLumiArray[0][iterator] not in range(1,9):
470 grchrist 1.11 RefLumiArray[0][iterator]=0
471 grchrist 1.15
472 grchrist 1.61 if not UsePSCol==-1:
473     if not RefLumiArray[0][iterator]==UsePSCol:
474     print "skipping LS",iterator
475     continue
476    
477 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):
478 abrinke1 1.5 if not JSON or RefRunNum in JSON:
479     if not JSON or iterator in JSON[RefRunNum]:
480     RefLumiRange.append(iterator)
481 grchrist 1.9
482 abrinke1 1.5 try:
483     nls = RefLumiRange[0]
484     LSRange = {}
485     except:
486     print "Run "+str(RefRunNum)+" has no good LS"
487     continue
488     if num_ls > len(RefLumiRange):
489 abrinke1 1.1 print "Run "+str(RefRunNum)+" is too short: from "+str(nls)+" to "+str(RefLumiRange[-1])+", while num_ls = "+str(num_ls)
490     continue
491     while nls < RefLumiRange[-1]-num_ls:
492 abrinke1 1.5 LSRange[nls] = []
493     counter = 0
494     for iterator in RefLumiRange:
495     if iterator >= nls and counter < num_ls:
496     LSRange[nls].append(iterator)
497     counter += 1
498 abrinke1 1.1 nls = LSRange[nls][-1]+1
499 grchrist 1.73 checkL1seedChangeALLPScols(trig_list,HLTL1PS) #for L1prescale changes
500 grchrist 1.17 #print "Run "+str(RefRunNum)+" contains LS from "+str(min(LSRange))+" to "+str(max(LSRange))
501 grchrist 1.8 for nls in sorted(LSRange.iterkeys()):
502 abrinke1 1.1 TriggerRates = RefParser.GetHLTRates(LSRange[nls])
503 grchrist 1.59 #L1Rate=RefParser.GetDeadTimeBeamActive(LSRange[nls])
504 awoodard 1.47
505     ## Clumsy way to append Stream A. Should choose correct method for calculating stream a based on ps column used in data taking.
506 grchrist 1.59
507 awoodard 1.57 if ('HLT_Stream_A' in trig_list) or all_triggers:
508     config = RateMonConfig(os.path.abspath(os.path.dirname(sys.argv[0])))
509     config.ReadCFG()
510     stream_mon = StreamMonitor()
511     core_a_rates = stream_mon.getStreamACoreRatesByLS(RefParser,LSRange[nls],config).values()
512     avg_core_a_rate = sum(core_a_rates)/len(LSRange[nls])
513     TriggerRates['HLT_Stream_A'] = [1,1,avg_core_a_rate,avg_core_a_rate]
514 grchrist 1.59
515     if DoL1:
516     L1RatesALL=RefParser.GetL1RatesALL(LSRange[nls])
517    
518     for L1seed in L1RatesALL.iterkeys():
519     TriggerRates[L1seed]=L1RatesALL[L1seed]
520 abrinke1 1.5
521     [inst, live, delivered, dead, pscols] = RefParser.GetAvLumiInfo(LSRange[nls])
522 grchrist 1.17 deadtimebeamactive=RefParser.GetDeadTimeBeamActive(LSRange[nls])
523 grchrist 1.59
524 abrinke1 1.2 physics = 1
525     active = 1
526 abrinke1 1.5 psi = 99
527     for iterator in LSRange[nls]: ##Gets lowest value of physics, active, and psi in the set of lumisections
528 abrinke1 1.2 if RefLumiArray[5][iterator] == 0:
529     physics = 0
530     if RefLumiArray[6][iterator] == 0:
531     active = 0
532 abrinke1 1.5 if RefLumiArray[0][iterator] < psi:
533     psi = RefLumiArray[0][iterator]
534 abrinke1 1.2
535 grchrist 1.17 MoreLumiMulti=LumiRangeGreens(RefMoreLumiArray,LSRange,nls,RefRunNum,deadtimebeamactive)
536 grchrist 1.10
537 abrinke1 1.5 if inst < 0 or live < 0 or delivered < 0:
538     print "Run "+str(RefRunNum)+" LS "+str(nls)+" inst lumi = "+str(inst)+" live lumi = "+str(live)+", delivered = "+str(delivered)+", physics = "+str(physics)+", active = "+str(active)
539 grchrist 1.10
540     LumiPageInfo[nls]=MoreLumiMulti
541    
542 abrinke1 1.1 for key in TriggerRates:
543 grchrist 1.59
544 grchrist 1.34 if NoVersion:
545     name = StripVersion(key)
546     else:
547     name=key
548 grchrist 1.59
549 grchrist 1.26 if not name in trig_list:
550 awoodard 1.50 if all_triggers:
551     trig_list.append(name)
552     else:
553     continue
554 grchrist 1.59
555 abrinke1 1.1 if not Rates.has_key(name):
556     Rates[name] = {}
557     Rates[name]["run"] = []
558     Rates[name]["ls"] = []
559     Rates[name]["ps"] = []
560     Rates[name]["inst_lumi"] = []
561     Rates[name]["live_lumi"] = []
562     Rates[name]["delivered_lumi"] = []
563     Rates[name]["deadtime"] = []
564     Rates[name]["rawrate"] = []
565     Rates[name]["rate"] = []
566     Rates[name]["rawxsec"] = []
567     Rates[name]["xsec"] = []
568 abrinke1 1.2 Rates[name]["physics"] = []
569     Rates[name]["active"] = []
570 abrinke1 1.5 Rates[name]["psi"] = []
571 grchrist 1.9
572 abrinke1 1.1 [avps, ps, rate, psrate] = TriggerRates[key]
573     Rates[name]["run"].append(RefRunNum)
574     Rates[name]["ls"].append(nls)
575     Rates[name]["ps"].append(ps)
576     Rates[name]["inst_lumi"].append(inst)
577     Rates[name]["live_lumi"].append(live)
578     Rates[name]["delivered_lumi"].append(delivered)
579 grchrist 1.17 Rates[name]["deadtime"].append(deadtimebeamactive)
580 abrinke1 1.1 Rates[name]["rawrate"].append(rate)
581 abrinke1 1.2 if live == 0:
582 awoodard 1.56 Rates[name]["rate"].append(0.0)
583 abrinke1 1.2 Rates[name]["rawxsec"].append(0.0)
584     Rates[name]["xsec"].append(0.0)
585     else:
586 awoodard 1.56 try:
587     Rates[name]["rate"].append(psrate/(1.0-deadtimebeamactive))
588     except:
589     Rates[name]["rate"].append(0.0)
590 abrinke1 1.2 Rates[name]["rawxsec"].append(rate/live)
591     Rates[name]["xsec"].append(psrate/live)
592     Rates[name]["physics"].append(physics)
593     Rates[name]["active"].append(active)
594 abrinke1 1.5 Rates[name]["psi"].append(psi)
595 grchrist 1.9
596 grchrist 1.10 #for keys, values in MoreLumiMulti.iteritems():
597     # Rates[name][keys].append(values)
598 grchrist 1.9 #print nls, name, keys, values
599 awoodard 1.47
600 grchrist 1.10 RateOutput = open(RefRunFile, 'wb') ##Save new Rates[] to RefRuns
601     pickle.dump(Rates, RateOutput, 2)
602     RateOutput.close()
603     LumiOutput = open(RefLumiFile,'wb')
604     pickle.dump(LumiPageInfo,LumiOutput, 2)
605     LumiOutput.close()
606    
607 grchrist 1.59 #for trig in Rates.iterkeys():
608     # print trig, Rates[trig]
609     return [Rates,LumiPageInfo,trig_list]
610 abrinke1 1.1
611 awoodard 1.50 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,NoVersion, linear,do_inst, TMDerr,wp_bool,all_triggers):
612 abrinke1 1.1 min_run = min(run_list)
613     max_run = max(run_list)
614 muell149 1.46
615     priot.has_been_called=False
616 grchrist 1.11
617 abrinke1 1.5 InputFit = {}
618     OutputFit = {}
619 awoodard 1.72 failed_paths = []
620 grchrist 1.14 first_trigger=True
621 grchrist 1.37
622     [[varX, varY, do_fit, save_root, save_png, fit_file]] = plot_properties
623    
624     RootNameTemplate = "%s_%sLS_%s_vs_%s_Run%s-%s.root"
625     RootFile = RootNameTemplate % (trig_name, num_ls, varX, varY, min_run, max_run)
626 abrinke1 1.1
627 amott 1.20 if not do_fit:
628     try:
629 abrinke1 1.1 pkl_file = open(fit_file, 'rb')
630 abrinke1 1.5 InputFit = pickle.load(pkl_file)
631 grchrist 1.40
632 amott 1.20 except:
633     print "ERROR: could not open fit file: %s" % (fit_file,)
634     if save_root:
635     try:
636     os.remove(RootFile)
637     except:
638     pass
639    
640 grchrist 1.38 trig_list_noV=[]
641     for trigs in trig_list:
642     trig_list_noV.append(StripVersion(trigs))
643     if NoVersion:
644     trig_list=trig_list_noV
645 grchrist 1.40
646 amott 1.20 ## check that all the triggers we ask to plot are in the input fit
647     if not save_fits:
648     goodtrig_list = []
649 grchrist 1.40 FitInputNoV={}
650 amott 1.20 for trig in trig_list:
651 grchrist 1.40
652     if NoVersion:
653     for trigger in InputFit.iterkeys():
654     FitInputNoV[StripVersion(trigger)]=InputFit[trigger]
655     InputFit=FitInputNoV
656 awoodard 1.56
657 amott 1.20 else:
658 grchrist 1.40 if not InputFit.has_key(trig):
659     print "WARNING: No Fit Prediction for Trigger %s, SKIPPING" % (trig,)
660     else:
661     goodtrig_list.append(trig)
662     trig_list = goodtrig_list
663 amott 1.18
664 grchrist 1.26 for print_trigger in sorted(Rates):
665 abrinke1 1.7 ##Limits Rates[] to runs in run_list
666     NewTrigger = {}
667 grchrist 1.34
668 grchrist 1.8 if not print_trigger in trig_list:
669 awoodard 1.50 if all_triggers:
670     trig_list.append(print_trigger)
671     else:
672     print "not in trig_list:",print_trigger, trig_list
673     continue
674 grchrist 1.34
675 abrinke1 1.7 for key in Rates[print_trigger]:
676     NewTrigger[key] = []
677     for iterator in range (len(Rates[print_trigger]["run"])):
678     if Rates[print_trigger]["run"][iterator] in run_list:
679     for key in Rates[print_trigger]:
680     NewTrigger[key].append(Rates[print_trigger][key][iterator])
681     Rates[print_trigger] = NewTrigger
682    
683 grchrist 1.59 #print print_trigger, Rates[print_trigger]
684 abrinke1 1.1 meanrawrate = sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])
685 grchrist 1.40
686 abrinke1 1.1 if not trig_name in print_trigger:
687 grchrist 1.59 #print "failed does not contain HLT",trig_name, print_trigger
688     pass
689 grchrist 1.40
690 abrinke1 1.1 if meanrawrate < min_rate:
691     continue
692     masked_trig = False
693     for mask in masked_triggers:
694     if str(mask) in print_trigger:
695     masked_trig = True
696     if masked_trig:
697     continue
698 grchrist 1.34
699 abrinke1 1.5 OutputFit[print_trigger] = {}
700 abrinke1 1.1
701     lowlumi = 0
702 abrinke1 1.7 meanlumi_init = median(Rates[print_trigger]["live_lumi"])
703 abrinke1 1.1 meanlumi = 0
704     highlumi = 0
705 abrinke1 1.67 lowrate = 0
706     meanrate = 0
707     highrate = 0
708 abrinke1 1.1 lowxsec = 0
709     meanxsec = 0
710     highxsec = 0
711     nlow = 0
712     nhigh = 0
713 grchrist 1.15
714 abrinke1 1.1 for iterator in range(len(Rates[print_trigger]["rate"])):
715     if Rates[print_trigger]["live_lumi"][iterator] <= meanlumi_init:
716 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):
717 abrinke1 1.67 meanrate+=Rates[print_trigger]["rate"][iterator]
718     lowrate+=Rates[print_trigger]["rate"][iterator]
719 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
720     lowxsec+=Rates[print_trigger]["xsec"][iterator]
721     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
722     lowlumi+=Rates[print_trigger]["live_lumi"][iterator]
723     nlow+=1
724     if Rates[print_trigger]["live_lumi"][iterator] > meanlumi_init:
725 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):
726 abrinke1 1.67 meanrate+=Rates[print_trigger]["rate"][iterator]
727     highrate+=Rates[print_trigger]["rate"][iterator]
728 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
729     highxsec+=Rates[print_trigger]["xsec"][iterator]
730     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
731     highlumi+=Rates[print_trigger]["live_lumi"][iterator]
732     nhigh+=1
733 abrinke1 1.7 try:
734 abrinke1 1.67 meanrate = meanrate/(nlow+nhigh)
735 abrinke1 1.7 meanxsec = meanxsec/(nlow+nhigh)
736     meanlumi = meanlumi/(nlow+nhigh)
737 abrinke1 1.67 if (nlow==0):
738     sloperate = (highrate/nhigh) / (highlumi/nhigh)
739     slopexsec = (highxsec/nhigh) / (highlumi/nhigh)
740     elif (nhigh==0):
741     sloperate = (lowrate/nlow) / (lowlumi/nlow)
742     slopexsec = (lowxsec/nlow) / (lowlumi/nlow)
743     else:
744     sloperate = ( (highrate/nhigh) - (lowrate/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
745     slopexsec = ( (highxsec/nhigh) - (lowxsec/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
746 abrinke1 1.7 except:
747 awoodard 1.72 # print str(print_trigger)+" has no good datapoints - setting initial xsec slope estimate to 0"
748 abrinke1 1.67 meanrate = median(Rates[print_trigger]["rate"])
749 abrinke1 1.7 meanxsec = median(Rates[print_trigger]["xsec"])
750     meanlumi = median(Rates[print_trigger]["live_lumi"])
751 abrinke1 1.67 sloperate = meanxsec
752 abrinke1 1.7 slopexsec = 0
753 abrinke1 1.1
754 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()
755    
756 amott 1.20 if not do_fit:
757 awoodard 1.62 try:
758     FitType = InputFit[print_trigger][0]
759 awoodard 1.72 except:
760     failed_paths.append([print_trigger,"This path did not exist in the monitorlist used to create the fit"])
761     continue
762     if FitType == "fit failed":
763     failure_comment = InputFit[print_trigger][1]
764     failed_paths.append([print_trigger, failure_comment])
765     continue
766     else:
767 awoodard 1.62 X0 = InputFit[print_trigger][1]
768     X1 = InputFit[print_trigger][2]
769     X2 = InputFit[print_trigger][3]
770     X3 = InputFit[print_trigger][4]
771     sigma = InputFit[print_trigger][5]/math.sqrt(num_ls)*3#Display 3 sigma band to show outliers more clearly
772     X0err= InputFit[print_trigger][7]
773    
774 grchrist 1.8 ## we are 2 lumis off when we start! -gets worse when we skip lumis
775 grchrist 1.17 it_offset=0
776 abrinke1 1.1 for iterator in range(len(Rates[print_trigger]["rate"])):
777     if not Rates[print_trigger]["run"][iterator] in run_list:
778     continue
779 abrinke1 1.67 ##Old Spike-killer: used average-xsec method, too clumsy, esp. for nonlinear triggers
780     #prediction = meanxsec + slopexsec * (Rates[print_trigger]["live_lumi"][iterator] - meanlumi)
781     #realvalue = Rates[print_trigger]["xsec"][iterator]
782    
783     ##New Spike-killer: gets average of nearest 4 LS
784 grchrist 1.68 try:
785     if (iterator > 2 and iterator+2 < len(Rates[print_trigger]["rate"])): ##2 LS before, 2 after
786     prediction = (Rates[print_trigger]["rate"][iterator-2]+Rates[print_trigger]["rate"][iterator-1]+Rates[print_trigger]["rate"][iterator+1]+Rates[print_trigger]["rate"][iterator+2])/4.0
787     elif (iterator > 2 and len(Rates[print_trigger]["rate"]) > 4 ): ##4 LS before
788     prediction = (Rates[print_trigger]["rate"][iterator-4]+Rates[print_trigger]["rate"][iterator-3]+Rates[print_trigger]["rate"][iterator-2]+Rates[print_trigger]["rate"][iterator-1])/4.0
789     elif (iterator+2 < len(Rates[print_trigger]["rate"]) and len(Rates[print_trigger]["rate"]) > 4 ): ##4 LS after
790     prediction = (Rates[print_trigger]["rate"][iterator+1]+Rates[print_trigger]["rate"][iterator+2]+Rates[print_trigger]["rate"][iterator+3]+Rates[print_trigger]["rate"][iterator+4])/4.0
791     else:
792     prediction = Rates[print_trigger]["rate"][iterator]
793     realvalue = Rates[print_trigger]["rate"][iterator]
794     except:
795     print "Error calculating prediction. Setting rates to defaults."
796 abrinke1 1.67 prediction = Rates[print_trigger]["rate"][iterator]
797 grchrist 1.68 #realvalue = Rates[print_trigger]["rate"][iterator]
798 grchrist 1.11
799     if pass_cuts(data_clean, realvalue, prediction, meanxsec, Rates, print_trigger, iterator, num_ls,LumiPageInfo,SubSystemOff,max_dt,print_info, trig_list):
800 grchrist 1.40
801 abrinke1 1.1 run_t.append(Rates[print_trigger]["run"][iterator])
802     ls_t.append(Rates[print_trigger]["ls"][iterator])
803     ps_t.append(Rates[print_trigger]["ps"][iterator])
804     inst_t.append(Rates[print_trigger]["inst_lumi"][iterator])
805 grchrist 1.17 live_t.append(Rates[print_trigger]["live_lumi"][iterator])
806     delivered_t.append(Rates[print_trigger]["delivered_lumi"][iterator])
807     deadtime_t.append(Rates[print_trigger]["deadtime"][iterator])
808 abrinke1 1.1 rawrate_t.append(Rates[print_trigger]["rawrate"][iterator])
809     rate_t.append(Rates[print_trigger]["rate"][iterator])
810     rawxsec_t.append(Rates[print_trigger]["rawxsec"][iterator])
811     xsec_t.append(Rates[print_trigger]["xsec"][iterator])
812 abrinke1 1.5 psi_t.append(Rates[print_trigger]["psi"][iterator])
813 abrinke1 1.1
814     e_run_t.append(0.0)
815     e_ls_t.append(0.0)
816     e_ps_t.append(0.0)
817     e_inst_t.append(14.14)
818     e_live_t.append(14.14)
819     e_delivered_t.append(14.14)
820 abrinke1 1.2 e_deadtime_t.append(0.01)
821 abrinke1 1.1 e_rawrate_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
822     e_rate_t.append(Rates[print_trigger]["ps"][iterator]*math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
823 abrinke1 1.5 e_psi_t.append(0.0)
824 abrinke1 1.2 if live_t[-1] == 0:
825     e_rawxsec_t.append(0)
826     e_xsec_t.append(0)
827     else:
828 grchrist 1.16 try:
829     e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["live_lumi"][iterator])
830     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])
831     except:
832     e_rawxsec_t.append(0.)
833     e_xsec_t.append(0.)
834 amott 1.20 if not do_fit:
835 grchrist 1.38 if not do_inst:
836     if FitType == "expo":
837 grchrist 1.43 rate_prediction = X0 + X1*math.exp(X2+X3*delivered_t[-1])
838 grchrist 1.38 else:
839     rate_prediction = X0 + X1*delivered_t[-1] + X2*delivered_t[-1]*delivered_t[-1] + X3*delivered_t[-1]*delivered_t[-1]*delivered_t[-1]
840 abrinke1 1.67
841 grchrist 1.38 else:
842     if FitType == "expo":
843 grchrist 1.43 rate_prediction = X0 + X1*math.exp(X2+X3*inst_t[-1])
844 grchrist 1.38 else:
845     rate_prediction = X0 + X1*inst_t[-1] + X2*inst_t[-1]*inst_t[-1] + X3*inst_t[-1]*inst_t[-1]*inst_t[-1]
846 abrinke1 1.5
847 abrinke1 1.2 if live_t[-1] == 0:
848 abrinke1 1.5 rawrate_fit_t.append(0)
849     rate_fit_t.append(0)
850 abrinke1 1.2 rawxsec_fit_t.append(0)
851     xsec_fit_t.append(0)
852 abrinke1 1.7 e_rawrate_fit_t.append(0)
853 awoodard 1.45 e_rate_fit_t.append(sigma)
854 abrinke1 1.7 e_rawxsec_fit_t.append(0)
855     e_xsec_fit_t.append(0)
856 abrinke1 1.67
857 abrinke1 1.2 else:
858 grchrist 1.11 if ps_t[-1]>0.0:
859     rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
860     else:
861     rawrate_fit_t.append(0.0)
862    
863 awoodard 1.52 rate_fit_t.append(rate_prediction)
864 awoodard 1.45 e_rate_fit_t.append(sigma)
865 abrinke1 1.5 rawxsec_fit_t.append(rawrate_fit_t[-1]/live_t[-1])
866     xsec_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/live_t[-1])
867 grchrist 1.38 try:
868    
869     if not TMDerr:
870 awoodard 1.45 e_rawrate_fit_t.append(sigma*rawrate_fit_t[-1]/rate_fit_t[-1])
871     e_rawxsec_fit_t.append(sigma*rawxsec_fit_t[-1]/rate_fit_t[-1])
872     e_xsec_fit_t.append(sigma*xsec_fit_t[-1]/rate_fit_t[-1])
873 grchrist 1.38 ###error from TMD predictions, calculated at 5e33
874     else:
875     e_rawrate_fit_t.append(X0err*inst_t[-1]/5000.)
876     e_rawxsec_fit_t.append(X0err/live_t[-1]*inst_t[-1]/5000.)
877     e_xsec_fit_t.append(X0err/live_t[-1]*inst_t[-1]/5000.)
878    
879     except:
880     print print_trigger, "has no fitted rate for LS", Rates[print_trigger]["ls"][iterator]
881 awoodard 1.45 e_rawrate_fit_t.append(sigma)
882     e_rawxsec_fit_t.append(sigma)
883     e_xsec_fit_t.append(sigma)
884 abrinke1 1.5
885 grchrist 1.16
886 grchrist 1.27 if (print_info and num_ls==1 and (fabs(rawrate_fit_t[-1]-rawrate_t[-1])>2.5*sqrt(sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])))):
887 grchrist 1.38 pass
888 grchrist 1.8
889 abrinke1 1.5 else: ##If the data point does not pass the data_clean filter
890 abrinke1 1.1 if debug_print:
891     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)
892    
893 abrinke1 1.5 ## End "for iterator in range(len(Rates[print_trigger]["rate"])):" loop
894 abrinke1 1.1
895 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]
896 grchrist 1.40
897 awoodard 1.63 [VX, VXE, x_label, VY, VYE, y_label, VF, VFE] = GetVXVY(plot_properties, fit_file, AllPlotArrays)
898 grchrist 1.40
899 abrinke1 1.5 if save_root or save_png:
900     c1 = TCanvas(str(varX),str(varY))
901     c1.SetName(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
902 grchrist 1.40
903     try:
904     gr1 = TGraphErrors(len(VX), VX, VY, VXE, VYE)
905     except:
906 awoodard 1.72 failure_comment = "In runs specified during creation of the fit file, there were no events for this path: probably due to high deadtime or low raw (prescaled) rate"
907     failed_paths.append([print_trigger,failure_comment])
908     if do_fit:
909     OutputFit[print_trigger] = ["fit failed",failure_comment]
910 abrinke1 1.67 continue
911 awoodard 1.72 if (len(VX)<20 and do_fit):
912     failure_comment = "In runs specified during creation of the fit file, there were less than 20 datapoints for this path: probably due to high deadtime or low raw (prescaled) rate"
913     failed_paths.append([print_trigger,failure_comment])
914     OutputFit[print_trigger] = ["fit failed",failure_comment]
915 grchrist 1.40 continue
916 awoodard 1.72
917 abrinke1 1.5 gr1.SetName("Graph_"+str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
918 awoodard 1.63 gr1.GetXaxis().SetTitle(x_label)
919     gr1.GetYaxis().SetTitle(y_label)
920 abrinke1 1.5 gr1.SetTitle(str(print_trigger))
921     gr1.SetMinimum(0)
922     gr1.SetMaximum(1.2*max(VY))
923     #gr1.GetXaxis().SetLimits(min(VX)-0.2*max(VX),1.2*max(VX))
924     gr1.GetXaxis().SetLimits(0,1.2*max(VX))
925     gr1.SetMarkerStyle(8)
926 grchrist 1.40
927 abrinke1 1.5 if fit_file:
928     gr1.SetMarkerSize(0.8)
929     else:
930     gr1.SetMarkerSize(0.5)
931     gr1.SetMarkerColor(2)
932    
933 amott 1.20 if not do_fit:
934 abrinke1 1.5 gr3 = TGraphErrors(len(VX), VX, VF, VXE, VFE)
935     gr3.SetMarkerStyle(8)
936     gr3.SetMarkerSize(0.4)
937     gr3.SetMarkerColor(4)
938     gr3.SetFillColor(4)
939     gr3.SetFillStyle(3003)
940 abrinke1 1.1
941 abrinke1 1.5 if do_fit:
942 grchrist 1.35 if "rate" in varY and not linear:
943 grchrist 1.61 f1d=0
944     f1d = TF1("f1d","pol1",0,8000)#linear
945 abrinke1 1.67 f1d.SetParameters(0.01,min(sum(VY)/sum(VX),sloperate)) ##Set Y-intercept near 0, slope either mean_rate/mean_lumi or est. slope (may be negative)
946 grchrist 1.61 f1d.SetLineColor(4)
947     f1d.SetLineWidth(2)
948 abrinke1 1.67 if nlow>0:
949     f1d.SetParLimits(0,0,1.5*lowrate/nlow) ##Keep Y-intercept in range of low-lumi rate points
950     else:
951     f1d.SetParLimits(0,0,1.5*sum(VY)/len(VY))
952     if (sloperate > 0):
953     if (sloperate > 0.5*sum(VY)/sum(VX)): ##Slope is substantially positive
954     f1d.SetParLimits(1,min(0.5*sloperate,0.5*sum(VY)/sum(VX)),1.5*sum(VY)/sum(VX))
955     else: ##Slope is somewhat positive or flat
956     f1d.SetParLimits(1,-0.1*sloperate,1.5*sum(VY)/sum(VX))
957     else: ##Slope is negative or flat
958     f1d.SetParLimits(1,1.5*sloperate,-0.1*sloperate)
959     gr1.Fit("f1d","QN","rob=0.90")
960     f1d_Chi2 = f1d.GetChisquare()/f1d.GetNDF()
961 grchrist 1.61
962 grchrist 1.41 f1a=0
963 abrinke1 1.67 f1a = TF1("f1a","pol2",0,8000)#quadratic
964     f1a.SetParameters(f1d.GetParameter(0),f1d.GetParameter(1),0) ##Initial values from linear fit
965 grchrist 1.61 f1a.SetLineColor(6)
966 abrinke1 1.5 f1a.SetLineWidth(2)
967 abrinke1 1.67 if nlow>0 and sloperate < 0.5*sum(VY)/sum(VX): ##Slope is not substantially positive
968     f1a.SetParLimits(0,0,1.5*lowrate/nlow) ##Keep Y-intercept in range of low-lumi rate points
969     else:
970     f1a.SetParLimits(0,0,max(min(VY),0.3*sum(VY)/len(VY))) ##Keep Y-intercept reasonably low
971     f1a.SetParLimits(1,-2.0*(max(VY)-min(VY))/(max(VX)-min(VX)),2.0*(max(VY)-min(VY))/(max(VX)-min(VX))) ##Reasonable bounds
972     f1a.SetParLimits(2,-2.0*max(VY)/(max(VX)*max(VX)),2.0*max(VY)/(max(VX)*max(VX))) ##Reasonable bounds
973 awoodard 1.48 gr1.Fit("f1a","QN","rob=0.90")
974 abrinke1 1.67 f1a_Chi2 = f1a.GetChisquare()/f1a.GetNDF()
975     f1a_BadMinimum = (f1a.GetMinimumX(5,7905,10)>2000 and f1a.GetMinimumX(5,7905,10)<7000) ##Don't allow minimum between 2000 and 7000
976 muell149 1.46
977 abrinke1 1.5 f1b = 0
978     f1c = 0
979 muell149 1.46 meanps = median(Rates[print_trigger]["ps"])
980     av_rte = mean(VY)
981    
982 grchrist 1.43 if True:
983 abrinke1 1.67 f1b = TF1("f1b","pol3",0,8000)#cubic
984     f1b.SetParameters(f1a.GetParameter(0),f1a.GetParameter(1),f1a.GetParameter(2),0) ##Initial values from quadratic fit
985 grchrist 1.43 f1b.SetLineColor(2)
986     f1b.SetLineWidth(2)
987 abrinke1 1.67 f1b.SetParLimits(0,0,max(min(VY),0.3*sum(VY)/len(VY))) ##Keep Y-intercept reasonably low
988     f1b.SetParLimits(1,-2.0*(max(VY)-min(VY))/(max(VX)-min(VX)),2.0*(max(VY)-min(VY))/(max(VX)-min(VX))) ##Reasonable bounds
989     f1b.SetParLimits(2,-2.0*max(VY)/(max(VX)*max(VX)),2.0*max(VY)/(max(VX)*max(VX))) ##Reasonable bounds
990     f1b.SetParLimits(3,0,2.0*max(VY)/(max(VX)*max(VX)*max(VX))) ##Reasonable bounds
991 awoodard 1.48 gr1.Fit("f1b","QN","rob=0.90")
992 grchrist 1.70 try:
993     f1b_Chi2 = f1b.GetChisquare()/f1b.GetNDF()
994     except ZeroDivisionError:
995     f1b_Chi2=10.0
996     print "Zero ndf for ",print_trigger
997 awoodard 1.72
998     f1b_BadMinimum = (f1b.GetMinimumX(5,7905,10)>2000 and f1b.GetMinimumX(5,7905,10)<7000)
999 grchrist 1.43 f1c = TF1("f1c","[0]+[1]*expo(2)",0,8000)
1000     f1c.SetLineColor(3)
1001     f1c.SetLineWidth(2)
1002 abrinke1 1.67 #f1c.SetParLimits(0,0,max(min(VY),0.3*sum(VY)/len(VY)))
1003     f1c.SetParLimits(0,0,max(min(VY),0.01*sum(VY)/len(VY))) ##Exponential fits should start low
1004     f1c.SetParLimits(1,max(VY)/math.exp(15.0),max(VY)/math.exp(2.0))
1005 grchrist 1.43 f1c.SetParLimits(2,0.0,0.0000000001)
1006 abrinke1 1.67 f1c.SetParLimits(3,2.0/max(VX),15.0/max(VX))
1007 awoodard 1.48 gr1.Fit("f1c","QN","rob=0.90")
1008 abrinke1 1.67 f1c_Chi2 = f1c.GetChisquare()/f1c.GetNDF()
1009     f1c_BadMinimum = ((f1c.GetMinimumX(5,7905,10)>2000 and f1c.GetMinimumX(5,7905,10)<7000)) or f1c.GetMaximum(min(VX),max(VX),10)/max(VY) > 2.0
1010     ##Some fits are so exponential, the graph ends early and returns a false low Chi2 value
1011 grchrist 1.15
1012 abrinke1 1.5 else: ##If this is not a rate plot
1013     f1a = TF1("f1a","pol1",0,8000)
1014     f1a.SetLineColor(4)
1015     f1a.SetLineWidth(2)
1016     if "xsec" in varY:
1017     f1a.SetParLimits(0,0,meanxsec*1.5)
1018     if slopexsec > 0:
1019     f1a.SetParLimits(1,0,max(VY)/max(VX))
1020     else:
1021     f1a.SetParLimits(1,2*slopexsec,-2*slopexsec)
1022 abrinke1 1.1 else:
1023 abrinke1 1.5 f1a.SetParLimits(0,-1000,1000)
1024     gr1.Fit("f1a","Q","rob=0.80")
1025 abrinke1 1.1
1026 grchrist 1.35 if (first_trigger):
1027 grchrist 1.59 print '%-50s %4s x0 x1 x2 x3 chi2 ndf chi2/ndf' % ('trigger', 'type')
1028 grchrist 1.35 first_trigger=False
1029 grchrist 1.40 try:
1030 abrinke1 1.67 print '%-50s | line | % .2f | +/-%.2f | % .2e | +/-%.1e | % .2e | +/-%.1e | % .2e | +/-%.1e | %7.0f | %4.0f | %5.2f | ' % (print_trigger, f1a.GetParameter(0), f1a.GetParError(0), f1a.GetParameter(1), f1a.GetParError(1), 0 , 0 , 0 , 0 , f1a.GetChisquare(), f1a.GetNDF(), f1a_Chi2)
1031 grchrist 1.40 except:
1032     pass
1033 awoodard 1.72
1034     chioffset=1.0 ##chioffset now a fraction; must be 10% better to use expo rather than quad, quad rather than line
1035     width = max([len(trigger_name) for trigger_name in trig_list])
1036 awoodard 1.48 if print_table or save_fits:
1037     if not do_fit:
1038     print "Can't have save_fits = True and do_fit = False"
1039     continue
1040 awoodard 1.72 if min([f1a_Chi2,f1b_Chi2,f1c_Chi2,f1d_Chi2]) > 500:#require a minimum chi^2/nDOF of 500
1041     failure_comment = "There were events for this path in the runs specified during the creation of the fit file, but the fit failed to converge"
1042     failed_paths.append([print_trigger,failure_comment])
1043     OutputFit[print_trigger] = ["fit failed",failure_comment]
1044     continue
1045     if "rate" in varY and not linear:
1046     if first_trigger:
1047     print '\n%-*s | TYPE | %-8s | %-11s | %-7s | %-10s | %-7s | %-10s | %-8s | %-10s | %-6s | %-4s |%-7s|' % (width,"TRIGGER", "X0","X0 ERROR","X1","X1 ERROR","X2","X2 ERROR","X3","X3 ERROR","CHI^2","DOF","CHI2/DOF")
1048     first_trigger = False
1049    
1050 grchrist 1.70 if ((f1c_Chi2 < (f1a_Chi2*chioffset) or f1a_BadMinimum) and ((f1c_Chi2 < f1b_Chi2) or f1b_BadMinimum) and f1c_Chi2 < (f1d_Chi2*chioffset) and not f1c_BadMinimum and len(VX)>1):
1051 awoodard 1.72 print '%-*s | expo | %-8.1f | +/-%-8.1f | %8.1e | +/-%.1e | %8.1e | +/-%.1e | %-8.1e | +/-%.1e | %6.0f | %4.0f | %5.1f | ' % (width,print_trigger, f1c.GetParameter(0) , f1c.GetParError(0) , f1c.GetParameter(1) , f1c.GetParError(1) , f1c.GetParameter(2), f1c.GetParError(2) ,f1c.GetParameter(3), f1c.GetParError(3) ,f1c.GetChisquare() , f1c.GetNDF() , f1c_Chi2)
1052 awoodard 1.48 f1c.SetLineColor(1)
1053     priot(wp_bool,print_trigger,meanps,f1d,f1c,"expo",av_rte)
1054     sigma = CalcSigma(VX, VY, f1c)*math.sqrt(num_ls)
1055     OutputFit[print_trigger] = ["expo", f1c.GetParameter(0) , f1c.GetParameter(1) , f1c.GetParameter(2) , f1c.GetParameter(3) , sigma , meanrawrate, f1c.GetParError(0) , f1c.GetParError(1) , f1c.GetParError(2) , f1c.GetParError(3)]
1056    
1057 grchrist 1.70 elif ((f1b_Chi2 < (f1a_Chi2*chioffset) or f1a_BadMinimum) and f1b_Chi2 < (f1d_Chi2*chioffset) and not f1b_BadMinimum and len(VX)>1):
1058 awoodard 1.72 print '%-*s | cube | %-8.1f | +/-%-8.1f | %8.1e | +/-%.1e | %8.1e | +/-%.1e | %-8.1e | +/-%.1e | %6.0f | %4.0f | %5.1f | ' % (width,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_Chi2)
1059 awoodard 1.48 f1b.SetLineColor(1)
1060     priot(wp_bool,print_trigger,meanps,f1d,f1b,"cubic",av_rte)
1061     sigma = CalcSigma(VX, VY, f1b)*math.sqrt(num_ls)
1062     OutputFit[print_trigger] = ["poly", f1b.GetParameter(0) , f1b.GetParameter(1) , f1b.GetParameter(2) , f1b.GetParameter(3) , sigma , meanrawrate, f1b.GetParError(0) , f1b.GetParError(1) , f1b.GetParError(2) , f1b.GetParError(3)]
1063    
1064 abrinke1 1.67 elif (f1a_Chi2 < (f1d_Chi2*chioffset)):
1065 awoodard 1.72 print '%-*s | quad | %-8.1f | +/-%-8.1f | %8.1e | +/-%.1e | %8.1e | +/-%.1e | %-8.1e | +/-%.1e | %6.0f | %4.0f | %5.1f | ' % (width,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_Chi2)
1066 awoodard 1.48 f1a.SetLineColor(1)
1067     priot(wp_bool,print_trigger,meanps,f1d,f1a,"quad",av_rte)
1068 awoodard 1.53 sigma = CalcSigma(VX, VY, f1a)*math.sqrt(num_ls)
1069 awoodard 1.48 OutputFit[print_trigger] = ["poly", f1a.GetParameter(0) , f1a.GetParameter(1) , f1a.GetParameter(2) , 0.0 , sigma , meanrawrate, f1a.GetParError(0) , f1a.GetParError(1) , f1a.GetParError(2) , 0.0]
1070 awoodard 1.72
1071 abrinke1 1.67 else:
1072 awoodard 1.72 print '%-*s | line | %-8.1f | +/-%-8.1f | %8.1e | +/-%.1e | %8.1e | +/-%.1e | %-8.1e | +/-%.1e | %6.0f | %4.0f | %5.1f | ' % (width,print_trigger, f1d.GetParameter(0) , f1d.GetParError(0) , f1d.GetParameter(1) , f1d.GetParError(1) , 0 , 0 , 0 , 0 , f1d.GetChisquare() , f1d.GetNDF() , f1d_Chi2)
1073 abrinke1 1.67 f1d.SetLineColor(1)
1074     priot(wp_bool,print_trigger,meanps,f1d,f1d,"line",av_rte)
1075     sigma = CalcSigma(VX, VY, f1d)*math.sqrt(num_ls)
1076     OutputFit[print_trigger] = ["poly", f1d.GetParameter(0) , f1d.GetParameter(1) , 0.0 , 0.0 , sigma , meanrawrate, f1d.GetParError(0) , f1d.GetParError(1) , 0.0 , 0.0]
1077 grchrist 1.65 else:
1078 awoodard 1.72 print '%-*s | quad | %-8.1f | +/-%-8.1f | %8.1e | +/-%.1e | %8.1e | +/-%.1e | %-8.1e | +/-%.1e | %6.0f | %4.0f | %5.1f | ' % (width,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_Chi2)
1079 grchrist 1.65 f1a.SetLineColor(1)
1080     #priot(wp_bool,print_trigger,meanps,f1d,f1a,"quad",av_rte)
1081     sigma = CalcSigma(VX, VY, f1a)*math.sqrt(num_ls)
1082     OutputFit[print_trigger] = ["poly", f1a.GetParameter(0) , f1a.GetParameter(1) , f1a.GetParameter(2) , 0.0 , sigma , meanrawrate, f1a.GetParError(0) , f1a.GetParError(1) , f1a.GetParError(2) , 0.0]
1083 awoodard 1.48
1084 abrinke1 1.5 if save_root or save_png:
1085     gr1.Draw("APZ")
1086 amott 1.20 if not do_fit:
1087 abrinke1 1.5 gr3.Draw("P3")
1088 grchrist 1.58
1089 abrinke1 1.5 if do_fit:
1090 grchrist 1.59 try:
1091 abrinke1 1.67 if ((f1c_Chi2 < (f1a_Chi2*chioffset) or f1a_BadMinimum ) and (f1c_Chi2 < f1b_Chi2 or f1b_BadMinimum ) and not f1c_BadMinimum ):
1092 grchrist 1.59 f1c.Draw("same")
1093 abrinke1 1.67 elif ( (f1b_Chi2 < (f1a_Chi2*chioffset) or f1a_BadMinimum) and not f1b_BadMinimum):
1094 grchrist 1.59 f1b.Draw("same")
1095     else:
1096     f1a.Draw("same")
1097    
1098     f1d.Draw("same")
1099     except:
1100     True
1101 grchrist 1.58
1102 abrinke1 1.5 c1.Update()
1103     if save_root:
1104     myfile = TFile( RootFile, 'UPDATE' )
1105     c1.Write()
1106     myfile.Close()
1107     if save_png:
1108     c1.SaveAs(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX)+".png")
1109 awoodard 1.72
1110     if len(failed_paths) > 0:
1111     if save_fits:
1112     print "\n***************NO FIT RECORDED FOR THE FOLLOWING PATHS***************"
1113     else:
1114     print "\n***************THE FOLLOWING PATHS HAVE BEEN SKIPPED BECAUSE THE FIT WAS MISSING***************"
1115     sorted_failed_paths = sorted(failed_paths, key=itemgetter(1))
1116     for error_comment, entries in groupby(sorted_failed_paths, key=itemgetter(1)):
1117     error_comment = error_comment.replace('this path','these paths')
1118     print '\n'+error_comment+'.'
1119     for entry in entries:
1120     print entry[0]
1121    
1122 abrinke1 1.1 if save_root:
1123 awoodard 1.72 print "\nOutput root file is "+str(RootFile)
1124 abrinke1 1.1
1125     if save_fits:
1126 amott 1.20 if os.path.exists(fit_file):
1127     os.remove(fit_file)
1128     FitOutputFile = open(fit_file, 'wb')
1129 abrinke1 1.5 pickle.dump(OutputFit, FitOutputFile, 2)
1130     FitOutputFile.close()
1131 amott 1.20 print "Output fit file is "+str(fit_file)
1132 abrinke1 1.1
1133     if print_table:
1134 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."
1135 grchrist 1.59 print '%50s%10s%10s%10s%10s%10s%10s%10s' % ("Trig", "fit", "p0", "p1", "p2", "p3", "Chi2", "Av raw")
1136 abrinke1 1.5 for print_trigger in OutputFit:
1137     _trigger = (print_trigger[:56] + '...') if len(print_trigger) > 59 else print_trigger
1138 abrinke1 1.1 try:
1139 abrinke1 1.5 if OutputFit[print_trigger][0] == "poly":
1140 grchrist 1.59 print '%50s%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))
1141 abrinke1 1.5 else:
1142 grchrist 1.59 print '%50s%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))
1143 abrinke1 1.1 except:
1144     print str(print_trigger)+" is somehow broken"
1145 amott 1.18 return RootFile
1146 abrinke1 1.5
1147     ############# SUPPORTING FUNCTIONS ################
1148    
1149    
1150 awoodard 1.47 def CalcSigma(var_x, var_y, func):
1151     residuals = []
1152     for x, y in zip(var_x,var_y):
1153 awoodard 1.53 residuals.append(y - func.Eval(x,0,0))
1154 awoodard 1.47
1155     res_squared = [i*i for i in residuals]
1156 awoodard 1.53 if len(res_squared) > 2:
1157     sigma = math.sqrt(sum(res_squared)/(len(res_squared)-2))
1158     else:
1159     sigma = 0
1160 awoodard 1.47 return sigma
1161    
1162 abrinke1 1.5 def GetJSON(json_file):
1163    
1164     input_file = open(json_file)
1165     file_content = input_file.read()
1166     inputRange = selectionParser(file_content)
1167     JSON = inputRange.runsandls()
1168     return JSON
1169     ##JSON is an array: JSON[run_number] = [1st ls, 2nd ls, 3rd ls ... nth ls]
1170    
1171     def MakePlotArrays():
1172     run_t = array.array('f')
1173     ls_t = array.array('f')
1174     ps_t = array.array('f')
1175     inst_t = array.array('f')
1176     live_t = array.array('f')
1177     delivered_t = array.array('f')
1178     deadtime_t = array.array('f')
1179     rawrate_t = array.array('f')
1180     rate_t = array.array('f')
1181     rawxsec_t = array.array('f')
1182     xsec_t = array.array('f')
1183     psi_t = array.array('f')
1184    
1185     e_run_t = array.array('f')
1186     e_ls_t = array.array('f')
1187     e_ps_t = array.array('f')
1188     e_inst_t = array.array('f')
1189     e_live_t = array.array('f')
1190     e_delivered_t = array.array('f')
1191     e_deadtime_t = array.array('f')
1192     e_rawrate_t = array.array('f')
1193     e_rate_t = array.array('f')
1194     e_rawxsec_t = array.array('f')
1195     e_xsec_t = array.array('f')
1196     e_psi_t = array.array('f')
1197    
1198     rawrate_fit_t = array.array('f')
1199     rate_fit_t = array.array('f')
1200     rawxsec_fit_t = array.array('f')
1201     xsec_fit_t = array.array('f')
1202     e_rawrate_fit_t = array.array('f')
1203     e_rate_fit_t = array.array('f')
1204     e_rawxsec_fit_t = array.array('f')
1205     e_xsec_fit_t = array.array('f')
1206    
1207     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]
1208    
1209    
1210     def GetVXVY(plot_properties, fit_file, AllPlotArrays):
1211    
1212     VF = "0"
1213     VFE = "0"
1214    
1215     [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
1216     for varX, varY, do_fit, save_root, save_png, fit_file in plot_properties:
1217     if varX == "run":
1218     VX = run_t
1219     VXE = run_t_e
1220 awoodard 1.63 x_label = "Run Number"
1221 abrinke1 1.5 elif varX == "ls":
1222     VX = ls_t
1223     VXE = e_ls_t
1224 awoodard 1.63 x_label = "Lumisection"
1225 abrinke1 1.5 elif varX == "ps":
1226     VX = ps_t
1227     VXE = e_ps_t
1228 awoodard 1.63 x_label = "Prescale"
1229 abrinke1 1.5 elif varX == "inst":
1230     VX = inst_t
1231     VXE = e_inst_t
1232 abrinke1 1.67 x_label = "Instantaneous Luminosity [10^{30} Hz/cm^{2}]"
1233    
1234 abrinke1 1.5 elif varX == "live":
1235     VX = live_t
1236     VXE = e_live_t
1237 abrinke1 1.67 x_label = "Instantaneous Luminosity [10^{30} Hz/cm^{2}]"
1238    
1239 abrinke1 1.5 elif varX == "delivered":
1240     VX = delivered_t
1241     VXE = e_delivered_t
1242 abrinke1 1.67 x_label = "Instantaneous Luminosity [10^{30} Hz/cm^{2}]"
1243    
1244 abrinke1 1.5 elif varX == "deadtime":
1245     VX = deadtime_t
1246     VXE = e_deadtime_t
1247 awoodard 1.63 x_label = "Deadtime"
1248 abrinke1 1.5 elif varX == "rawrate":
1249     VX = rawrate_t
1250     VXE = e_rawrate_t
1251 awoodard 1.63 x_label = "Raw Rate [Hz]"
1252 abrinke1 1.5 elif varX == "rate":
1253     VX = rate_t
1254     VXE = e_rate_t
1255 awoodard 1.63 x_label = "Rate [Hz]"
1256 abrinke1 1.5 elif varX == "rawxsec":
1257     VX = rawxsec_t
1258     VXE = e_rawxsec_t
1259 grchrist 1.65 x_label = "Cross Section"
1260 abrinke1 1.5 elif varX == "xsec":
1261     VX = xsec_t
1262     VXE = e_xsec_t
1263 grchrist 1.65 x_label = "Cross Section"
1264 abrinke1 1.5 elif varX == "psi":
1265     VX = psi_t
1266     VXE = e_psi_t
1267 awoodard 1.63 x_label = "Prescale Index"
1268 abrinke1 1.5 else:
1269     print "No valid variable entered for X"
1270     continue
1271     if varY == "run":
1272     VY = run_t
1273     VYE = run_t_e
1274 awoodard 1.63 y_label = "Run Number"
1275 abrinke1 1.5 elif varY == "ls":
1276     VY = ls_t
1277     VYE = e_ls_t
1278 awoodard 1.63 y_label = "Lumisection"
1279 abrinke1 1.5 elif varY == "ps":
1280     VY = ps_t
1281     VYE = e_ps_t
1282 awoodard 1.63 y_label = "Prescale"
1283 abrinke1 1.5 elif varY == "inst":
1284     VY = inst_t
1285     VYE = e_inst_t
1286 grchrist 1.65 y_label = "Instantaneous Luminosity"
1287 abrinke1 1.5 elif varY == "live":
1288     VY = live_t
1289     VYE = e_live_t
1290 grchrist 1.65 y_label = "Instantaneous Luminosity"
1291 abrinke1 1.5 elif varY == "delivered":
1292     VY = delivered_t
1293     VYE = e_delivered_t
1294 grchrist 1.65 y_label = "Instantaneous Luminosity"
1295 abrinke1 1.5 elif varY == "deadtime":
1296     VY = deadtime_t
1297     VYE = e_deadtime_t
1298 awoodard 1.63 y_label = "Deadtime"
1299 abrinke1 1.5 elif varY == "rawrate":
1300     VY = rawrate_t
1301     VYE = e_rawrate_t
1302 awoodard 1.63 y_label = "Raw Rate [Hz]"
1303 abrinke1 1.5 if fit_file:
1304     VF = rawrate_fit_t
1305     VFE = e_rawrate_fit_t
1306     elif varY == "rate":
1307     VY = rate_t
1308     VYE = e_rate_t
1309 awoodard 1.63 y_label = "Rate [Hz]"
1310 abrinke1 1.5 if fit_file:
1311     VF = rate_fit_t
1312     VFE = e_rate_fit_t
1313     elif varY == "rawxsec":
1314     VY = rawxsec_t
1315     VYE = e_rawxsec_t
1316 grchrist 1.65 y_label = "Cross Section"
1317 abrinke1 1.5 if fit_file:
1318     VF = rawxsec_fit_t
1319     VFE = e_rawxsec_fit_t
1320     elif varY == "xsec":
1321     VY = xsec_t
1322     VYE = e_xsec_t
1323 grchrist 1.65 y_label = "Cross Section"
1324 abrinke1 1.5 if fit_file:
1325     VF = xsec_fit_t
1326     VFE = e_xsec_fit_t
1327     elif varY == "psi":
1328     VY = psi_t
1329     VYE = e_psi_t
1330 awoodard 1.63 y_label = "Prescale Index"
1331 abrinke1 1.5 else:
1332     print "No valid variable entered for Y"
1333     continue
1334    
1335 awoodard 1.63 return [VX, VXE, x_label, VY, VYE, y_label, VF, VFE]
1336 abrinke1 1.5
1337 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):
1338 grchrist 1.17 it_offset=0
1339 grchrist 1.15 Passed=True
1340     subsystemfailed=[]
1341 grchrist 1.11
1342 grchrist 1.8 if num_ls==1:
1343     ##fit is 2 ls ahead of real rate
1344 grchrist 1.10 LS=Rates[print_trigger]["ls"][iterator]
1345     #print "ls=",LS,
1346     LSRange=LumiPageInfo[LS]["LSRange"]
1347     #print LSRange,
1348     LS2=LSRange[-1]
1349     #LS2=LSRange.pop()
1350     #print "LS2=",LS2
1351     #print LumiPageInfo[LS]
1352     lumidict={}
1353     lumidict=LumiPageInfo[LS]
1354 grchrist 1.11
1355     if print_info:
1356     if (iterator==0 and print_trigger==trig_list[0]):
1357 grchrist 1.69 print '%10s%10s%10s%10s%10s%10s%10s%15s%20s%15s' % ("Status", "Run", "LS", "Physics", "Active", "Deadtime", " MaxDeadTime", " Passed all subsystems?", " List of Subsystems", " Spike killing")
1358 grchrist 1.11
1359     ## if SubSystemOff["All"]:
1360     ## for keys in LumiPageInfo[LS]:
1361     ## #print LS, keys, LumiPageInfo[LS][keys]
1362     ## if not LumiPageInfo[LS][keys]:
1363     ## Passed=False
1364     ## subsystemfailed.append(keys)
1365     ## break
1366     ## else:
1367     if SubSystemOff["Mu"] or SubSystemOff["All"]:
1368 grchrist 1.38 if not (LumiPageInfo[LS]["rpc"] and LumiPageInfo[LS]["dt0"] and LumiPageInfo[LS]["dtp"] and LumiPageInfo[LS]["dtm"] and LumiPageInfo[LS]["cscp"] and LumiPageInfo[LS]["cscm"]):
1369 grchrist 1.11 Passed=False
1370     subsystemfailed.append("Mu")
1371     if SubSystemOff["HCal"] or SubSystemOff["All"]:
1372     if not (LumiPageInfo[LS]["hbhea"] and LumiPageInfo[LS]["hbheb"] and LumiPageInfo[LS]["hbhec"]):
1373     Passed=False
1374     subsystemfailed.append("HCal")
1375     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["hf"]):
1376     Passed=False
1377     subsystemfailed.append("HCal-EndCap")
1378     if SubSystemOff["ECal"] or SubSystemOff["All"]:
1379     if not (LumiPageInfo[LS]["ebp"] and LumiPageInfo[LS]["ebm"]):
1380     Passed=False
1381     subsystemfailed.append("ECal")
1382     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["eep"] and LumiPageInfo[LS]["eem"] and LumiPageInfo[LS]["esp"] or LumiPageInfo[LS]["esm"]):
1383     Passed=False
1384     subsystemfailed.append("ECal-EndCap")
1385     if SubSystemOff["Tracker"] or SubSystemOff["All"]:
1386     if not (LumiPageInfo[LS]["tob"] and LumiPageInfo[LS]["tibtid"] and LumiPageInfo[LS]["bpix"] and LumiPageInfo[LS]["fpix"]):
1387     Passed=False
1388     subsystemfailed.append("Tracker")
1389     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["tecp"] and LumiPageInfo[LS]["tecm"]):
1390     Passed=False
1391     subsystemfailed.append("Tracker-EndCap")
1392     if SubSystemOff["Beam"] or SubSystemOff["All"]:
1393     if not(LumiPageInfo[LS]["b1pres"] and LumiPageInfo[LS]["b2pres"] and LumiPageInfo[LS]["b1stab"] and LumiPageInfo[LS]["b2stab"]):
1394     Passed=False
1395     subsystemfailed.append("Beam")
1396 grchrist 1.12 else:
1397     Passed=True
1398 grchrist 1.10
1399 awoodard 1.72 if not data_clean or (
1400 grchrist 1.11 Rates[print_trigger]["physics"][iterator] == 1
1401 grchrist 1.8 and Rates[print_trigger]["active"][iterator] == 1
1402 grchrist 1.17 and Rates[print_trigger]["deadtime"][iterator] < max_dt
1403 grchrist 1.11 #and Rates[print_trigger]["psi"][iterator] > 0
1404 grchrist 1.10 and Passed
1405 abrinke1 1.67 and (realvalue >0.6*prediction and realvalue<1.5*prediction)
1406     and Rates[print_trigger]["rawrate"][iterator] > 0.04
1407 grchrist 1.8 ):
1408 grchrist 1.26 if (print_info and num_ls==1 and (realvalue <0.4*prediction or realvalue>2.5*prediction)):
1409     pass
1410 grchrist 1.59 ##print '%-50s%10s%10s%10s%10s%10s%10s%10s%15s%20s' % (print_trigger,"Passed", 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)
1411 grchrist 1.8 return True
1412     else:
1413 grchrist 1.13 if (print_info and print_trigger==trig_list[0] and num_ls==1):
1414 grchrist 1.69 prediction_match=True
1415     if (realvalue >0.6*prediction and realvalue<1.5*prediction):
1416     prediction_match=False
1417     print '%10s%10s%10s%10s%10s%10s%10s%15s%20s%15s' % ("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, prediction_match )
1418 grchrist 1.26
1419 grchrist 1.8 return False
1420 abrinke1 1.5
1421 grchrist 1.10
1422     #### LumiRangeGreens ####
1423     ####inputs: RefMoreLumiArray --dict with lumi page info in LS by LS blocks,
1424     #### LRange --list range over lumis,
1425     #### nls --number of lumisections
1426     #### RefRunNum --run number
1427     ####
1428     ####outputs RangeMoreLumi --lumi page info in dict LSRange blocks with lumi, added items Run and LSRange
1429 grchrist 1.17 def LumiRangeGreens(RefMoreLumiArray,LSRange,nls,RefRunNum,deadtimebeamactive):
1430 grchrist 1.9
1431     RangeMoreLumi={}
1432     for keys,values in RefMoreLumiArray.iteritems():
1433     RangeMoreLumi[keys]=1
1434 grchrist 1.10
1435 grchrist 1.9 for iterator in LSRange[nls]:
1436     for keys, values in RefMoreLumiArray.iteritems():
1437     if RefMoreLumiArray[keys][iterator]==0:
1438     RangeMoreLumi[keys]=0
1439 grchrist 1.10 RangeMoreLumi['LSRange']=LSRange[nls]
1440     RangeMoreLumi['Run']=RefRunNum
1441 grchrist 1.17 RangeMoreLumi['DeadTimeBeamActive']=deadtimebeamactive
1442 grchrist 1.9 return RangeMoreLumi
1443    
1444 grchrist 1.10 #### CheckLumis ####
1445     ####inputs:
1446     #### PageLumiInfo --dict of LS with dict of some lumipage info
1447     #### Rates --dict of triggernames with dict of info
1448     def checkLS(Rates, PageLumiInfo,trig_list):
1449     rateslumis=Rates[trig_list[-1]]["ls"]
1450     keys=PageLumiInfo.keys()
1451     print "lumi run=",PageLumiInfo[keys[-1]]["Run"]
1452     ll=0
1453     for ls in keys:
1454     print ls,rateslumis[ll]
1455     ll=ll+1
1456     return False
1457    
1458    
1459 grchrist 1.71 def checkL1seedChangeALLPScols(trig_list,HLTL1PS):
1460     L1PSchangedic={}
1461     nps=0
1462     for HLTkey in trig_list:
1463     if HLTkey=='HLT_Stream_A':
1464     continue
1465     #print HLTkey
1466     try:
1467     dict=HLTL1PS[StripVersion(HLTkey)]
1468     #print "dict=",dict
1469     except:
1470     #print HLTkey, StripVersion(HLTkey)
1471     exit(2)
1472     HLTL1dummy={}
1473     for L1seed in dict.iterkeys():
1474     #print L1seed
1475     dummyL1seedlist=[]
1476     #print dict[L1seed]
1477     dummy=dict[L1seed]
1478     L1seedchangedummy=[]
1479     L1fulldummy=[]
1480     nps=len(dict[L1seed])
1481     #print "nps=",nps
1482     for PScol in range(0,len(dict[L1seed])):
1483     PScoldummy=PScol+1
1484     if PScoldummy>(len(dict)-1):
1485     PScoldummy=len(dict)-1
1486     #print PScol, PScoldummy, dummy[PScol]
1487    
1488     if dummy[PScol]==dummy[PScoldummy]:
1489     L1seedchangedummy.append(PScol)
1490     else:
1491     L1seedchangedummy.append(PScol)
1492     for ps in L1seedchangedummy:
1493     L1fulldummy.append(L1seedchangedummy)
1494     #print "L1seed change ", L1seedchangedummy, "full=",L1fulldummy
1495     L1seedchangedummy=[]
1496     for ps in L1seedchangedummy:
1497     L1fulldummy.append(L1seedchangedummy)
1498     #print "L1full=",L1fulldummy
1499     HLTL1dummy[L1seed]=L1fulldummy
1500     #print HLTL1dummy
1501     HLTL1_seedchanges=commonL1PS(HLTL1dummy,nps)
1502     print HLTkey, HLTL1_seedchanges
1503    
1504     def commonL1PS(HLTL1dummy, nps):
1505     ### find commmon elements in L1 seeds
1506     HLTL1_seedchanges=[]
1507     for PScol in range(0,nps):
1508    
1509     L1seedslist=HLTL1dummy.keys()
1510     L1tupletmp=set(tuple(HLTL1dummy[L1seedslist.pop()][PScol]))
1511     while len(L1seedslist)>0:
1512     L1tupletmp2=set(tuple(HLTL1dummy[L1seedslist.pop()][PScol]))
1513     L1tupletmp=L1tupletmp & L1tupletmp2
1514     if list(tuple(L1tupletmp)) not in HLTL1_seedchanges:
1515     HLTL1_seedchanges.append(list(tuple(L1tupletmp)))
1516     #print HLTL1_seedchanges
1517     return HLTL1_seedchanges
1518 grchrist 1.9
1519 grchrist 1.29
1520 abrinke1 1.1 if __name__=='__main__':
1521 grchrist 1.29 global thisyear
1522 abrinke1 1.1 main()