test
CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
estimatePileup2.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 VERSION='2.00'
3 import os, sys
4 import coral
5 import optparse
6 from RecoLuminosity.LumiDB import sessionManager,csvSelectionParser,selectionParser,lumiCorrections,lumiCalcAPI
7 
8 beamChoices=['PROTPHYS','IONPHYS']
9 
10 class pileupParameters(object):
11  def __init__(self):
12  self.minBiasXsec=71300 #unit microbarn
13  self.rotationRate=11245.613 # for 3.5 TeV Beam energy
14  self.NBX=3564
15  self.rotationTime=1.0/self.rotationRate
16  self.lumiSectionLen=2**18*self.rotationTime
17 
18 def fillPileupHistogram (bxlumiinfo,pileupHistName,maxPileupBin,
19  runNumber=0, hist = None, debug = False):
20  '''
21  bxlumiinfo:[[cmslsnum(0),avgdelivered(1),avgrecorded(2),bxlumiarray[3]]]
22 
23  Given luminfo , deadfrac info and run number, will (create if necessary
24  and) fill histogram with expected pileup distribution. If a
25  histogram is created, it is owned by the user and is his/her
26  responsibility to clean up the memory.'''
27  if hist:
28  maxBin = hist.GetNbinsX()
29  upper = int( hist.GetBinLowEdge(maxBin) + \
30  hist.GetBinWidth(maxBin) + 0.25 )
31  else:
32  histname = '%s_%s' % (pileupHistName, runNumber)
33  hist = ROOT.TH1D (histname, histname, maxPileupBin + 1,
34  -0.5,maxPileupBin + 0.5)
35  upper = maxPileupBin
36  p=pileupParameters()
37  for perlsinfo in bxlumiinfo:
38  cmslsnum=perlsinfo[0]
39  avgdelivered=perlsinfo[1]
40  avgrecorded=perlsinfo[2]
41  bxdata=perlsinfo[3]
42  bxidx=bxdata[0]
43  bxvaluelist=bxdata[1]
44  #calculate livefrac
45  livetime=1
46  if avgrecorded<0:
47  avgrecorded=0
48  if avgdelivered:
49  livetime=avgrecorded/avgdelivered
50  else:
51  livetime=0
52  for idx,bxvalue in enumerate(bxvaluelist):
53  xingIntLumi=bxvalue * p.lumiSectionLen * livetime
54  if options.minBiasXsec:
55  mean = bxvalue * options.minBiasXsec * p.rotationTime
56  else:
57  mean = bxvalue * p.minBiasXsec * p.rotationTime
58  if mean > 100:
59  if runNumber:
60  print "mean number of pileup events > 100 for run %d, lum %d : m %f l %f" % \
61  (runNumber, lumiSection, mean, bxvalue)
62  else:
63  print "mean number of pileup events > 100 for lum %d: m %f l %f" % \
64  (cmslsnum, mean, bxvalue)
65  totalProb = 0
66  for obs in range (upper):
67  prob = ROOT.TMath.Poisson (obs, mean)
68  totalProb += prob
69  hist.Fill (obs, prob * xingIntLumi)
70  if debug:
71  xing=bxidx[idx]
72  print "ls", lumiSection, "xing", xing, "inst", bxvalue, \
73  "mean", mean, "totalProb", totalProb, 1 - totalProb
74  print " hist mean", hist.GetMean()
75  if totalProb < 1:
76  hist.Fill (obs, (1 - totalProb) * xingIntLumi)
77  return hist
78 
79 ##############################
80 ## ######################## ##
81 ## ## ################## ## ##
82 ## ## ## Main Program ## ## ##
83 ## ## ################## ## ##
84 ## ######################## ##
85 ##############################
86 
87 if __name__ == '__main__':
88  beamModeChoices = [ "","stable", "quiet", "either"]
89  amodetagChoices = [ "PROTPHYS","IONPHYS" ]
90  xingAlgoChoices =[ "OCC1","OCC2","ET"]
91  parser = optparse.OptionParser ("Usage: %prog [--options] output.root",
92  description = "Script to estimate pileup distribution using xing instantaneous luminosity information and minimum bias cross section. Output is TH1D stored in root file")
93  dbGroup = optparse.OptionGroup (parser, "Database Options")
94  inputGroup = optparse.OptionGroup (parser, "Input Options")
95  pileupGroup = optparse.OptionGroup (parser, "Pileup Options")
96  dbGroup.add_option ('-c', dest = 'connect', action = 'store',
97  default='frontier://LumiCalc/CMS_LUMI_PROD',
98  help = 'connect string to lumiDB ,default %default')
99  dbGroup.add_option ('-P', dest = 'authpath', action = 'store',
100  help = 'path to authentication file')
101  dbGroup.add_option ('--siteconfpath', dest = 'siteconfpath', action = 'store',
102  help = 'specific path to site-local-config.xml file, default to $CMS_PATH/SITECONF/local/JobConfig, ' \
103  'if path undefined, fallback to cern proxy&server')
104  dbGroup.add_option ('--debug', dest = 'debug', action = 'store_true',
105  help = 'debug')
106 
107  inputGroup.add_option ('-r', dest = 'runnumber', action = 'store',
108  help = 'run number')
109  inputGroup.add_option ('-i', dest = 'inputfile', action = 'store',
110  help = 'lumi range selection file')
111  inputGroup.add_option ('-b', dest = 'beamstatus', default='', choices=beamModeChoices,
112  help = "beam mode, optional for delivered action, default ('%%default' out of %s)" % beamModeChoices)
113  inputGroup.add_option ('--hltpath', dest = 'hltpath', action = 'store',
114  help = 'specific hltpath to calculate the recorded luminosity, default to all')
115  inputGroup.add_option ('--csvInput', dest = 'csvInput', type='string', default='',
116  help = 'Use CSV file from lumiCalc2.py instead of lumiDB')
117  inputGroup.add_option ('--without-correction', dest = 'withoutFineCorrection', action='store_true',
118  help = 'not apply fine correction on calibration')
119  pileupGroup.add_option('--algoname',dest='algoname',default='OCC1',
120  help = 'lumi algorithm , default %default')
121  pileupGroup.add_option ('--xingMinLum', dest = 'xingMinLum', type='float', default = 1.0e-3,
122  help = 'Minimum luminosity considered for "lsbylsXing" action, default %default')
123  pileupGroup.add_option ('--minBiasXsec', dest = 'minBiasXsec', type='float', default = 71300,
124  help = 'Minimum bias cross section assumed (in mb), default %default mb')
125  pileupGroup.add_option ('--histName', dest='pileupHistName', type='string', default = 'pileup',
126  help = 'Histrogram name (default %default)')
127  pileupGroup.add_option ('--maxPileupBin', dest='maxPileupBin', type='int', default = 10,
128  help = 'Maximum bin in pileup histogram (default %default)')
129  pileupGroup.add_option ('--saveRuns', dest='saveRuns', action='store_true',
130  help = 'Save pileup histograms for individual runs')
131  pileupGroup.add_option ('--debugLumi', dest='debugLumi',
132  action='store_true',
133  help = 'Print out debug info for individual lumis')
134  pileupGroup.add_option ('--nowarning', dest = 'nowarning',
135  action = 'store_true', default = False,
136  help = 'suppress bad for lumi warnings' )
137  parser.add_option_group (dbGroup)
138  parser.add_option_group (inputGroup)
139  parser.add_option_group (pileupGroup)
140  # parse arguments
141  try:
142  (options, args) = parser.parse_args()
143  except Exception as e:
144  print e
145  if not args:
146  parser.print_usage()
147  sys.exit()
148  if len (args) != 1:
149  parser.print_usage()
150  raise RuntimeError("Exactly one output file must be given")
151  output = args[0]
152  finecorrections=None
153  # get database session hooked up
154  if options.authpath:
155  os.environ['CORAL_AUTH_PATH'] = options.authpath
156 
157  svc=sessionManager.sessionManager(options.connect,authpath=options.authpath,debugON=options.debug)
158 
159  if not options.inputfile and not options.runnumber and not options.csvInput:
160  raise "must specify either a run (-r), an input run selection file (-i), or an input CSV file (--csvInput)"
161 
162  runDict = {}
163  if options.csvInput:
164  import re
165  # we're going to read in the CSV file and use this as not only
166  # the selection of which run/events to use, but also the
167  # source of the lumi information.
168  sepRE = re.compile (r'[\s,;:]+')
169  events = open (options.csvInput, 'r')
170  for line in events:
171  pieces = sepRE.split (line.strip())
172  if len (pieces) < 6:
173  continue
174  if len (pieces) % 2:
175  # not an even number
176  continue
177  try:
178  run, cmslsnum = int ( pieces[0] ), int ( pieces[1] )
179  delivered, recorded = float( pieces[2] ), float( pieces[3] )
180  xingIdx = [int(myidx) for myidx in pieces[4::2] ]
181  xingVal = [float(myval) for myval in pieces[5::2] ]
182 
183  #xingInstLumiArray = [( int(orbit), float(lum) ) \
184  # for orbit, lum in zip( pieces[4::2],
185  # pieces[5::2] ) ]
186  except:
187  continue
188  runDict.setdefault(run,[]).append([cmslsnum,delivered,recorded,(xingIdx,xingVal)])
189  #{run:[[cmslsnum,delivered,recorded,xingInstlumiArray]..]}
190  events.close()
191  else:
192  session=svc.openSession(isReadOnly=True,cpp2sqltype=[('unsigned int','NUMBER(10)'),('unsigned long long','NUMBER(20)')])
193  finecorrections=None
194  if options.runnumber:
195  inputRange = {int(options.runnumber):None}
196  else:
197  basename, extension = os.path.splitext (options.inputfile)
198  if extension == '.csv': # if file ends with .csv, use csv
199  # parser, else parse as json file
200  fileparsingResult = csvSelectionParser.csvSelectionParser (options.inputfile)
201  else:
202  f = open (options.inputfile, 'r')
203  inputfilecontent = f.read()
204  inputRange = selectionParser.selectionParser (inputfilecontent).runsandls()
205  if not inputRange:
206  print 'failed to parse the input file', options.inputfile
207  raise
208  if not options.withoutFineCorrection:
209  rruns=inputRange.keys()
210  schema=session.nominalSchema()
211  session.transaction().start(True)
212  finecorrections=lumiCorrections.correctionsForRange(schema,rruns)
213  session.transaction().commit()
214  session.transaction().start(True)
215  schema=session.nominalSchema()
216  lumiData=lumiCalcAPI.lumiForRange(schema,inputRange,beamstatus=options.beamstatus,withBXInfo=True,bxAlgo=options.algoname,xingMinLum=options.xingMinLum,withBeamIntensity=False,datatag=None,finecorrections=finecorrections)
217  session.transaction().commit()
218  # {run:[lumilsnum(0),cmslsnum(1),timestamp(2),beamstatus(3),beamenergy(4),deliveredlumi(5),recordedlumi(6),calibratedlumierror(7),(bxidx,bxvalues,bxerrs)(8),None]}}
219  ##convert lumiData to lumiDict format #{run:[[cmslsnum,avg]]}
220  for runnum,perrundata in lumiData.items():
221  bxlumiinfo=[]
222  for perlsdata in perrundata:
223  cmslsnum=perlsdata[1]
224  deliveredlumi=perlsdata[5]
225  recordedlumi=perlsdata[6]
226  bxlist=perlsdata[8]
227  bxlumiinfo.append([cmslsnum,deliveredlumi,recordedlumi,bxlist])
228  runDict.setdefault(runnum,[]).append([cmslsnum,deliveredlumi,recordedlumi,bxlist])
229  #print 'runDict ',runDict
230 
231  import ROOT
232  pileupHist = ROOT.TH1D (options.pileupHistName,options.pileupHistName,
233  options.maxPileupBin + 1,
234  -0.5, options.maxPileupBin + 0.5)
235  histList = []
236  for runNumber, lumiList in sorted( runDict.iteritems() ):
237  if options.saveRuns:
238  hist = fillPileupHistogram (lumiList,options.pileupHistName,options.maxPileupBin,
239  runNumber = runNumber,
240  debug = options.debugLumi)
241  pileupHist.Add (hist)
242  histList.append (hist)
243  else:
244  fillPileupHistogram (lumiList, options.pileupHistName,options.maxPileupBin,
245  hist = pileupHist,
246  debug = options.debugLumi)
247  histFile = ROOT.TFile.Open (output, 'recreate')
248  if not histFile:
249  raise RuntimeError("Could not open '%s' as an output root file" % output)
250  pileupHist.Write()
251  for hist in histList:
252  hist.Write()
253  histFile.Close()
254  sys.exit()
start
Definition: start.py:1
bitset_utilities::append
boost::dynamic_bitset append(const boost::dynamic_bitset<> &bs1, const boost::dynamic_bitset<> &bs2)
this method takes two bitsets bs1 and bs2 and returns result of bs2 appended to the end of bs1 ...
Definition: bitset_append.cc:10
Generated for CMSSW Reference Manual by doxygen 1.8.5