ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.72
Committed: Thu Nov 15 14:29:05 2012 UTC (12 years, 5 months ago) by awoodard
Content type: text/x-python
Branch: MAIN
Changes since 1.71: +57 -49 lines
Log Message:
improving treatment of paths without fits; adding min chi^2/ndof

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.70 DoL1=True
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 awoodard 1.50 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.10
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.59 if DoL1:
453     L1HLTseeds=RefParser.GetL1HLTseeds()
454     for HLTkey in trig_list:
455     #print name
456     if "L1" in HLTkey:
457     continue
458     else:
459     try:
460     for L1seed in L1HLTseeds[HLTkey]:
461     if L1seed not in trig_list:
462     trig_list.append(L1seed)
463     except:
464     pass
465 abrinke1 1.5 for iterator in RefLumiArray[0]: ##Makes array of LS with proper PAP and JSON properties
466 grchrist 1.11 ##cheap way of getting PSCol None-->0
467 amott 1.18 if RefLumiArray[0][iterator] not in range(1,9):
468 grchrist 1.11 RefLumiArray[0][iterator]=0
469 grchrist 1.15
470 grchrist 1.61 if not UsePSCol==-1:
471     if not RefLumiArray[0][iterator]==UsePSCol:
472     print "skipping LS",iterator
473     continue
474    
475 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):
476 abrinke1 1.5 if not JSON or RefRunNum in JSON:
477     if not JSON or iterator in JSON[RefRunNum]:
478     RefLumiRange.append(iterator)
479 grchrist 1.9
480 abrinke1 1.5 try:
481     nls = RefLumiRange[0]
482     LSRange = {}
483     except:
484     print "Run "+str(RefRunNum)+" has no good LS"
485     continue
486     if num_ls > len(RefLumiRange):
487 abrinke1 1.1 print "Run "+str(RefRunNum)+" is too short: from "+str(nls)+" to "+str(RefLumiRange[-1])+", while num_ls = "+str(num_ls)
488     continue
489     while nls < RefLumiRange[-1]-num_ls:
490 abrinke1 1.5 LSRange[nls] = []
491     counter = 0
492     for iterator in RefLumiRange:
493     if iterator >= nls and counter < num_ls:
494     LSRange[nls].append(iterator)
495     counter += 1
496 abrinke1 1.1 nls = LSRange[nls][-1]+1
497 grchrist 1.71 ###checkL1seedChangeALLPScols(trig_list,HLTL1PS) #for L1prescale changes
498 grchrist 1.17 #print "Run "+str(RefRunNum)+" contains LS from "+str(min(LSRange))+" to "+str(max(LSRange))
499 grchrist 1.8 for nls in sorted(LSRange.iterkeys()):
500 abrinke1 1.1 TriggerRates = RefParser.GetHLTRates(LSRange[nls])
501 grchrist 1.59 #L1Rate=RefParser.GetDeadTimeBeamActive(LSRange[nls])
502 awoodard 1.47
503     ## Clumsy way to append Stream A. Should choose correct method for calculating stream a based on ps column used in data taking.
504 grchrist 1.59
505 awoodard 1.57 if ('HLT_Stream_A' in trig_list) or all_triggers:
506     config = RateMonConfig(os.path.abspath(os.path.dirname(sys.argv[0])))
507     config.ReadCFG()
508     stream_mon = StreamMonitor()
509     core_a_rates = stream_mon.getStreamACoreRatesByLS(RefParser,LSRange[nls],config).values()
510     avg_core_a_rate = sum(core_a_rates)/len(LSRange[nls])
511     TriggerRates['HLT_Stream_A'] = [1,1,avg_core_a_rate,avg_core_a_rate]
512 grchrist 1.59
513     if DoL1:
514     L1RatesALL=RefParser.GetL1RatesALL(LSRange[nls])
515    
516     for L1seed in L1RatesALL.iterkeys():
517     TriggerRates[L1seed]=L1RatesALL[L1seed]
518 abrinke1 1.5
519     [inst, live, delivered, dead, pscols] = RefParser.GetAvLumiInfo(LSRange[nls])
520 grchrist 1.17 deadtimebeamactive=RefParser.GetDeadTimeBeamActive(LSRange[nls])
521 grchrist 1.59
522 abrinke1 1.2 physics = 1
523     active = 1
524 abrinke1 1.5 psi = 99
525     for iterator in LSRange[nls]: ##Gets lowest value of physics, active, and psi in the set of lumisections
526 abrinke1 1.2 if RefLumiArray[5][iterator] == 0:
527     physics = 0
528     if RefLumiArray[6][iterator] == 0:
529     active = 0
530 abrinke1 1.5 if RefLumiArray[0][iterator] < psi:
531     psi = RefLumiArray[0][iterator]
532 abrinke1 1.2
533 grchrist 1.17 MoreLumiMulti=LumiRangeGreens(RefMoreLumiArray,LSRange,nls,RefRunNum,deadtimebeamactive)
534 grchrist 1.10
535 abrinke1 1.5 if inst < 0 or live < 0 or delivered < 0:
536     print "Run "+str(RefRunNum)+" LS "+str(nls)+" inst lumi = "+str(inst)+" live lumi = "+str(live)+", delivered = "+str(delivered)+", physics = "+str(physics)+", active = "+str(active)
537 grchrist 1.10
538     LumiPageInfo[nls]=MoreLumiMulti
539    
540 abrinke1 1.1 for key in TriggerRates:
541 grchrist 1.59
542 grchrist 1.34 if NoVersion:
543     name = StripVersion(key)
544     else:
545     name=key
546 grchrist 1.59
547 grchrist 1.26 if not name in trig_list:
548 awoodard 1.50 if all_triggers:
549     trig_list.append(name)
550     else:
551     continue
552 grchrist 1.59
553 abrinke1 1.1 if not Rates.has_key(name):
554     Rates[name] = {}
555     Rates[name]["run"] = []
556     Rates[name]["ls"] = []
557     Rates[name]["ps"] = []
558     Rates[name]["inst_lumi"] = []
559     Rates[name]["live_lumi"] = []
560     Rates[name]["delivered_lumi"] = []
561     Rates[name]["deadtime"] = []
562     Rates[name]["rawrate"] = []
563     Rates[name]["rate"] = []
564     Rates[name]["rawxsec"] = []
565     Rates[name]["xsec"] = []
566 abrinke1 1.2 Rates[name]["physics"] = []
567     Rates[name]["active"] = []
568 abrinke1 1.5 Rates[name]["psi"] = []
569 grchrist 1.9
570 abrinke1 1.1 [avps, ps, rate, psrate] = TriggerRates[key]
571     Rates[name]["run"].append(RefRunNum)
572     Rates[name]["ls"].append(nls)
573     Rates[name]["ps"].append(ps)
574     Rates[name]["inst_lumi"].append(inst)
575     Rates[name]["live_lumi"].append(live)
576     Rates[name]["delivered_lumi"].append(delivered)
577 grchrist 1.17 Rates[name]["deadtime"].append(deadtimebeamactive)
578 abrinke1 1.1 Rates[name]["rawrate"].append(rate)
579 abrinke1 1.2 if live == 0:
580 awoodard 1.56 Rates[name]["rate"].append(0.0)
581 abrinke1 1.2 Rates[name]["rawxsec"].append(0.0)
582     Rates[name]["xsec"].append(0.0)
583     else:
584 awoodard 1.56 try:
585     Rates[name]["rate"].append(psrate/(1.0-deadtimebeamactive))
586     except:
587     Rates[name]["rate"].append(0.0)
588 abrinke1 1.2 Rates[name]["rawxsec"].append(rate/live)
589     Rates[name]["xsec"].append(psrate/live)
590     Rates[name]["physics"].append(physics)
591     Rates[name]["active"].append(active)
592 abrinke1 1.5 Rates[name]["psi"].append(psi)
593 grchrist 1.9
594 grchrist 1.10 #for keys, values in MoreLumiMulti.iteritems():
595     # Rates[name][keys].append(values)
596 grchrist 1.9 #print nls, name, keys, values
597 awoodard 1.47
598 grchrist 1.10 RateOutput = open(RefRunFile, 'wb') ##Save new Rates[] to RefRuns
599     pickle.dump(Rates, RateOutput, 2)
600     RateOutput.close()
601     LumiOutput = open(RefLumiFile,'wb')
602     pickle.dump(LumiPageInfo,LumiOutput, 2)
603     LumiOutput.close()
604    
605 grchrist 1.59 #for trig in Rates.iterkeys():
606     # print trig, Rates[trig]
607     return [Rates,LumiPageInfo,trig_list]
608 abrinke1 1.1
609 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):
610 abrinke1 1.1 min_run = min(run_list)
611     max_run = max(run_list)
612 muell149 1.46
613     priot.has_been_called=False
614 grchrist 1.11
615 abrinke1 1.5 InputFit = {}
616     OutputFit = {}
617 awoodard 1.72 failed_paths = []
618 grchrist 1.14 first_trigger=True
619 grchrist 1.37
620     [[varX, varY, do_fit, save_root, save_png, fit_file]] = plot_properties
621    
622     RootNameTemplate = "%s_%sLS_%s_vs_%s_Run%s-%s.root"
623     RootFile = RootNameTemplate % (trig_name, num_ls, varX, varY, min_run, max_run)
624 abrinke1 1.1
625 amott 1.20 if not do_fit:
626     try:
627 abrinke1 1.1 pkl_file = open(fit_file, 'rb')
628 abrinke1 1.5 InputFit = pickle.load(pkl_file)
629 grchrist 1.40
630 amott 1.20 except:
631     print "ERROR: could not open fit file: %s" % (fit_file,)
632     if save_root:
633     try:
634     os.remove(RootFile)
635     except:
636     pass
637    
638 grchrist 1.38 trig_list_noV=[]
639     for trigs in trig_list:
640     trig_list_noV.append(StripVersion(trigs))
641     if NoVersion:
642     trig_list=trig_list_noV
643 grchrist 1.40
644 amott 1.20 ## check that all the triggers we ask to plot are in the input fit
645     if not save_fits:
646     goodtrig_list = []
647 grchrist 1.40 FitInputNoV={}
648 amott 1.20 for trig in trig_list:
649 grchrist 1.40
650     if NoVersion:
651     for trigger in InputFit.iterkeys():
652     FitInputNoV[StripVersion(trigger)]=InputFit[trigger]
653     InputFit=FitInputNoV
654 awoodard 1.56
655 amott 1.20 else:
656 grchrist 1.40 if not InputFit.has_key(trig):
657     print "WARNING: No Fit Prediction for Trigger %s, SKIPPING" % (trig,)
658     else:
659     goodtrig_list.append(trig)
660     trig_list = goodtrig_list
661 amott 1.18
662 grchrist 1.26 for print_trigger in sorted(Rates):
663 abrinke1 1.7 ##Limits Rates[] to runs in run_list
664     NewTrigger = {}
665 grchrist 1.34
666 grchrist 1.8 if not print_trigger in trig_list:
667 awoodard 1.50 if all_triggers:
668     trig_list.append(print_trigger)
669     else:
670     print "not in trig_list:",print_trigger, trig_list
671     continue
672 grchrist 1.34
673 abrinke1 1.7 for key in Rates[print_trigger]:
674     NewTrigger[key] = []
675     for iterator in range (len(Rates[print_trigger]["run"])):
676     if Rates[print_trigger]["run"][iterator] in run_list:
677     for key in Rates[print_trigger]:
678     NewTrigger[key].append(Rates[print_trigger][key][iterator])
679     Rates[print_trigger] = NewTrigger
680    
681 grchrist 1.59 #print print_trigger, Rates[print_trigger]
682 abrinke1 1.1 meanrawrate = sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])
683 grchrist 1.40
684 abrinke1 1.1 if not trig_name in print_trigger:
685 grchrist 1.59 #print "failed does not contain HLT",trig_name, print_trigger
686     pass
687 grchrist 1.40
688 abrinke1 1.1 if meanrawrate < min_rate:
689     continue
690     masked_trig = False
691     for mask in masked_triggers:
692     if str(mask) in print_trigger:
693     masked_trig = True
694     if masked_trig:
695     continue
696 grchrist 1.34
697 abrinke1 1.5 OutputFit[print_trigger] = {}
698 abrinke1 1.1
699     lowlumi = 0
700 abrinke1 1.7 meanlumi_init = median(Rates[print_trigger]["live_lumi"])
701 abrinke1 1.1 meanlumi = 0
702     highlumi = 0
703 abrinke1 1.67 lowrate = 0
704     meanrate = 0
705     highrate = 0
706 abrinke1 1.1 lowxsec = 0
707     meanxsec = 0
708     highxsec = 0
709     nlow = 0
710     nhigh = 0
711 grchrist 1.15
712 abrinke1 1.1 for iterator in range(len(Rates[print_trigger]["rate"])):
713     if Rates[print_trigger]["live_lumi"][iterator] <= meanlumi_init:
714 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):
715 abrinke1 1.67 meanrate+=Rates[print_trigger]["rate"][iterator]
716     lowrate+=Rates[print_trigger]["rate"][iterator]
717 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
718     lowxsec+=Rates[print_trigger]["xsec"][iterator]
719     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
720     lowlumi+=Rates[print_trigger]["live_lumi"][iterator]
721     nlow+=1
722     if Rates[print_trigger]["live_lumi"][iterator] > meanlumi_init:
723 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):
724 abrinke1 1.67 meanrate+=Rates[print_trigger]["rate"][iterator]
725     highrate+=Rates[print_trigger]["rate"][iterator]
726 abrinke1 1.1 meanxsec+=Rates[print_trigger]["xsec"][iterator]
727     highxsec+=Rates[print_trigger]["xsec"][iterator]
728     meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
729     highlumi+=Rates[print_trigger]["live_lumi"][iterator]
730     nhigh+=1
731 abrinke1 1.7 try:
732 abrinke1 1.67 meanrate = meanrate/(nlow+nhigh)
733 abrinke1 1.7 meanxsec = meanxsec/(nlow+nhigh)
734     meanlumi = meanlumi/(nlow+nhigh)
735 abrinke1 1.67 if (nlow==0):
736     sloperate = (highrate/nhigh) / (highlumi/nhigh)
737     slopexsec = (highxsec/nhigh) / (highlumi/nhigh)
738     elif (nhigh==0):
739     sloperate = (lowrate/nlow) / (lowlumi/nlow)
740     slopexsec = (lowxsec/nlow) / (lowlumi/nlow)
741     else:
742     sloperate = ( (highrate/nhigh) - (lowrate/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
743     slopexsec = ( (highxsec/nhigh) - (lowxsec/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
744 abrinke1 1.7 except:
745 awoodard 1.72 # print str(print_trigger)+" has no good datapoints - setting initial xsec slope estimate to 0"
746 abrinke1 1.67 meanrate = median(Rates[print_trigger]["rate"])
747 abrinke1 1.7 meanxsec = median(Rates[print_trigger]["xsec"])
748     meanlumi = median(Rates[print_trigger]["live_lumi"])
749 abrinke1 1.67 sloperate = meanxsec
750 abrinke1 1.7 slopexsec = 0
751 abrinke1 1.1
752 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()
753    
754 amott 1.20 if not do_fit:
755 awoodard 1.62 try:
756     FitType = InputFit[print_trigger][0]
757 awoodard 1.72 except:
758     failed_paths.append([print_trigger,"This path did not exist in the monitorlist used to create the fit"])
759     continue
760     if FitType == "fit failed":
761     failure_comment = InputFit[print_trigger][1]
762     failed_paths.append([print_trigger, failure_comment])
763     continue
764     else:
765 awoodard 1.62 X0 = InputFit[print_trigger][1]
766     X1 = InputFit[print_trigger][2]
767     X2 = InputFit[print_trigger][3]
768     X3 = InputFit[print_trigger][4]
769     sigma = InputFit[print_trigger][5]/math.sqrt(num_ls)*3#Display 3 sigma band to show outliers more clearly
770     X0err= InputFit[print_trigger][7]
771    
772 grchrist 1.8 ## we are 2 lumis off when we start! -gets worse when we skip lumis
773 grchrist 1.17 it_offset=0
774 abrinke1 1.1 for iterator in range(len(Rates[print_trigger]["rate"])):
775     if not Rates[print_trigger]["run"][iterator] in run_list:
776     continue
777 abrinke1 1.67 ##Old Spike-killer: used average-xsec method, too clumsy, esp. for nonlinear triggers
778     #prediction = meanxsec + slopexsec * (Rates[print_trigger]["live_lumi"][iterator] - meanlumi)
779     #realvalue = Rates[print_trigger]["xsec"][iterator]
780    
781     ##New Spike-killer: gets average of nearest 4 LS
782 grchrist 1.68 try:
783     if (iterator > 2 and iterator+2 < len(Rates[print_trigger]["rate"])): ##2 LS before, 2 after
784     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
785     elif (iterator > 2 and len(Rates[print_trigger]["rate"]) > 4 ): ##4 LS before
786     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
787     elif (iterator+2 < len(Rates[print_trigger]["rate"]) and len(Rates[print_trigger]["rate"]) > 4 ): ##4 LS after
788     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
789     else:
790     prediction = Rates[print_trigger]["rate"][iterator]
791     realvalue = Rates[print_trigger]["rate"][iterator]
792     except:
793     print "Error calculating prediction. Setting rates to defaults."
794 abrinke1 1.67 prediction = Rates[print_trigger]["rate"][iterator]
795 grchrist 1.68 #realvalue = Rates[print_trigger]["rate"][iterator]
796 grchrist 1.11
797     if pass_cuts(data_clean, realvalue, prediction, meanxsec, Rates, print_trigger, iterator, num_ls,LumiPageInfo,SubSystemOff,max_dt,print_info, trig_list):
798 grchrist 1.40
799 abrinke1 1.1 run_t.append(Rates[print_trigger]["run"][iterator])
800     ls_t.append(Rates[print_trigger]["ls"][iterator])
801     ps_t.append(Rates[print_trigger]["ps"][iterator])
802     inst_t.append(Rates[print_trigger]["inst_lumi"][iterator])
803 grchrist 1.17 live_t.append(Rates[print_trigger]["live_lumi"][iterator])
804     delivered_t.append(Rates[print_trigger]["delivered_lumi"][iterator])
805     deadtime_t.append(Rates[print_trigger]["deadtime"][iterator])
806 abrinke1 1.1 rawrate_t.append(Rates[print_trigger]["rawrate"][iterator])
807     rate_t.append(Rates[print_trigger]["rate"][iterator])
808     rawxsec_t.append(Rates[print_trigger]["rawxsec"][iterator])
809     xsec_t.append(Rates[print_trigger]["xsec"][iterator])
810 abrinke1 1.5 psi_t.append(Rates[print_trigger]["psi"][iterator])
811 abrinke1 1.1
812     e_run_t.append(0.0)
813     e_ls_t.append(0.0)
814     e_ps_t.append(0.0)
815     e_inst_t.append(14.14)
816     e_live_t.append(14.14)
817     e_delivered_t.append(14.14)
818 abrinke1 1.2 e_deadtime_t.append(0.01)
819 abrinke1 1.1 e_rawrate_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
820     e_rate_t.append(Rates[print_trigger]["ps"][iterator]*math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
821 abrinke1 1.5 e_psi_t.append(0.0)
822 abrinke1 1.2 if live_t[-1] == 0:
823     e_rawxsec_t.append(0)
824     e_xsec_t.append(0)
825     else:
826 grchrist 1.16 try:
827     e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["live_lumi"][iterator])
828     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])
829     except:
830     e_rawxsec_t.append(0.)
831     e_xsec_t.append(0.)
832 amott 1.20 if not do_fit:
833 grchrist 1.38 if not do_inst:
834     if FitType == "expo":
835 grchrist 1.43 rate_prediction = X0 + X1*math.exp(X2+X3*delivered_t[-1])
836 grchrist 1.38 else:
837     rate_prediction = X0 + X1*delivered_t[-1] + X2*delivered_t[-1]*delivered_t[-1] + X3*delivered_t[-1]*delivered_t[-1]*delivered_t[-1]
838 abrinke1 1.67
839 grchrist 1.38 else:
840     if FitType == "expo":
841 grchrist 1.43 rate_prediction = X0 + X1*math.exp(X2+X3*inst_t[-1])
842 grchrist 1.38 else:
843     rate_prediction = X0 + X1*inst_t[-1] + X2*inst_t[-1]*inst_t[-1] + X3*inst_t[-1]*inst_t[-1]*inst_t[-1]
844 abrinke1 1.5
845 abrinke1 1.2 if live_t[-1] == 0:
846 abrinke1 1.5 rawrate_fit_t.append(0)
847     rate_fit_t.append(0)
848 abrinke1 1.2 rawxsec_fit_t.append(0)
849     xsec_fit_t.append(0)
850 abrinke1 1.7 e_rawrate_fit_t.append(0)
851 awoodard 1.45 e_rate_fit_t.append(sigma)
852 abrinke1 1.7 e_rawxsec_fit_t.append(0)
853     e_xsec_fit_t.append(0)
854 abrinke1 1.67
855 abrinke1 1.2 else:
856 grchrist 1.11 if ps_t[-1]>0.0:
857     rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
858     else:
859     rawrate_fit_t.append(0.0)
860    
861 awoodard 1.52 rate_fit_t.append(rate_prediction)
862 awoodard 1.45 e_rate_fit_t.append(sigma)
863 abrinke1 1.5 rawxsec_fit_t.append(rawrate_fit_t[-1]/live_t[-1])
864     xsec_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/live_t[-1])
865 grchrist 1.38 try:
866    
867     if not TMDerr:
868 awoodard 1.45 e_rawrate_fit_t.append(sigma*rawrate_fit_t[-1]/rate_fit_t[-1])
869     e_rawxsec_fit_t.append(sigma*rawxsec_fit_t[-1]/rate_fit_t[-1])
870     e_xsec_fit_t.append(sigma*xsec_fit_t[-1]/rate_fit_t[-1])
871 grchrist 1.38 ###error from TMD predictions, calculated at 5e33
872     else:
873     e_rawrate_fit_t.append(X0err*inst_t[-1]/5000.)
874     e_rawxsec_fit_t.append(X0err/live_t[-1]*inst_t[-1]/5000.)
875     e_xsec_fit_t.append(X0err/live_t[-1]*inst_t[-1]/5000.)
876    
877     except:
878     print print_trigger, "has no fitted rate for LS", Rates[print_trigger]["ls"][iterator]
879 awoodard 1.45 e_rawrate_fit_t.append(sigma)
880     e_rawxsec_fit_t.append(sigma)
881     e_xsec_fit_t.append(sigma)
882 abrinke1 1.5
883 grchrist 1.16
884 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"])))):
885 grchrist 1.38 pass
886 grchrist 1.8
887 abrinke1 1.5 else: ##If the data point does not pass the data_clean filter
888 abrinke1 1.1 if debug_print:
889     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)
890    
891 abrinke1 1.5 ## End "for iterator in range(len(Rates[print_trigger]["rate"])):" loop
892 abrinke1 1.1
893 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]
894 grchrist 1.40
895 awoodard 1.63 [VX, VXE, x_label, VY, VYE, y_label, VF, VFE] = GetVXVY(plot_properties, fit_file, AllPlotArrays)
896 grchrist 1.40
897 abrinke1 1.5 if save_root or save_png:
898     c1 = TCanvas(str(varX),str(varY))
899     c1.SetName(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
900 grchrist 1.40
901     try:
902     gr1 = TGraphErrors(len(VX), VX, VY, VXE, VYE)
903     except:
904 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"
905     failed_paths.append([print_trigger,failure_comment])
906     if do_fit:
907     OutputFit[print_trigger] = ["fit failed",failure_comment]
908 abrinke1 1.67 continue
909 awoodard 1.72 if (len(VX)<20 and do_fit):
910     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"
911     failed_paths.append([print_trigger,failure_comment])
912     OutputFit[print_trigger] = ["fit failed",failure_comment]
913 grchrist 1.40 continue
914 awoodard 1.72
915 abrinke1 1.5 gr1.SetName("Graph_"+str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
916 awoodard 1.63 gr1.GetXaxis().SetTitle(x_label)
917     gr1.GetYaxis().SetTitle(y_label)
918 abrinke1 1.5 gr1.SetTitle(str(print_trigger))
919     gr1.SetMinimum(0)
920     gr1.SetMaximum(1.2*max(VY))
921     #gr1.GetXaxis().SetLimits(min(VX)-0.2*max(VX),1.2*max(VX))
922     gr1.GetXaxis().SetLimits(0,1.2*max(VX))
923     gr1.SetMarkerStyle(8)
924 grchrist 1.40
925 abrinke1 1.5 if fit_file:
926     gr1.SetMarkerSize(0.8)
927     else:
928     gr1.SetMarkerSize(0.5)
929     gr1.SetMarkerColor(2)
930    
931 amott 1.20 if not do_fit:
932 abrinke1 1.5 gr3 = TGraphErrors(len(VX), VX, VF, VXE, VFE)
933     gr3.SetMarkerStyle(8)
934     gr3.SetMarkerSize(0.4)
935     gr3.SetMarkerColor(4)
936     gr3.SetFillColor(4)
937     gr3.SetFillStyle(3003)
938 abrinke1 1.1
939 abrinke1 1.5 if do_fit:
940 grchrist 1.35 if "rate" in varY and not linear:
941 grchrist 1.61 f1d=0
942     f1d = TF1("f1d","pol1",0,8000)#linear
943 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)
944 grchrist 1.61 f1d.SetLineColor(4)
945     f1d.SetLineWidth(2)
946 abrinke1 1.67 if nlow>0:
947     f1d.SetParLimits(0,0,1.5*lowrate/nlow) ##Keep Y-intercept in range of low-lumi rate points
948     else:
949     f1d.SetParLimits(0,0,1.5*sum(VY)/len(VY))
950     if (sloperate > 0):
951     if (sloperate > 0.5*sum(VY)/sum(VX)): ##Slope is substantially positive
952     f1d.SetParLimits(1,min(0.5*sloperate,0.5*sum(VY)/sum(VX)),1.5*sum(VY)/sum(VX))
953     else: ##Slope is somewhat positive or flat
954     f1d.SetParLimits(1,-0.1*sloperate,1.5*sum(VY)/sum(VX))
955     else: ##Slope is negative or flat
956     f1d.SetParLimits(1,1.5*sloperate,-0.1*sloperate)
957     gr1.Fit("f1d","QN","rob=0.90")
958     f1d_Chi2 = f1d.GetChisquare()/f1d.GetNDF()
959 grchrist 1.61
960 grchrist 1.41 f1a=0
961 abrinke1 1.67 f1a = TF1("f1a","pol2",0,8000)#quadratic
962     f1a.SetParameters(f1d.GetParameter(0),f1d.GetParameter(1),0) ##Initial values from linear fit
963 grchrist 1.61 f1a.SetLineColor(6)
964 abrinke1 1.5 f1a.SetLineWidth(2)
965 abrinke1 1.67 if nlow>0 and sloperate < 0.5*sum(VY)/sum(VX): ##Slope is not substantially positive
966     f1a.SetParLimits(0,0,1.5*lowrate/nlow) ##Keep Y-intercept in range of low-lumi rate points
967     else:
968     f1a.SetParLimits(0,0,max(min(VY),0.3*sum(VY)/len(VY))) ##Keep Y-intercept reasonably low
969     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
970     f1a.SetParLimits(2,-2.0*max(VY)/(max(VX)*max(VX)),2.0*max(VY)/(max(VX)*max(VX))) ##Reasonable bounds
971 awoodard 1.48 gr1.Fit("f1a","QN","rob=0.90")
972 abrinke1 1.67 f1a_Chi2 = f1a.GetChisquare()/f1a.GetNDF()
973     f1a_BadMinimum = (f1a.GetMinimumX(5,7905,10)>2000 and f1a.GetMinimumX(5,7905,10)<7000) ##Don't allow minimum between 2000 and 7000
974 muell149 1.46
975 abrinke1 1.5 f1b = 0
976     f1c = 0
977 muell149 1.46 meanps = median(Rates[print_trigger]["ps"])
978     av_rte = mean(VY)
979    
980 grchrist 1.43 if True:
981 abrinke1 1.67 f1b = TF1("f1b","pol3",0,8000)#cubic
982     f1b.SetParameters(f1a.GetParameter(0),f1a.GetParameter(1),f1a.GetParameter(2),0) ##Initial values from quadratic fit
983 grchrist 1.43 f1b.SetLineColor(2)
984     f1b.SetLineWidth(2)
985 abrinke1 1.67 f1b.SetParLimits(0,0,max(min(VY),0.3*sum(VY)/len(VY))) ##Keep Y-intercept reasonably low
986     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
987     f1b.SetParLimits(2,-2.0*max(VY)/(max(VX)*max(VX)),2.0*max(VY)/(max(VX)*max(VX))) ##Reasonable bounds
988     f1b.SetParLimits(3,0,2.0*max(VY)/(max(VX)*max(VX)*max(VX))) ##Reasonable bounds
989 awoodard 1.48 gr1.Fit("f1b","QN","rob=0.90")
990 grchrist 1.70 try:
991     f1b_Chi2 = f1b.GetChisquare()/f1b.GetNDF()
992     except ZeroDivisionError:
993     f1b_Chi2=10.0
994     print "Zero ndf for ",print_trigger
995 awoodard 1.72
996     f1b_BadMinimum = (f1b.GetMinimumX(5,7905,10)>2000 and f1b.GetMinimumX(5,7905,10)<7000)
997 grchrist 1.43 f1c = TF1("f1c","[0]+[1]*expo(2)",0,8000)
998     f1c.SetLineColor(3)
999     f1c.SetLineWidth(2)
1000 abrinke1 1.67 #f1c.SetParLimits(0,0,max(min(VY),0.3*sum(VY)/len(VY)))
1001     f1c.SetParLimits(0,0,max(min(VY),0.01*sum(VY)/len(VY))) ##Exponential fits should start low
1002     f1c.SetParLimits(1,max(VY)/math.exp(15.0),max(VY)/math.exp(2.0))
1003 grchrist 1.43 f1c.SetParLimits(2,0.0,0.0000000001)
1004 abrinke1 1.67 f1c.SetParLimits(3,2.0/max(VX),15.0/max(VX))
1005 awoodard 1.48 gr1.Fit("f1c","QN","rob=0.90")
1006 abrinke1 1.67 f1c_Chi2 = f1c.GetChisquare()/f1c.GetNDF()
1007     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
1008     ##Some fits are so exponential, the graph ends early and returns a false low Chi2 value
1009 grchrist 1.15
1010 abrinke1 1.5 else: ##If this is not a rate plot
1011     f1a = TF1("f1a","pol1",0,8000)
1012     f1a.SetLineColor(4)
1013     f1a.SetLineWidth(2)
1014     if "xsec" in varY:
1015     f1a.SetParLimits(0,0,meanxsec*1.5)
1016     if slopexsec > 0:
1017     f1a.SetParLimits(1,0,max(VY)/max(VX))
1018     else:
1019     f1a.SetParLimits(1,2*slopexsec,-2*slopexsec)
1020 abrinke1 1.1 else:
1021 abrinke1 1.5 f1a.SetParLimits(0,-1000,1000)
1022     gr1.Fit("f1a","Q","rob=0.80")
1023 abrinke1 1.1
1024 grchrist 1.35 if (first_trigger):
1025 grchrist 1.59 print '%-50s %4s x0 x1 x2 x3 chi2 ndf chi2/ndf' % ('trigger', 'type')
1026 grchrist 1.35 first_trigger=False
1027 grchrist 1.40 try:
1028 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)
1029 grchrist 1.40 except:
1030     pass
1031 awoodard 1.72
1032     chioffset=1.0 ##chioffset now a fraction; must be 10% better to use expo rather than quad, quad rather than line
1033     width = max([len(trigger_name) for trigger_name in trig_list])
1034 awoodard 1.48 if print_table or save_fits:
1035     if not do_fit:
1036     print "Can't have save_fits = True and do_fit = False"
1037     continue
1038 awoodard 1.72 if min([f1a_Chi2,f1b_Chi2,f1c_Chi2,f1d_Chi2]) > 500:#require a minimum chi^2/nDOF of 500
1039     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"
1040     failed_paths.append([print_trigger,failure_comment])
1041     OutputFit[print_trigger] = ["fit failed",failure_comment]
1042     continue
1043     if "rate" in varY and not linear:
1044     if first_trigger:
1045     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")
1046     first_trigger = False
1047    
1048 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):
1049 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)
1050 awoodard 1.48 f1c.SetLineColor(1)
1051     priot(wp_bool,print_trigger,meanps,f1d,f1c,"expo",av_rte)
1052     sigma = CalcSigma(VX, VY, f1c)*math.sqrt(num_ls)
1053     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)]
1054    
1055 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):
1056 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)
1057 awoodard 1.48 f1b.SetLineColor(1)
1058     priot(wp_bool,print_trigger,meanps,f1d,f1b,"cubic",av_rte)
1059     sigma = CalcSigma(VX, VY, f1b)*math.sqrt(num_ls)
1060     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)]
1061    
1062 abrinke1 1.67 elif (f1a_Chi2 < (f1d_Chi2*chioffset)):
1063 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)
1064 awoodard 1.48 f1a.SetLineColor(1)
1065     priot(wp_bool,print_trigger,meanps,f1d,f1a,"quad",av_rte)
1066 awoodard 1.53 sigma = CalcSigma(VX, VY, f1a)*math.sqrt(num_ls)
1067 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]
1068 awoodard 1.72
1069 abrinke1 1.67 else:
1070 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)
1071 abrinke1 1.67 f1d.SetLineColor(1)
1072     priot(wp_bool,print_trigger,meanps,f1d,f1d,"line",av_rte)
1073     sigma = CalcSigma(VX, VY, f1d)*math.sqrt(num_ls)
1074     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]
1075 grchrist 1.65 else:
1076 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)
1077 grchrist 1.65 f1a.SetLineColor(1)
1078     #priot(wp_bool,print_trigger,meanps,f1d,f1a,"quad",av_rte)
1079     sigma = CalcSigma(VX, VY, f1a)*math.sqrt(num_ls)
1080     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]
1081 awoodard 1.48
1082 abrinke1 1.5 if save_root or save_png:
1083     gr1.Draw("APZ")
1084 amott 1.20 if not do_fit:
1085 abrinke1 1.5 gr3.Draw("P3")
1086 grchrist 1.58
1087 abrinke1 1.5 if do_fit:
1088 grchrist 1.59 try:
1089 abrinke1 1.67 if ((f1c_Chi2 < (f1a_Chi2*chioffset) or f1a_BadMinimum ) and (f1c_Chi2 < f1b_Chi2 or f1b_BadMinimum ) and not f1c_BadMinimum ):
1090 grchrist 1.59 f1c.Draw("same")
1091 abrinke1 1.67 elif ( (f1b_Chi2 < (f1a_Chi2*chioffset) or f1a_BadMinimum) and not f1b_BadMinimum):
1092 grchrist 1.59 f1b.Draw("same")
1093     else:
1094     f1a.Draw("same")
1095    
1096     f1d.Draw("same")
1097     except:
1098     True
1099 grchrist 1.58
1100 abrinke1 1.5 c1.Update()
1101     if save_root:
1102     myfile = TFile( RootFile, 'UPDATE' )
1103     c1.Write()
1104     myfile.Close()
1105     if save_png:
1106     c1.SaveAs(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX)+".png")
1107 awoodard 1.72
1108     if len(failed_paths) > 0:
1109     if save_fits:
1110     print "\n***************NO FIT RECORDED FOR THE FOLLOWING PATHS***************"
1111     else:
1112     print "\n***************THE FOLLOWING PATHS HAVE BEEN SKIPPED BECAUSE THE FIT WAS MISSING***************"
1113     sorted_failed_paths = sorted(failed_paths, key=itemgetter(1))
1114     for error_comment, entries in groupby(sorted_failed_paths, key=itemgetter(1)):
1115     error_comment = error_comment.replace('this path','these paths')
1116     print '\n'+error_comment+'.'
1117     for entry in entries:
1118     print entry[0]
1119    
1120 abrinke1 1.1 if save_root:
1121 awoodard 1.72 print "\nOutput root file is "+str(RootFile)
1122 abrinke1 1.1
1123     if save_fits:
1124 amott 1.20 if os.path.exists(fit_file):
1125     os.remove(fit_file)
1126     FitOutputFile = open(fit_file, 'wb')
1127 abrinke1 1.5 pickle.dump(OutputFit, FitOutputFile, 2)
1128     FitOutputFile.close()
1129 amott 1.20 print "Output fit file is "+str(fit_file)
1130 abrinke1 1.1
1131     if print_table:
1132 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."
1133 grchrist 1.59 print '%50s%10s%10s%10s%10s%10s%10s%10s' % ("Trig", "fit", "p0", "p1", "p2", "p3", "Chi2", "Av raw")
1134 abrinke1 1.5 for print_trigger in OutputFit:
1135     _trigger = (print_trigger[:56] + '...') if len(print_trigger) > 59 else print_trigger
1136 abrinke1 1.1 try:
1137 abrinke1 1.5 if OutputFit[print_trigger][0] == "poly":
1138 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))
1139 abrinke1 1.5 else:
1140 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))
1141 abrinke1 1.1 except:
1142     print str(print_trigger)+" is somehow broken"
1143 amott 1.18 return RootFile
1144 abrinke1 1.5
1145     ############# SUPPORTING FUNCTIONS ################
1146    
1147    
1148 awoodard 1.47 def CalcSigma(var_x, var_y, func):
1149     residuals = []
1150     for x, y in zip(var_x,var_y):
1151 awoodard 1.53 residuals.append(y - func.Eval(x,0,0))
1152 awoodard 1.47
1153     res_squared = [i*i for i in residuals]
1154 awoodard 1.53 if len(res_squared) > 2:
1155     sigma = math.sqrt(sum(res_squared)/(len(res_squared)-2))
1156     else:
1157     sigma = 0
1158 awoodard 1.47 return sigma
1159    
1160 abrinke1 1.5 def GetJSON(json_file):
1161    
1162     input_file = open(json_file)
1163     file_content = input_file.read()
1164     inputRange = selectionParser(file_content)
1165     JSON = inputRange.runsandls()
1166     return JSON
1167     ##JSON is an array: JSON[run_number] = [1st ls, 2nd ls, 3rd ls ... nth ls]
1168    
1169     def MakePlotArrays():
1170     run_t = array.array('f')
1171     ls_t = array.array('f')
1172     ps_t = array.array('f')
1173     inst_t = array.array('f')
1174     live_t = array.array('f')
1175     delivered_t = array.array('f')
1176     deadtime_t = array.array('f')
1177     rawrate_t = array.array('f')
1178     rate_t = array.array('f')
1179     rawxsec_t = array.array('f')
1180     xsec_t = array.array('f')
1181     psi_t = array.array('f')
1182    
1183     e_run_t = array.array('f')
1184     e_ls_t = array.array('f')
1185     e_ps_t = array.array('f')
1186     e_inst_t = array.array('f')
1187     e_live_t = array.array('f')
1188     e_delivered_t = array.array('f')
1189     e_deadtime_t = array.array('f')
1190     e_rawrate_t = array.array('f')
1191     e_rate_t = array.array('f')
1192     e_rawxsec_t = array.array('f')
1193     e_xsec_t = array.array('f')
1194     e_psi_t = array.array('f')
1195    
1196     rawrate_fit_t = array.array('f')
1197     rate_fit_t = array.array('f')
1198     rawxsec_fit_t = array.array('f')
1199     xsec_fit_t = array.array('f')
1200     e_rawrate_fit_t = array.array('f')
1201     e_rate_fit_t = array.array('f')
1202     e_rawxsec_fit_t = array.array('f')
1203     e_xsec_fit_t = array.array('f')
1204    
1205     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]
1206    
1207    
1208     def GetVXVY(plot_properties, fit_file, AllPlotArrays):
1209    
1210     VF = "0"
1211     VFE = "0"
1212    
1213     [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
1214     for varX, varY, do_fit, save_root, save_png, fit_file in plot_properties:
1215     if varX == "run":
1216     VX = run_t
1217     VXE = run_t_e
1218 awoodard 1.63 x_label = "Run Number"
1219 abrinke1 1.5 elif varX == "ls":
1220     VX = ls_t
1221     VXE = e_ls_t
1222 awoodard 1.63 x_label = "Lumisection"
1223 abrinke1 1.5 elif varX == "ps":
1224     VX = ps_t
1225     VXE = e_ps_t
1226 awoodard 1.63 x_label = "Prescale"
1227 abrinke1 1.5 elif varX == "inst":
1228     VX = inst_t
1229     VXE = e_inst_t
1230 abrinke1 1.67 x_label = "Instantaneous Luminosity [10^{30} Hz/cm^{2}]"
1231    
1232 abrinke1 1.5 elif varX == "live":
1233     VX = live_t
1234     VXE = e_live_t
1235 abrinke1 1.67 x_label = "Instantaneous Luminosity [10^{30} Hz/cm^{2}]"
1236    
1237 abrinke1 1.5 elif varX == "delivered":
1238     VX = delivered_t
1239     VXE = e_delivered_t
1240 abrinke1 1.67 x_label = "Instantaneous Luminosity [10^{30} Hz/cm^{2}]"
1241    
1242 abrinke1 1.5 elif varX == "deadtime":
1243     VX = deadtime_t
1244     VXE = e_deadtime_t
1245 awoodard 1.63 x_label = "Deadtime"
1246 abrinke1 1.5 elif varX == "rawrate":
1247     VX = rawrate_t
1248     VXE = e_rawrate_t
1249 awoodard 1.63 x_label = "Raw Rate [Hz]"
1250 abrinke1 1.5 elif varX == "rate":
1251     VX = rate_t
1252     VXE = e_rate_t
1253 awoodard 1.63 x_label = "Rate [Hz]"
1254 abrinke1 1.5 elif varX == "rawxsec":
1255     VX = rawxsec_t
1256     VXE = e_rawxsec_t
1257 grchrist 1.65 x_label = "Cross Section"
1258 abrinke1 1.5 elif varX == "xsec":
1259     VX = xsec_t
1260     VXE = e_xsec_t
1261 grchrist 1.65 x_label = "Cross Section"
1262 abrinke1 1.5 elif varX == "psi":
1263     VX = psi_t
1264     VXE = e_psi_t
1265 awoodard 1.63 x_label = "Prescale Index"
1266 abrinke1 1.5 else:
1267     print "No valid variable entered for X"
1268     continue
1269     if varY == "run":
1270     VY = run_t
1271     VYE = run_t_e
1272 awoodard 1.63 y_label = "Run Number"
1273 abrinke1 1.5 elif varY == "ls":
1274     VY = ls_t
1275     VYE = e_ls_t
1276 awoodard 1.63 y_label = "Lumisection"
1277 abrinke1 1.5 elif varY == "ps":
1278     VY = ps_t
1279     VYE = e_ps_t
1280 awoodard 1.63 y_label = "Prescale"
1281 abrinke1 1.5 elif varY == "inst":
1282     VY = inst_t
1283     VYE = e_inst_t
1284 grchrist 1.65 y_label = "Instantaneous Luminosity"
1285 abrinke1 1.5 elif varY == "live":
1286     VY = live_t
1287     VYE = e_live_t
1288 grchrist 1.65 y_label = "Instantaneous Luminosity"
1289 abrinke1 1.5 elif varY == "delivered":
1290     VY = delivered_t
1291     VYE = e_delivered_t
1292 grchrist 1.65 y_label = "Instantaneous Luminosity"
1293 abrinke1 1.5 elif varY == "deadtime":
1294     VY = deadtime_t
1295     VYE = e_deadtime_t
1296 awoodard 1.63 y_label = "Deadtime"
1297 abrinke1 1.5 elif varY == "rawrate":
1298     VY = rawrate_t
1299     VYE = e_rawrate_t
1300 awoodard 1.63 y_label = "Raw Rate [Hz]"
1301 abrinke1 1.5 if fit_file:
1302     VF = rawrate_fit_t
1303     VFE = e_rawrate_fit_t
1304     elif varY == "rate":
1305     VY = rate_t
1306     VYE = e_rate_t
1307 awoodard 1.63 y_label = "Rate [Hz]"
1308 abrinke1 1.5 if fit_file:
1309     VF = rate_fit_t
1310     VFE = e_rate_fit_t
1311     elif varY == "rawxsec":
1312     VY = rawxsec_t
1313     VYE = e_rawxsec_t
1314 grchrist 1.65 y_label = "Cross Section"
1315 abrinke1 1.5 if fit_file:
1316     VF = rawxsec_fit_t
1317     VFE = e_rawxsec_fit_t
1318     elif varY == "xsec":
1319     VY = xsec_t
1320     VYE = e_xsec_t
1321 grchrist 1.65 y_label = "Cross Section"
1322 abrinke1 1.5 if fit_file:
1323     VF = xsec_fit_t
1324     VFE = e_xsec_fit_t
1325     elif varY == "psi":
1326     VY = psi_t
1327     VYE = e_psi_t
1328 awoodard 1.63 y_label = "Prescale Index"
1329 abrinke1 1.5 else:
1330     print "No valid variable entered for Y"
1331     continue
1332    
1333 awoodard 1.63 return [VX, VXE, x_label, VY, VYE, y_label, VF, VFE]
1334 abrinke1 1.5
1335 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):
1336 grchrist 1.17 it_offset=0
1337 grchrist 1.15 Passed=True
1338     subsystemfailed=[]
1339 grchrist 1.11
1340 grchrist 1.8 if num_ls==1:
1341     ##fit is 2 ls ahead of real rate
1342 grchrist 1.10 LS=Rates[print_trigger]["ls"][iterator]
1343     #print "ls=",LS,
1344     LSRange=LumiPageInfo[LS]["LSRange"]
1345     #print LSRange,
1346     LS2=LSRange[-1]
1347     #LS2=LSRange.pop()
1348     #print "LS2=",LS2
1349     #print LumiPageInfo[LS]
1350     lumidict={}
1351     lumidict=LumiPageInfo[LS]
1352 grchrist 1.11
1353     if print_info:
1354     if (iterator==0 and print_trigger==trig_list[0]):
1355 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")
1356 grchrist 1.11
1357     ## if SubSystemOff["All"]:
1358     ## for keys in LumiPageInfo[LS]:
1359     ## #print LS, keys, LumiPageInfo[LS][keys]
1360     ## if not LumiPageInfo[LS][keys]:
1361     ## Passed=False
1362     ## subsystemfailed.append(keys)
1363     ## break
1364     ## else:
1365     if SubSystemOff["Mu"] or SubSystemOff["All"]:
1366 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"]):
1367 grchrist 1.11 Passed=False
1368     subsystemfailed.append("Mu")
1369     if SubSystemOff["HCal"] or SubSystemOff["All"]:
1370     if not (LumiPageInfo[LS]["hbhea"] and LumiPageInfo[LS]["hbheb"] and LumiPageInfo[LS]["hbhec"]):
1371     Passed=False
1372     subsystemfailed.append("HCal")
1373     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["hf"]):
1374     Passed=False
1375     subsystemfailed.append("HCal-EndCap")
1376     if SubSystemOff["ECal"] or SubSystemOff["All"]:
1377     if not (LumiPageInfo[LS]["ebp"] and LumiPageInfo[LS]["ebm"]):
1378     Passed=False
1379     subsystemfailed.append("ECal")
1380     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["eep"] and LumiPageInfo[LS]["eem"] and LumiPageInfo[LS]["esp"] or LumiPageInfo[LS]["esm"]):
1381     Passed=False
1382     subsystemfailed.append("ECal-EndCap")
1383     if SubSystemOff["Tracker"] or SubSystemOff["All"]:
1384     if not (LumiPageInfo[LS]["tob"] and LumiPageInfo[LS]["tibtid"] and LumiPageInfo[LS]["bpix"] and LumiPageInfo[LS]["fpix"]):
1385     Passed=False
1386     subsystemfailed.append("Tracker")
1387     if (SubSystemOff["EndCap"] or SubSystemOff["All"]) and not (LumiPageInfo[LS]["tecp"] and LumiPageInfo[LS]["tecm"]):
1388     Passed=False
1389     subsystemfailed.append("Tracker-EndCap")
1390     if SubSystemOff["Beam"] or SubSystemOff["All"]:
1391     if not(LumiPageInfo[LS]["b1pres"] and LumiPageInfo[LS]["b2pres"] and LumiPageInfo[LS]["b1stab"] and LumiPageInfo[LS]["b2stab"]):
1392     Passed=False
1393     subsystemfailed.append("Beam")
1394 grchrist 1.12 else:
1395     Passed=True
1396 grchrist 1.10
1397 awoodard 1.72 if not data_clean or (
1398 grchrist 1.11 Rates[print_trigger]["physics"][iterator] == 1
1399 grchrist 1.8 and Rates[print_trigger]["active"][iterator] == 1
1400 grchrist 1.17 and Rates[print_trigger]["deadtime"][iterator] < max_dt
1401 grchrist 1.11 #and Rates[print_trigger]["psi"][iterator] > 0
1402 grchrist 1.10 and Passed
1403 abrinke1 1.67 and (realvalue >0.6*prediction and realvalue<1.5*prediction)
1404     and Rates[print_trigger]["rawrate"][iterator] > 0.04
1405 grchrist 1.8 ):
1406 grchrist 1.26 if (print_info and num_ls==1 and (realvalue <0.4*prediction or realvalue>2.5*prediction)):
1407     pass
1408 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)
1409 grchrist 1.8 return True
1410     else:
1411 grchrist 1.13 if (print_info and print_trigger==trig_list[0] and num_ls==1):
1412 grchrist 1.69 prediction_match=True
1413     if (realvalue >0.6*prediction and realvalue<1.5*prediction):
1414     prediction_match=False
1415     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 )
1416 grchrist 1.26
1417 grchrist 1.8 return False
1418 abrinke1 1.5
1419 grchrist 1.10
1420     #### LumiRangeGreens ####
1421     ####inputs: RefMoreLumiArray --dict with lumi page info in LS by LS blocks,
1422     #### LRange --list range over lumis,
1423     #### nls --number of lumisections
1424     #### RefRunNum --run number
1425     ####
1426     ####outputs RangeMoreLumi --lumi page info in dict LSRange blocks with lumi, added items Run and LSRange
1427 grchrist 1.17 def LumiRangeGreens(RefMoreLumiArray,LSRange,nls,RefRunNum,deadtimebeamactive):
1428 grchrist 1.9
1429     RangeMoreLumi={}
1430     for keys,values in RefMoreLumiArray.iteritems():
1431     RangeMoreLumi[keys]=1
1432 grchrist 1.10
1433 grchrist 1.9 for iterator in LSRange[nls]:
1434     for keys, values in RefMoreLumiArray.iteritems():
1435     if RefMoreLumiArray[keys][iterator]==0:
1436     RangeMoreLumi[keys]=0
1437 grchrist 1.10 RangeMoreLumi['LSRange']=LSRange[nls]
1438     RangeMoreLumi['Run']=RefRunNum
1439 grchrist 1.17 RangeMoreLumi['DeadTimeBeamActive']=deadtimebeamactive
1440 grchrist 1.9 return RangeMoreLumi
1441    
1442 grchrist 1.10 #### CheckLumis ####
1443     ####inputs:
1444     #### PageLumiInfo --dict of LS with dict of some lumipage info
1445     #### Rates --dict of triggernames with dict of info
1446     def checkLS(Rates, PageLumiInfo,trig_list):
1447     rateslumis=Rates[trig_list[-1]]["ls"]
1448     keys=PageLumiInfo.keys()
1449     print "lumi run=",PageLumiInfo[keys[-1]]["Run"]
1450     ll=0
1451     for ls in keys:
1452     print ls,rateslumis[ll]
1453     ll=ll+1
1454     return False
1455    
1456    
1457 grchrist 1.71 def checkL1seedChangeALLPScols(trig_list,HLTL1PS):
1458     L1PSchangedic={}
1459     nps=0
1460     for HLTkey in trig_list:
1461     if HLTkey=='HLT_Stream_A':
1462     continue
1463     #print HLTkey
1464     try:
1465     dict=HLTL1PS[StripVersion(HLTkey)]
1466     #print "dict=",dict
1467     except:
1468     #print HLTkey, StripVersion(HLTkey)
1469     exit(2)
1470     HLTL1dummy={}
1471     for L1seed in dict.iterkeys():
1472     #print L1seed
1473     dummyL1seedlist=[]
1474     #print dict[L1seed]
1475     dummy=dict[L1seed]
1476     L1seedchangedummy=[]
1477     L1fulldummy=[]
1478     nps=len(dict[L1seed])
1479     #print "nps=",nps
1480     for PScol in range(0,len(dict[L1seed])):
1481     PScoldummy=PScol+1
1482     if PScoldummy>(len(dict)-1):
1483     PScoldummy=len(dict)-1
1484     #print PScol, PScoldummy, dummy[PScol]
1485    
1486     if dummy[PScol]==dummy[PScoldummy]:
1487     L1seedchangedummy.append(PScol)
1488     else:
1489     L1seedchangedummy.append(PScol)
1490     for ps in L1seedchangedummy:
1491     L1fulldummy.append(L1seedchangedummy)
1492     #print "L1seed change ", L1seedchangedummy, "full=",L1fulldummy
1493     L1seedchangedummy=[]
1494     for ps in L1seedchangedummy:
1495     L1fulldummy.append(L1seedchangedummy)
1496     #print "L1full=",L1fulldummy
1497     HLTL1dummy[L1seed]=L1fulldummy
1498     #print HLTL1dummy
1499     HLTL1_seedchanges=commonL1PS(HLTL1dummy,nps)
1500     print HLTkey, HLTL1_seedchanges
1501    
1502     def commonL1PS(HLTL1dummy, nps):
1503     ### find commmon elements in L1 seeds
1504     HLTL1_seedchanges=[]
1505     for PScol in range(0,nps):
1506    
1507     L1seedslist=HLTL1dummy.keys()
1508     L1tupletmp=set(tuple(HLTL1dummy[L1seedslist.pop()][PScol]))
1509     while len(L1seedslist)>0:
1510     L1tupletmp2=set(tuple(HLTL1dummy[L1seedslist.pop()][PScol]))
1511     L1tupletmp=L1tupletmp & L1tupletmp2
1512     if list(tuple(L1tupletmp)) not in HLTL1_seedchanges:
1513     HLTL1_seedchanges.append(list(tuple(L1tupletmp)))
1514     #print HLTL1_seedchanges
1515     return HLTL1_seedchanges
1516 grchrist 1.9
1517 grchrist 1.29
1518 abrinke1 1.1 if __name__=='__main__':
1519 grchrist 1.29 global thisyear
1520 abrinke1 1.1 main()