| 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 |
|
| 10 |
from ROOT import gROOT, TCanvas, TF1, TGraph, TGraphErrors, TPaveStats, gPad, gStyle
|
| 11 |
from ROOT import TFile, TPaveText
|
| 12 |
from ROOT import gBenchmark
|
| 13 |
import array
|
| 14 |
import math
|
| 15 |
|
| 16 |
from selectionParser import selectionParser
|
| 17 |
|
| 18 |
def main():
|
| 19 |
|
| 20 |
## Can use any combination of LowestRunNumber, HighestRunNumber, and NumberOfRuns -
|
| 21 |
## just modify "ExistingRuns.sort" and for "run in ExistingRuns" accordingly
|
| 22 |
|
| 23 |
LowestRunNumber = 176000
|
| 24 |
HighestRunNumber = 180252
|
| 25 |
NumberOfRuns = 1
|
| 26 |
|
| 27 |
run_list = []
|
| 28 |
ExistingRuns = GetLatestRunNumber(LowestRunNumber)
|
| 29 |
#ExistingRuns.sort(reverse = True) ##Allows you to count down from last run
|
| 30 |
|
| 31 |
for run in ExistingRuns:
|
| 32 |
#if NumberOfRuns > 0:
|
| 33 |
if run <= HighestRunNumber:
|
| 34 |
run_list.append(run)
|
| 35 |
NumberOfRuns-=1
|
| 36 |
|
| 37 |
## ###### TO CREATE FITS #########
|
| 38 |
## #run_list = [179497,179547,179558,179563,179889,179959,179977,180072,180076,180093,180241,180250,180252]
|
| 39 |
## #run_list = [180241]
|
| 40 |
## trig_name = "IsoMu"
|
| 41 |
## num_ls = 10
|
| 42 |
## physics_active_psi = True ##Requires that physics and active be on, and that the prescale column is not 0
|
| 43 |
## #JSON = [] ##To not use a JSON file, just leave the array empty
|
| 44 |
## JSON = GetJSON("Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt") ##Returns array JSON[runs][ls_list]
|
| 45 |
|
| 46 |
## debug_print = False
|
| 47 |
|
| 48 |
## min_rate = 0.1
|
| 49 |
## print_table = False
|
| 50 |
## data_clean = True ##Gets rid of anomalous rate points, reqires physics_active_psi (PAP) and deadtime < 20%
|
| 51 |
## ##plot_properties = [varX, varY, do_fit, save_root, save_png, fit_file]
|
| 52 |
## plot_properties = [["delivered", "rate", True, False, False, ""]]
|
| 53 |
|
| 54 |
## masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
|
| 55 |
## save_fits = True
|
| 56 |
|
| 57 |
|
| 58 |
###### TO SEE RATE VS PREDICTION ########
|
| 59 |
run_list = [180252]
|
| 60 |
|
| 61 |
trig_name = "Mu"
|
| 62 |
num_ls = 1
|
| 63 |
physics_active_psi = True
|
| 64 |
JSON = []
|
| 65 |
debug_print = False
|
| 66 |
|
| 67 |
min_rate = 1.0
|
| 68 |
print_table = False
|
| 69 |
data_clean = False
|
| 70 |
##plot_properties = [varX, varY, do_fit, save_root, save_png, fit_file]
|
| 71 |
plot_properties = [["ls", "rawrate", False, True, False, "Fits/2011/Fit_HLT_10LS_Run176023to180252.pkl"]]
|
| 72 |
masked_triggers = ["AlCa_", "DST_", "HLT_L1", "HLT_L2", "HLT_Zero"]
|
| 73 |
save_fits = False
|
| 74 |
|
| 75 |
|
| 76 |
######## END PARAMETERS - CALL FUNCTIONS ##########
|
| 77 |
Rates = GetDBRates(run_list, trig_name, num_ls, physics_active_psi, JSON, debug_print)
|
| 78 |
MakePlots(Rates, run_list, trig_name, num_ls, min_rate, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print)
|
| 79 |
|
| 80 |
|
| 81 |
def GetDBRates(run_list,trig_name,num_ls,physics_active_psi,JSON,debug_print):
|
| 82 |
|
| 83 |
Rates = {}
|
| 84 |
## Save in RefRuns with name dependent on trig_name, num_ls, JSON, and physics_active_psi
|
| 85 |
if JSON:
|
| 86 |
if physics_active_psi:
|
| 87 |
RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_JPAP.pkl"
|
| 88 |
else:
|
| 89 |
RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_JSON.pkl"
|
| 90 |
else:
|
| 91 |
if physics_active_psi:
|
| 92 |
RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS_PAP.pkl"
|
| 93 |
else:
|
| 94 |
RefRunNameTemplate = "RefRuns/2011/Rates_%s_%sLS.pkl"
|
| 95 |
|
| 96 |
RefRunFile = RefRunNameTemplate % (trig_name,num_ls)
|
| 97 |
RefRunFileHLT = RefRunNameTemplate % ("HLT",num_ls)
|
| 98 |
try: ##Open an existing RefRun file with the same parameters and trigger name
|
| 99 |
pkl_file = open(RefRunFile, 'rb')
|
| 100 |
Rates = pickle.load(pkl_file)
|
| 101 |
pkl_file.close()
|
| 102 |
os.remove(RefRunFile)
|
| 103 |
except:
|
| 104 |
try: ##Open an existing RefRun file with the same parameters and HLT for trigger name
|
| 105 |
pkl_file = open(RefRunFileHLT)
|
| 106 |
HLTRates = pickle.load(pkl_file)
|
| 107 |
for key in HLTRates:
|
| 108 |
if trig_name in str(key):
|
| 109 |
Rates[key] = HLTRates[key]
|
| 110 |
print str(RefRunFile)+" does not exist. Creating ..."
|
| 111 |
except:
|
| 112 |
print str(RefRunFile)+" does not exist. Creating ..."
|
| 113 |
|
| 114 |
for RefRunNum in run_list:
|
| 115 |
|
| 116 |
if JSON:
|
| 117 |
if not RefRunNum in JSON:
|
| 118 |
continue
|
| 119 |
|
| 120 |
try:
|
| 121 |
ExistsAlready = False
|
| 122 |
for key in Rates:
|
| 123 |
if RefRunNum in Rates[key]["run"]:
|
| 124 |
ExistsAlready = True
|
| 125 |
break
|
| 126 |
if ExistsAlready:
|
| 127 |
continue
|
| 128 |
except:
|
| 129 |
print "Getting info for run "+str(RefRunNum)
|
| 130 |
|
| 131 |
if RefRunNum < 1:
|
| 132 |
continue
|
| 133 |
|
| 134 |
if True: ##Placeholder
|
| 135 |
if True: #May replace with "try" - for now it's good to know when problems happen
|
| 136 |
RefParser = DatabaseParser()
|
| 137 |
RefParser.RunNumber = RefRunNum
|
| 138 |
RefParser.ParseRunSetup()
|
| 139 |
RefLumiRangePhysicsActive = RefParser.GetLSRange(1,9999) ##Gets array of all LS with physics and active on
|
| 140 |
RefLumiArray = RefParser.GetLumiInfo() ##Gets array of all existing LS and their lumi info
|
| 141 |
RefLumiRange = []
|
| 142 |
|
| 143 |
for iterator in RefLumiArray[0]: ##Makes array of LS with proper PAP and JSON properties
|
| 144 |
if not physics_active_psi or (RefLumiArray[5][iterator] == 1 and RefLumiArray[6][iterator] == 1 and RefLumiArray[0][iterator] > 0):
|
| 145 |
if not JSON or RefRunNum in JSON:
|
| 146 |
if not JSON or iterator in JSON[RefRunNum]:
|
| 147 |
RefLumiRange.append(iterator)
|
| 148 |
|
| 149 |
try:
|
| 150 |
nls = RefLumiRange[0]
|
| 151 |
LSRange = {}
|
| 152 |
except:
|
| 153 |
print "Run "+str(RefRunNum)+" has no good LS"
|
| 154 |
continue
|
| 155 |
if num_ls > len(RefLumiRange):
|
| 156 |
print "Run "+str(RefRunNum)+" is too short: from "+str(nls)+" to "+str(RefLumiRange[-1])+", while num_ls = "+str(num_ls)
|
| 157 |
continue
|
| 158 |
while nls < RefLumiRange[-1]-num_ls:
|
| 159 |
LSRange[nls] = []
|
| 160 |
counter = 0
|
| 161 |
for iterator in RefLumiRange:
|
| 162 |
if iterator >= nls and counter < num_ls:
|
| 163 |
LSRange[nls].append(iterator)
|
| 164 |
counter += 1
|
| 165 |
nls = LSRange[nls][-1]+1
|
| 166 |
|
| 167 |
print "Run "+str(RefRunNum)+" contains LS from "+str(min(LSRange))+" to "+str(max(LSRange))
|
| 168 |
for nls in LSRange:
|
| 169 |
TriggerRates = RefParser.GetHLTRates(LSRange[nls])
|
| 170 |
|
| 171 |
[inst, live, delivered, dead, pscols] = RefParser.GetAvLumiInfo(LSRange[nls])
|
| 172 |
|
| 173 |
physics = 1
|
| 174 |
active = 1
|
| 175 |
psi = 99
|
| 176 |
for iterator in LSRange[nls]: ##Gets lowest value of physics, active, and psi in the set of lumisections
|
| 177 |
if RefLumiArray[5][iterator] == 0:
|
| 178 |
physics = 0
|
| 179 |
if RefLumiArray[6][iterator] == 0:
|
| 180 |
active = 0
|
| 181 |
if RefLumiArray[0][iterator] < psi:
|
| 182 |
psi = RefLumiArray[0][iterator]
|
| 183 |
|
| 184 |
if inst < 0 or live < 0 or delivered < 0:
|
| 185 |
print "Run "+str(RefRunNum)+" LS "+str(nls)+" inst lumi = "+str(inst)+" live lumi = "+str(live)+", delivered = "+str(delivered)+", physics = "+str(physics)+", active = "+str(active)
|
| 186 |
|
| 187 |
for key in TriggerRates:
|
| 188 |
if not trig_name in key:
|
| 189 |
continue
|
| 190 |
name = key
|
| 191 |
if re.match('.*_v[0-9]+',name): ##Removes _v#
|
| 192 |
name = name[:name.rfind('_')]
|
| 193 |
|
| 194 |
if not Rates.has_key(name):
|
| 195 |
Rates[name] = {}
|
| 196 |
Rates[name]["run"] = []
|
| 197 |
Rates[name]["ls"] = []
|
| 198 |
Rates[name]["ps"] = []
|
| 199 |
Rates[name]["inst_lumi"] = []
|
| 200 |
Rates[name]["live_lumi"] = []
|
| 201 |
Rates[name]["delivered_lumi"] = []
|
| 202 |
Rates[name]["deadtime"] = []
|
| 203 |
Rates[name]["rawrate"] = []
|
| 204 |
Rates[name]["rate"] = []
|
| 205 |
Rates[name]["rawxsec"] = []
|
| 206 |
Rates[name]["xsec"] = []
|
| 207 |
Rates[name]["physics"] = []
|
| 208 |
Rates[name]["active"] = []
|
| 209 |
Rates[name]["psi"] = []
|
| 210 |
[avps, ps, rate, psrate] = TriggerRates[key]
|
| 211 |
Rates[name]["run"].append(RefRunNum)
|
| 212 |
Rates[name]["ls"].append(nls)
|
| 213 |
Rates[name]["ps"].append(ps)
|
| 214 |
Rates[name]["inst_lumi"].append(inst)
|
| 215 |
Rates[name]["live_lumi"].append(live)
|
| 216 |
Rates[name]["delivered_lumi"].append(delivered)
|
| 217 |
Rates[name]["deadtime"].append(dead)
|
| 218 |
Rates[name]["rawrate"].append(rate)
|
| 219 |
if live == 0:
|
| 220 |
Rates[name]["rate"].append(0)
|
| 221 |
Rates[name]["rawxsec"].append(0.0)
|
| 222 |
Rates[name]["xsec"].append(0.0)
|
| 223 |
else:
|
| 224 |
Rates[name]["rate"].append(psrate/(1.0-dead))
|
| 225 |
Rates[name]["rawxsec"].append(rate/live)
|
| 226 |
Rates[name]["xsec"].append(psrate/live)
|
| 227 |
Rates[name]["physics"].append(physics)
|
| 228 |
Rates[name]["active"].append(active)
|
| 229 |
Rates[name]["psi"].append(psi)
|
| 230 |
#except: ##If we replace "if True:" with "try:"
|
| 231 |
#print "Failed to parse run "+str(RefRunNum)
|
| 232 |
|
| 233 |
RateOutput = open(RefRunFile, 'wb') ##Save new Rates[] to RefRuns
|
| 234 |
pickle.dump(Rates, RateOutput, 2)
|
| 235 |
RateOutput.close()
|
| 236 |
return Rates
|
| 237 |
|
| 238 |
def MakePlots(Rates, run_list, trig_name, num_ls, min_rate, print_table, data_clean, plot_properties, masked_triggers, save_fits, debug_print):
|
| 239 |
min_run = min(run_list)
|
| 240 |
max_run = max(run_list)
|
| 241 |
|
| 242 |
InputFit = {}
|
| 243 |
OutputFit = {}
|
| 244 |
|
| 245 |
RootNameTemplate = "%s_%sLS_Run%sto%s.root"
|
| 246 |
RootFile = RootNameTemplate % (trig_name, num_ls, min_run, max_run)
|
| 247 |
|
| 248 |
for varX, varY, do_fit, save_root, save_png, fit_file in plot_properties:
|
| 249 |
if fit_file:
|
| 250 |
pkl_file = open(fit_file, 'rb')
|
| 251 |
InputFit = pickle.load(pkl_file)
|
| 252 |
pkl_file.close()
|
| 253 |
if save_root:
|
| 254 |
try:
|
| 255 |
os.remove(RootFile)
|
| 256 |
except:
|
| 257 |
break
|
| 258 |
|
| 259 |
for print_trigger in Rates:
|
| 260 |
##Limits Rates[] to runs in run_list
|
| 261 |
NewTrigger = {}
|
| 262 |
for key in Rates[print_trigger]:
|
| 263 |
NewTrigger[key] = []
|
| 264 |
for iterator in range (len(Rates[print_trigger]["run"])):
|
| 265 |
if Rates[print_trigger]["run"][iterator] in run_list:
|
| 266 |
for key in Rates[print_trigger]:
|
| 267 |
NewTrigger[key].append(Rates[print_trigger][key][iterator])
|
| 268 |
Rates[print_trigger] = NewTrigger
|
| 269 |
|
| 270 |
meanrawrate = sum(Rates[print_trigger]["rawrate"])/len(Rates[print_trigger]["rawrate"])
|
| 271 |
if not trig_name in print_trigger:
|
| 272 |
continue
|
| 273 |
if meanrawrate < min_rate:
|
| 274 |
continue
|
| 275 |
masked_trig = False
|
| 276 |
for mask in masked_triggers:
|
| 277 |
if str(mask) in print_trigger:
|
| 278 |
masked_trig = True
|
| 279 |
if masked_trig:
|
| 280 |
continue
|
| 281 |
|
| 282 |
OutputFit[print_trigger] = {}
|
| 283 |
|
| 284 |
lowlumi = 0
|
| 285 |
meanlumi_init = median(Rates[print_trigger]["live_lumi"])
|
| 286 |
meanlumi = 0
|
| 287 |
highlumi = 0
|
| 288 |
lowxsec = 0
|
| 289 |
meanxsec = 0
|
| 290 |
highxsec = 0
|
| 291 |
nlow = 0
|
| 292 |
nhigh = 0
|
| 293 |
for iterator in range(len(Rates[print_trigger]["rate"])):
|
| 294 |
if Rates[print_trigger]["live_lumi"][iterator] <= meanlumi_init:
|
| 295 |
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] < 0.20 and Rates[print_trigger]["psi"][iterator] > 0 and Rates[print_trigger]["live_lumi"] > 500):
|
| 296 |
meanxsec+=Rates[print_trigger]["xsec"][iterator]
|
| 297 |
lowxsec+=Rates[print_trigger]["xsec"][iterator]
|
| 298 |
meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
|
| 299 |
lowlumi+=Rates[print_trigger]["live_lumi"][iterator]
|
| 300 |
nlow+=1
|
| 301 |
if Rates[print_trigger]["live_lumi"][iterator] > meanlumi_init:
|
| 302 |
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] < 0.20 and Rates[print_trigger]["psi"][iterator] > 0 and Rates[print_trigger]["live_lumi"] > 500):
|
| 303 |
meanxsec+=Rates[print_trigger]["xsec"][iterator]
|
| 304 |
highxsec+=Rates[print_trigger]["xsec"][iterator]
|
| 305 |
meanlumi+=Rates[print_trigger]["live_lumi"][iterator]
|
| 306 |
highlumi+=Rates[print_trigger]["live_lumi"][iterator]
|
| 307 |
nhigh+=1
|
| 308 |
try:
|
| 309 |
meanxsec = meanxsec/(nlow+nhigh)
|
| 310 |
meanlumi = meanlumi/(nlow+nhigh)
|
| 311 |
slopexsec = ( (highxsec/nhigh) - (lowxsec/nlow) ) / ( (highlumi/nhigh) - (lowlumi/nlow) )
|
| 312 |
except:
|
| 313 |
print str(print_trigger)+" has no good datapoints - setting initial xsec slope estimate to 0"
|
| 314 |
meanxsec = median(Rates[print_trigger]["xsec"])
|
| 315 |
meanlumi = median(Rates[print_trigger]["live_lumi"])
|
| 316 |
slopexsec = 0
|
| 317 |
|
| 318 |
[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()
|
| 319 |
|
| 320 |
if fit_file:
|
| 321 |
FitType = InputFit[print_trigger][0]
|
| 322 |
X0 = InputFit[print_trigger][1]
|
| 323 |
X1 = InputFit[print_trigger][2]
|
| 324 |
X2 = InputFit[print_trigger][3]
|
| 325 |
X3 = InputFit[print_trigger][4]
|
| 326 |
Chi2 = InputFit[print_trigger][5]
|
| 327 |
##print str(print_trigger)+" "+str(FitType)+" "+str(X0)+" "+str(X1)+" "+str(X2)+" "+str(X3)
|
| 328 |
|
| 329 |
for iterator in range(len(Rates[print_trigger]["rate"])):
|
| 330 |
if not Rates[print_trigger]["run"][iterator] in run_list:
|
| 331 |
continue
|
| 332 |
prediction = meanxsec + slopexsec * (Rates[print_trigger]["live_lumi"][iterator] - meanlumi)
|
| 333 |
realvalue = Rates[print_trigger]["xsec"][iterator]
|
| 334 |
if not data_clean or ( ((realvalue > 0.4*prediction and realvalue < 2.5*prediction) or (realvalue > 0.4*meanxsec and realvalue < 2.5*meanxsec) or prediction < 0 ) and Rates[print_trigger]["physics"][iterator] == 1 and Rates[print_trigger]["active"][iterator] == 1 and Rates[print_trigger]["deadtime"][iterator] < 0.20 and Rates[print_trigger]["psi"][iterator] > 0):
|
| 335 |
run_t.append(Rates[print_trigger]["run"][iterator])
|
| 336 |
ls_t.append(Rates[print_trigger]["ls"][iterator])
|
| 337 |
ps_t.append(Rates[print_trigger]["ps"][iterator])
|
| 338 |
inst_t.append(Rates[print_trigger]["inst_lumi"][iterator])
|
| 339 |
live_t.append(Rates[print_trigger]["live_lumi"][iterator])
|
| 340 |
delivered_t.append(Rates[print_trigger]["delivered_lumi"][iterator])
|
| 341 |
deadtime_t.append(Rates[print_trigger]["deadtime"][iterator])
|
| 342 |
rawrate_t.append(Rates[print_trigger]["rawrate"][iterator])
|
| 343 |
rate_t.append(Rates[print_trigger]["rate"][iterator])
|
| 344 |
rawxsec_t.append(Rates[print_trigger]["rawxsec"][iterator])
|
| 345 |
xsec_t.append(Rates[print_trigger]["xsec"][iterator])
|
| 346 |
psi_t.append(Rates[print_trigger]["psi"][iterator])
|
| 347 |
|
| 348 |
e_run_t.append(0.0)
|
| 349 |
e_ls_t.append(0.0)
|
| 350 |
e_ps_t.append(0.0)
|
| 351 |
e_inst_t.append(14.14)
|
| 352 |
e_live_t.append(14.14)
|
| 353 |
e_delivered_t.append(14.14)
|
| 354 |
e_deadtime_t.append(0.01)
|
| 355 |
e_rawrate_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
|
| 356 |
e_rate_t.append(Rates[print_trigger]["ps"][iterator]*math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3)))
|
| 357 |
e_psi_t.append(0.0)
|
| 358 |
if live_t[-1] == 0:
|
| 359 |
e_rawxsec_t.append(0)
|
| 360 |
e_xsec_t.append(0)
|
| 361 |
else:
|
| 362 |
e_rawxsec_t.append(math.sqrt(Rates[print_trigger]["rawrate"][iterator]/(num_ls*23.3))/Rates[print_trigger]["live_lumi"][iterator])
|
| 363 |
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])
|
| 364 |
|
| 365 |
if fit_file:
|
| 366 |
if FitType == "expo":
|
| 367 |
rate_prediction = X0 + X1*math.exp(X2*delivered_t[-1])
|
| 368 |
else:
|
| 369 |
rate_prediction = X0 + X1*delivered_t[-1] + X2*delivered_t[-1]*delivered_t[-1] + X3*delivered_t[-1]*delivered_t[-1]*delivered_t[-1]
|
| 370 |
## if rate_t[-1] < 0.7 * rate_prediction or rate_t[-1] > 1.4 * rate_prediction:
|
| 371 |
## 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])
|
| 372 |
|
| 373 |
if live_t[-1] == 0:
|
| 374 |
rawrate_fit_t.append(0)
|
| 375 |
rate_fit_t.append(0)
|
| 376 |
rawxsec_fit_t.append(0)
|
| 377 |
xsec_fit_t.append(0)
|
| 378 |
e_rawrate_fit_t.append(0)
|
| 379 |
e_rate_fit_t.append(math.sqrt(Chi2))
|
| 380 |
e_rawxsec_fit_t.append(0)
|
| 381 |
e_xsec_fit_t.append(0)
|
| 382 |
else:
|
| 383 |
rawrate_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/(ps_t[-1]))
|
| 384 |
rate_fit_t.append(rate_prediction)
|
| 385 |
rawxsec_fit_t.append(rawrate_fit_t[-1]/live_t[-1])
|
| 386 |
xsec_fit_t.append(rate_prediction*(1.0-deadtime_t[-1])/live_t[-1])
|
| 387 |
e_rawrate_fit_t.append(math.sqrt(Chi2)*rawrate_fit_t[-1]/rate_fit_t[-1])
|
| 388 |
e_rate_fit_t.append(math.sqrt(Chi2))
|
| 389 |
e_rawxsec_fit_t.append(math.sqrt(Chi2)*rawxsec_fit_t[-1]/rate_fit_t[-1])
|
| 390 |
e_xsec_fit_t.append(math.sqrt(Chi2)*xsec_fit_t[-1]/rate_fit_t[-1])
|
| 391 |
|
| 392 |
else: ##If the data point does not pass the data_clean filter
|
| 393 |
if debug_print:
|
| 394 |
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)
|
| 395 |
|
| 396 |
## End "for iterator in range(len(Rates[print_trigger]["rate"])):" loop
|
| 397 |
|
| 398 |
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]
|
| 399 |
[VX, VXE, VY, VYE, VF, VFE] = GetVXVY(plot_properties, fit_file, AllPlotArrays)
|
| 400 |
|
| 401 |
if save_root or save_png:
|
| 402 |
c1 = TCanvas(str(varX),str(varY))
|
| 403 |
c1.SetName(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
|
| 404 |
|
| 405 |
gr1 = TGraphErrors(len(VX), VX, VY, VXE, VYE)
|
| 406 |
gr1.SetName("Graph_"+str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX))
|
| 407 |
gr1.GetXaxis().SetTitle(varX)
|
| 408 |
gr1.GetYaxis().SetTitle(varY)
|
| 409 |
gr1.SetTitle(str(print_trigger))
|
| 410 |
gr1.SetMinimum(0)
|
| 411 |
gr1.SetMaximum(1.2*max(VY))
|
| 412 |
#gr1.GetXaxis().SetLimits(min(VX)-0.2*max(VX),1.2*max(VX))
|
| 413 |
gr1.GetXaxis().SetLimits(0,1.2*max(VX))
|
| 414 |
gr1.SetMarkerStyle(8)
|
| 415 |
if fit_file:
|
| 416 |
gr1.SetMarkerSize(0.8)
|
| 417 |
else:
|
| 418 |
gr1.SetMarkerSize(0.5)
|
| 419 |
gr1.SetMarkerColor(2)
|
| 420 |
|
| 421 |
if fit_file:
|
| 422 |
gr3 = TGraphErrors(len(VX), VX, VF, VXE, VFE)
|
| 423 |
gr3.SetMarkerStyle(8)
|
| 424 |
gr3.SetMarkerSize(0.4)
|
| 425 |
gr3.SetMarkerColor(4)
|
| 426 |
gr3.SetFillColor(4)
|
| 427 |
gr3.SetFillStyle(3003)
|
| 428 |
|
| 429 |
if do_fit:
|
| 430 |
if "rate" in varY:
|
| 431 |
f1a = TF1("f1a","pol2",0,8000)
|
| 432 |
f1a.SetLineColor(4)
|
| 433 |
f1a.SetLineWidth(2)
|
| 434 |
f1a.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
|
| 435 |
f1a.SetParLimits(1,0,2.0*max(VY)/(max(VX)*max(VX)))
|
| 436 |
#gr1.Fit("f1a","B","Q")
|
| 437 |
gr1.Fit("f1a","Q","rob=0.90")
|
| 438 |
|
| 439 |
f1b = 0
|
| 440 |
f1c = 0
|
| 441 |
if True:
|
| 442 |
f1b = TF1("f1b","pol3",0,8000)
|
| 443 |
f1b.SetLineColor(2)
|
| 444 |
f1b.SetLineWidth(2)
|
| 445 |
f1b.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
|
| 446 |
f1b.SetParLimits(1,0,f1a.GetParameter(1)+0.0000001)
|
| 447 |
f1b.SetParLimits(2,0,f1a.GetParameter(2)+0.0000000001)
|
| 448 |
f1b.SetParLimits(3,0,2.0*max(VY)/(max(VX)*max(VX)*max(VX)))
|
| 449 |
gr1.Fit("f1b","Q","rob=0.90")
|
| 450 |
#if f1b.GetChisquare()/f1b.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
|
| 451 |
print "X0 = "+str(f1a.GetParameter(0))+" X1 = "+str(f1a.GetParameter(1))+" X2 = "+str(f1a.GetParameter(2))
|
| 452 |
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()))
|
| 453 |
print "X0 = "+str(f1b.GetParameter(0))+" X1 = "+str(f1b.GetParameter(1))+" X2 = "+str(f1b.GetParameter(2))+" X3 = "+str(f1b.GetParameter(3))
|
| 454 |
|
| 455 |
f1c = TF1("f1c","[0]+[1]*expo(2)",0,8000)
|
| 456 |
f1c.SetLineColor(3)
|
| 457 |
f1c.SetLineWidth(2)
|
| 458 |
f1c.SetParLimits(0,0,0.2*(sum(VY)/len(VY))+0.8*min(VY))
|
| 459 |
f1c.SetParLimits(1,max(VY)/math.exp(10.0),max(VY)/math.exp(2.0))
|
| 460 |
f1c.SetParLimits(2,0.0,0.0000000001)
|
| 461 |
f1c.SetParLimits(3,2.0/max(VX),10.0/max(VX))
|
| 462 |
print str(max(VY)/math.exp(2.0))+" "+str(10.0/max(VX))
|
| 463 |
gr1.Fit("f1c","Q","rob=0.90")
|
| 464 |
#if f1c.GetChisquare()/f1c.GetNDF() < f1a.GetChisquare()/f1a.GetNDF():
|
| 465 |
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()))
|
| 466 |
print "X0 = "+str(f1c.GetParameter(0))+" X1 = "+str(f1c.GetParameter(1))+" X2 = "+str(f1c.GetParameter(2))+" X3 = "+str(f1c.GetParameter(3))
|
| 467 |
|
| 468 |
else: ##If this is not a rate plot
|
| 469 |
f1a = TF1("f1a","pol1",0,8000)
|
| 470 |
f1a.SetLineColor(4)
|
| 471 |
f1a.SetLineWidth(2)
|
| 472 |
if "xsec" in varY:
|
| 473 |
f1a.SetParLimits(0,0,meanxsec*1.5)
|
| 474 |
if slopexsec > 0:
|
| 475 |
f1a.SetParLimits(1,0,max(VY)/max(VX))
|
| 476 |
else:
|
| 477 |
f1a.SetParLimits(1,2*slopexsec,-2*slopexsec)
|
| 478 |
else:
|
| 479 |
f1a.SetParLimits(0,-1000,1000)
|
| 480 |
gr1.Fit("f1a","Q","rob=0.80")
|
| 481 |
|
| 482 |
if save_root or save_png:
|
| 483 |
gr1.Draw("APZ")
|
| 484 |
## ##Option to draw stats box
|
| 485 |
## p1 = TPaveStats()
|
| 486 |
## p1 = gr1.GetListOfFunctions().FindObject("stats")
|
| 487 |
## print p1
|
| 488 |
## gr1.PaintStats(f1b).Draw("same")
|
| 489 |
if fit_file:
|
| 490 |
gr3.Draw("P3")
|
| 491 |
if do_fit:
|
| 492 |
f1a.Draw("same")
|
| 493 |
try:
|
| 494 |
f1b.Draw("same")
|
| 495 |
f1c.Draw("same")
|
| 496 |
except:
|
| 497 |
True
|
| 498 |
c1.Update()
|
| 499 |
if save_root:
|
| 500 |
myfile = TFile( RootFile, 'UPDATE' )
|
| 501 |
c1.Write()
|
| 502 |
myfile.Close()
|
| 503 |
if save_png:
|
| 504 |
c1.SaveAs(str(print_trigger)+"_"+str(varY)+"_vs_"+str(varX)+".png")
|
| 505 |
|
| 506 |
|
| 507 |
if print_table or save_fits:
|
| 508 |
if not do_fit:
|
| 509 |
print "Can't have save_fits = True and do_fit = False"
|
| 510 |
continue
|
| 511 |
if f1c.GetChisquare()/f1c.GetNDF() < 0.95*f1a.GetChisquare()/f1a.GetNDF() and f1c.GetChisquare()/f1c.GetNDF() < 0.95*f1b.GetChisquare()/f1b.GetNDF():
|
| 512 |
OutputFit[print_trigger] = ["expo", f1c.GetParameter(0), f1c.GetParameter(1), f1c.GetParameter(3), 0.0, f1c.GetChisquare()/f1c.GetNDF(), meanrawrate]
|
| 513 |
elif f1b.GetChisquare()/f1b.GetNDF() < 0.95*f1a.GetChisquare()/f1a.GetNDF():
|
| 514 |
OutputFit[print_trigger] = ["poly", f1b.GetParameter(0), f1b.GetParameter(1), f1b.GetParameter(2), f1b.GetParameter(3), f1b.GetChisquare()/f1b.GetNDF(), meanrawrate]
|
| 515 |
else:
|
| 516 |
OutputFit[print_trigger] = ["poly", f1a.GetParameter(0), f1a.GetParameter(1), f1a.GetParameter(2), 0.0, f1a.GetChisquare()/f1a.GetNDF(), meanrawrate]
|
| 517 |
|
| 518 |
if save_root:
|
| 519 |
print "Output root file is "+str(RootFile)
|
| 520 |
|
| 521 |
if save_fits:
|
| 522 |
FitNameTemplate = "Fits/2011/Fit_%s_%sLS_Run%sto%s.pkl"
|
| 523 |
FitFile = FitNameTemplate % (trig_name, num_ls, min_run, max_run)
|
| 524 |
if os.path.exists(FitFile):
|
| 525 |
os.remove(FitFile)
|
| 526 |
FitOutputFile = open(FitFile, 'wb')
|
| 527 |
pickle.dump(OutputFit, FitOutputFile, 2)
|
| 528 |
FitOutputFile.close()
|
| 529 |
print "Output fit file is "+str(FitFile)
|
| 530 |
|
| 531 |
if print_table:
|
| 532 |
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."
|
| 533 |
print '%60s%10s%10s%10s%10s%10s%10s%10s' % ("Trig", "fit", "p0", "p1", "p2", "p3", "Chi2", "Av raw")
|
| 534 |
for print_trigger in OutputFit:
|
| 535 |
_trigger = (print_trigger[:56] + '...') if len(print_trigger) > 59 else print_trigger
|
| 536 |
try:
|
| 537 |
if OutputFit[print_trigger][0] == "poly":
|
| 538 |
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))
|
| 539 |
else:
|
| 540 |
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))
|
| 541 |
except:
|
| 542 |
print str(print_trigger)+" is somehow broken"
|
| 543 |
|
| 544 |
|
| 545 |
############# SUPPORTING FUNCTIONS ################
|
| 546 |
|
| 547 |
|
| 548 |
def GetJSON(json_file):
|
| 549 |
|
| 550 |
input_file = open(json_file)
|
| 551 |
file_content = input_file.read()
|
| 552 |
inputRange = selectionParser(file_content)
|
| 553 |
JSON = inputRange.runsandls()
|
| 554 |
return JSON
|
| 555 |
##JSON is an array: JSON[run_number] = [1st ls, 2nd ls, 3rd ls ... nth ls]
|
| 556 |
|
| 557 |
def MakePlotArrays():
|
| 558 |
run_t = array.array('f')
|
| 559 |
ls_t = array.array('f')
|
| 560 |
ps_t = array.array('f')
|
| 561 |
inst_t = array.array('f')
|
| 562 |
live_t = array.array('f')
|
| 563 |
delivered_t = array.array('f')
|
| 564 |
deadtime_t = array.array('f')
|
| 565 |
rawrate_t = array.array('f')
|
| 566 |
rate_t = array.array('f')
|
| 567 |
rawxsec_t = array.array('f')
|
| 568 |
xsec_t = array.array('f')
|
| 569 |
psi_t = array.array('f')
|
| 570 |
|
| 571 |
e_run_t = array.array('f')
|
| 572 |
e_ls_t = array.array('f')
|
| 573 |
e_ps_t = array.array('f')
|
| 574 |
e_inst_t = array.array('f')
|
| 575 |
e_live_t = array.array('f')
|
| 576 |
e_delivered_t = array.array('f')
|
| 577 |
e_deadtime_t = array.array('f')
|
| 578 |
e_rawrate_t = array.array('f')
|
| 579 |
e_rate_t = array.array('f')
|
| 580 |
e_rawxsec_t = array.array('f')
|
| 581 |
e_xsec_t = array.array('f')
|
| 582 |
e_psi_t = array.array('f')
|
| 583 |
|
| 584 |
rawrate_fit_t = array.array('f')
|
| 585 |
rate_fit_t = array.array('f')
|
| 586 |
rawxsec_fit_t = array.array('f')
|
| 587 |
xsec_fit_t = array.array('f')
|
| 588 |
e_rawrate_fit_t = array.array('f')
|
| 589 |
e_rate_fit_t = array.array('f')
|
| 590 |
e_rawxsec_fit_t = array.array('f')
|
| 591 |
e_xsec_fit_t = array.array('f')
|
| 592 |
|
| 593 |
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]
|
| 594 |
|
| 595 |
|
| 596 |
def GetVXVY(plot_properties, fit_file, AllPlotArrays):
|
| 597 |
|
| 598 |
VF = "0"
|
| 599 |
VFE = "0"
|
| 600 |
|
| 601 |
[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
|
| 602 |
for varX, varY, do_fit, save_root, save_png, fit_file in plot_properties:
|
| 603 |
if varX == "run":
|
| 604 |
VX = run_t
|
| 605 |
VXE = run_t_e
|
| 606 |
elif varX == "ls":
|
| 607 |
VX = ls_t
|
| 608 |
VXE = e_ls_t
|
| 609 |
elif varX == "ps":
|
| 610 |
VX = ps_t
|
| 611 |
VXE = e_ps_t
|
| 612 |
elif varX == "inst":
|
| 613 |
VX = inst_t
|
| 614 |
VXE = e_inst_t
|
| 615 |
elif varX == "live":
|
| 616 |
VX = live_t
|
| 617 |
VXE = e_live_t
|
| 618 |
elif varX == "delivered":
|
| 619 |
VX = delivered_t
|
| 620 |
VXE = e_delivered_t
|
| 621 |
elif varX == "deadtime":
|
| 622 |
VX = deadtime_t
|
| 623 |
VXE = e_deadtime_t
|
| 624 |
elif varX == "rawrate":
|
| 625 |
VX = rawrate_t
|
| 626 |
VXE = e_rawrate_t
|
| 627 |
elif varX == "rate":
|
| 628 |
VX = rate_t
|
| 629 |
VXE = e_rate_t
|
| 630 |
elif varX == "rawxsec":
|
| 631 |
VX = rawxsec_t
|
| 632 |
VXE = e_rawxsec_t
|
| 633 |
elif varX == "xsec":
|
| 634 |
VX = xsec_t
|
| 635 |
VXE = e_xsec_t
|
| 636 |
elif varX == "psi":
|
| 637 |
VX = psi_t
|
| 638 |
VXE = e_psi_t
|
| 639 |
else:
|
| 640 |
print "No valid variable entered for X"
|
| 641 |
continue
|
| 642 |
if varY == "run":
|
| 643 |
VY = run_t
|
| 644 |
VYE = run_t_e
|
| 645 |
elif varY == "ls":
|
| 646 |
VY = ls_t
|
| 647 |
VYE = e_ls_t
|
| 648 |
elif varY == "ps":
|
| 649 |
VY = ps_t
|
| 650 |
VYE = e_ps_t
|
| 651 |
elif varY == "inst":
|
| 652 |
VY = inst_t
|
| 653 |
VYE = e_inst_t
|
| 654 |
elif varY == "live":
|
| 655 |
VY = live_t
|
| 656 |
VYE = e_live_t
|
| 657 |
elif varY == "delivered":
|
| 658 |
VY = delivered_t
|
| 659 |
VYE = e_delivered_t
|
| 660 |
elif varY == "deadtime":
|
| 661 |
VY = deadtime_t
|
| 662 |
VYE = e_deadtime_t
|
| 663 |
elif varY == "rawrate":
|
| 664 |
VY = rawrate_t
|
| 665 |
VYE = e_rawrate_t
|
| 666 |
if fit_file:
|
| 667 |
VF = rawrate_fit_t
|
| 668 |
VFE = e_rawrate_fit_t
|
| 669 |
elif varY == "rate":
|
| 670 |
VY = rate_t
|
| 671 |
VYE = e_rate_t
|
| 672 |
if fit_file:
|
| 673 |
VF = rate_fit_t
|
| 674 |
VFE = e_rate_fit_t
|
| 675 |
elif varY == "rawxsec":
|
| 676 |
VY = rawxsec_t
|
| 677 |
VYE = e_rawxsec_t
|
| 678 |
if fit_file:
|
| 679 |
VF = rawxsec_fit_t
|
| 680 |
VFE = e_rawxsec_fit_t
|
| 681 |
elif varY == "xsec":
|
| 682 |
VY = xsec_t
|
| 683 |
VYE = e_xsec_t
|
| 684 |
if fit_file:
|
| 685 |
VF = xsec_fit_t
|
| 686 |
VFE = e_xsec_fit_t
|
| 687 |
elif varY == "psi":
|
| 688 |
VY = psi_t
|
| 689 |
VYE = e_psi_t
|
| 690 |
else:
|
| 691 |
print "No valid variable entered for Y"
|
| 692 |
continue
|
| 693 |
|
| 694 |
return [VX, VXE, VY, VYE, VF, VFE]
|
| 695 |
|
| 696 |
|
| 697 |
if __name__=='__main__':
|
| 698 |
main()
|
