test/miniAOD/example_fwlite_checkJets.py

0001 #! /usr/bin/env python3
0002
0003 from __future__ import print_function
0004 import ROOT
0005 import sys
0006 from DataFormats.FWLite import Events, Handle
0007 from math import *
0008
0009 def deltaPhi(a,b) :
0010  r = a-b
0011  while r>2*pi : r-=2*pi
0012  while r<-2*pi : r+=2*pi
0013  return r
0014
0015 def deltaR(a,b) :
0016   dphi=deltaPhi(a.phi(),b.phi())
0017   return sqrt(dphi*dphi+(a.eta()-b.eta())**2)
0018
0019 events = Events (['test.root'])
0020
0021 handleGJ1  = Handle ("std::vector<pat::Jet>")
0022 handleGJ2  = Handle ("std::vector<pat::Jet>")
0023 handleGP  = Handle ("std::vector<reco::GenParticle>")
0024
0025 # for now, label is just a tuple of strings that is initialized just
0026 # like and edm::InputTag
0027 labelGJ1 = ("slimmedJets","","PAT")
0028 labelGJ2 = ("patJetsAK5PFCHS","","S2")
0029 labelGP = ("prunedGenParticles")
0030
0031 ROOT.gROOT.SetBatch()        # don't pop up canvases
0032 ROOT.gROOT.SetStyle('Plain') # white background
0033 zptHist = ROOT.TH1F ("zpt", "Z Pt", 50, 0, 500)
0034
0035 # loop over events
0036 count= 0
0037 for event in events:
0038     count+=1
0039     if count % 1000 == 0 :
0040     print(count)
0041     event.getByLabel (labelGJ1, handleGJ1)
0042     event.getByLabel (labelGJ2, handleGJ2)
0043     event.getByLabel (labelGP, handleGP)
0044     # get the product
0045     jets1 = handleGJ1.product()
0046     jets2 = handleGJ2.product()
0047
0048     for j1,j2 in zip(jets1,jets2)  :
0049             if abs(j1.eta()) < 2.5 and j1.pt() > 20 and j1.chargedHadronEnergyFraction() > 0.05 :
0050         if abs(j1.pt()-j2.pt())/(j1.pt()+j2.pt()) >0.05 :
0051             print("Mismatch at record ", count)
0052             print("Bad match is : pt %s vs %s, b-tag %s vs %s, MC flavour %s vs %s " %(j1.pt(),j2.pt(),j1.bDiscriminator("combinedSecondaryVertexBJetTags"),j2.bDiscriminator("combinedSecondaryVertexBJetTags"),j1.partonFlavour(),j2.partonFlavour()))
0053             print("Jet eta and phi" ,j1.eta(),j1.phi())
0054             print(" alljets ")
0055                 for j11,j22 in zip(jets1,jets2)  :
0056                             print("  Jet pt %s vs %s, b-tag %s vs %s, MC flavour %s vs %s " % (j11.pt(),j22.pt(),j11.bDiscriminator("combinedSecondaryVertexBJetTags"),j22.bDiscriminator("combinedSecondaryVertexBJetTags"),j11.partonFlavour(),j22.partonFlavour()))
0057     #       print "gen parts"
0058     #       genparts = handleGP.product()
0059     #       for gp in genparts :
0060     #           if abs(gp.eta()-j1.eta()) < 0.3 :
0061     #               print deltaR(j1.p4(),gp.p4()),gp.pdgId(),gp.status(),gp.pt(),gp.eta(),gp.phi()
0062
0063