Reference Manual

Classes
class	LumiInfo
	LumiInfo Class More...

class	LumiInfoCont
	LumiInfoCont Class More...

Functions
def	loadEvents
	General Functions More...

def	makeEDFplot

Variables
string	action = 'store_true'

list	allowedEDF = ['time', 'instLum', 'instIntLum']
	## Main More...

tuple	cont = LumiInfoCont(args[0], **options.__dict__)
	load Luminosity info ## More...

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 ## More...

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.

References python.multivaluedict.append().

 
 def loadEvents (filename, cont, options):
     eventsDict = {}
     print "loading events from '%s'" % filename
     events = open (filename, 'r')
     runIndex, lumiIndex, eventIndex, weightIndex = 0, 1, 2, 3
     if options.relOrder:
         lumiIndex, eventIndex = 2,1
     minPieces = 3
     totalWeight = 0.
     if options.weights:
         minPieces = 4
     for line in events:
         pieces = sepRE.split (line.strip())
         if len (pieces) < minPieces:
             if nonSpaceRE.search (line):
                 print "skipping", line
             continue
         try:
             run, lumi, event = int( pieces[runIndex]   ), \
                                int( pieces[lumiIndex]  ), \
                                int( pieces[eventIndex] )
         except:
             continue
         key = (run, lumi)
         if not cont.has_key (key):
             if options.ignore:
                 print "Warning, %s is not found in the lumi information" \
                       % key.__str__()
                 continue
             else:
                 raise RuntimeError, "%s is not found in lumi information.  Use '--ignoreNoLumiEvents' option to ignore these events and continue." \
                       % key.__str__()
         if options.edfMode != 'time' and not cont[key].xingInfo:
             if options.ignore:
                 print "Warning, %s does not have Xing information" \
                       % key.__str__()
                 continue
             else:
                 raise RuntimeError, "%s  does not have Xing information.  Use '--ignoreNoLumiEvents' option to ignore these events and continue." \
                       % key.__str__()            
         if options.weights:
             weight = float (pieces[weightIndex])
         else:
             weight = 1
         eventsDict.setdefault( key, []).append( (event, weight) )
         totalWeight += weight
     events.close()
     return eventsDict, totalWeight
 

python.multivaluedict.append

def append

Definition: multivaluedict.py:73

generateEDF.loadEvents

def loadEvents

General Functions

Definition: generateEDF.py:234

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

Definition at line 284 of file generateEDF.py.

