ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.17
Committed: Tue Mar 20 12:15:43 2012 UTC (13 years, 1 month ago) by grchrist
Content type: text/x-python
Branch: MAIN
Changes since 1.16: +69 -80 lines
Log Message:
added beam veto on fitted ls and deadtime beam active now parsed directly from db

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