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