ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.18
Committed: Wed Mar 21 15:04:17 2012 UTC (13 years, 1 month ago) by amott
Content type: text/x-python
Branch: MAIN
CVS Tags: V00-00-15
Changes since 1.17: +178 -80 lines
Log Message:
Lots of updates for the secondary shifter

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