Classes
class	LumiInfo
	#################### ## ## LumiInfo Class ## ## #################### ## More...
class	LumiInfoCont
	######################## ## ## LumiInfoCont Class ## ## ######################## ## More...
Functions
def	loadEvents
	####################### ## ## General Functions ## ## ####################### ##
def	makeEDFplot
Variables
string	action = 'store_true'
list	allowedEDF = ['time', 'instLum', 'instIntLum']
	################ ## ## ########## ## ## ## ## Main ## ## ## ## ########## ## ## ################ ##
tuple	cont = LumiInfoCont(args[0], **options.__dict__)
	load Luminosity info ##
string	default = 'Empirical Distribution Function'
	done = False
string	help = 'title of plot (default %default)'
tuple	inputGroup = optparse.OptionGroup(parser, "Input Options")
tuple	modeGroup = optparse.OptionGroup(parser, "Mode Options")
tuple	nonSpaceRE = re.compile(r'\S')
tuple	parser = optparse.OptionParser("Usage: %prog [options] lumi.csv events.txt output.png", description='Script for generating EDF curves. See https://twiki.cern.ch/twiki/bin/viewauth/CMS/SWGuideGenerateEDF for more details.')
tuple	pieces = sepRE.split(line)
tuple	plotGroup = optparse.OptionGroup(parser, "Plot Options")
int	prevRecLumi = 0
tuple	rangeGroup = optparse.OptionGroup(parser, "Range Options")
int	recLumIndex = 0
list	recLumis = []
	look for which runs correspond to what total ## recorded integrated luminosity ##
tuple	recLumValue = float(piece)
tuple	sepRE = re.compile(r'[\s,;:]+')
string	type = 'string'

Function Documentation

def generateEDF::loadEvents	(	filename,
		cont,
		options
	)

####################### ## ## General Functions ## ## ####################### ##

Definition at line 234 of file generateEDF.py.

00235                                         :
00236     eventsDict = {}
00237     print "loading events from '%s'" % filename
00238     events = open (filename, 'r')
00239     runIndex, lumiIndex, eventIndex, weightIndex = 0, 1, 2, 3
00240     if options.relOrder:
00241         lumiIndex, eventIndex = 2,1
00242     minPieces = 3
00243     totalWeight = 0.
00244     if options.weights:
00245         minPieces = 4
00246     for line in events:
00247         pieces = sepRE.split (line.strip())
00248         if len (pieces) < minPieces:
00249             if nonSpaceRE.search (line):
00250                 print "skipping", line
00251             continue
00252         try:
00253             run, lumi, event = int( pieces[runIndex]   ), \
00254                                int( pieces[lumiIndex]  ), \
00255                                int( pieces[eventIndex] )
00256         except:
00257             continue
00258         key = (run, lumi)
00259         if not cont.has_key (key):
00260             if options.ignore:
00261                 print "Warning, %s is not found in the lumi information" \
00262                       % key.__str__()
00263                 continue
00264             else:
00265                 raise RuntimeError, "%s is not found in lumi information.  Use '--ignoreNoLumiEvents' option to ignore these events and continue." \
00266                       % key.__str__()
00267         if options.edfMode != 'time' and not cont[key].xingInfo:
00268             if options.ignore:
00269                 print "Warning, %s does not have Xing information" \
00270                       % key.__str__()
00271                 continue
00272             else:
00273                 raise RuntimeError, "%s  does not have Xing information.  Use '--ignoreNoLumiEvents' option to ignore these events and continue." \
00274                       % key.__str__()            
00275         if options.weights:
00276             weight = float (pieces[weightIndex])
00277         else:
00278             weight = 1
00279         eventsDict.setdefault( key, []).append( (event, weight) )
00280         totalWeight += weight
00281     events.close()
00282     return eventsDict, totalWeight
00283

def generateEDF::makeEDFplot	(	lumiCont,
		eventsDict,
		totalWeight,
		outputFile,
		options
	)

Definition at line 284 of file generateEDF.py.

