PatUtils/src/LeptonJetIsolationAngle.cc

0001 #include "PhysicsTools/PatUtils/interface/LeptonJetIsolationAngle.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include "DataFormats/Math/interface/deltaR.h"
0006
0007 #include <vector>
0008
0009 using namespace pat;
0010
0011 // constructor
0012 LeptonJetIsolationAngle::LeptonJetIsolationAngle(edm::ConsumesCollector&& iC)
0013     : jetToken_(iC.consumes<reco::CaloJetCollection>(edm::InputTag("iterativeCone5CaloJets"))),
0014       electronsToken_(iC.consumes<std::vector<reco::GsfElectron> >(edm::InputTag("pixelMatchGsfElectrons"))) {}
0015
0016 // destructor
0017 LeptonJetIsolationAngle::~LeptonJetIsolationAngle() {}
0018
0019 // calculate the JetIsoA for the lepton object
0020 float LeptonJetIsolationAngle::calculate(const Electron& theElectron,
0021                                          const edm::Handle<edm::View<reco::Track> >& trackHandle,
0022                                          const edm::Event& iEvent) {
0023   CLHEP::HepLorentzVector theElectronHLV(theElectron.px(), theElectron.py(), theElectron.pz(), theElectron.energy());
0024   return this->calculate(theElectronHLV, trackHandle, iEvent);
0025 }
0026 float LeptonJetIsolationAngle::calculate(const Muon& theMuon,
0027                                          const edm::Handle<edm::View<reco::Track> >& trackHandle,
0028                                          const edm::Event& iEvent) {
0029   CLHEP::HepLorentzVector theMuonHLV(theMuon.px(), theMuon.py(), theMuon.pz(), theMuon.energy());
0030   return this->calculate(theMuonHLV, trackHandle, iEvent);
0031 }
0032
0033 // calculate the JetIsoA for the lepton's HLV
0034 float LeptonJetIsolationAngle::calculate(const CLHEP::HepLorentzVector& aLepton,
0035                                          const edm::Handle<edm::View<reco::Track> >& trackHandle,
0036                                          const edm::Event& iEvent) {
0037   // FIXME: this is an ugly temporary workaround, JetMET+egamma should come up with a better tool
0038   // retrieve the jets
0039   edm::Handle<reco::CaloJetCollection> jetHandle;
0040   iEvent.getByToken(jetToken_, jetHandle);
0041   reco::CaloJetCollection jetColl = *(jetHandle.product());
0042   // retrieve the electrons which might be in the jet list
0043   edm::Handle<std::vector<reco::GsfElectron> > electronsHandle;
0044   iEvent.getByToken(electronsToken_, electronsHandle);
0045   std::vector<reco::GsfElectron> electrons = *electronsHandle;
0046   // determine the set of isolated electrons
0047   std::vector<Electron> isoElectrons;
0048   for (size_t ie = 0; ie < electrons.size(); ie++) {
0049     Electron anElectron(electrons[ie]);
0050     if (anElectron.pt() > 10 && trkIsolator_.calculate(anElectron, *trackHandle) < 3.0) {
0051       isoElectrons.push_back(electrons[ie]);
0052     }
0053   }
0054   // determine the collections of jets, cleaned from electrons
0055   std::vector<reco::CaloJet> theJets;
0056   for (reco::CaloJetCollection::const_iterator itJet = jetColl.begin(); itJet != jetColl.end(); itJet++) {
0057     float mindr2 = 9999.;
0058     for (size_t ie = 0; ie < isoElectrons.size(); ie++) {
0059       float dr2 = ::deltaR2(*itJet, isoElectrons[ie]);
0060       if (dr2 < mindr2)
0061         mindr2 = dr2;
0062     }
0063     float mindr = std::sqrt(mindr2);
0064     // yes, all cuts hardcoded buts, but it's a second-order effect
0065     if (itJet->et() > 15 && mindr > 0.3)
0066       theJets.push_back(reco::CaloJet(*itJet));
0067   }
0068   // calculate finally the isolation angle
0069   float isoAngle = 1000;  // default to some craze impossible number to inhibit compiler warnings
0070   for (std::vector<reco::CaloJet>::const_iterator itJet = theJets.begin(); itJet != theJets.end(); itJet++) {
0071     float curDR = this->spaceAngle(aLepton, *itJet);
0072     if (curDR < isoAngle)
0073       isoAngle = curDR;
0074   }
0075   return isoAngle;
0076 }
0077
0078 // calculate the angle between two vectors in 3d eucledian space
0079 float LeptonJetIsolationAngle::spaceAngle(const CLHEP::HepLorentzVector& aLepton, const reco::CaloJet& aJet) {
0080   return acos(sin(aJet.theta()) * cos(aJet.phi()) * sin(aLepton.theta()) * cos(aLepton.phi()) +
0081               sin(aJet.theta()) * sin(aJet.phi()) * sin(aLepton.theta()) * sin(aLepton.phi()) +
0082               cos(aJet.theta()) * cos(aLepton.theta()));
0083 }