print_tbl_nvtx_metxy_mean_var_eta.py

print_tbl_nvtx_metxy_mean_var_eta.py

#!/usr/bin/env python
# Tai Sakuma <sakuma@fnal.gov>
import ROOT
import sys
import math
import json
import re
import signal
from optparse import OptionParser

sys.path.insert(1, './libs')
import Binning
from funcCMSSW import *

ROOT.gROOT.SetBatch(1)

##____________________________________________________________________________||
GLOBAL_LAST = False

##____________________________________________________________________________||
parser = OptionParser()
parser.add_option('-i', '--inputPath', default = None, action = 'store', type = 'string')
parser.add_option('--inputPathListFile', default = None, action = 'store', type = 'string')
parser.add_option('-j', '--jsonPath', default = None, action = 'store', type = 'string')
parser.add_option("-n", "--nevents", action = "store", default = -1, type = 'long', help = "maximum number of events to process")
parser.add_option("-v", "--verbose", action="store_true", default=False)
(options, args) = parser.parse_args(sys.argv)

##____________________________________________________________________________||
absEta_cuts = [3, 3.5, 4, 4.5, 5, 5.5]

##____________________________________________________________________________||
def main():
    printHeader()
    inputPathList = mkInputPathList(options.inputPath, options.inputPathListFile)
    if noEvents(inputPathList): return
    counts, stats = count(inputPathList)
    printCounts(counts, stats)

##____________________________________________________________________________||
def printHeader():
    print '%4s' % 'nvtx',
    print '%12s' % 'abs.eta.max',
    print '%5s' % 'n',
    print '%10s' % 'met.x',
    print '%13s' % 'met.x.var',
    print '%10s' % 'met.y',
    print '%13s' % 'met.y.var',
    print 

##____________________________________________________________________________||
def printCounts(counts, stats):
    keys = counts.keys()
    keys.sort()
    for k in keys:
        n = counts[k]
        metx = stats[k]["metx"]
        metx2 = stats[k]["metx2"]
        metx_mean = metx/n
        metx_var = (metx2 - (metx**2)/n)/(n-1) if n > 1 else 0

        mety = stats[k]["mety"]
        mety2 = stats[k]["mety2"]
        mety_mean = mety/n
        mety_var = (mety2 - (mety**2)/n)/(n-1) if n > 1 else 0

        print '%4d'    % k[0],
        print '%12.2f' % k[1],
        print '%5d'    % counts[k],
        print '%10.3f' % metx_mean,
        print '%13.3f' % metx_var,
        print '%10.3f' % mety_mean,
        print '%13.3f' % mety_var,
        print 

##____________________________________________________________________________||
def count(inputPathList):

    signal.signal(signal.SIGINT, handler)
    events = Events(inputPathList, maxEvents = options.nevents)

    handleVertex = Handle("std::vector<reco::Vertex>")
    handlePFMET = Handle("std::vector<reco::PFMET>")
    handlePFCandidate = Handle("std::vector<reco::PFCandidate>") 

    eventSelector = EventSelector()
    lumiBlockCertification = LuminosityBlockCertification(options.jsonPath)
    printFileNameIfChange = PrintFileNameIfChange(options.verbose)

    counts = { }
    stats = { }

    for event in events:

        if GLOBAL_LAST: break

        printFileNameIfChange(event)

        if not lumiBlockCertification(event): continue
        
        if not eventSelector(event): continue

        event.getByLabel(("primaryVertexFilter", "", ""), handleVertex)
        vertices = handleVertex.product()
        nvtx = len(vertices)

        event.getByLabel(("particleFlow", "", ""), handlePFCandidate)
        pfcands = handlePFCandidate.product()

        p4list = [createLorentzVector(c) for c in pfcands]

        for eta_cut in absEta_cuts:
            p4list_eta = [p for p in p4list if abs(p.Eta()) < eta_cut]
            sump4 = sum(p4list_eta, ROOT.TLorentzVector())
            met = -sump4

            key = (nvtx, eta_cut)
            if key not in counts:
                counts[key] = 0
                stats[key] = {"metx": 0.0, "metx2": 0.0, "mety": 0.0, "mety2": 0.0 }
            counts[key] += 1
            stats[key]["metx"] += met.Px()
            stats[key]["metx2"] += met.Px()**2
            stats[key]["mety"] += met.Py()
            stats[key]["mety2"] += met.Py()**2

    return counts, stats

##____________________________________________________________________________||
class EventSelector(object):
    def __init__(self):
        self.handleTriggerPaths = Handle("std::vector<pat::TriggerPath>") 
    def __call__(self, event):

        run = event.eventAuxiliary().run()

        event.getByLabel(('patTrigger', '', ''), self.handleTriggerPaths)
        triggerPaths = self.handleTriggerPaths.product()

        triggerPaths_ = [e for e in triggerPaths if re.match('^HLT_ZeroBias($|_v.*$)$', e.name())]

        if not any([e.wasAccept() for e in triggerPaths_]): return False

        return True

##____________________________________________________________________________||
def createLorentzVector(cand):
    return ROOT.TLorentzVector(cand.px(), cand.py(), cand.pz(), cand.energy())

##____________________________________________________________________________||
def DeltaR(cand1, cand2):
    return DeltaR_etaPhi(cand1.eta(), cand1.phi(), cand2.eta(), cand2.phi())

##____________________________________________________________________________||
def DeltaR_etaPhi(eta1, phi1, eta2, phi2):
    deltaEta = eta1 - eta2
    deltaPhi = abs(normalizedPhi(phi1 - phi2))
    return math.sqrt(deltaEta**2 + deltaPhi**2)

##____________________________________________________________________________||
def normalizedPhi(phi):
    while phi < - math.pi: phi += 2*math.pi;
    while phi > math.pi: phi -= 2*math.pi;
    return phi;

##____________________________________________________________________________||
def handler( signum, frame ):
    global GLOBAL_LAST
    GLOBAL_LAST = True

##____________________________________________________________________________||
loadLibraries()
from DataFormats.FWLite import Events, Handle

##____________________________________________________________________________||
if __name__ == '__main__':
    main()