pileupCalc Namespace Reference

Classes
class	EquidistantBinning

Functions
def	fillPileupHistogram
def	MyErf
def	parseInputFile
def	poisson

Variables
tuple	binning = EquidistantBinning(options.numPileupBins, 0., options.maxPileupBin)
list	choices = ['true', 'observed']
string	default = 'pileup'
tuple	hContents = np.zeros(binning.centers.shape)
string	help = 'input Run/LS file for your analysis in JSON format'
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]
tuple	options = parser.parse_args()
	output = options.outputfile
tuple	parser = argparse.ArgumentParser(description = "Script to estimate pileup distribution using bunch instantaneous luminosity information and minimum bias cross section. Outputs a TH1D of the pileup distribution stored in a ROOT file.")
	## Main Program More...
tuple	pileupHist
tuple	req_group = parser.add_argument_group('required arguments')
tuple	Sqrt2 = np.sqrt(2)
float	type = 69200.0
string	VERSION = '1.00'

Function Documentation

def pileupCalc.fillPileupHistogram	(		lumiInfo,
			calcOption,
			binning,
			hContents,
			minbXsec,
			run,
			ls
	)

lumiinfo:[intlumi per LS, mean interactions ]

intlumi is the deadtime corrected average integrated lumi per lumisection

Definition at line 69 of file pileupCalc.py.

References MyErf(), poisson(), and print().

def fillPileupHistogram(lumiInfo, calcOption, binning, hContents, minbXsec, run, ls):
    '''
    lumiinfo:[intlumi per LS, mean interactions ]

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

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

    # First, re-constitute lumi distribution for this LS from RMS:
    if RMSInt > 0:
        areaAbove = MyErf((AveNumInt-binning.edges)/Sqrt2/RMSInt)
        ## area above edge, so areaAbove[i]-areaAbove[i+1] = area in bin
        ProbFromRMS = .5*(areaAbove[:-1]-areaAbove[1:])
    else:
        ProbFromRMS = np.zeros(hContents.shape)
        obs = binning.find(AveNumInt)
        if ( obs < binning.num ) and ( AveNumInt >= 1.0E-5 ): # just ignore zero values
            ProbFromRMS[obs] = 1.0
        
    if calcOption == 'true':  # Just put distribution into histogram
        if RMSInt > 0:
            hContents += ProbFromRMS*LSintLumi
            totalProb = np.sum(ProbFromRMS)

            if 1.0-totalProb > 0.01:
                print("Run %d, LS %d: Significant probability density outside of your histogram (mean %.2f," % (run, ls, AveNumInt))
                print("rms %.2f, integrated probability %.3f). Consider using a higher value of --maxPileupBin." % (RMSInt, totalProb))
        else:
            hContents[obs] += LSintLumi ## obs = FindBin(AveNumInt), -1 because hContents has no overflows
    else: # have to convolute with a poisson distribution to get observed Nint
        if not hasattr(binning, "poissConv"): ## only depends on binning, cache
            ## poissConv[i,j] = TMath.Poisson(e[i], c[j])
            binning.poissConv = poisson(
                binning.edges[:-1,np.newaxis], ## e'[i,] = e[i]
                binning.centers[np.newaxis,:]) ## c'[,j] = c[j]
        # prob[i] = sum_j ProbFromRMS[j]*TMath.Poisson(e[i], c[j])
        prob = np.sum(binning.poissConv * ProbFromRMS[np.newaxis,:], axis=1)
        hContents += prob*LSintLumi
        #if ( not np.all(ProbFromRMS == 0) ) and 1.0-np.sum(prob) > 0.01:
        #    print("Run %d, LS %d: significant probability density outside of your histogram, %f" % (run, ls, np.sum(prob)))
        #    print("Consider using a higher value of --maxPileupBin")