generateEDF Namespace Reference

Classes
class	LumiInfo
	LumiInfo Class More...
class	LumiInfoCont
	LumiInfoCont Class More...

Functions
def	loadEvents (filename, cont, options)
	General Functions More...
def	makeEDFplot (lumiCont, eventsDict, totalWeight, outputFile, options)

Function Documentation

def generateEDF.loadEvents	(		filename,
			cont,
			options
	)

General Functions

Definition at line 236 of file generateEDF.py.

References mps_setup.append, createfilelist.int, and edm.print().

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 key not in cont:
            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

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

Definition at line 286 of file generateEDF.py.

References objects.autophobj.float, list(), SiStripPI.max, min(), and edm.print().

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( six.iteritems(eventsDict) ):
            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 = list(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( six.iteritems(eventsDict) ):
            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]), end=' ')
            if index > step:
                denom = xVals[index] - xVals[index - step]
                numer = yVals[index] - yVals[index - step]
                if denom:
                    print(" %8f" % (numer / denom), end=' ')
            if 0 == index % step:
                print(" **", end=' ') ## 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)