PythonUtilities/scripts/generateEDF.py

0001 #! /usr/bin/env python3
0002
0003 from __future__ import print_function
0004 from __future__ import division
0005 from builtins import zip
0006 from builtins import object
0007 from past.utils import old_div
0008 from builtins import range
0009 import sys
0010 import re
0011 import optparse
0012 from pprint import pprint
0013 import array
0014 import ROOT
0015 import math
0016
0017 sepRE      = re.compile (r'[\s,;:]+')
0018 nonSpaceRE = re.compile (r'\S')
0019
0020 ##########################
0021 ## #################### ##
0022 ## ## LumiInfo Class ## ##
0023 ## #################### ##
0024 ##########################
0025
0026 class LumiInfo (object):
0027
0028     lastSingleXingRun = 136175
0029     lumiSectionLength = 23.310779
0030
0031     def __init__ (self, line):
0032         self.totInstLum  = 0.
0033         self.aveInstLum  = 0.
0034         self.numXings    = 0
0035         self.instLums    = []
0036         self.events      = []
0037         self.xingInfo    = False
0038         self.badXingInfo = False
0039         pieces = sepRE.split (line.strip())
0040         size = len (pieces)
0041         if size % 2:
0042             raise RuntimeError("Odd number of pieces")
0043         if size < 4:
0044             raise RuntimeError("Not enough pieces")
0045         try:
0046             self.run       = int   (pieces[0])
0047             self.lumi      = int   (pieces[1])
0048             self.delivered = float (pieces[2])
0049             self.recorded  = float (pieces[3])
0050         except:
0051             raise RuntimeError("Pieces not right format")
0052         if size > 4:
0053             try:
0054                 for xing, lum in zip (pieces[4::2],pieces[5::2]):
0055                     xing = int   (xing)
0056                     lum  = float (lum)
0057                     self.instLums.append( (xing, lum) )
0058                     self.totInstLum += lum
0059                     self.numXings += 1
0060             except:
0061                 raise RuntimeError("Inst Lumi Info malformed")
0062             self.aveInstLum = old_div(self.totInstLum, (self.numXings))
0063             self.xingInfo   = True
0064         self.key       = (self.run, self.lumi)
0065         self.keyString = self.key.__str__()
0066
0067
0068     def fixXingInfo (self):
0069         if self.numXings:
0070             # You shouldn't try and fix an event if it already has
0071             # xing information.
0072             raise RuntimeError("This event %s already has Xing information" \
0073                   % self.keyString)
0074         if self.run > LumiInfo.lastSingleXingRun:
0075             # this run may have more than one crossing.  I don't know
0076             # how to fix this.
0077             self.badXingInfo = True
0078             return False
0079         self.numXings = 1
0080         xing = 1
0081         self.aveInstLum = self.totInstLum = lum = \
0082                           old_div(self.delivered, LumiInfo.lumiSectionLength)
0083         self.instLums.append( (xing, lum) )
0084         self.xingInfo = True
0085         return True
0086
0087
0088     def deadtime (self):
0089         if not self.delivered:
0090             return 1
0091         return 1 - (old_div(self.recorded, self.delivered))
0092
0093
0094     def __str__ (self):
0095         return "%6d, %4d: %6.1f (%4.1f%%) %4.2f (%3d)" % \
0096                (self.run,
0097                 self.lumi,
0098                 self.delivered,
0099                 self.deadtime(),
0100                 self.totInstLum,
0101                 self.numXings)
0102
0103
0104 ##############################
0105 ## ######################## ##
0106 ## ## LumiInfoCont Class ## ##
0107 ## ######################## ##
0108 ##############################
0109
0110 class LumiInfoCont (dict):
0111
0112     def __init__ (self, filename, **kwargs):
0113         print("loading luminosity information from '%s'." % filename)
0114         source = open (filename, 'r')
0115         self.minMaxKeys = ['totInstLum', 'aveInstLum', 'numXings',
0116                      'delivered', 'recorded']
0117         self._min = {}
0118         self._max = {}
0119         self.totalRecLum = 0.
0120         self.xingInfo    = False
0121         self.allowNoXing = kwargs.get ('ignore')
0122         self.noWarnings  = kwargs.get ('noWarnings')
0123         self.minRun      = 0
0124         self.maxRun      = 0
0125         self.minIntLum   = 0
0126         self.maxIntLum   = 0
0127
0128         for key in self.minMaxKeys:
0129             self._min[key] = -1
0130             self._max[key] =  0
0131         for line in source:
0132             try:
0133                 lumi = LumiInfo (line)
0134             except:
0135                 continue
0136             self[lumi.key] = lumi
0137             self.totalRecLum += lumi.recorded
0138             if not self.xingInfo and lumi.xingInfo:
0139                 self.xingInfo = True
0140             if lumi.xingInfo:
0141                 #print "yes", lumi.keyString
0142                 if not self.xingInfo:
0143                     print("huh?")
0144             for key in self.minMaxKeys:
0145                 val = getattr (lumi, key)
0146                 if val < self._min[key] or self._min[key] < 0:
0147                     self._min[key] = val
0148                 if val > self._max[key] or not self._max[key]:
0149                     self._max[key] = val
0150         source.close()
0151         ######################################################
0152         ## Now that everything is setup, switch integrated  ##
0153         ## luminosity to more reasonable units.             ##
0154         ######################################################
0155         # the default is '1/mb', but that's just silly.
0156         self.invunits = 'nb'
0157         lumFactor = 1e3
0158         if   self.totalRecLum > 1e9:
0159             lumFactor = 1e9
0160             self.invunits = 'fb'
0161         elif self.totalRecLum > 1e6:
0162             lumFactor = 1e6
0163             self.invunits = 'pb'
0164         # use lumFactor to make everything consistent
0165         #print "units", self.invunits, "factor", lumFactor
0166         self.totalRecLum /= lumFactor
0167         for lumis in self.values():
0168             lumis.delivered /= lumFactor
0169             lumis.recorded  /= lumFactor
0170         # Probably want to rename this next subroutine, but I'll leave
0171         # it alone for now...
0172         self._integrateContainer()
0173
0174
0175
0176     def __str__ (self):
0177         retval = 'run,     lum     del ( dt  ) inst (#xng)\n'
0178         for key, value in sorted (self.items()):
0179             retval += "%s\n" % value
0180         return retval
0181
0182
0183     def min (self, key):
0184         return self._min[key]
0185
0186
0187     def max (self, key):
0188         return self._max[key]
0189
0190
0191     def keys (self):
0192         return sorted (dict.keys (self))
0193
0194
0195     def iteritems (self):
0196         return sorted (dict.iteritems (self))
0197
0198
0199     def _integrateContainer (self):
0200         # calculate numbers for recorded integrated luminosity
0201         total = 0.
0202         for key, lumi in self.items():
0203             total += lumi.recorded
0204             lumi.totalRecorded = total
0205             lumi.fracRecorded  = old_div(total, self.totalRecLum)
0206         # calculate numbers for average xing instantaneous luminosity
0207         if not self.xingInfo:
0208             # nothing to do here
0209             return
0210         xingKeyList = []
0211         maxAveInstLum = 0.
0212         for key, lumi in self.items():
0213             if not lumi.xingInfo and not lumi.fixXingInfo():
0214                 if not self.noWarnings:
0215                     print("Do not have lumi xing info for %s" % lumi.keyString)
0216                 if not self.allowNoXing:
0217                     print("Setting no Xing info flag")
0218                     self.xingInfo = False
0219                     return
0220                 continue
0221             xingKeyList.append( (lumi.aveInstLum, key) )
0222             if lumi.aveInstLum > maxAveInstLum:
0223                 maxAveInstLum = lumi.aveInstLum
0224         xingKeyList.sort()
0225         total = 0.
0226         for tup in xingKeyList:
0227             lumi = self[tup[1]]
0228             total += lumi.recorded
0229             lumi.totalAXILrecorded = total
0230             lumi.fracAXILrecorded  = old_div(total, self.totalRecLum)
0231             lumi.fracAXIL = old_div(lumi.aveInstLum, maxAveInstLum)
0232
0233
0234 #############################
0235 ## ####################### ##
0236 ## ## General Functions ## ##
0237 ## ####################### ##
0238 #############################
0239
0240 def loadEvents (filename, cont, options):
0241     eventsDict = {}
0242     print("loading events from '%s'" % filename)
0243     events = open (filename, 'r')
0244     runIndex, lumiIndex, eventIndex, weightIndex = 0, 1, 2, 3
0245     if options.relOrder:
0246         lumiIndex, eventIndex = 2,1
0247     minPieces = 3
0248     totalWeight = 0.
0249     if options.weights:
0250         minPieces = 4
0251     for line in events:
0252         pieces = sepRE.split (line.strip())
0253         if len (pieces) < minPieces:
0254             if nonSpaceRE.search (line):
0255                 print("skipping", line)
0256             continue
0257         try:
0258             run, lumi, event = int( pieces[runIndex]   ), \
0259                                int( pieces[lumiIndex]  ), \
0260                                int( pieces[eventIndex] )
0261         except:
0262             continue
0263         key = (run, lumi)
0264         if key not in cont:
0265             if options.ignore:
0266                 print("Warning, %s is not found in the lumi information" \
0267                       % key.__str__())
0268                 continue
0269             else:
0270                 raise RuntimeError("%s is not found in lumi information.  Use '--ignoreNoLumiEvents' option to ignore these events and continue." \
0271                       % key.__str__())
0272         if options.edfMode != 'time' and not cont[key].xingInfo:
0273             if options.ignore:
0274                 print("Warning, %s does not have Xing information" \
0275                       % key.__str__())
0276                 continue
0277             else:
0278                 raise RuntimeError("%s  does not have Xing information.  Use '--ignoreNoLumiEvents' option to ignore these events and continue." \
0279                       % key.__str__())
0280         if options.weights:
0281             weight = float (pieces[weightIndex])
0282         else:
0283             weight = 1
0284         eventsDict.setdefault( key, []).append( (event, weight) )
0285         totalWeight += weight
0286     events.close()
0287     return eventsDict, totalWeight
0288
0289
0290 def makeEDFplot (lumiCont, eventsDict, totalWeight, outputFile, options):
0291     # make TGraph
0292     xVals      = [0]
0293     yVals      = [0]
0294     expectedVals = [0]
0295     predVals   = [0]
0296     weight = 0
0297     expectedChunks = []
0298     ########################
0299     ## Time Ordering Mode ##
0300     ########################
0301     if 'time' == options.edfMode:
0302         # if we have a minimum run number, clear the lists
0303         if lumiCont.minRun or lumiCont.minIntLum:
0304             xVals      = []
0305             yVals      = []
0306             expectedVals = []
0307             predVals   = []
0308         # loop over events
0309         for key, eventList in sorted( eventsDict.items() ):
0310             usePoints = True
0311             # should we add this point?
0312             if lumiCont.minRun and lumiCont.minRun > key[0] or \
0313                lumiCont.maxRun and lumiCont.maxRun < key[0]:
0314                 usePoints = False
0315             for event in eventList:
0316                 weight += event[1]
0317                 if not usePoints:
0318                     continue
0319                 factor = old_div(weight, totalWeight)
0320                 try:
0321                     intLum = lumiCont[key].totalRecorded
0322                 except:
0323                     raise RuntimeError("key %s not found in lumi information" \
0324                           % key.__str__())
0325                 if lumiCont.minIntLum and lumiCont.minIntLum > intLum or \
0326                    lumiCont.maxIntLum and lumiCont.maxIntLum < intLum:
0327                     continue
0328                 lumFrac = old_div(intLum, lumiCont.totalRecLum)
0329                 xVals.append( lumiCont[key].totalRecorded)
0330                 yVals.append (factor)
0331                 expectedVals.append (lumFrac)
0332                 predVals.append   (lumFrac * options.pred)
0333         # put on the last point if we aren't giving a maximum run
0334         if not lumiCont.maxRun and not lumiCont.maxIntLum:
0335             xVals.append (lumiCont.totalRecLum)
0336             yVals.append (1)
0337             expectedVals.append (1)
0338             predVals.append (options.pred)
0339         ####################
0340         ## Reset Expected ##
0341         ####################
0342         if options.resetExpected:
0343             slope = old_div((yVals[-1] - yVals[0]), (xVals[-1] - xVals[0]))
0344             print("slope", slope)
0345             for index, old in enumerate (expectedVals):
0346                 expectedVals[index] = yVals[0] + \
0347                                     slope * (xVals[index] - xVals[0])
0348         #############################################
0349         ## Break Expected by Integrated Luminosity ##
0350         #############################################
0351         if options.breakExpectedIntLum:
0352             breakExpectedIntLum = []
0353             for chunk in options.breakExpectedIntLum:
0354                 pieces = sepRE.split (chunk)
0355                 try:
0356                     for piece in pieces:
0357                         breakExpectedIntLum.append( float(piece) )
0358                 except:
0359                     raise RuntimeError("'%s' from '%s' is not a valid float" \
0360                           % (piece, chunk))
0361             breakExpectedIntLum.sort()
0362             boundaries = []
0363             breakIndex = 0
0364             done = False
0365             for index, xPos in enumerate (xVals):
0366                 if xPos > breakExpectedIntLum[breakIndex]:
0367                     boundaries.append (index)
0368                     while breakIndex < len (breakExpectedIntLum):
0369                         breakIndex += 1
0370                         if breakIndex >= len (breakExpectedIntLum):
0371                             done = True
0372                             break
0373                         # If this next position is different, than
0374                         # we're golden.  Otherwise, let it go through
0375                         # the loop again.
0376                         if xPos <= breakExpectedIntLum[breakIndex]:
0377                             break
0378                     if done:
0379                         break
0380             # do we have any boundaries?
0381             if not boundaries:
0382                 raise RuntimeError("No values of 'breakExpectedIntLum' are in current range.")
0383             # is the first boundary at 0?  If not, add 0
0384             if boundaries[0]:
0385                 boundaries.insert (0, 0)
0386             # is the last boundary at the end?  If not, make the end a
0387             # boundary
0388             if boundaries[-1] != len (xVals) - 1:
0389                 boundaries.append( len (xVals) - 1 )
0390             rangeList = list(zip (boundaries, boundaries[1:]))
0391             for thisRange in rangeList:
0392                 upper = thisRange[1]
0393                 lower = thisRange[0]
0394                 slope = old_div((yVals[upper] - yVals[lower]), \
0395                         (xVals[upper] - xVals[lower]))
0396                 print("slope", slope)
0397                 # now go over the range inclusively
0398                 pairList = []
0399                 for index in range (lower, upper + 1):
0400                     newExpected = yVals[lower] + \
0401                                 slope * (xVals[index] - xVals[lower])
0402                     pairList.append( (xVals[index], newExpected) )
0403                     expectedVals[index] = newExpected
0404                 expectedChunks.append (pairList)
0405     ###########################################
0406     ## Instantanous Luminosity Ordering Mode ##
0407     ###########################################
0408     elif 'instLum' == options.edfMode or 'instIntLum' == options.edfMode:
0409         eventTupList = []
0410         if not lumiCont.xingInfo:
0411             raise RuntimeError("Luminosity Xing information missing.")
0412         for key, eventList in sorted( eventsDict.items() ):
0413             try:
0414                 lumi =  lumiCont[key]
0415                 instLum   = lumi.aveInstLum
0416                 fracAXIL  = lumi.fracAXILrecorded
0417                 totalAXIL = lumi.totalAXILrecorded
0418             except:
0419                 raise RuntimeError("key %s not found in lumi information" \
0420                       % key.__str__())
0421             for event in eventList:
0422                 eventTupList.append( (instLum, fracAXIL, totalAXIL, key,
0423                                       event[0], event[1], ) )
0424         eventTupList.sort()
0425         for eventTup in eventTupList:
0426             weight += eventTup[5]
0427             factor = old_div(weight, totalWeight)
0428             if 'instLum' == options.edfMode:
0429                 xVals.append (eventTup[0])
0430             else:
0431                 xVals.append (eventTup[2])
0432             yVals.append (factor)
0433             expectedVals.append (eventTup[1])
0434             predVals.append   (eventTup[1] * options.pred)
0435     else:
0436         raise RuntimeError("It looks like Charles screwed up if you are seeing this.")
0437
0438     size = len (xVals)
0439     step = int (old_div(math.sqrt(size), 2) + 1)
0440     if options.printValues:
0441         for index in range (size):
0442             print("%8f %8f %8f" % (xVals[index], yVals[index], expectedVals[index]), end=' ')
0443             if index > step:
0444                 denom = xVals[index] - xVals[index - step]
0445                 numer = yVals[index] - yVals[index - step]
0446                 if denom:
0447                     print(" %8f" % (old_div(numer, denom)), end=' ')
0448             if 0 == index % step:
0449                 print(" **", end=' ') ## indicates statistically independent
0450                              ## slope measurement
0451             print()
0452         print()
0453
0454     xArray = array.array ('d', xVals)
0455     yArray = array.array ('d', yVals)
0456     expected = array.array ('d', expectedVals)
0457     graph = ROOT.TGraph( size, xArray, yArray)
0458     graph.SetTitle (options.title)
0459     graph.SetMarkerStyle (20)
0460     expectedGraph = ROOT.TGraph( size, xArray, expected)
0461     expectedGraph.SetLineColor (ROOT.kRed)
0462     expectedGraph.SetLineWidth (3)
0463     if options.noDataPoints:
0464         expectedGraph.SetLineStyle (2)
0465
0466     # run statistical tests
0467     if options.weights:
0468         print("average weight per event:", old_div(weight, ( size - 1)))
0469     maxDistance = ROOT.TMath.KolmogorovTest (size, yArray,
0470                                              size, expected,
0471                                              "M")
0472     prob = ROOT.TMath.KolmogorovProb( maxDistance * math.sqrt( size ) )
0473
0474     # display everything
0475     ROOT.gROOT.SetStyle('Plain')
0476     ROOT.gROOT.SetBatch()
0477     c1 = ROOT.TCanvas()
0478     graph.GetXaxis().SetRangeUser (min (xVals), max (xVals))
0479     minValue = min (min(yVals), min(expected))
0480     if options.pred:
0481         minValue = min (minValue, min (predVals))
0482     graph.GetYaxis().SetRangeUser (minValue,
0483                                    max (max(yVals), max(expected), max(predVals)))
0484     graph.SetLineWidth (3)
0485     if options.noDataPoints:
0486         graph.Draw ("AL")
0487     else:
0488         graph.Draw ("ALP")
0489     if 'instLum' == options.edfMode:
0490         graph.GetXaxis().SetTitle ("Average Xing Inst. Luminosity (1/ub/s)")
0491         graph.GetXaxis().SetRangeUser (0., lumiCont.max('aveInstLum'))
0492     else:
0493         if 'instIntLum' == options.edfMode:
0494             graph.GetXaxis().SetTitle ("Integrated Luminosity - Inst. Lum. Ordered (1/%s)" \
0495                                        % lumiCont.invunits)
0496         else:
0497             graph.GetXaxis().SetTitle ("Integrated Luminosity (1/%s)" \
0498                                        % lumiCont.invunits)
0499     graph.GetYaxis().SetTitle ("Fraction of Events Seen")
0500     expectedGraphs = []
0501     if expectedChunks:
0502         for index, chunk in enumerate (expectedChunks):
0503             expectedXarray = array.array ('d', [item[0] for item in chunk])
0504             expectedYarray = array.array ('d', [item[1] for item in chunk])
0505             expectedGraph = ROOT.TGraph( len(chunk),
0506                                          expectedXarray,
0507                                          expectedYarray )
0508             expectedGraph.SetLineWidth (3)
0509             if options.noDataPoints:
0510                 expectedGraph.SetLineStyle (2)
0511             if index % 2:
0512                 expectedGraph.SetLineColor (ROOT.kBlue)
0513             else:
0514                 expectedGraph.SetLineColor (ROOT.kRed)
0515             expectedGraph.Draw("L")
0516             expectedGraphs.append (expectedGraph)
0517         exptectedGraph = expectedGraphs[0]
0518     else:
0519         expectedGraph.Draw ("L")
0520     green = 0
0521     if options.pred:
0522         predArray = array.array ('d', predVals)
0523         green = ROOT.TGraph (size, xArray, predArray)
0524         green.SetLineWidth (3)
0525         green.SetLineColor (8)
0526         green.Draw ('l')
0527     legend = ROOT.TLegend(0.15, 0.65, 0.50, 0.85)
0528     legend.SetFillStyle (0)
0529     legend.SetLineColor(ROOT.kWhite)
0530     observed = 'Observed'
0531     if options.weights:
0532         observed += ' (weighted)'
0533     legend.AddEntry(graph, observed,"PL")
0534     if options.resetExpected:
0535         legend.AddEntry(expectedGraph,  "Expected from partial yield","L")
0536     else:
0537         legend.AddEntry(expectedGraph,  "Expected from total yield","L")
0538     if options.pred:
0539         legend.AddEntry(green, options.predLabel,"L")
0540     legend.AddEntry("","D_{stat}=%.3f, N=%d" % (maxDistance, size),"")
0541     legend.AddEntry("","P_{KS}=%.3f" % prob,"")
0542     legend.Draw()
0543
0544     # save file
0545     c1.Print (outputFile)
0546
0547
0548 ######################
0549 ## ################ ##
0550 ## ## ########## ## ##
0551 ## ## ## Main ## ## ##
0552 ## ## ########## ## ##
0553 ## ################ ##
0554 ######################
0555
0556 if __name__ == '__main__':
0557     ##########################
0558     ## command line options ##
0559     ##########################
0560     allowedEDF = ['time', 'instLum', 'instIntLum']
0561     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.')
0562     plotGroup  = optparse.OptionGroup (parser, "Plot Options")
0563     rangeGroup = optparse.OptionGroup (parser, "Range Options")
0564     inputGroup = optparse.OptionGroup (parser, "Input Options")
0565     modeGroup  = optparse.OptionGroup (parser, "Mode Options")
0566     plotGroup.add_option ('--title', dest='title', type='string',
0567                           default = 'Empirical Distribution Function',
0568                           help = 'title of plot (default %default)')
0569     plotGroup.add_option ('--predicted', dest='pred', type='float',
0570                           default = 0,
0571                           help = 'factor by which predicted curve is greater than observed')
0572     plotGroup.add_option ('--predLabel', dest='predLabel', type='string',
0573                           default = 'Predicted',
0574                           help = 'label of predicted in legend')
0575     plotGroup.add_option ('--noDataPoints', dest='noDataPoints',
0576                           action='store_true',
0577                           help="Draw lines but no points for data")
0578     rangeGroup.add_option ('--minRun', dest='minRun', type='int', default=0,
0579                            help='Minimum run number to consider')
0580     rangeGroup.add_option ('--maxRun', dest='maxRun', type='int', default=0,
0581                            help='Maximum run number to consider')
0582     rangeGroup.add_option ('--minIntLum', dest='minIntLum', type='float', default=0,
0583                            help='Minimum integrated luminosity to consider')
0584     rangeGroup.add_option ('--maxIntLum', dest='maxIntLum', type='float', default=0,
0585                            help='Maximum integrated luminosity to consider')
0586     rangeGroup.add_option ('--resetExpected', dest='resetExpected',
0587                            action='store_true',
0588                            help='Reset expected from total yield to highest point considered')
0589     rangeGroup.add_option ('--breakExpectedIntLum', dest='breakExpectedIntLum',
0590                            type='string', action='append', default=[],
0591                            help='Break expected curve into pieces at integrated luminosity boundaries')
0592     inputGroup.add_option ('--ignoreNoLumiEvents', dest='ignore',
0593                            action='store_true',
0594                            help = 'Ignore (with a warning) events that do not have a lumi section')
0595     inputGroup.add_option ('--noWarnings', dest='noWarnings',
0596                            action='store_true',
0597                            help = 'Do not print warnings about missing luminosity information')
0598     inputGroup.add_option ('--runEventLumi', dest='relOrder',
0599                            action='store_true',
0600                            help = 'Parse event list assuming Run, Event #, Lumi# order')
0601     inputGroup.add_option ('--weights', dest='weights', action='store_true',
0602                            help = 'Read fourth column as a weight')
0603     modeGroup.add_option ('--print', dest='printValues', action='store_true',
0604                           help = 'Print X and Y values of EDF plot')
0605     modeGroup.add_option ('--runsWithLumis', dest='runsWithLumis',
0606                           type='string',action='append', default=[],
0607                           help='Print out run and lumi sections corresponding to integrated luminosities provided and then exits')
0608     modeGroup.add_option ('--edfMode', dest='edfMode', type='string',
0609                           default='time',
0610                           help="EDF Mode %s (default '%%default')" % allowedEDF)
0611     parser.add_option_group (plotGroup)
0612     parser.add_option_group (rangeGroup)
0613     parser.add_option_group (inputGroup)
0614     parser.add_option_group (modeGroup)
0615     (options, args) = parser.parse_args()
0616
0617     if options.edfMode not in allowedEDF:
0618         raise RuntimeError("edfMode (currently '%s') must be one of %s" \
0619               % (options.edfMode, allowedEDF))
0620
0621     if len (args) != 3 and not (options.runsWithLumis and len(args) >= 1):
0622         raise RuntimeError("Must provide lumi.csv, events.txt, and output.png")
0623
0624
0625     ##########################
0626     ## load Luminosity info ##
0627     ##########################
0628     cont = LumiInfoCont (args[0], **options.__dict__)
0629     cont.minRun    = options.minRun
0630     cont.maxRun    = options.maxRun
0631     cont.minIntLum = options.minIntLum
0632     cont.maxIntLum = options.maxIntLum
0633
0634     ##################################################
0635     ## look for which runs correspond to what total ##
0636     ## recorded integrated luminosity               ##
0637     ##################################################
0638     if options.runsWithLumis:
0639         recLumis = []
0640         for line in options.runsWithLumis:
0641             pieces = sepRE.split (line)
0642             for piece in pieces:
0643                 try:
0644                     recLumValue = float (piece)
0645                 except:
0646                     raise RuntimeError("'%s' in '%s' is not a float" % \
0647                           (piece, line))
0648                 if recLumValue <= 0:
0649                     raise RuntimeError("You must provide positive values for -runsWithLumis ('%f' given)" % recLumValue)
0650                 recLumis.append (recLumValue)
0651         if not recLumis:
0652             raise RuntimeError("What did Charles do now?")
0653         recLumis.sort()
0654         recLumIndex = 0
0655         recLumValue = recLumis [recLumIndex]
0656         prevRecLumi = 0.
0657         done = False
0658         for key, lumi in cont.items():
0659             if prevRecLumi >= recLumValue and recLumValue < lumi.totalRecorded:
0660                 # found it
0661                 print("%s contains total recorded lumi %f" % \
0662                       (key.__str__(), recLumValue))
0663                 while True:
0664                     recLumIndex += 1
0665                     if recLumIndex == len (recLumis):
0666                         done = True
0667                         break
0668                     recLumValue = recLumis [recLumIndex]
0669                     if prevRecLumi >= recLumValue and recLumValue < lumi.totalRecorded:
0670                         # found it
0671                         print("%s contains total recorded lumi %f" % \
0672                               (key.__str__(), recLumValue))
0673                     else:
0674                         break
0675                 if done:
0676                     break
0677             prevRecLumi = lumi.totalRecorded
0678         if recLumIndex < len (recLumis):
0679             print("Theses lumis not found: %s" % recLumis[recLumIndex:])
0680         sys.exit()
0681
0682     ####################
0683     ## make EDF plots ##
0684     ####################
0685     if options.edfMode != 'time' and not cont.xingInfo:
0686         raise RuntimeError("'%s' does not have Xing info" % args[0])
0687     eventsDict, totalWeight = loadEvents (args[1], cont, options)
0688     makeEDFplot (cont, eventsDict, totalWeight, args[2], options)