pileupCalc Namespace Reference

Functions
def	fillPileupHistogram
def	MyErf
def	parseInputFile

Variables
list	CalculationModeChoices = ['truth', 'observed']
string	default = 'PileupCalc.root'
string	help = 'output root file'
tuple	histFile = ROOT.TFile.Open(output, 'recreate')
tuple	inpf = open(options.inputfile, 'r')
tuple	inputfilecontent = inpf.read()
tuple	inputPileupRange = parseInputFile(options.inputLumiJSON)
tuple	inputRange = selectionParser.selectionParser(inputfilecontent)
list	LSPUlist = inputPileupRange[run]
list	lumiInfo = LSPUlist[LSnumber]
	nbins = options.numPileupBins
list	output = args[0]
tuple	parser
	## Main Program More...
tuple	pileupHist
	type = float,
	upper = options.maxPileupBin
string	VERSION = '1.00'

Function Documentation

def pileupCalc.fillPileupHistogram	(		lumiInfo,
			calcOption,
			hist,
			minbXsec,
			Nbins
	)

lumiinfo:[intlumi per LS, mean interactions ]

intlumi is the deadtime corrected average integraged lumi per lumisection

Definition at line 54 of file pileupCalc.py.

References MyErf(), and mathSSE.sqrt().

def fillPileupHistogram (lumiInfo, calcOption, hist, minbXsec, Nbins):
    '''
    lumiinfo:[intlumi per LS, mean interactions ]

    intlumi is the deadtime corrected average integraged lumi per lumisection
    '''

    LSintLumi = lumiInfo[0]
    RMSInt = lumiInfo[1]*minbXsec
    AveNumInt = lumiInfo[2]*minbXsec

    #coeff = 0

    #if RMSInt > 0:
    #    coeff = 1.0/RMSInt/sqrt(6.283185)

    #expon = 2.0*RMSInt*RMSInt

    Sqrt2 = sqrt(2)

    ##Nbins = hist.GetXaxis().GetNbins()

    ProbFromRMS = []
    BinWidth = hist.GetBinWidth(1)

    # First, re-constitute lumi distribution for this LS from RMS:
    if RMSInt > 0:

        AreaLnew = -10.
        AreaL = 0

        for obs in range (Nbins):
            #Old Gaussian normalization; inaccurate for small rms and large bins
            #val = hist.GetBinCenter(obs+1)
            #prob = coeff*exp(-1.0*(val-AveNumInt)*(val-AveNumInt)/expon)
            #ProbFromRMS.append(prob)
            
            left = hist.GetBinLowEdge(obs+1)
            right = left+BinWidth

            argR = (AveNumInt-right)/Sqrt2/RMSInt
            AreaR = MyErf(argR)

            if AreaLnew<-5.:
                argL = (AveNumInt-left)/Sqrt2/RMSInt
                AreaL = MyErf(argL)
            else:
                AreaL = AreaLnew
                AreaLnew = AreaR  # save R bin value for L next time

            NewProb = (AreaL-AreaR)*0.5

            ProbFromRMS.append(NewProb)

            #print left, right, argL, argR, AreaL, AreaR, NewProb

    else:
        obs = hist.FindBin(AveNumInt)
        for bin in range (Nbins):
            ProbFromRMS.append(0.0)
        if obs<Nbins+1:            
            ProbFromRMS[obs] = 1.0
        if AveNumInt < 1.0E-5:
           ProbFromRMS[obs] = 0.  # just ignore zero values
        
    if calcOption == 'true':  # Just put distribution into histogram
        if RMSInt > 0:
            totalProb = 0
            for obs in range (Nbins):
                prob = ProbFromRMS[obs]
                val = hist.GetBinCenter(obs+1)
                #print obs, val, RMSInt,coeff,expon,prob
                totalProb += prob
                hist.Fill (val, prob * LSintLumi)
                
            if 1.0-totalProb > 0.01:
                print "Significant probability density outside of your histogram"
                print "Consider using a higher value of --maxPileupBin"
                print "Mean %f, RMS %f, Integrated probability %f" % (AveNumInt,RMSInt,totalProb)
            #    hist.Fill (val, (1 - totalProb) * LSintLumi)
        else:
            hist.Fill(AveNumInt,LSintLumi)
    else: # have to convolute with a poisson distribution to get observed Nint
        totalProb = 0
        Peak = 0
        BinWidth = hist.GetBinWidth(1)
        for obs in range (Nbins):
            Peak = hist.GetBinCenter(obs+1)
            RMSWeight = ProbFromRMS[obs]
            for bin in range (Nbins):
                val = hist.GetBinCenter(bin+1)-0.5*BinWidth
                prob = ROOT.TMath.Poisson (val, Peak)
                totalProb += prob
                hist.Fill (val, prob * LSintLumi * RMSWeight)

        if 1.0-totalProb > 0.01:
            print "Significant probability density outside of your histogram"
            print "Consider using a higher value of --maxPileupBin"


    return hist