test/miniAOD/example_fwlite_checkJets.py

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