00285                                                                         :
00286     # make TGraph
00287     xVals      = [0]
00288     yVals      = [0]
00289     expectedVals = [0]
00290     predVals   = [0]
00291     weight = 0
00292     expectedChunks = []
00293     ########################
00294     ## Time Ordering Mode ##
00295     ########################
00296     if 'time' == options.edfMode:
00297         # if we have a minimum run number, clear the lists
00298         if lumiCont.minRun or lumiCont.minIntLum:
00299             xVals      = []
00300             yVals      = []
00301             expectedVals = []
00302             predVals   = []
00303         # loop over events
00304         for key, eventList in sorted( eventsDict.iteritems() ):
00305             usePoints = True
00306             # should we add this point?
00307             if lumiCont.minRun and lumiCont.minRun > key[0] or \
00308                lumiCont.maxRun and lumiCont.maxRun < key[0]:
00309                 usePoints = False
00310             for event in eventList:
00311                 weight += event[1]
00312                 if not usePoints:
00313                     continue
00314                 factor = weight / totalWeight
00315                 try:
00316                     intLum = lumiCont[key].totalRecorded
00317                 except:
00318                     raise RuntimeError, "key %s not found in lumi information" \
00319                           % key.__str__()
00320                 if lumiCont.minIntLum and lumiCont.minIntLum > intLum or \
00321                    lumiCont.maxIntLum and lumiCont.maxIntLum < intLum:
00322                     continue
00323                 lumFrac = intLum / lumiCont.totalRecLum
00324                 xVals.append( lumiCont[key].totalRecorded)
00325                 yVals.append (factor)
00326                 expectedVals.append (lumFrac)
00327                 predVals.append   (lumFrac * options.pred)
00328         # put on the last point if we aren't giving a maximum run
00329         if not lumiCont.maxRun and not lumiCont.maxIntLum:
00330             xVals.append (lumiCont.totalRecLum)
00331             yVals.append (1)
00332             expectedVals.append (1)
00333             predVals.append (options.pred)
00334         ####################
00335         ## Reset Expected ##
00336         ####################
00337         if options.resetExpected:
00338             slope = (yVals[-1] - yVals[0]) / (xVals[-1] - xVals[0])
00339             print "slope", slope
00340             for index, old in enumerate (expectedVals):
00341                 expectedVals[index] = yVals[0] + \
00342                                     slope * (xVals[index] - xVals[0])
00343         #############################################
00344         ## Break Expected by Integrated Luminosity ##
00345         #############################################
00346         if options.breakExpectedIntLum:
00347             breakExpectedIntLum = []
00348             for chunk in options.breakExpectedIntLum:
00349                 pieces = sepRE.split (chunk)
00350                 try:
00351                     for piece in pieces:
00352                         breakExpectedIntLum.append( float(piece) )
00353                 except:
00354                     raise RuntimeError, "'%s' from '%s' is not a valid float" \
00355                           % (piece, chunk)
00356             breakExpectedIntLum.sort()
00357             boundaries = []
00358             breakIndex = 0
00359             done = False
00360             for index, xPos in enumerate (xVals):
00361                 if xPos > breakExpectedIntLum[breakIndex]:
00362                     boundaries.append (index)
00363                     while breakIndex < len (breakExpectedIntLum):
00364                         breakIndex += 1
00365                         if breakIndex >= len (breakExpectedIntLum):
00366                             done = True
00367                             break
00368                         # If this next position is different, than
00369                         # we're golden.  Otherwise, let it go through
00370                         # the loop again.
00371                         if xPos <= breakExpectedIntLum[breakIndex]:
00372                             break
00373                     if done:
00374                         break
00375             # do we have any boundaries?
00376             if not boundaries:
00377                 raise RuntimeError, "No values of 'breakExpectedIntLum' are in current range."
00378             # is the first boundary at 0?  If not, add 0
00379             if boundaries[0]:
00380                 boundaries.insert (0, 0)
00381             # is the last boundary at the end?  If not, make the end a
00382             # boundary
00383             if boundaries[-1] != len (xVals) - 1:
00384                 boundaries.append( len (xVals) - 1 )
00385             rangeList = zip (boundaries, boundaries[1:])
00386             for thisRange in rangeList:
00387                 upper = thisRange[1]
00388                 lower = thisRange[0]
00389                 slope = (yVals[upper] - yVals[lower]) / \
00390                         (xVals[upper] - xVals[lower])
00391                 print "slope", slope
00392                 # now go over the range inclusively
00393                 pairList = []
00394                 for index in range (lower, upper + 1):
00395                     newExpected = yVals[lower] + \
00396                                 slope * (xVals[index] - xVals[lower])
00397                     pairList.append( (xVals[index], newExpected) )
00398                     expectedVals[index] = newExpected
00399                 expectedChunks.append (pairList)
00400     ###########################################
00401     ## Instantanous Luminosity Ordering Mode ##
00402     ###########################################
00403     elif 'instLum' == options.edfMode or 'instIntLum' == options.edfMode:
00404         eventTupList = []
00405         if not lumiCont.xingInfo:
00406             raise RuntimeError, "Luminosity Xing information missing."
00407         for key, eventList in sorted( eventsDict.iteritems() ):
00408             try:
00409                 lumi =  lumiCont[key]
00410                 instLum   = lumi.aveInstLum
00411                 fracAXIL  = lumi.fracAXILrecorded
00412                 totalAXIL = lumi.totalAXILrecorded
00413             except:
00414                 raise RuntimeError, "key %s not found in lumi information" \
00415                       % key.__str__()
00416             for event in eventList:
00417                 eventTupList.append( (instLum, fracAXIL, totalAXIL, key,
00418                                       event[0], event[1], ) )
00419         eventTupList.sort()
00420         for eventTup in eventTupList:
00421             weight += eventTup[5]
00422             factor = weight / totalWeight
00423             if 'instLum' == options.edfMode:
00424                 xVals.append (eventTup[0])
00425             else:
00426                 xVals.append (eventTup[2])
00427             yVals.append (factor)
00428             expectedVals.append (eventTup[1])
00429             predVals.append   (eventTup[1] * options.pred)
00430     else:
00431         raise RuntimeError, "It looks like Charles screwed up if you are seeing this."
00432 
00433     size = len (xVals)
00434     step = int (math.sqrt(size) / 2 + 1)
00435     if options.printValues:
00436         for index in range (size):
00437             print "%8f %8f %8f" % (xVals[index], yVals[index], expectedVals[index]),
00438             if index > step:
00439                 denom = xVals[index] - xVals[index - step]
00440                 numer = yVals[index] - yVals[index - step]
00441                 if denom:
00442                     print " %8f" % (numer / denom),
00443             if 0 == index % step:
00444                 print " **", ## indicates statistically independent
00445                              ## slope measurement
00446             print
00447         print
00448 
00449     xArray = array.array ('d', xVals)
00450     yArray = array.array ('d', yVals)
00451     expected = array.array ('d', expectedVals)
00452     graph = ROOT.TGraph( size, xArray, yArray)
00453     graph.SetTitle (options.title)
00454     graph.SetMarkerStyle (20)
00455     expectedGraph = ROOT.TGraph( size, xArray, expected)
00456     expectedGraph.SetLineColor (ROOT.kRed)
00457     expectedGraph.SetLineWidth (3)
00458     if options.noDataPoints:
00459         expectedGraph.SetLineStyle (2)
00460 
00461     # run statistical tests
00462     if options.weights:
00463         print "average weight per event:", weight / ( size - 1)
00464     maxDistance = ROOT.TMath.KolmogorovTest (size, yArray,
00465                                              size, expected,
00466                                              "M")
00467     prob = ROOT.TMath.KolmogorovProb( maxDistance * math.sqrt( size ) )
00468 
00469     # display everything
00470     ROOT.gROOT.SetStyle('Plain')
00471     ROOT.gROOT.SetBatch()
00472     c1 = ROOT.TCanvas()
00473     graph.GetXaxis().SetRangeUser (min (xVals), max (xVals))
00474     minValue = min (min(yVals), min(expected))
00475     if options.pred:
00476         minValue = min (minValue, min (predVals))
00477     graph.GetYaxis().SetRangeUser (minValue,
00478                                    max (max(yVals), max(expected), max(predVals)))
00479     graph.SetLineWidth (3)
00480     if options.noDataPoints:
00481         graph.Draw ("AL")
00482     else:
00483         graph.Draw ("ALP")
00484     if 'instLum' == options.edfMode:
00485         graph.GetXaxis().SetTitle ("Average Xing Inst. Luminosity (1/ub/s)")
00486         graph.GetXaxis().SetRangeUser (0., lumiCont.max('aveInstLum'))
00487     else:
00488         if 'instIntLum' == options.edfMode:
00489             graph.GetXaxis().SetTitle ("Integrated Luminosity - Inst. Lum. Ordered (1/%s)" \
00490                                        % lumiCont.invunits)
00491         else:
00492             graph.GetXaxis().SetTitle ("Integrated Luminosity (1/%s)" \
00493                                        % lumiCont.invunits)
00494     graph.GetYaxis().SetTitle ("Fraction of Events Seen")
00495     expectedGraphs = []
00496     if expectedChunks:
00497         for index, chunk in enumerate (expectedChunks):
00498             expectedXarray = array.array ('d', [item[0] for item in chunk])
00499             expectedYarray = array.array ('d', [item[1] for item in chunk])
00500             expectedGraph = ROOT.TGraph( len(chunk),
00501                                          expectedXarray,
00502                                          expectedYarray )
00503             expectedGraph.SetLineWidth (3)
00504             if options.noDataPoints:
00505                 expectedGraph.SetLineStyle (2)
00506             if index % 2:
00507                 expectedGraph.SetLineColor (ROOT.kBlue)
00508             else:
00509                 expectedGraph.SetLineColor (ROOT.kRed)
00510             expectedGraph.Draw("L")
00511             expectedGraphs.append (expectedGraph)
00512         exptectedGraph = expectedGraphs[0]
00513     else:
00514         expectedGraph.Draw ("L")
00515     green = 0
00516     if options.pred:
00517         predArray = array.array ('d', predVals)
00518         green = ROOT.TGraph (size, xArray, predArray)
00519         green.SetLineWidth (3)
00520         green.SetLineColor (8)
00521         green.Draw ('l')
00522     legend = ROOT.TLegend(0.15, 0.65, 0.50, 0.85)
00523     legend.SetFillStyle (0)
00524     legend.SetLineColor(ROOT.kWhite)
00525     observed = 'Observed'
00526     if options.weights:
00527         observed += ' (weighted)'
00528     legend.AddEntry(graph, observed,"PL")
00529     if options.resetExpected:
00530         legend.AddEntry(expectedGraph,  "Expected from partial yield","L")
00531     else:
00532         legend.AddEntry(expectedGraph,  "Expected from total yield","L")
00533     if options.pred:
00534         legend.AddEntry(green, options.predLabel,"L")
00535     legend.AddEntry("","D_{stat}=%.3f, N=%d" % (maxDistance, size),"")
00536     legend.AddEntry("","P_{KS}=%.3f" % prob,"")
00537     legend.Draw()
00538 
00539     # save file
00540     c1.Print (outputFile)
00541

