Reference Manual

Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 import os, sys
00003 import coral
00004 import array
00005 import optparse
00006 from RecoLuminosity.LumiDB import csvSelectionParser, selectionParser
00007 import RecoLuminosity.LumiDB.lumiQueryAPI as LumiQueryAPI
00008 import re
00009 
00010 from pprint import pprint
00011 
00012 def fillPileupHistogram (deadTable, parameters,
00013                          runNumber = 0, hist = None, debug = False,
00014                          mode='deadtable'):
00015     '''Given a deadtable and run number, will (create if necessary
00016     and) fill histogram with expected pileup distribution.  If a
00017     histogram is created, it is owned by the user and is his/her
00018     responsibility to clean up the memory.'''
00019     if hist:
00020         maxBin = hist.GetNbinsX()
00021         upper = int( hist.GetBinLowEdge(maxBin) + \
00022                      hist.GetBinWidth(maxBin) + 0.25 )
00023     else:
00024         histname = '%s_%s' % (parameters.pileupHistName, runNumber)
00025         hist = ROOT.TH1F (histname, histname, parameters.maxPileupBin + 1,
00026                           -0.5, parameters.maxPileupBin + 0.5)
00027         upper  = parameters.maxPileupBin
00028     for lumiSection, deadArray in sorted (deadTable.iteritems()):
00029         if mode == 'csv':
00030             numerator     = float (deadArray[1])
00031             denominator   = float (deadArray[0])
00032             instLumiArray =        deadArray[2]
00033             livetime = 1
00034             if numerator < 0:
00035                 numerator = 0
00036             if denominator:
00037                 livetime = numerator / denominator
00038         else:
00039             # we got everything from lumiDB
00040             if len(deadArray) <= parameters.xingIndex:
00041                 # for some reason the xing instantaneous luminosity
00042                 # information isn't there.  Print a warning and then skip
00043                 # it:
00044                 if parameters.noWarnings:
00045                     continue
00046                 if runNumber:
00047                     print "No Xing Instantaneous luminosity information for run %d, lumi section %d" \
00048                           % (runNumber, lumiSection)
00049                 else:
00050                     print "No Xing Instantaneous luminosity information for lumi section %d" \
00051                           % lumiSection
00052                 continue
00053             numerator   = float (deadArray[0])
00054             denominator = float (deadArray[2] * deadArray[4])
00055             xingInstLumiArray = deadArray[parameters.xingIndex]
00056             # here we only want the instantaneous luminosities and don't
00057             # care which crosings they fall in.  So we only want the odd 
00058             instLumiArray = [(xingInstLumiArray[index], xingInstLumiArray[index + 1]) \
00059                              for index in xrange( 0, len (xingInstLumiArray), 2 ) ]
00060             livetime = 1
00061             if numerator < 0:
00062                 numerator = 0
00063             if denominator:
00064                 livetime = 1 - numerator / denominator
00065         # totalInstLumi = reduce(lambda x, y: x+y, instLumiArray) # not needed
00066         for xing, xingInstLumi in instLumiArray:
00067             xingIntLumi = xingInstLumi * parameters.lumiSectionLen * livetime
00068             mean = xingInstLumi * parameters.minBiasXsec * \
00069                    parameters.rotationTime
00070             if mean > 100:
00071                 if runNumber:
00072                     print "mean number of pileup events > 100 for run %d, lum %d : m %f l %f" % \
00073                           (runNumber, lumiSection, mean, xingInstLumi)
00074                 else:
00075                     print "mean number of pileup events > 100 for lum %d: m %f l %f" % \
00076                           (lumiSection, mean, xingInstLumi)
00077             totalProb = 0
00078             for obs in range (upper):
00079                 prob = ROOT.TMath.Poisson (obs, mean)
00080                 totalProb += prob
00081                 hist.Fill (obs, prob * xingIntLumi)
00082             if debug:
00083                 print "ls", lumiSection, "xing", xing, "inst", xingInstLumi, \
00084                       "mean", mean, "totalProb", totalProb, 1 - totalProb
00085                 print "  hist mean", hist.GetMean()
00086             if totalProb < 1:
00087                 hist.Fill (obs, (1 - totalProb) * xingIntLumi)
00088     return hist
00089     
00090 
00091 
00092 ##############################
00093 ## ######################## ##
00094 ## ## ################## ## ##
00095 ## ## ## Main Program ## ## ##
00096 ## ## ################## ## ##
00097 ## ######################## ##
00098 ##############################
00099 
00100 if __name__ == '__main__':
00101     parameters = LumiQueryAPI.ParametersObject()
00102     beamModeChoices = ["","stable", "quiet", "either"]
00103     parser = optparse.OptionParser ("Usage: %prog [--options] output.root",
00104                                     description = "Script to estimate pileup distribution using xing instantaneous luminosity information and minimum bias cross section.  Output is TH1F stored in root file")
00105     dbGroup     = optparse.OptionGroup (parser, "Database Options")
00106     inputGroup  = optparse.OptionGroup (parser, "Input Options")
00107     pileupGroup = optparse.OptionGroup (parser, "Pileup Options")
00108     dbGroup.add_option     ('--parameters', dest = 'connect', action = 'store',
00109                            help = 'connect string to lumiDB (optional, default to frontier://LumiCalc/CMS_LUMI_PROD)')
00110     dbGroup.add_option     ('-P', dest = 'authpath', action = 'store',
00111                              help = 'path to authentication file')
00112     dbGroup.add_option     ('--siteconfpath', dest = 'siteconfpath', action = 'store',
00113                              help = 'specific path to site-local-config.xml file, default to $CMS_PATH/SITECONF/local/JobConfig, ' \
00114                              'if path undefined, fallback to cern proxy&server')
00115     dbGroup.add_option     ('--debug', dest = 'debug', action = 'store_true',
00116                             help = 'debug')
00117     inputGroup.add_option  ('-r', dest = 'runnumber', action = 'store',
00118                             help = 'run number')
00119     inputGroup.add_option  ('-i', dest = 'inputfile', action = 'store',
00120                             help = 'lumi range selection file')
00121     inputGroup.add_option  ('-b', dest = 'beammode', default='', choices=beamModeChoices,
00122                             help = "beam mode, optional for delivered action, default ('%%default' out of %s)" % beamModeChoices)
00123     inputGroup.add_option  ('--lumiversion', dest = 'lumiversion', type='string', default='0001',
00124                             help = 'lumi data version, optional for all, default %default')
00125     inputGroup.add_option  ('--hltpath', dest = 'hltpath', action = 'store',
00126                            help = 'specific hltpath to calculate the recorded luminosity, default to all')
00127     inputGroup.add_option  ('--csvInput', dest = 'csvInput', type='string', default='',
00128                             help = 'Use CSV file from lumiCalc.py instead of lumiDB')
00129     pileupGroup.add_option ('--xingMinLum', dest = 'xingMinLum', type='float', default = 1e-3,
00130                             help = 'Minimum luminosity considered for "lsbylsXing" action, default %default')
00131     pileupGroup.add_option ('--minBiasXsec', dest = 'minBiasXsec', type='float', default = parameters.minBiasXsec,
00132                             help = 'Minimum bias cross section assumed (in mb), default %default mb')
00133     pileupGroup.add_option ('--histName', dest='histName', type='string', default = parameters.pileupHistName,
00134                             help = 'Histrogram name (default %default)')
00135     pileupGroup.add_option ('--maxPileupBin', dest='maxPileupBin', type='int', default = parameters.maxPileupBin,
00136                             help = 'Maximum bin in pileup histogram (default %default)')
00137     pileupGroup.add_option ('--saveRuns', dest='saveRuns', action='store_true',
00138                             help = 'Save pileup histograms for individual runs')
00139     pileupGroup.add_option ('--debugLumi', dest='debugLumi',
00140                             action='store_true',
00141                             help = 'Print out debug info for individual lumis')
00142     pileupGroup.add_option ('--nowarning', dest = 'nowarning',
00143                             action = 'store_true', default = False,
00144                             help = 'suppress bad for lumi warnings' )
00145     parser.add_option_group (dbGroup)
00146     parser.add_option_group (inputGroup)
00147     parser.add_option_group (pileupGroup)
00148     # parse arguments
00149     (options, args) = parser.parse_args()
00150     import ROOT 
00151     if not args:
00152         parser.print_usage()
00153         sys.exit()
00154     if len (args) != 1:
00155         parser.print_usage()
00156         raise RuntimeError, "Exactly one output file must be given"
00157     output = args[0]
00158 
00159     # get database session hooked up
00160     if options.authpath:
00161         os.environ['CORAL_AUTH_PATH'] = options.authpath
00162 
00163     ## Save what we need in the parameters object
00164     parameters.verbose        = True
00165     parameters.noWarnings     = options.nowarning
00166     parameters.lumiXing       = True
00167     parameters.lumiversion    = options.lumiversion
00168     if options.beammode=='stable':
00169         parameters.beammode    = 'STABLE BEAMS'
00170     parameters.xingMinLum     = options.xingMinLum
00171     parameters.minBiasXsec    = options.minBiasXsec
00172     parameters.pileupHistName = options.histName
00173     parameters.maxPileupBin   = options.maxPileupBin
00174 
00175     session, svc =  \
00176              LumiQueryAPI.setupSession (options.connect or \
00177                                         'frontier://LumiCalc/CMS_LUMI_PROD',
00178                                         options.siteconfpath, parameters,options.debug)
00179 
00180     ## Let's start the fun
00181     if not options.inputfile and not options.runnumber and not options.csvInput:
00182         raise "must specify either a run (-r), an input run selection file (-i), or an input CSV file (--csvInput)"
00183     pileupHist = ROOT.TH1F (parameters.pileupHistName, parameters.pileupHistName,
00184                       parameters.maxPileupBin + 1,
00185                       -0.5, parameters.maxPileupBin + 0.5)
00186     histList = []
00187     if options.csvInput:
00188         # we're going to read in the CSV file and use this as not only
00189         # the selection of which run/events to use, but also the
00190         # source of the lumi information.
00191         sepRE = re.compile (r'[\s,;:]+')
00192         runLumiDict = {}    
00193         events = open (options.csvInput, 'r')
00194         csvDict = {}
00195         for line in events:
00196             pieces = sepRE.split (line.strip())
00197             if len (pieces) < 6:
00198                 continue
00199             if len (pieces) % 2:
00200                 # not an even number
00201                 continue
00202             try:
00203                 run,       lumi     = int  ( pieces[0] ), int  ( pieces[1] )
00204                 delivered, recorded = float( pieces[2] ), float( pieces[3] )
00205                 xingInstLumiArray = [( int(orbit), float(lum) ) \
00206                                      for orbit, lum in zip( pieces[4::2],
00207                                                             pieces[5::2] ) ]
00208             except:
00209                 continue
00210             csvDict.setdefault (run, {})[lumi] = \
00211                                ( delivered, recorded, xingInstLumiArray )
00212         events.close()
00213         for runNumber, lumiDict in sorted( csvDict.iteritems() ):
00214             if options.saveRuns:
00215                 hist = fillPileupHistogram (lumiDict, parameters,
00216                                             runNumber = runNumber,
00217                                             debug = options.debugLumi,
00218                                             mode='csv')
00219                 pileupHist.Add (hist)
00220                 histList.append (hist)
00221             else:
00222                 fillPileupHistogram (lumiDict, parameters,
00223                                      hist = pileupHist,
00224                                      debug = options.debugLumi,
00225                                      mode='csv')
00226             
00227         histFile = ROOT.TFile.Open (output, 'recreate')
00228         if not histFile:
00229             raise RuntimeError, \
00230                   "Could not open '%s' as an output root file" % output
00231         pileupHist.Write()
00232         for hist in histList:
00233             hist.Write()
00234         histFile.Close()
00235         #pprint (csvDict)
00236         sys.exit()
00237 
00238         
00239     ## Get input source
00240     if options.runnumber:
00241         inputRange = options.runnumber
00242     else:
00243         basename, extension = os.path.splitext (options.inputfile)
00244         if extension == '.csv': # if file ends with .csv, use csv
00245                                 # parser, else parse as json file
00246             fileparsingResult = csvSelectionParser.csvSelectionParser (options.inputfile)
00247         else:
00248             f = open (options.inputfile, 'r')
00249             inputfilecontent = f.read()
00250             inputRange =  selectionParser.selectionParser (inputfilecontent)
00251         if not inputRange:
00252             print 'failed to parse the input file', options.inputfile
00253             raise 
00254 
00255     recordedData  = LumiQueryAPI.recordedLumiForRange  (session, parameters, inputRange)
00256     ## pprint (recordedData)
00257     for runDTarray in recordedData:
00258         runNumber = runDTarray[0]
00259         deadTable = runDTarray[2]
00260         if options.saveRuns:
00261             hist = fillPileupHistogram (deadTable, parameters,
00262                                         runNumber = runNumber,
00263                                         debug = options.debugLumi)
00264             pileupHist.Add (hist)
00265             histList.append (hist)
00266         else:
00267             fillPileupHistogram (deadTable, parameters,
00268                                  hist = pileupHist,
00269                                  debug = options.debugLumi)
00270     histFile = ROOT.TFile.Open (output, 'recreate')
00271     if not histFile:
00272         raise RuntimeError, \
00273               "Could not open '%s' as an output root file" % output
00274     pileupHist.Write()
00275     for hist in histList:
00276         hist.Write()
00277     histFile.Close()
00278
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/RecoLuminosity/LumiDB/scripts/estimatePileup.py