References bookConverter.max, and min().

 
 def makeEDFplot (lumiCont, eventsDict, totalWeight, outputFile, options):
     # make TGraph
     xVals      = [0]
     yVals      = [0]
     expectedVals = [0]
     predVals   = [0]
     weight = 0
     expectedChunks = []
     ########################
     ## Time Ordering Mode ##
     ########################
     if 'time' == options.edfMode:
         # if we have a minimum run number, clear the lists
         if lumiCont.minRun or lumiCont.minIntLum:
             xVals      = []
             yVals      = []
             expectedVals = []
             predVals   = []
         # loop over events
         for key, eventList in sorted( eventsDict.iteritems() ):
             usePoints = True
             # should we add this point?
             if lumiCont.minRun and lumiCont.minRun > key[0] or \
                lumiCont.maxRun and lumiCont.maxRun < key[0]:
                 usePoints = False
             for event in eventList:
                 weight += event[1]
                 if not usePoints:
                     continue
                 factor = weight / totalWeight
                 try:
                     intLum = lumiCont[key].totalRecorded
                 except:
                     raise RuntimeError, "key %s not found in lumi information" \
                           % key.__str__()
                 if lumiCont.minIntLum and lumiCont.minIntLum > intLum or \
                    lumiCont.maxIntLum and lumiCont.maxIntLum < intLum:
                     continue
                 lumFrac = intLum / lumiCont.totalRecLum
                 xVals.append( lumiCont[key].totalRecorded)
                 yVals.append (factor)
                 expectedVals.append (lumFrac)
                 predVals.append   (lumFrac * options.pred)
         # put on the last point if we aren't giving a maximum run
         if not lumiCont.maxRun and not lumiCont.maxIntLum:
             xVals.append (lumiCont.totalRecLum)
             yVals.append (1)
             expectedVals.append (1)
             predVals.append (options.pred)
         ####################
         ## Reset Expected ##
         ####################
         if options.resetExpected:
             slope = (yVals[-1] - yVals[0]) / (xVals[-1] - xVals[0])
             print "slope", slope
             for index, old in enumerate (expectedVals):
                 expectedVals[index] = yVals[0] + \
                                     slope * (xVals[index] - xVals[0])
         #############################################
         ## Break Expected by Integrated Luminosity ##
         #############################################
         if options.breakExpectedIntLum:
             breakExpectedIntLum = []
             for chunk in options.breakExpectedIntLum:
                 pieces = sepRE.split (chunk)
                 try:
                     for piece in pieces:
                         breakExpectedIntLum.append( float(piece) )
                 except:
                     raise RuntimeError, "'%s' from '%s' is not a valid float" \
                           % (piece, chunk)
             breakExpectedIntLum.sort()
             boundaries = []
             breakIndex = 0
             done = False
             for index, xPos in enumerate (xVals):
                 if xPos > breakExpectedIntLum[breakIndex]:
                     boundaries.append (index)
                     while breakIndex < len (breakExpectedIntLum):
                         breakIndex += 1
                         if breakIndex >= len (breakExpectedIntLum):
                             done = True
                             break
                         # If this next position is different, than
                         # we're golden.  Otherwise, let it go through
                         # the loop again.
                         if xPos <= breakExpectedIntLum[breakIndex]:
                             break
                     if done:
                         break
             # do we have any boundaries?
             if not boundaries:
                 raise RuntimeError, "No values of 'breakExpectedIntLum' are in current range."
             # is the first boundary at 0?  If not, add 0
             if boundaries[0]:
                 boundaries.insert (0, 0)
             # is the last boundary at the end?  If not, make the end a
             # boundary
             if boundaries[-1] != len (xVals) - 1:
                 boundaries.append( len (xVals) - 1 )
             rangeList = zip (boundaries, boundaries[1:])
             for thisRange in rangeList:
                 upper = thisRange[1]
                 lower = thisRange[0]
                 slope = (yVals[upper] - yVals[lower]) / \
                         (xVals[upper] - xVals[lower])
                 print "slope", slope
                 # now go over the range inclusively
                 pairList = []
                 for index in range (lower, upper + 1):
                     newExpected = yVals[lower] + \
                                 slope * (xVals[index] - xVals[lower])
                     pairList.append( (xVals[index], newExpected) )
                     expectedVals[index] = newExpected
                 expectedChunks.append (pairList)
     ###########################################
     ## Instantanous Luminosity Ordering Mode ##
     ###########################################
     elif 'instLum' == options.edfMode or 'instIntLum' == options.edfMode:
         eventTupList = []
         if not lumiCont.xingInfo:
             raise RuntimeError, "Luminosity Xing information missing."
         for key, eventList in sorted( eventsDict.iteritems() ):
             try:
                 lumi =  lumiCont[key]
                 instLum   = lumi.aveInstLum
                 fracAXIL  = lumi.fracAXILrecorded
                 totalAXIL = lumi.totalAXILrecorded
             except:
                 raise RuntimeError, "key %s not found in lumi information" \
                       % key.__str__()
             for event in eventList:
                 eventTupList.append( (instLum, fracAXIL, totalAXIL, key,
                                       event[0], event[1], ) )
         eventTupList.sort()
         for eventTup in eventTupList:
             weight += eventTup[5]
             factor = weight / totalWeight
             if 'instLum' == options.edfMode:
                 xVals.append (eventTup[0])
             else:
                 xVals.append (eventTup[2])
             yVals.append (factor)
             expectedVals.append (eventTup[1])
             predVals.append   (eventTup[1] * options.pred)
     else:
         raise RuntimeError, "It looks like Charles screwed up if you are seeing this."
 
     size = len (xVals)
     step = int (math.sqrt(size) / 2 + 1)
     if options.printValues:
         for index in range (size):
             print "%8f %8f %8f" % (xVals[index], yVals[index], expectedVals[index]),
             if index > step:
                 denom = xVals[index] - xVals[index - step]
                 numer = yVals[index] - yVals[index - step]
                 if denom:
                     print " %8f" % (numer / denom),
             if 0 == index % step:
                 print " **", ## indicates statistically independent
                              ## slope measurement
             print
         print
 
     xArray = array.array ('d', xVals)
     yArray = array.array ('d', yVals)
     expected = array.array ('d', expectedVals)
     graph = ROOT.TGraph( size, xArray, yArray)
     graph.SetTitle (options.title)
     graph.SetMarkerStyle (20)
     expectedGraph = ROOT.TGraph( size, xArray, expected)
     expectedGraph.SetLineColor (ROOT.kRed)
     expectedGraph.SetLineWidth (3)
     if options.noDataPoints:
         expectedGraph.SetLineStyle (2)
 
     # run statistical tests
     if options.weights:
         print "average weight per event:", weight / ( size - 1)
     maxDistance = ROOT.TMath.KolmogorovTest (size, yArray,
                                              size, expected,
                                              "M")
     prob = ROOT.TMath.KolmogorovProb( maxDistance * math.sqrt( size ) )
 
     # display everything
     ROOT.gROOT.SetStyle('Plain')
     ROOT.gROOT.SetBatch()
     c1 = ROOT.TCanvas()
     graph.GetXaxis().SetRangeUser (min (xVals), max (xVals))
     minValue = min (min(yVals), min(expected))
     if options.pred:
         minValue = min (minValue, min (predVals))
     graph.GetYaxis().SetRangeUser (minValue,
                                    max (max(yVals), max(expected), max(predVals)))
     graph.SetLineWidth (3)
     if options.noDataPoints:
         graph.Draw ("AL")
     else:
         graph.Draw ("ALP")
     if 'instLum' == options.edfMode:
         graph.GetXaxis().SetTitle ("Average Xing Inst. Luminosity (1/ub/s)")
         graph.GetXaxis().SetRangeUser (0., lumiCont.max('aveInstLum'))
     else:
         if 'instIntLum' == options.edfMode:
             graph.GetXaxis().SetTitle ("Integrated Luminosity - Inst. Lum. Ordered (1/%s)" \
                                        % lumiCont.invunits)
         else:
             graph.GetXaxis().SetTitle ("Integrated Luminosity (1/%s)" \
                                        % lumiCont.invunits)
     graph.GetYaxis().SetTitle ("Fraction of Events Seen")
     expectedGraphs = []
     if expectedChunks:
         for index, chunk in enumerate (expectedChunks):
             expectedXarray = array.array ('d', [item[0] for item in chunk])
             expectedYarray = array.array ('d', [item[1] for item in chunk])
             expectedGraph = ROOT.TGraph( len(chunk),
                                          expectedXarray,
                                          expectedYarray )
             expectedGraph.SetLineWidth (3)
             if options.noDataPoints:
                 expectedGraph.SetLineStyle (2)
             if index % 2:
                 expectedGraph.SetLineColor (ROOT.kBlue)
             else:
                 expectedGraph.SetLineColor (ROOT.kRed)
             expectedGraph.Draw("L")
             expectedGraphs.append (expectedGraph)
         exptectedGraph = expectedGraphs[0]
     else:
         expectedGraph.Draw ("L")
     green = 0
     if options.pred:
         predArray = array.array ('d', predVals)
         green = ROOT.TGraph (size, xArray, predArray)
         green.SetLineWidth (3)
         green.SetLineColor (8)
         green.Draw ('l')
     legend = ROOT.TLegend(0.15, 0.65, 0.50, 0.85)
     legend.SetFillStyle (0)
     legend.SetLineColor(ROOT.kWhite)
     observed = 'Observed'
     if options.weights:
         observed += ' (weighted)'
     legend.AddEntry(graph, observed,"PL")
     if options.resetExpected:
         legend.AddEntry(expectedGraph,  "Expected from partial yield","L")
     else:
         legend.AddEntry(expectedGraph,  "Expected from total yield","L")
     if options.pred:
         legend.AddEntry(green, options.predLabel,"L")
     legend.AddEntry("","D_{stat}=%.3f, N=%d" % (maxDistance, size),"")
     legend.AddEntry("","P_{KS}=%.3f" % prob,"")
     legend.Draw()
 
     # save file
     c1.Print (outputFile)
 