ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/RateMonShiftTool_dev/DatabaseRatePredictor.py
Revision: 1.14
Committed: Fri Mar 16 10:11:38 2012 UTC (13 years, 1 month ago) by grchrist
Content type: text/x-python
Branch: MAIN
Changes since 1.13: +22 -14 lines
Log Message:
added versions and output for fitting

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