ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.20
Committed: Wed Mar 21 16:11:02 2012 UTC (13 years, 1 month ago) by amott
Content type: text/x-python
Branch: MAIN
CVS Tags: V00-00-20, V00-00-19, V00-00-18, V00-00-17
Changes since 1.19: +35 -25 lines
Log Message:
Fixed lots of bugs

File Contents

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