PhotonIDValueMapProducer.cc

Go to the documentation of this file.

#include "DataFormats/Common/interface/ValueMap.h"
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "DataFormats/PatCandidates/interface/Photon.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
#include "RecoEgamma/EgammaTools/interface/MultiToken.h"
#include "RecoEgamma/EgammaTools/interface/Utils.h"
#include "FWCore/Utilities/interface/isFinite.h"

namespace {

  // This template function finds whether theCandidate is in thefootprint
  // collection. It is templated to be able to handle both reco and pat
  // photons (from AOD and miniAOD, respectively).
  template <class T>
  bool isInFootprint(const T& footprint,
                     const edm::Ptr<reco::Candidate>& candidate)
  {
      for (auto& it : footprint) {
          if (it.key() == candidate.key())
              return true;
      }
      return false;
  }

  struct CachingPtrCandidate {
    CachingPtrCandidate(const reco::Candidate* cPtr, bool isAOD) : 
      candidate(cPtr), 
      track(isAOD? &*static_cast<const reco::PFCandidate* >(cPtr)->trackRef() : nullptr),
      packed(isAOD? nullptr : static_cast<const pat::PackedCandidate*>(cPtr))
    {}
    
    const reco::Candidate* candidate;
    const reco::Track* track;
    const pat::PackedCandidate* packed;
  };

  void getImpactParameters(const CachingPtrCandidate& candidate, const reco::Vertex& pv, float& dxy, float& dz)
  {
    if (candidate.track != nullptr) {
      dxy = candidate.track->dxy(pv.position());
      dz = candidate.track->dz(pv.position());
    } else {
      dxy = candidate.packed->dxy(pv.position());
      dz = candidate.packed->dz(pv.position());
    }
  }

  // Some helper functions that are needed to access info in
  // AOD vs miniAOD
  reco::PFCandidate::ParticleType getCandidatePdgId(const reco::Candidate* candidate, bool isAOD)
  {
      if (isAOD)
        return static_cast<const reco::PFCandidate*>(candidate)->particleId();

      // the neutral hadrons and charged hadrons can be of pdgId types
      // only 130 (K0L) and +-211 (pi+-) in packed candidates
      const int pdgId = static_cast<const pat::PackedCandidate*>(candidate)->pdgId();
      if (pdgId == 22)
          return reco::PFCandidate::gamma;
      else if (abs(pdgId) == 130) // PDG ID for K0L
          return reco::PFCandidate::h0;
      else if (abs(pdgId) == 211) // PDG ID for pi+-
          return reco::PFCandidate::h;
      else
          return reco::PFCandidate::X;
  }

};

class PhotonIDValueMapProducer : public edm::global::EDProducer<> {

public:
    explicit PhotonIDValueMapProducer(const edm::ParameterSet&);
    ~PhotonIDValueMapProducer() override {}

    static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
    void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;

    // This function computes charged hadron isolation with respect to multiple
    // PVs and returns the worst of the found isolation values. The function
    // implements the computation method taken directly from Run 1 code of
    // H->gamma gamma, specifically from the class CiCPhotonID of the
    // HiggsTo2photons anaysis code. Template is introduced to handle reco/pat
    // photons and aod/miniAOD PF candidates collections
    float computeWorstPFChargedIsolation(
            const reco::Photon& photon,
            const edm::View<reco::Candidate>& pfCands,
            const reco::VertexCollection& vertices,
            const reco::Vertex& pv,
            unsigned char options, bool isAOD) const;

    // check whether a non-null preshower is there
    const bool usesES_;

    // Dual Tokens for AOD and MiniAOD case
    const MultiTokenT<edm::View<reco::Photon>>    src_;
    const MultiTokenT<EcalRecHitCollection>       ebRecHits_;
    const MultiTokenT<EcalRecHitCollection>       eeRecHits_;
    const MultiTokenT<EcalRecHitCollection>       esRecHits_;
    const MultiTokenT<reco::VertexCollection>     vtxToken_;
    const MultiTokenT<edm::View<reco::Candidate>> pfCandsToken_;
    const edm::EDGetToken                         particleBasedIsolationToken_;

};

constexpr int nVars_ = 19;

const std::string names[nVars_] = {
    // Cluster shapes
    "phoFull5x5SigmaIEtaIEta", // 0
    "phoFull5x5SigmaIEtaIPhi",
    "phoFull5x5E1x3",
    "phoFull5x5E2x2",
    "phoFull5x5E2x5Max",
    "phoFull5x5E5x5", // 5
    "phoESEffSigmaRR",
    // Cluster shape ratios
    "phoFull5x5E1x3byE5x5",
    "phoFull5x5E2x2byE5x5",
    "phoFull5x5E2x5byE5x5",
    // Isolations
    "phoChargedIsolation", // 10
    "phoNeutralHadronIsolation",
    "phoPhotonIsolation",
    "phoWorstChargedIsolation",
    "phoWorstChargedIsolationConeVeto",
    "phoWorstChargedIsolationConeVetoPVConstr", // 15
    // PFCluster Isolation
    "phoTrkIsolation",
    "phoHcalPFClIsolation",
    "phoEcalPFClIsolation"
};

