test/miniAOD/example_packedCandidates.py

0001 from __future__ import print_function
0002 # import ROOT in batch mode
0003 from builtins import range
0004 import sys
0005 oldargv = sys.argv[:]
0006 sys.argv = [ '-b-' ]
0007 import ROOT
0008 ROOT.gROOT.SetBatch(True)
0009 sys.argv = oldargv
0010
0011 # load FWLite C++ libraries
0012 ROOT.gSystem.Load("libFWCoreFWLite.so");
0013 ROOT.gSystem.Load("libDataFormatsFWLite.so");
0014 ROOT.FWLiteEnabler.enable()
0015
0016 # load FWlite python libraries
0017 from DataFormats.FWLite import Handle, Events
0018
0019 # define deltaR
0020 from math import hypot, pi
0021 def deltaR(a,b):
0022     dphi = abs(a.phi()-b.phi());
0023     if dphi < pi: dphi = 2*pi-dphi
0024     return hypot(a.eta()-b.eta(),dphi)
0025
0026 muons, muonLabel = Handle("std::vector<pat::Muon>"), "slimmedMuons"
0027 electrons, electronLabel = Handle("std::vector<pat::Electron>"), "slimmedElectrons"
0028 jets, jetLabel = Handle("std::vector<pat::Jet>"), "slimmedJets"
0029 pfs, pfLabel = Handle("std::vector<pat::PackedCandidate>"), "packedPFCandidates"
0030
0031 # open file (you can use 'edmFileUtil -d /store/whatever.root' to get the physical file name)
0032 events = Events("patMiniAOD_standard.root")
0033
0034 for iev,event in enumerate(events):
0035     event.getByLabel(muonLabel, muons)
0036     event.getByLabel(electronLabel, electrons)
0037     event.getByLabel(pfLabel, pfs)
0038     event.getByLabel(jetLabel, jets)
0039
0040     print("\nEvent %d: run %6d, lumi %4d, event %12d" % (iev,event.eventAuxiliary().run(), event.eventAuxiliary().luminosityBlock(),event.eventAuxiliary().event()))
0041
0042     # Let's compute lepton PF Isolation with R=0.2, 0.5 GeV threshold on neutrals, and deltaBeta corrections
0043     leps  = [ p for p in muons.product() ] + [ p for p in electrons.product() ]
0044     for lep in leps:
0045         # skip those below 5 GeV, which we don't care about
0046         if lep.pt() < 5: continue
0047         # initialize sums
0048         charged = 0
0049         neutral = 0
0050         pileup  = 0
0051         # now get a list of the PF candidates used to build this lepton, so to exclude them
0052         footprint = set()
0053         for i in range(lep.numberOfSourceCandidatePtrs()):
0054             footprint.add(lep.sourceCandidatePtr(i).key()) # the key is the index in the pf collection
0055         # now loop on pf candidates
0056         for ipf,pf in enumerate(pfs.product()):
0057             if deltaR(pf,lep) < 0.2:
0058                 # pfcandidate-based footprint removal
0059                 if ipf in footprint: continue
0060                 # add up
0061                 if (pf.charge() == 0):
0062                     if pf.pt() > 0.5: neutral += pf.pt()
0063                 elif pf.fromPV() >= 2:
0064                     charged += pf.pt()
0065                 else:
0066                     if pf.pt() > 0.5: pileup += pf.pt()
0067         # do deltaBeta
0068         iso = charged + max(0, neutral-0.5*pileup)
0069         print("%-8s of pt %6.1f, eta %+4.2f: relIso = %5.2f" % (
0070                     "muon" if abs(lep.pdgId())==13 else "electron",
0071                     lep.pt(), lep.eta(), iso/lep.pt()))
0072
0073     # Let's compute the fraction of charged pt from particles with dz < 0.1 cm
0074     for i,j in enumerate(jets.product()):
0075         if j.pt() < 40 or abs(j.eta()) > 2.4: continue
0076         sums = [0,0]
0077         for id in range(j.numberOfDaughters()):
0078             dau = j.daughter(id)
0079             if (dau.charge() == 0): continue
0080             sums[ abs(dau.dz())<0.1 ] += dau.pt()
0081         sum = sums[0]+sums[1]
0082         print("Jet with pt %6.1f, eta %+4.2f, beta(0.1) = %+5.3f, pileup mva disc %+.2f" % (
0083                 j.pt(),j.eta(), sums[1]/sum if sum else 0, j.userFloat("pileupJetId:fullDiscriminant")))
0084
0085     # Let's check the calorimeter response for hadrons (after PF hadron calibration)
0086     for i,j in enumerate(pfs.product()):
0087         if not j.isIsolatedChargedHadron(): continue
0088         print("Isolated charged hadron candidate with pt %6.1f, eta %+4.2f, calo/track energy = %+5.3f, hcal/calo energy %+5.3f" % (
0089                 j.pt(),j.eta(), j.rawCaloFraction(), j.hcalFraction()))
0090
0091     if iev > 10: break