test/miniAOD/example_packedCandidates.py

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