// options and bitflags
constexpr float coneSizeDR2 = 0.3*0.3;
constexpr float dxyMax = 0.1;
constexpr float dzMax = 0.2;
constexpr float dRveto2Barrel = 0.02*0.02;
constexpr float dRveto2Endcap = 0.02*0.02;
constexpr float ptMin = 0.1;

const unsigned char PV_CONSTRAINT  = 0x1;
const unsigned char DR_VETO = 0x2;
const unsigned char PT_MIN_THRESH = 0x8;

PhotonIDValueMapProducer::PhotonIDValueMapProducer(const edm::ParameterSet& cfg)
    : usesES_(!cfg.getParameter<edm::InputTag>("esReducedRecHitCollection").label().empty()
          || !cfg.getParameter<edm::InputTag>("esReducedRecHitCollectionMiniAOD").label().empty())
    , src_(               consumesCollector(), cfg, "src", "srcMiniAOD")
    , ebRecHits_   (src_, consumesCollector(), cfg, "ebReducedRecHitCollection", "ebReducedRecHitCollectionMiniAOD")
    , eeRecHits_   (src_, consumesCollector(), cfg, "eeReducedRecHitCollection", "eeReducedRecHitCollectionMiniAOD")
    , esRecHits_   (src_, consumesCollector(), cfg, "esReducedRecHitCollection", "esReducedRecHitCollectionMiniAOD")
    , vtxToken_    (src_, consumesCollector(), cfg, "vertices", "verticesMiniAOD")
    , pfCandsToken_(src_, consumesCollector(), cfg, "pfCandidates", "pfCandidatesMiniAOD")
    , particleBasedIsolationToken_(mayConsume<edm::ValueMap<std::vector<reco::PFCandidateRef>>>(
          cfg.getParameter<edm::InputTag>("particleBasedIsolation")) /* ...only for AOD... */ )
{

    // Declare producibles
    for (int i = 0; i < nVars_; ++i)
        produces<edm::ValueMap<float>>(names[i]);
}

void PhotonIDValueMapProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const
{
    // Get the handles
    auto src           = src_.getValidHandle(iEvent);
    auto vertices      = vtxToken_.getValidHandle(iEvent);
    auto pfCandsHandle = pfCandsToken_.getValidHandle(iEvent);

    bool isAOD = src_.getGoodTokenIndex() == 0;
    edm::Handle<edm::ValueMap<std::vector<reco::PFCandidateRef>>> particleBasedIsolationMap;
    if (isAOD) { // this exists only in AOD
        iEvent.getByToken(particleBasedIsolationToken_, particleBasedIsolationMap);
    } else if (!src->empty()) {
        edm::Ptr<pat::Photon> test(src->ptrAt(0));
        if (test.isNull() || !test.isAvailable()) {
            throw cms::Exception("InvalidConfiguration")
                << "DataFormat is detected as miniAOD but cannot cast to pat::Photon!";
        }
    }

    // Configure Lazy Tools, which will compute 5x5 quantities
    auto lazyToolnoZS = usesES_ ? noZS::EcalClusterLazyTools(iEvent, iSetup,
                                        ebRecHits_.get(iEvent), eeRecHits_.get(iEvent), esRecHits_.get(iEvent))
                                : noZS::EcalClusterLazyTools(iEvent, iSetup,
                                        ebRecHits_.get(iEvent), eeRecHits_.get(iEvent));

    // Get PV
    if (vertices->empty())
        return; // skip the event if no PV found
    const reco::Vertex& pv = vertices->front();

    std::vector<float> vars[nVars_];

    // reco::Photon::superCluster() is virtual so we can exploit polymorphism
    for(auto const& iPho : src->ptrs())
    {
        //
        // Compute full 5x5 quantities
        //
        const auto& seed = *(iPho->superCluster()->seed());

        // For full5x5_sigmaIetaIeta, for 720 we use: lazy tools for AOD,
        // and userFloats or lazy tools for miniAOD. From some point in 72X and on, one can
        // retrieve the full5x5 directly from the object with ->full5x5_sigmaIetaIeta()
        // for both formats.
        std::vector<float> vCov = lazyToolnoZS.localCovariances(seed);
        vars[0].push_back(edm::isNotFinite(vCov[0]) ? 0. : sqrt(vCov[0]));
        vars[1].push_back(vCov[1]);
        vars[2].push_back(lazyToolnoZS.e1x3(seed));
        vars[3].push_back(lazyToolnoZS.e2x2(seed));
        vars[4].push_back(lazyToolnoZS.e2x5Max(seed));
        vars[5].push_back(lazyToolnoZS.e5x5(seed));
        vars[6].push_back(lazyToolnoZS.eseffsirir(*(iPho->superCluster())));
        vars[7].push_back(vars[2].back() / vars[5].back());
        vars[8].push_back(vars[3].back() / vars[5].back());
        vars[9].push_back(vars[4].back() / vars[5].back());

        //
        // Compute absolute uncorrected isolations with footprint removal
        //

        // First, find photon direction with respect to the good PV
        math::XYZVector phoWrtVtx(
            iPho->superCluster()->x() - pv.x(), iPho->superCluster()->y() - pv.y(), iPho->superCluster()->z() - pv.z());

        // isolation sums
        float chargedIsoSum = 0.;
        float neutralHadronIsoSum = 0.;
        float photonIsoSum = 0.;

        // Loop over all PF candidates
        for(auto const& iCand : pfCandsHandle->ptrs())
        {
            // Here, the type will be a simple reco::Candidate. We cast it
            // for full PFCandidate or PackedCandidate below as necessary

            // One would think that we should check that this iCand from the
            // generic PF collection is not identical to the iPho photon for
            // which we are computing the isolations. However, it turns out to
            // be unnecessary. Below, in the function isInFootprint(), we drop
            // this iCand if it is in the footprint, and this always removes
            // the iCand if it matches the iPho. The explicit check at this
            // point is not totally trivial because of non-triviality of
            // implementation of this check for miniAOD (PackedCandidates of
            // the PF collection do not contain the supercluser link, so can't
            // use that).
            //
            // if( isAOD ) {
            //     if( ((const edm::Ptr<reco::PFCandidate>)iCand)->superClusterRef() == iPho->superCluster() )
            //     continue;
            // }

            // Check if this candidate is within the isolation cone
            float dR2 = deltaR2(phoWrtVtx.Eta(), phoWrtVtx.Phi(), iCand->eta(), iCand->phi());
            if (dR2 > coneSizeDR2)
                continue;

            // Check if this candidate is not in the footprint
            if (isAOD) {
                if(isInFootprint((*particleBasedIsolationMap)[iPho], iCand))
                    continue;
            } else {
                edm::Ptr<pat::Photon> patPhotonPtr(iPho);
                if(isInFootprint(patPhotonPtr->associatedPackedPFCandidates(), iCand))
                    continue;
            }

            // Find candidate type
            reco::PFCandidate::ParticleType thisCandidateType = getCandidatePdgId(&*iCand, isAOD);

            // Increment the appropriate isolation sum
            if (thisCandidateType == reco::PFCandidate::h) {
                // for charged hadrons, additionally check consistency
                // with the PV
                float dxy = -999;
                float dz = -999;

                getImpactParameters(CachingPtrCandidate(&*iCand, isAOD), pv, dxy, dz);

                if (fabs(dxy) > dxyMax || fabs(dz) > dzMax)
                    continue;

                // The candidate is eligible, increment the isolaiton
                chargedIsoSum += iCand->pt();
            }

            if (thisCandidateType == reco::PFCandidate::h0)
                neutralHadronIsoSum += iCand->pt();

            if (thisCandidateType == reco::PFCandidate::gamma)
                photonIsoSum += iCand->pt();
        }

        vars[10].push_back(chargedIsoSum);
        vars[11].push_back(neutralHadronIsoSum);
        vars[12].push_back(photonIsoSum);

        // Worst isolation computed with no vetos or ptMin cut, as in Run 1 Hgg code.
        unsigned char options = 0;
        vars[13].push_back(computeWorstPFChargedIsolation(*iPho, *pfCandsHandle, *vertices, pv, options, isAOD));

        // Worst isolation computed with cone vetos and a ptMin cut, as in Run 2 Hgg code.
        options |= PT_MIN_THRESH | DR_VETO;
        vars[14].push_back(computeWorstPFChargedIsolation(*iPho, *pfCandsHandle, *vertices, pv, options, isAOD));

        // Like before, but adding primary vertex constraint
        options |= PV_CONSTRAINT;
        vars[15].push_back(computeWorstPFChargedIsolation(*iPho, *pfCandsHandle, *vertices, pv, options, isAOD));

        // PFCluster Isolations
        vars[16].push_back(iPho->trkSumPtSolidConeDR04());
        if (isAOD) {
            vars[17].push_back(0.f);
            vars[18].push_back(0.f);
        } else {
            edm::Ptr<pat::Photon> patPhotonPtr{ iPho };
            vars[17].push_back(patPhotonPtr->hcalPFClusterIso());
            vars[18].push_back(patPhotonPtr->ecalPFClusterIso());
        }
    }

    // write the value maps
    for (int i = 0; i < nVars_; ++i) {
        writeValueMap(iEvent, src, vars[i], names[i]);
    }
}

void PhotonIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  // photonIDValueMapProducer
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("particleBasedIsolation",           edm::InputTag("particleBasedIsolation","gedPhotons"));
  desc.add<edm::InputTag>("src",                              edm::InputTag("gedPhotons"));
  desc.add<edm::InputTag>("srcMiniAOD",                       edm::InputTag("slimmedPhotons","","@skipCurrentProcess"));
  desc.add<edm::InputTag>("esReducedRecHitCollectionMiniAOD", edm::InputTag("reducedEgamma","reducedESRecHits"));
  desc.add<edm::InputTag>("eeReducedRecHitCollection",        edm::InputTag("reducedEcalRecHitsEE"));
  desc.add<edm::InputTag>("pfCandidates",                     edm::InputTag("particleFlow"));
  desc.add<edm::InputTag>("vertices",                         edm::InputTag("offlinePrimaryVertices"));
  desc.add<edm::InputTag>("ebReducedRecHitCollectionMiniAOD", edm::InputTag("reducedEgamma","reducedEBRecHits"));
  desc.add<edm::InputTag>("eeReducedRecHitCollectionMiniAOD", edm::InputTag("reducedEgamma","reducedEERecHits"));
  desc.add<edm::InputTag>("esReducedRecHitCollection",        edm::InputTag("reducedEcalRecHitsES"));
  desc.add<edm::InputTag>("pfCandidatesMiniAOD",              edm::InputTag("packedPFCandidates"));
  desc.add<edm::InputTag>("verticesMiniAOD",                  edm::InputTag("offlineSlimmedPrimaryVertices"));
  desc.add<edm::InputTag>("ebReducedRecHitCollection",        edm::InputTag("reducedEcalRecHitsEB"));
  descriptions.add("photonIDValueMapProducer", desc);
}


// Charged isolation with respect to the worst vertex. See more
// comments above at the function declaration.
float PhotonIDValueMapProducer::computeWorstPFChargedIsolation(const reco::Photon& photon, const edm::View<reco::Candidate>& pfCands,
    const reco::VertexCollection& vertices, const reco::Vertex& pv, unsigned char options, bool isAOD) const
{
    float worstIsolation = 0.0;

    const float dRveto2 = photon.isEB() ? dRveto2Barrel : dRveto2Endcap;

    std::vector<CachingPtrCandidate> chargedCands;
    chargedCands.reserve(pfCands.size());
    for (auto const& aCand : pfCands){
      
      // require that PFCandidate is a charged hadron
      reco::PFCandidate::ParticleType thisCandidateType = getCandidatePdgId(&aCand, isAOD);
      if (thisCandidateType != reco::PFCandidate::h)
        continue;
      
      if ((options & PT_MIN_THRESH) && aCand.pt() < ptMin)
        continue;

      chargedCands.emplace_back(&aCand, isAOD);
    }

    // Calculate isolation sum separately for each vertex
    for (unsigned int ivtx = 0; ivtx < vertices.size(); ++ivtx) {

        // Shift the photon according to the vertex
        const reco::VertexRef vtx(&vertices, ivtx);
        math::XYZVector phoWrtVtx(photon.superCluster()->x() - vtx->x(),
            photon.superCluster()->y() - vtx->y(), photon.superCluster()->z() - vtx->z());

    const float phoWrtVtxPhi = phoWrtVtx.phi();
    const float phoWrtVtxEta = phoWrtVtx.eta();

        float sum = 0;
        // Loop over the PFCandidates
        for (auto const& aCCand : chargedCands) {

            auto iCand = aCCand.candidate;
            float dR2 = deltaR2(phoWrtVtxEta, phoWrtVtxPhi, iCand->eta(), iCand->phi());
            if (dR2 > coneSizeDR2  ||
                    (options & DR_VETO && dR2 < dRveto2))
                continue;

        float dxy = -999;
            float dz = -999;
            if (options & PV_CONSTRAINT)
                getImpactParameters(aCCand, pv, dxy, dz);
            else
                getImpactParameters(aCCand, *vtx, dxy, dz);

            if (fabs(dxy) > dxyMax || fabs(dz) > dzMax)
                continue;


            sum += iCand->pt();
        }

        worstIsolation = std::max(sum, worstIsolation);
    }

    return worstIsolation;
}

DEFINE_FWK_MODULE(PhotonIDValueMapProducer);