ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.36
Committed: Mon Apr 9 10:53:10 2012 UTC (13 years ago) by grchrist
Content type: text/x-python
Branch: MAIN
CVS Tags: V00-00-26
Changes since 1.35: +1 -1 lines
Log Message:
latest and greated fit file

File Contents

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