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 from math import sqrt
00010 
00011 from pprint import pprint
00012 
00013 def CalcPileup (deadTable, parameters, mode='deadtable'):
00014     '''Given a deadtable, will calculate parameters of pileup distribution. Return formatted
00015     string with LumiSection, LS integrated lumi, RMS of bunch to bunch lumi and pileup.'''
00016 
00017     LumiString = ""
00018     LumiArray = []
00019 
00020     for lumiSection, deadArray in sorted (deadTable.iteritems()):
00021         numerator = 0
00022         if mode == 'csv':
00023             numerator     = float (deadArray[1])
00024             denominator   = float (deadArray[0])
00025             instLumiArray =        deadArray[2]
00026             livetime = 1
00027             if numerator < 0:
00028                 numerator = 0
00029             if denominator:
00030                 livetime = numerator / denominator
00031 
00032         else:
00033             print "no csv input! Doh!"
00034             return
00035         # totalInstLumi = reduce(lambda x, y: x+y, instLumiArray) # not needed
00036         if lumiSection > 0:
00037              TotalLumi = 0 
00038              TotalInt = 0
00039              TotalInt2 = 0
00040              TotalWeight = 0
00041              TotalWeight2 = 0
00042              FilledXings = 0
00043              for xing, xingInstLumi in instLumiArray:
00044                 xingIntLumi = xingInstLumi * parameters.lumiSectionLen * livetime
00045                 mean = xingInstLumi * parameters.rotationTime
00046                 if mean > 100:
00047                     if runNumber:
00048                         print "mean number of pileup events > 100 for run %d, lum %d : m %f l %f" % \
00049                           (runNumber, lumiSection, mean, xingInstLumi)
00050                     else:
00051                         print "mean number of pileup events > 100 for lum %d: m %f l %f" % \
00052                           (lumiSection, mean, xingInstLumi)
00053 #                print "mean number of pileup events for lum %d: m %f idx %d l %f" % (lumiSection, mean, xing, xingIntLumi)
00054 
00055                 if xingInstLumi > 0.1:
00056 
00057                     TotalLumi = TotalLumi+xingIntLumi
00058                     TotalInt+= mean*xingIntLumi
00059                     FilledXings = FilledXings+1
00060                     #print "xing inst lumi %f %f %d" % (xingIntLumi,TotalLumi,FilledXings)
00061 
00062              #compute weighted mean, then loop again to get weighted RMS       
00063              MeanInt = 0
00064              if TotalLumi >0:
00065                  MeanInt = TotalInt/TotalLumi
00066              for xing, xingInstLumi in instLumiArray:
00067                  if xingInstLumi > 0.1:
00068                      xingIntLumi = xingInstLumi * parameters.lumiSectionLen * livetime
00069                      mean = xingInstLumi * parameters.rotationTime
00070                      TotalInt2+= xingIntLumi*(mean-MeanInt)*(mean-MeanInt)
00071                      TotalWeight+= xingIntLumi
00072                      TotalWeight2+= xingIntLumi*xingIntLumi
00073 
00074 
00075 
00076         if ((lumiSection > 0)):
00077             #print " LS, Total lumi, filled xings %d, %f, %d" %(lumiSection,TotalLumi,FilledXings)
00078             if FilledXings > 0:
00079                 AveLumi = TotalLumi/FilledXings
00080             else:
00081                 AveLumi = 0
00082             RMSLumi = 0
00083             Denom = TotalWeight*TotalWeight-TotalWeight2
00084             if TotalLumi > 0 and Denom > 0:
00085                 RMSLumi = sqrt(TotalWeight/(TotalWeight*TotalWeight-TotalWeight2)*TotalInt2)
00086             LumiString = "[%d,%2.4e,%2.4e,%2.4e]," % (lumiSection, TotalLumi, RMSLumi, MeanInt)
00087             LumiArray.append(lumiSection)
00088             LumiArray.append(TotalLumi)  # should really weight by total lumi in LS
00089             LumiArray.append(RMSLumi)
00090             LumiArray.append(MeanInt)
00091             
00092     return LumiArray
00093     
00094 
00095 
00096 ##############################
00097 ## ######################## ##
00098 ## ## ################## ## ##
00099 ## ## ## Main Program ## ## ##
00100 ## ## ################## ## ##
00101 ## ######################## ##
00102 ##############################
00103 
00104 # modified from the estimatePileup.py script in RecoLuminosity/LumiDB
00105 # 5 Jan, 2012  Mike Hildreth
00106 # As of now, .csv file is the only allowed input.  See the old script for a variety
00107 # of different ways to access the data.
00108 
00109 if __name__ == '__main__':
00110     parameters = LumiQueryAPI.ParametersObject()
00111     parser = optparse.OptionParser ("Usage: %prog [--options] output.root",
00112                                     description = "Script to estimate average instantaneous bunch crossing luminosity using xing instantaneous luminosity information. Output is JSON format file with one entry per LumiSection")
00113     inputGroup  = optparse.OptionGroup (parser, "Input Options")
00114     pileupGroup = optparse.OptionGroup (parser, "Pileup Options")
00115     inputGroup.add_option  ('--csvInput', dest = 'csvInput', type='string', default='',
00116                             help = 'Use CSV file from lumiCalc.py instead of lumiDB')
00117     parser.add_option_group (inputGroup)
00118     parser.add_option_group (pileupGroup)
00119     # parse arguments
00120     (options, args) = parser.parse_args()
00121 
00122     if not args:
00123         parser.print_usage()
00124         sys.exit()
00125     if len (args) != 1:
00126         parser.print_usage()
00127         raise RuntimeError, "Please specify an output file as your last argument"
00128     output = args[0]
00129 
00130     ## Let's start the fun
00131     if not options.csvInput:
00132         raise "you must specify an input CSV file with (--csvInput)"
00133 
00134     OUTPUTLINE = ""
00135     if options.csvInput:
00136         # we're going to read in the CSV file and use this as not only
00137         # the selection of which run/events to use, but also the
00138         # source of the lumi information.
00139         sepRE = re.compile (r'[\s,;:]+')
00140         events = open (options.csvInput, 'r')
00141         OldRun = -1
00142 
00143         InGap = 0;
00144         GapDict = {}
00145         LastValidLumi = []
00146         LastDelivered = 0
00147 
00148 
00149         OUTPUTLINE+='{'
00150         for line in events:
00151             runLumiDict = {}    
00152             csvDict = {}
00153             pieces = sepRE.split (line.strip())
00154             if len (pieces) < 8: # means we are missing data; keep track of LS, lumi
00155                 InGap = 1
00156                 try:
00157                     run,       lumi     = int  ( pieces[0] ), int  ( pieces[2] )
00158                     delivered, recorded = float( pieces[7] ), float( pieces[8] )
00159                 except:
00160                     if pieces[0] != 'run':
00161                         print " cannot parse csv file "
00162                     InGap = 0
00163                     continue
00164                 GapDict[lumi] = [delivered, recorded]
00165                 continue
00166             #if len (pieces) % 2:
00167                 # not an even number
00168             #    continue
00169             try:
00170                 run,       lumi     = int  ( pieces[0] ), int  ( pieces[2] )
00171                 delivered, recorded = float( pieces[7] ), float( pieces[8] )
00172                 xingInstLumiArray = [( int(orbit), float(lum) ) \
00173                                      for orbit, lum in zip( pieces[9::2],
00174                                                             pieces[10::2] ) ]
00175             except:
00176                 print " Bad Parsing: Check if the input format has changed"
00177                 print pieces[0],pieces[1],pieces[2],pieces[3],pieces[4],pieces[5],pieces[6]
00178                 continue
00179 
00180             csvDict.setdefault (run, {})[lumi] = \
00181                                ( delivered, recorded, xingInstLumiArray )
00182 
00183             if run != OldRun:
00184                 if OldRun>0:
00185                     if InGap == 1:  # We have some LS's at the end with no data
00186                         lastLumiS = 0
00187                         for lumiS, lumiInfo in sorted ( GapDict.iteritems() ):
00188                             record = lumiInfo[1]
00189                             lastLumiS = lumiS
00190                             if record > 0.01:
00191 
00192                                 peakratio = lumiInfo[0]/LastDelivered # roughly, ratio of inst lumi
00193                                 pileup = LastValidLumi[3]*peakratio     # scale for this LS
00194                                 aveLumi = 0
00195                                 if lumiInfo[0] >0:
00196                                     aveLumi = LastValidLumi[1]*peakratio*lumiInfo[1]/lumiInfo[0]  # scale by rec/del
00197                                 LumiString = "[%d,%2.4e,%2.4e,%2.4e]," % (lumiS, aveLumi, LastValidLumi[2],pileup)
00198                                 OUTPUTLINE += LumiString
00199                             else:
00200                                 LumiString = "[%d,0.0,0.0,0.0]," % (lumiS)
00201                                 OUTPUTLINE+=LumiString
00202                         #add one empty one at the end
00203                         LumiString = "[%d,0.0,0.0,0.0]," % (lastLumiS+1)
00204                         OUTPUTLINE+=LumiString
00205 
00206                         InGap = 0
00207                         GapDict.clear()
00208 
00209                     # add one empty lumi section at the end of each run, to mesh with JSON files
00210                     LumiString = "[%d,0.0,0.0,0.0]," % (LastValidLumi[0]+1)
00211                     OUTPUTLINE+=LumiString
00212 
00213                     lastindex=len(OUTPUTLINE)-1
00214                     trunc = OUTPUTLINE[0:lastindex]
00215                     OUTPUTLINE = trunc
00216                     OUTPUTLINE += '], '
00217                 OldRun = run
00218                 OUTPUTLINE+= ('"%d":' % run )
00219                 OUTPUTLINE+= ' ['
00220 
00221                 if lumi == 2:  # there is a missing LS=1 for this run
00222                     OUTPUTLINE+= '[1,0.0,0.0,0.0],'
00223 
00224             for runNumber, lumiDict in sorted( csvDict.iteritems() ):
00225 
00226                 LumiArray = CalcPileup (lumiDict, parameters,
00227                                      mode='csv')
00228 
00229                 LastValidLumi = LumiArray
00230                 LastDelivered = lumiDict[LumiArray[0]][0] 
00231 
00232                 if InGap == 1:  # We have some gap before this in entry in this run
00233                     for lumiS, lumiInfo in sorted ( GapDict.iteritems() ):
00234                         peakratio = lumiInfo[0]/LastDelivered # roughly, ratio of inst lumi
00235                         pileup = LumiArray[3]*peakratio     # scale for this LS
00236                         aveLumi = 0
00237                         if lumiInfo[0] > 0:
00238                             aveLumi = LumiArray[1]*peakratio*lumiInfo[1]/lumiInfo[0]  # scale by rec/del
00239                         LumiString = "[%d,%2.4e,%2.4e,%2.4e]," % (lumiS, aveLumi, LumiArray[2],pileup)
00240                         OUTPUTLINE += LumiString
00241                     InGap = 0
00242                 LumiString = "[%d,%2.4e,%2.4e,%2.4e]," % (LumiArray[0], LumiArray[1], LumiArray[2], LumiArray[3])
00243                 OUTPUTLINE += LumiString
00244 
00245 
00246         lastindex=len(OUTPUTLINE)-1
00247         trunc = OUTPUTLINE[0:lastindex]
00248         OUTPUTLINE = trunc
00249         OUTPUTLINE += ']}'
00250         events.close()
00251 
00252         outputfile = open(output,'w')
00253         if not outputfile:
00254             raise RuntimeError, \
00255                "Could not open '%s' as an output JSON file" % output
00256                     
00257         outputfile.write(OUTPUTLINE)
00258         outputfile.close()
00259 
00260         sys.exit()
00261 
00262
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoLuminosity/LumiDB/scripts/estimatePileup_makeJSON.py
