CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
estimatePileup.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 import os, sys
3 import coral
4 import array
5 import optparse
6 from RecoLuminosity.LumiDB import csvSelectionParser, selectionParser
7 import RecoLuminosity.LumiDB.lumiQueryAPI as LumiQueryAPI
8 import re
9 
10 from pprint import pprint
11 
12 def fillPileupHistogram (deadTable, parameters,
13  runNumber = 0, hist = None, debug = False,
14  mode='deadtable'):
15  '''Given a deadtable and run number, will (create if necessary
16  and) fill histogram with expected pileup distribution. If a
17  histogram is created, it is owned by the user and is his/her
18  responsibility to clean up the memory.'''
19  if hist:
20  maxBin = hist.GetNbinsX()
21  upper = int( hist.GetBinLowEdge(maxBin) + \
22  hist.GetBinWidth(maxBin) + 0.25 )
23  else:
24  histname = '%s_%s' % (parameters.pileupHistName, runNumber)
25  hist = ROOT.TH1D (histname, histname, parameters.maxPileupBin + 1,
26  -0.5, parameters.maxPileupBin + 0.5)
27  upper = parameters.maxPileupBin
28  for lumiSection, deadArray in sorted (deadTable.iteritems()):
29  if mode == 'csv':
30  numerator = float (deadArray[1])
31  denominator = float (deadArray[0])
32  instLumiArray = deadArray[2]
33  livetime = 1
34  if numerator < 0:
35  numerator = 0
36  if denominator:
37  livetime = numerator / denominator
38  else:
39  # we got everything from lumiDB
40  if len(deadArray) <= parameters.xingIndex:
41  # for some reason the xing instantaneous luminosity
42  # information isn't there. Print a warning and then skip
43  # it:
44  if parameters.noWarnings:
45  continue
46  if runNumber:
47  print "No Xing Instantaneous luminosity information for run %d, lumi section %d" \
48  % (runNumber, lumiSection)
49  else:
50  print "No Xing Instantaneous luminosity information for lumi section %d" \
51  % lumiSection
52  continue
53  numerator = float (deadArray[0])
54  denominator = float (deadArray[2] * deadArray[4])
55  xingInstLumiArray = deadArray[parameters.xingIndex]
56  # here we only want the instantaneous luminosities and don't
57  # care which crosings they fall in. So we only want the odd
58  instLumiArray = [(xingInstLumiArray[index], xingInstLumiArray[index + 1]) \
59  for index in xrange( 0, len (xingInstLumiArray), 2 ) ]
60  livetime = 1
61  if numerator < 0:
62  numerator = 0
63  if denominator:
64  livetime = 1 - numerator / denominator
65  # totalInstLumi = reduce(lambda x, y: x+y, instLumiArray) # not needed
66  for xing, xingInstLumi in instLumiArray:
67  xingIntLumi = xingInstLumi * parameters.lumiSectionLen * livetime
68  mean = xingInstLumi * parameters.minBiasXsec * \
69  parameters.rotationTime
70  if mean > 100:
71  if runNumber:
72  print "mean number of pileup events > 100 for run %d, lum %d : m %f l %f" % \
73  (runNumber, lumiSection, mean, xingInstLumi)
74  else:
75  print "mean number of pileup events > 100 for lum %d: m %f l %f" % \
76  (lumiSection, mean, xingInstLumi)
77 
78  for obs in range (upper):
79  prob = ROOT.TMath.Poisson (obs, mean)
80  totalProb += prob
81  hist.Fill (obs, prob * xingIntLumi)
82  if debug:
83  print "ls", lumiSection, "xing", xing, "inst", xingInstLumi, \
84  "mean", mean, "totalProb", totalProb, 1 - totalProb
85  print " hist mean", hist.GetMean()
86  if totalProb < 1:
87  hist.Fill (obs, (1 - totalProb) * xingIntLumi)
88  return hist
89 
90 
91 
92 ##############################
93 ## ######################## ##
94 ## ## ################## ## ##
95 ## ## ## Main Program ## ## ##
96 ## ## ################## ## ##
97 ## ######################## ##
98 ##############################
99 
100 if __name__ == '__main__':
101  parameters = LumiQueryAPI.ParametersObject()
102  beamModeChoices = ["","stable", "quiet", "either"]
103  parser = optparse.OptionParser ("Usage: %prog [--options] output.root",
104  description = "Script to estimate pileup distribution using xing instantaneous luminosity information and minimum bias cross section. Output is TH1D stored in root file")
105  dbGroup = optparse.OptionGroup (parser, "Database Options")
106  inputGroup = optparse.OptionGroup (parser, "Input Options")
107  pileupGroup = optparse.OptionGroup (parser, "Pileup Options")
108  dbGroup.add_option ('--parameters', dest = 'connect', action = 'store',
109  help = 'connect string to lumiDB (optional, default to frontier://LumiCalc/CMS_LUMI_PROD)')
110  dbGroup.add_option ('-P', dest = 'authpath', action = 'store',
111  help = 'path to authentication file')
112  dbGroup.add_option ('--siteconfpath', dest = 'siteconfpath', action = 'store',
113  help = 'specific path to site-local-config.xml file, default to $CMS_PATH/SITECONF/local/JobConfig, ' \
114  'if path undefined, fallback to cern proxy&server')
115  dbGroup.add_option ('--debug', dest = 'debug', action = 'store_true',
116  help = 'debug')
117  inputGroup.add_option ('-r', dest = 'runnumber', action = 'store',
118  help = 'run number')
119  inputGroup.add_option ('-i', dest = 'inputfile', action = 'store',
120  help = 'lumi range selection file')
121  inputGroup.add_option ('-b', dest = 'beammode', default='', choices=beamModeChoices,
122  help = "beam mode, optional for delivered action, default ('%%default' out of %s)" % beamModeChoices)
123  inputGroup.add_option ('--lumiversion', dest = 'lumiversion', type='string', default='0001',
124  help = 'lumi data version, optional for all, default %default')
125  inputGroup.add_option ('--hltpath', dest = 'hltpath', action = 'store',
126  help = 'specific hltpath to calculate the recorded luminosity, default to all')
127  inputGroup.add_option ('--csvInput', dest = 'csvInput', type='string', default='',
128  help = 'Use CSV file from lumiCalc.py instead of lumiDB')
129  pileupGroup.add_option ('--xingMinLum', dest = 'xingMinLum', type='float', default = 1e-3,
130  help = 'Minimum luminosity considered for "lsbylsXing" action, default %default')
131  pileupGroup.add_option ('--minBiasXsec', dest = 'minBiasXsec', type='float', default = parameters.minBiasXsec,
132  help = 'Minimum bias cross section assumed (in mb), default %default mb')
133  pileupGroup.add_option ('--histName', dest='histName', type='string', default = parameters.pileupHistName,
134  help = 'Histrogram name (default %default)')
135  pileupGroup.add_option ('--maxPileupBin', dest='maxPileupBin', type='int', default = parameters.maxPileupBin,
136  help = 'Maximum bin in pileup histogram (default %default)')
137  pileupGroup.add_option ('--saveRuns', dest='saveRuns', action='store_true',
138  help = 'Save pileup histograms for individual runs')
139  pileupGroup.add_option ('--debugLumi', dest='debugLumi',
140  action='store_true',
141  help = 'Print out debug info for individual lumis')
142  pileupGroup.add_option ('--nowarning', dest = 'nowarning',
143  action = 'store_true', default = False,
144  help = 'suppress bad for lumi warnings' )
145  parser.add_option_group (dbGroup)
146  parser.add_option_group (inputGroup)
147  parser.add_option_group (pileupGroup)
148  # parse arguments
149  (options, args) = parser.parse_args()
150  import ROOT
151  if not args:
152  parser.print_usage()
153  sys.exit()
154  if len (args) != 1:
155  parser.print_usage()
156  raise RuntimeError("Exactly one output file must be given")
157  output = args[0]
158 
159  # get database session hooked up
160  if options.authpath:
161  os.environ['CORAL_AUTH_PATH'] = options.authpath
162 
163  ## Save what we need in the parameters object
164  parameters.verbose = True
165  parameters.noWarnings = options.nowarning
166  parameters.lumiXing = True
167  parameters.lumiversion = options.lumiversion
168  if options.beammode=='stable':
169  parameters.beammode = 'STABLE BEAMS'
170  parameters.xingMinLum = options.xingMinLum
171  parameters.minBiasXsec = options.minBiasXsec
172  parameters.pileupHistName = options.histName
173  parameters.maxPileupBin = options.maxPileupBin
174 
175  session, svc = \
176  LumiQueryAPI.setupSession (options.connect or \
177  'frontier://LumiCalc/CMS_LUMI_PROD',
178  options.siteconfpath, parameters,options.debug)
179 
180  ## Let's start the fun
181  if not options.inputfile and not options.runnumber and not options.csvInput:
182  raise "must specify either a run (-r), an input run selection file (-i), or an input CSV file (--csvInput)"
183  pileupHist = ROOT.TH1D(parameters.pileupHistName, parameters.pileupHistName,
184  1001,
185  0.0, 25.)
186  histList = []
187  if options.csvInput:
188  # we're going to read in the CSV file and use this as not only
189  # the selection of which run/events to use, but also the
190  # source of the lumi information.
191  sepRE = re.compile (r'[\s,;:]+')
192  runLumiDict = {}
193  events = open (options.csvInput, 'r')
194  csvDict = {}
195  for line in events:
196  pieces = sepRE.split (line.strip())
197  if len (pieces) < 8:
198  continue
199  #if len (pieces) % 2:
200  # # not an even number
201  # continue
202  try:
203  run, lumi = int ( pieces[0] ), int ( pieces[2] )
204  delivered, recorded = float( pieces[7] ), float( pieces[8] )
205  xingInstLumiArray = [( int(orbit), float(lum) ) \
206  for orbit, lum in zip( pieces[9::2],
207  pieces[10::2] ) ]
208  except:
209  continue
210  csvDict.setdefault (run, {})[lumi] = \
211  ( delivered, recorded, xingInstLumiArray )
212  events.close()
213  for runNumber, lumiDict in sorted( csvDict.iteritems() ):
214  if options.saveRuns:
215  hist = fillPileupHistogram (lumiDict, parameters,
216  runNumber = runNumber,
217  debug = options.debugLumi,
218  mode='csv')
219  pileupHist.Add (hist)
220  histList.append (hist)
221  else:
222  fillPileupHistogram (lumiDict, parameters,
223  hist = pileupHist,
224  debug = options.debugLumi,
225  mode='csv')
226 
227  histFile = ROOT.TFile.Open (output, 'recreate')
228  if not histFile:
229  raise RuntimeError("Could not open '%s' as an output root file" % output)
230  pileupHist.Write()
231  for hist in histList:
232  hist.Write()
233  histFile.Close()
234  #pprint (csvDict)
235  sys.exit()
236 
237 
238  ## Get input source
239  if options.runnumber:
240  inputRange = options.runnumber
241  else:
242  basename, extension = os.path.splitext (options.inputfile)
243  if extension == '.csv': # if file ends with .csv, use csv
244  # parser, else parse as json file
245  fileparsingResult = csvSelectionParser.csvSelectionParser (options.inputfile)
246  else:
247  f = open (options.inputfile, 'r')
248  inputfilecontent = f.read()
249  inputRange = selectionParser.selectionParser (inputfilecontent)
250  if not inputRange:
251  print 'failed to parse the input file', options.inputfile
252  raise
253 
254  recordedData = LumiQueryAPI.recordedLumiForRange (session, parameters, inputRange)
255  ## pprint (recordedData)
256  for runDTarray in recordedData:
257  runNumber = runDTarray[0]
258  deadTable = runDTarray[2]
259  if options.saveRuns:
260  hist = fillPileupHistogram (deadTable, parameters,
261  runNumber = runNumber,
262  debug = options.debugLumi)
263  pileupHist.Add (hist)
264  histList.append (hist)
265  else:
266  fillPileupHistogram (deadTable, parameters,
267  hist = pileupHist,
268  debug = options.debugLumi)
269  histFile = ROOT.TFile.Open (output, 'recreate')
270  if not histFile:
271  raise RuntimeError("Could not open '%s' as an output root file" % output)
272  pileupHist.Write()
273  for hist in histList:
274  hist.Write()
275  histFile.Close()
276 
Generated for CMSSW Reference Manual by doxygen 1.8.5