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