Reference Manual

Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 VERSION='2.00'
00003 import os, sys
00004 import coral
00005 import optparse
00006 from RecoLuminosity.LumiDB import sessionManager,csvSelectionParser,selectionParser,lumiCorrections,lumiCalcAPI
00007 
00008 beamChoices=['PROTPHYS','IONPHYS']
00009 
00010 class pileupParameters(object):
00011     def __init__(self):
00012         self.minBiasXsec=71300 #unit microbarn
00013         self.rotationRate=11245.613 # for 3.5 TeV Beam energy
00014         self.NBX=3564
00015         self.rotationTime=1.0/self.rotationRate
00016         self.lumiSectionLen=2**18*self.rotationTime
00017 
00018 def fillPileupHistogram (bxlumiinfo,pileupHistName,maxPileupBin,
00019                          runNumber=0, hist = None, debug = False):
00020     '''
00021     bxlumiinfo:[[cmslsnum(0),avgdelivered(1),avgrecorded(2),bxlumiarray[3]]]
00022 
00023     Given luminfo , deadfrac info and run number, will (create if necessary
00024     and) fill histogram with expected pileup distribution.  If a
00025     histogram is created, it is owned by the user and is his/her
00026     responsibility to clean up the memory.'''
00027     if hist:
00028         maxBin = hist.GetNbinsX()
00029         upper = int( hist.GetBinLowEdge(maxBin) + \
00030                      hist.GetBinWidth(maxBin) + 0.25 )
00031     else:
00032         histname = '%s_%s' % (pileupHistName, runNumber)
00033         hist = ROOT.TH1D (histname, histname, maxPileupBin + 1,
00034                           -0.5,maxPileupBin + 0.5)
00035         upper = maxPileupBin
00036     p=pileupParameters()
00037     for perlsinfo in bxlumiinfo:
00038         cmslsnum=perlsinfo[0]
00039         avgdelivered=perlsinfo[1]
00040         avgrecorded=perlsinfo[2]
00041         bxdata=perlsinfo[3]
00042         bxidx=bxdata[0]
00043         bxvaluelist=bxdata[1]
00044         #calculate livefrac
00045         livetime=1
00046         if avgrecorded<0:
00047             avgrecorded=0
00048         if avgdelivered:
00049             livetime=avgrecorded/avgdelivered
00050         else:
00051             livetime=0
00052         for idx,bxvalue in enumerate(bxvaluelist):
00053             xingIntLumi=bxvalue * p.lumiSectionLen * livetime
00054             if options.minBiasXsec:
00055                 mean = bxvalue * options.minBiasXsec * p.rotationTime
00056             else:
00057                 mean = bxvalue * p.minBiasXsec * p.rotationTime
00058             if mean > 100:
00059                 if runNumber:
00060                     print "mean number of pileup events > 100 for run %d, lum %d : m %f l %f" % \
00061                           (runNumber, lumiSection, mean, bxvalue)
00062                 else:
00063                     print "mean number of pileup events > 100 for lum %d: m %f l %f" % \
00064                           (cmslsnum, mean, bxvalue)
00065             totalProb = 0
00066             for obs in range (upper):
00067                 prob = ROOT.TMath.Poisson (obs, mean)
00068                 totalProb += prob
00069                 hist.Fill (obs, prob * xingIntLumi)
00070             if debug:
00071                 xing=bxidx[idx]
00072                 print "ls", lumiSection, "xing", xing, "inst", bxvalue, \
00073                       "mean", mean, "totalProb", totalProb, 1 - totalProb
00074                 print "  hist mean", hist.GetMean()
00075             if totalProb < 1:
00076                 hist.Fill (obs, (1 - totalProb) * xingIntLumi)
00077     return hist
00078     
00079 ##############################
00080 ## ######################## ##
00081 ## ## ################## ## ##
00082 ## ## ## Main Program ## ## ##
00083 ## ## ################## ## ##
00084 ## ######################## ##
00085 ##############################
00086 
00087 if __name__ == '__main__':
00088     beamModeChoices = [ "","stable", "quiet", "either"]
00089     amodetagChoices = [ "PROTPHYS","IONPHYS" ]
00090     xingAlgoChoices =[ "OCC1","OCC2","ET"]
00091     parser = optparse.OptionParser ("Usage: %prog [--options] output.root",
00092                                     description = "Script to estimate pileup distribution using xing instantaneous luminosity information and minimum bias cross section.  Output is TH1D stored in root file")
00093     dbGroup     = optparse.OptionGroup (parser, "Database Options")
00094     inputGroup  = optparse.OptionGroup (parser, "Input Options")
00095     pileupGroup = optparse.OptionGroup (parser, "Pileup Options")
00096     dbGroup.add_option     ('-c', dest = 'connect', action = 'store',
00097                             default='frontier://LumiCalc/CMS_LUMI_PROD',
00098                             help = 'connect string to lumiDB ,default %default')
00099     dbGroup.add_option     ('-P', dest = 'authpath', action = 'store',
00100                              help = 'path to authentication file')
00101     dbGroup.add_option     ('--siteconfpath', dest = 'siteconfpath', action = 'store',
00102                              help = 'specific path to site-local-config.xml file, default to $CMS_PATH/SITECONF/local/JobConfig, ' \
00103                              'if path undefined, fallback to cern proxy&server')
00104     dbGroup.add_option     ('--debug', dest = 'debug', action = 'store_true',
00105                             help = 'debug')
00106     
00107     inputGroup.add_option  ('-r', dest = 'runnumber', action = 'store',
00108                             help = 'run number')
00109     inputGroup.add_option  ('-i', dest = 'inputfile', action = 'store',
00110                             help = 'lumi range selection file')
00111     inputGroup.add_option  ('-b', dest = 'beamstatus', default='', choices=beamModeChoices,
00112                             help = "beam mode, optional for delivered action, default ('%%default' out of %s)" % beamModeChoices)
00113     inputGroup.add_option  ('--hltpath', dest = 'hltpath', action = 'store',
00114                             help = 'specific hltpath to calculate the recorded luminosity, default to all')
00115     inputGroup.add_option  ('--csvInput', dest = 'csvInput', type='string', default='',
00116                             help = 'Use CSV file from lumiCalc2.py instead of lumiDB')
00117     inputGroup.add_option  ('--without-correction', dest = 'withoutFineCorrection', action='store_true',
00118                             help = 'not apply fine correction on calibration')
00119     pileupGroup.add_option('--algoname',dest='algoname',default='OCC1',
00120                             help = 'lumi algorithm , default %default')
00121     pileupGroup.add_option ('--xingMinLum', dest = 'xingMinLum', type='float', default = 1.0e-3,
00122                             help = 'Minimum luminosity considered for "lsbylsXing" action, default %default')
00123     pileupGroup.add_option ('--minBiasXsec', dest = 'minBiasXsec', type='float', default = 71300,
00124                             help = 'Minimum bias cross section assumed (in mb), default %default mb')
00125     pileupGroup.add_option ('--histName', dest='pileupHistName', type='string', default = 'pileup',
00126                             help = 'Histrogram name (default %default)')
00127     pileupGroup.add_option ('--maxPileupBin', dest='maxPileupBin', type='int', default = 10,
00128                             help = 'Maximum bin in pileup histogram (default %default)')
00129     pileupGroup.add_option ('--saveRuns', dest='saveRuns', action='store_true',
00130                             help = 'Save pileup histograms for individual runs')
00131     pileupGroup.add_option ('--debugLumi', dest='debugLumi',
00132                             action='store_true',
00133                             help = 'Print out debug info for individual lumis')
00134     pileupGroup.add_option ('--nowarning', dest = 'nowarning',
00135                             action = 'store_true', default = False,
00136                             help = 'suppress bad for lumi warnings' )
00137     parser.add_option_group (dbGroup)
00138     parser.add_option_group (inputGroup)
00139     parser.add_option_group (pileupGroup)
00140     # parse arguments
00141     try:
00142         (options, args) = parser.parse_args()
00143     except Exception , e:
00144         print e
00145     if not args:
00146         parser.print_usage()
00147         sys.exit()
00148     if len (args) != 1:
00149         parser.print_usage()
00150         raise RuntimeError, "Exactly one output file must be given"
00151     output = args[0]
00152     finecorrections=None
00153     # get database session hooked up
00154     if options.authpath:
00155         os.environ['CORAL_AUTH_PATH'] = options.authpath
00156 
00157     svc=sessionManager.sessionManager(options.connect,authpath=options.authpath,debugON=options.debug)
00158 
00159     if not options.inputfile and not options.runnumber and not options.csvInput:
00160         raise "must specify either a run (-r), an input run selection file (-i), or an input CSV file (--csvInput)"
00161 
00162     runDict = {}
00163     if options.csvInput:
00164         import re
00165         # we're going to read in the CSV file and use this as not only
00166         # the selection of which run/events to use, but also the
00167         # source of the lumi information.
00168         sepRE = re.compile (r'[\s,;:]+')
00169         events = open (options.csvInput, 'r')
00170         for line in events:
00171             pieces = sepRE.split (line.strip())
00172             if len (pieces) < 6:
00173                 continue
00174             if len (pieces) % 2:
00175                 # not an even number
00176                 continue
00177             try:
00178                 run,       cmslsnum = int  ( pieces[0] ), int  ( pieces[1] )
00179                 delivered, recorded = float( pieces[2] ), float( pieces[3] )
00180                 xingIdx = [int(myidx) for myidx in  pieces[4::2] ]
00181                 xingVal = [float(myval) for myval in pieces[5::2] ]
00182                 
00183                 #xingInstLumiArray = [( int(orbit), float(lum) ) \
00184                 #                     for orbit, lum in zip( pieces[4::2],
00185                 #                                            pieces[5::2] ) ]
00186             except:
00187                 continue
00188             runDict.setdefault(run,[]).append([cmslsnum,delivered,recorded,(xingIdx,xingVal)])
00189             #{run:[[cmslsnum,delivered,recorded,xingInstlumiArray]..]}
00190         events.close()
00191     else:
00192         session=svc.openSession(isReadOnly=True,cpp2sqltype=[('unsigned int','NUMBER(10)'),('unsigned long long','NUMBER(20)')])
00193         finecorrections=None
00194         if options.runnumber:
00195             inputRange = {int(options.runnumber):None}
00196         else:
00197             basename, extension = os.path.splitext (options.inputfile)
00198             if extension == '.csv': # if file ends with .csv, use csv
00199                 # parser, else parse as json file
00200                 fileparsingResult = csvSelectionParser.csvSelectionParser (options.inputfile)
00201             else:
00202                 f = open (options.inputfile, 'r')
00203                 inputfilecontent = f.read()
00204                 inputRange =  selectionParser.selectionParser (inputfilecontent).runsandls()
00205         if not inputRange:
00206             print 'failed to parse the input file', options.inputfile
00207             raise
00208         if not options.withoutFineCorrection:
00209             rruns=inputRange.keys()
00210             schema=session.nominalSchema()
00211             session.transaction().start(True)
00212             finecorrections=lumiCorrections.correctionsForRange(schema,rruns)
00213             session.transaction().commit()
00214         session.transaction().start(True)
00215         schema=session.nominalSchema()
00216         lumiData=lumiCalcAPI.lumiForRange(schema,inputRange,beamstatus=options.beamstatus,withBXInfo=True,bxAlgo=options.algoname,xingMinLum=options.xingMinLum,withBeamIntensity=False,datatag=None,finecorrections=finecorrections)
00217         session.transaction().commit()
00218         # {run:[lumilsnum(0),cmslsnum(1),timestamp(2),beamstatus(3),beamenergy(4),deliveredlumi(5),recordedlumi(6),calibratedlumierror(7),(bxidx,bxvalues,bxerrs)(8),None]}}
00219     ##convert lumiData to lumiDict format #{run:[[cmslsnum,avg]]}
00220         for runnum,perrundata in lumiData.items():
00221             bxlumiinfo=[]
00222             for perlsdata in perrundata:
00223                 cmslsnum=perlsdata[1]
00224                 deliveredlumi=perlsdata[5]
00225                 recordedlumi=perlsdata[6]
00226                 bxlist=perlsdata[8]
00227                 bxlumiinfo.append([cmslsnum,deliveredlumi,recordedlumi,bxlist])
00228                 runDict.setdefault(runnum,[]).append([cmslsnum,deliveredlumi,recordedlumi,bxlist])
00229     #print 'runDict ',runDict
00230     
00231     import ROOT 
00232     pileupHist = ROOT.TH1D (options.pileupHistName,options.pileupHistName,
00233                       options.maxPileupBin + 1,
00234                       -0.5, options.maxPileupBin + 0.5)
00235     histList = []
00236     for runNumber, lumiList in sorted( runDict.iteritems() ):
00237         if options.saveRuns:
00238             hist = fillPileupHistogram (lumiList,options.pileupHistName,options.maxPileupBin,
00239                                         runNumber = runNumber,
00240                                         debug = options.debugLumi)
00241             pileupHist.Add (hist)
00242             histList.append (hist)
00243         else:
00244             fillPileupHistogram (lumiList, options.pileupHistName,options.maxPileupBin,
00245                                  hist = pileupHist,
00246                                  debug = options.debugLumi)            
00247     histFile = ROOT.TFile.Open (output, 'recreate')
00248     if not histFile:
00249         raise RuntimeError, \
00250               "Could not open '%s' as an output root file" % output
00251     pileupHist.Write()
00252     for hist in histList:
00253         hist.Write()
00254     histFile.Close()
00255     sys.exit()
/data/refman/pasoursint/CMSSW_5_3_10/src/RecoLuminosity/LumiDB/scripts/estimatePileup2.py