Variable Documentation

string generateEDF::action = 'store_true'

Definition at line 570 of file generateEDF.py.

list generateEDF::allowedEDF = ['time', 'instLum', 'instIntLum']

################ ## ## ########## ## ## ## ## Main ## ## ## ## ########## ## ## ################ ##

command line options ##

Definition at line 554 of file generateEDF.py.

tuple generateEDF::cont = LumiInfoCont(args[0], **options.__dict__)

load Luminosity info ##

Definition at line 622 of file generateEDF.py.

string generateEDF::default = 'Empirical Distribution Function'

Definition at line 561 of file generateEDF.py.

generateEDF::done = False

Definition at line 651 of file generateEDF.py.

string generateEDF::help = 'title of plot (default %default)'

Definition at line 562 of file generateEDF.py.

tuple generateEDF::inputGroup = optparse.OptionGroup(parser, "Input Options")

Definition at line 558 of file generateEDF.py.

tuple generateEDF::modeGroup = optparse.OptionGroup(parser, "Mode Options")

Definition at line 559 of file generateEDF.py.

tuple generateEDF::nonSpaceRE = re.compile(r'\S')

Definition at line 12 of file generateEDF.py.

tuple generateEDF::parser = optparse.OptionParser("Usage: %prog [options] lumi.csv events.txt output.png", description='Script for generating EDF curves. See https://twiki.cern.ch/twiki/bin/viewauth/CMS/SWGuideGenerateEDF for more details.')

Definition at line 555 of file generateEDF.py.

tuple generateEDF::pieces = sepRE.split(line)

Definition at line 635 of file generateEDF.py.

tuple generateEDF::plotGroup = optparse.OptionGroup(parser, "Plot Options")

Definition at line 556 of file generateEDF.py.

generateEDF::prevRecLumi = 0

Definition at line 650 of file generateEDF.py.

tuple generateEDF::rangeGroup = optparse.OptionGroup(parser, "Range Options")

Definition at line 557 of file generateEDF.py.

int generateEDF::recLumIndex = 0

Definition at line 648 of file generateEDF.py.

list generateEDF::recLumis = []