EgammaIsolationAlgos/src/PFBlockBasedIsolation.cc

0001 #include "RecoEgamma/EgammaIsolationAlgos/interface/PFBlockBasedIsolation.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
0004 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0005 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0006 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
0007
0008 #include <cmath>
0009
0010 //--------------------------------------------------------------------------------------------------
0011
0012 PFBlockBasedIsolation::PFBlockBasedIsolation() {
0013   // Default Constructor.
0014 }
0015
0016 //--------------------------------------------------------------------------------------------------
0017 PFBlockBasedIsolation::~PFBlockBasedIsolation() {}
0018
0019 void PFBlockBasedIsolation::setup(const edm::ParameterSet& conf) { coneSize_ = conf.getParameter<double>("coneSize"); }
0020
0021 std::vector<reco::PFCandidateRef> PFBlockBasedIsolation::calculate(
0022     math::XYZTLorentzVectorD p4,
0023     const reco::PFCandidateRef pfEGCand,
0024     const edm::Handle<reco::PFCandidateCollection> pfCandidateHandle) {
0025   std::vector<reco::PFCandidateRef> myVec;
0026
0027   math::XYZVector candidateMomentum(p4.px(), p4.py(), p4.pz());
0028   math::XYZVector candidateDirection = candidateMomentum.Unit();
0029
0030   const reco::PFCandidate::ElementsInBlocks& theElementsInpfEGcand = (*pfEGCand).elementsInBlocks();
0031   reco::PFCandidate::ElementsInBlocks::const_iterator ieg = theElementsInpfEGcand.begin();
0032   const reco::PFBlockRef egblock = ieg->first;
0033
0034   unsigned nObj = pfCandidateHandle->size();
0035   for (unsigned int lCand = 0; lCand < nObj; lCand++) {
0036     reco::PFCandidateRef pfCandRef(reco::PFCandidateRef(pfCandidateHandle, lCand));
0037
0038     float dR = 0.0;
0039     if (coneSize_ < 10.0) {
0040       dR = deltaR(candidateDirection.Eta(), candidateDirection.Phi(), pfCandRef->eta(), pfCandRef->phi());
0041       if (dR > coneSize_)
0042         continue;
0043     }
0044
0045     const reco::PFCandidate::ElementsInBlocks& theElementsInPFcand = pfCandRef->elementsInBlocks();
0046
0047     bool elementFound = false;
0048     for (reco::PFCandidate::ElementsInBlocks::const_iterator ipf = theElementsInPFcand.begin();
0049          ipf < theElementsInPFcand.end();
0050          ++ipf) {
0051       if (ipf->first == egblock && !elementFound) {
0052         for (ieg = theElementsInpfEGcand.begin(); ieg < theElementsInpfEGcand.end(); ++ieg) {
0053           if (ipf->second == ieg->second && !elementFound) {
0054             if (elementPassesCleaning(pfCandRef, pfEGCand)) {
0055               myVec.push_back(pfCandRef);
0056               elementFound = true;
0057             }
0058           }
0059         }
0060       }
0061     }
0062   }
0063
0064   return myVec;
0065 }
0066
0067 bool PFBlockBasedIsolation::elementPassesCleaning(const reco::PFCandidateRef& pfCand,
0068                                                   const reco::PFCandidateRef& pfEGCand) {
0069   if (pfCand->particleId() == reco::PFCandidate::h)
0070     return passesCleaningChargedHadron(pfCand, pfEGCand);
0071   else if (pfCand->particleId() == reco::PFCandidate::h0)
0072     return passesCleaningNeutralHadron(pfCand, pfEGCand);
0073   else if (pfCand->particleId() == reco::PFCandidate::gamma)
0074     return passesCleaningPhoton(pfCand, pfEGCand);
0075   else
0076     return true;  //doesnt really matter here as if its not a photon,neutral or charged it wont be included in isolation
0077 }
0078
0079 //currently the record of which candidates came from the charged hadron is acceptable, no further cleaning is needed
0080 bool PFBlockBasedIsolation::passesCleaningChargedHadron(const reco::PFCandidateRef& pfCand,
0081                                                         const reco::PFCandidateRef& pfEGCand) {
0082   return true;
0083 }
0084
0085 //neutral hadrons are not part of the PF E/gamma reco, therefore they cant currently come from an electron/photon and so should be rejected
0086 //but we still think there may be some useful info here and given we can easily
0087 //fix this at AOD level, we will auto accept them for now and clean later
0088 bool PFBlockBasedIsolation::passesCleaningNeutralHadron(const reco::PFCandidateRef& pfCand,
0089                                                         const reco::PFCandidateRef& pfEGCand) {
0090   return true;
0091 }
0092
0093 //the highest et ECAL element of the photon must match to the electron superclusters or one of its sub clusters
0094 bool PFBlockBasedIsolation::passesCleaningPhoton(const reco::PFCandidateRef& pfCand,
0095                                                  const reco::PFCandidateRef& pfEGCand) {
0096   bool passesCleaning = false;
0097   const reco::PFBlockElementCluster* ecalClusWithMaxEt = getHighestEtECALCluster(*pfCand);
0098   if (ecalClusWithMaxEt) {
0099     if (ecalClusWithMaxEt->superClusterRef().isNonnull() &&
0100         ecalClusWithMaxEt->superClusterRef()->seed()->seed() ==
0101             pfEGCand->superClusterRef()
0102                 ->seed()
0103                 ->seed()) {  //being sure to match, some concerned about different collections, shouldnt be but to be safe
0104       passesCleaning = true;
0105     } else {
0106       for (auto cluster : pfEGCand->superClusterRef()->clusters()) {
0107         //the PF clusters there are in two different collections so cant reference match
0108         //but we can match on the seed id, no clusters can share a seed so if the seeds are
0109         //equal, it must be the same cluster
0110         if (ecalClusWithMaxEt->clusterRef()->seed() == cluster->seed()) {
0111           passesCleaning = true;
0112         }
0113       }  //end of loop over clusters
0114     }
0115   }  //end of null check for highest ecal cluster
0116   return passesCleaning;
0117 }
0118
0119 const reco::PFBlockElementCluster* PFBlockBasedIsolation::getHighestEtECALCluster(const reco::PFCandidate& pfCand) {
0120   float maxECALEt = -1;
0121   const reco::PFBlockElement* maxEtECALCluster = nullptr;
0122   const reco::PFCandidate::ElementsInBlocks& elementsInPFCand = pfCand.elementsInBlocks();
0123   for (auto& elemIndx : elementsInPFCand) {
0124     const reco::PFBlockElement* elem =
0125         elemIndx.second < elemIndx.first->elements().size() ? &elemIndx.first->elements()[elemIndx.second] : nullptr;
0126     if (elem && elem->type() == reco::PFBlockElement::ECAL && elem->clusterRef()->pt() > maxECALEt) {
0127       maxECALEt = elem->clusterRef()->pt();
0128       maxEtECALCluster = elem;
0129     }
0130   }
0131   return dynamic_cast<const reco::PFBlockElementCluster*>(maxEtECALCluster);
0132 }