PFRecoTauClusterVariables.cc

Go to the documentation of this file.

#include "RecoTauTag/RecoTau/interface/PFRecoTauClusterVariables.h"

namespace {
  struct PFTau_traits {
    typedef reco::PFTau Tau_t;
    typedef const std::vector<reco::CandidatePtr>& Ret_t;
  };
  struct PATTau_traits {
    typedef pat::Tau Tau_t;
    typedef reco::CandidatePtrVector Ret_t;
  };

  template <typename T>
  typename T::Ret_t getGammas_T(const typename T::Tau_t& tau, bool signal) {
    return typename T::Ret_t();
  }
  template <>
  const std::vector<reco::CandidatePtr>& getGammas_T<PFTau_traits>(const reco::PFTau& tau, bool signal) {
    if (signal) {
      return tau.signalGammaCands();
    }
    return tau.isolationGammaCands();
  }

  template <>
  reco::CandidatePtrVector getGammas_T<PATTau_traits>(const pat::Tau& tau, bool signal) {
    if (signal) {
      return tau.signalGammaCands();
    }
    return tau.isolationGammaCands();
  }

  bool isInside(float photon_pt, float deta, float dphi) {
    constexpr double stripEtaAssociationDistance_0p95_p0 = 0.197077;
    constexpr double stripEtaAssociationDistance_0p95_p1 = 0.658701;
    constexpr double stripPhiAssociationDistance_0p95_p0 = 0.352476;
    constexpr double stripPhiAssociationDistance_0p95_p1 = 0.707716;
    if (photon_pt == 0) {
      return false;
    }
    if ((dphi < 0.3 && dphi < std::max(0.05,
                                       stripPhiAssociationDistance_0p95_p0 *
                                           std::pow(photon_pt, -stripPhiAssociationDistance_0p95_p1))) &&
        (deta < 0.15 && deta < std::max(0.05,
                                        stripEtaAssociationDistance_0p95_p0 *
                                            std::pow(photon_pt, -stripEtaAssociationDistance_0p95_p1)))) {
      return true;
    }
    return false;
  }
}  // namespace

namespace reco {
  namespace tau {
    float lead_track_chi2(const reco::PFTau& tau) {
      float LeadingTracknormalizedChi2 = 0;
      const reco::CandidatePtr& leadingCharged = tau.leadChargedHadrCand();
      if (leadingCharged.isNonnull()) {
        const reco::PFCandidate* pfcand = dynamic_cast<const reco::PFCandidate*>(leadingCharged.get());
        if (pfcand != nullptr) {
          reco::TrackRef tref = pfcand->trackRef();
          if (tref.isNonnull()) {
            LeadingTracknormalizedChi2 = tref->normalizedChi2();
          }
        } else {
          const pat::PackedCandidate* patcand = dynamic_cast<const pat::PackedCandidate*>(leadingCharged.get());
          if (patcand != nullptr && patcand->hasTrackDetails()) {
            LeadingTracknormalizedChi2 = patcand->pseudoTrack().normalizedChi2();
          }
        }
      }
      return LeadingTracknormalizedChi2;
    }
    float eratio(const reco::PFTau& tau) {
      float ecal_en_in_signal_pf_cands = 0;
      float hcal_en_in_signal_pf_cands = 0;
      for (const auto& signal_cand : tau.signalCands()) {
        const reco::PFCandidate* signal_pfcand = dynamic_cast<const reco::PFCandidate*>(signal_cand.get());
        if (signal_pfcand != nullptr) {
          ecal_en_in_signal_pf_cands += signal_pfcand->ecalEnergy();
          hcal_en_in_signal_pf_cands += signal_pfcand->hcalEnergy();
        }
        // TauReco@MiniAOD: recalculate for PackedCandidate if added to MiniAOD event content
      }
      float total = ecal_en_in_signal_pf_cands + hcal_en_in_signal_pf_cands;
      if (total == 0.) {
        return -1.;
      }
      return ecal_en_in_signal_pf_cands / total;
    }
    float eratio(const pat::Tau& tau) {
      float ecal_en_in_signal_cands = tau.ecalEnergy();
      float hcal_en_in_signal_cands = tau.hcalEnergy();
      float total = ecal_en_in_signal_cands + hcal_en_in_signal_cands;
      if (total == 0.) {
        return -1.;
      }
      return ecal_en_in_signal_cands / total;
    }
    template <typename T>
    float pt_weighted_dx_T(const typename T::Tau_t& tau, int mode, int var, int decaymode) {
      float sum_pt = 0.;
      float sum_dx_pt = 0.;
      float signalrad = std::max(0.05, std::min(0.1, 3. / std::max(1., tau.pt())));
      int is3prong = (decaymode == 10);
      const auto& cands = getGammas_T<T>(tau, mode < 2);
      for (const auto& cand : cands) {
        // only look at electrons/photons with pT > 0.5
        if (cand->pt() < 0.5) {
          continue;
        }
        float dr = reco::deltaR(cand->eta(), cand->phi(), tau.eta(), tau.phi());
        float deta = std::abs(cand->eta() - tau.eta());
        float dphi = std::abs(reco::deltaPhi(cand->phi(), tau.phi()));
        float pt = cand->pt();
        bool flag = isInside(pt, deta, dphi);
        if (is3prong == 0) {
          if (mode == 2 || (mode == 0 && dr < signalrad) || (mode == 1 && dr > signalrad)) {
            sum_pt += pt;
            if (var == 0)
              sum_dx_pt += pt * dr;
            else if (var == 1)
              sum_dx_pt += pt * deta;
            else if (var == 2)
              sum_dx_pt += pt * dphi;
          }
        } else if (is3prong == 1) {
          if ((mode == 2 && flag == false) || (mode == 1 && flag == true) || mode == 0) {
            sum_pt += pt;
            if (var == 0)
              sum_dx_pt += pt * dr;
            else if (var == 1)
              sum_dx_pt += pt * deta;
            else if (var == 2)
              sum_dx_pt += pt * dphi;
          }
        }
      }
      if (sum_pt > 0.) {
        return sum_dx_pt / sum_pt;
      }
      return 0.;
    }
    float pt_weighted_dx(const reco::PFTau& tau, int mode, int var, int decaymode) {
      return pt_weighted_dx_T<PFTau_traits>(tau, mode, var, decaymode);
    }
    float pt_weighted_dx(const pat::Tau& tau, int mode, int var, int decaymode) {
      return pt_weighted_dx_T<PATTau_traits>(tau, mode, var, decaymode);
    }
    unsigned int n_photons_total(const reco::PFTau& tau) {
      unsigned int n_photons = 0;
      for (auto& cand : tau.signalGammaCands()) {
        if (cand->pt() > 0.5)
          ++n_photons;
      }
      for (auto& cand : tau.isolationGammaCands()) {
        if (cand->pt() > 0.5)
          ++n_photons;
      }
      return n_photons;
    }
    unsigned int n_photons_total(const pat::Tau& tau) {
      unsigned int n_photons = 0;
      for (auto& cand : tau.signalGammaCands()) {
        if (cand->pt() > 0.5)
          ++n_photons;
      }
      for (auto& cand : tau.isolationGammaCands()) {
        if (cand->pt() > 0.5)
          ++n_photons;
      }
      return n_photons;
    }

    bool fillIsoMVARun2Inputs(float* mvaInput,
                              const pat::Tau& tau,
                              int mvaOpt,
                              const std::string& nameCharged,
                              const std::string& nameNeutral,
                              const std::string& namePu,
                              const std::string& nameOutside,
                              const std::string& nameFootprint) {
      int tauDecayMode = tau.decayMode();

      if (((mvaOpt == kOldDMwoLT || mvaOpt == kOldDMwLT || mvaOpt == kDBoldDMwLT || mvaOpt == kPWoldDMwLT ||
            mvaOpt == kDBoldDMwLTwGJ) &&
           (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
          ((mvaOpt == kNewDMwoLT || mvaOpt == kNewDMwLT || mvaOpt == kDBnewDMwLT || mvaOpt == kPWnewDMwLT ||
            mvaOpt == kDBnewDMwLTwGJ) &&
           (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 ||
            tauDecayMode == 10 || tauDecayMode == 11))) {
        float chargedIsoPtSum = tau.tauID(nameCharged);
        float neutralIsoPtSum = tau.tauID(nameNeutral);
        float puCorrPtSum = tau.tauID(namePu);
        float photonPtSumOutsideSignalCone = tau.tauID(nameOutside);
        float footprintCorrection = tau.tauID(nameFootprint);

        float decayDistX = tau.flightLength().x();
        float decayDistY = tau.flightLength().y();
        float decayDistZ = tau.flightLength().z();
        float decayDistMag = std::sqrt(decayDistX * decayDistX + decayDistY * decayDistY + decayDistZ * decayDistZ);

        // --- The following 5 variables differ slightly between AOD & MiniAOD
        //     because they are recomputed using packedCandidates saved in the tau
        float nPhoton = reco::tau::n_photons_total(tau);
        float ptWeightedDetaStrip = reco::tau::pt_weighted_deta_strip(tau, tauDecayMode);
        float ptWeightedDphiStrip = reco::tau::pt_weighted_dphi_strip(tau, tauDecayMode);
        float ptWeightedDrSignal = reco::tau::pt_weighted_dr_signal(tau, tauDecayMode);
        float ptWeightedDrIsolation = reco::tau::pt_weighted_dr_iso(tau, tauDecayMode);
        // ---
        float leadingTrackChi2 = tau.leadingTrackNormChi2();
        float eRatio = reco::tau::eratio(tau);

        // Difference between measured and maximally allowed Gottfried-Jackson angle
        float gjAngleDiff = -999;
        if (tauDecayMode == 10) {
          double mTau = 1.77682;
          double mAOne = tau.p4().M();
          double pAOneMag = tau.p();
          double argumentThetaGJmax = (std::pow(mTau, 2) - std::pow(mAOne, 2)) / (2 * mTau * pAOneMag);
          double argumentThetaGJmeasured =
              (tau.p4().px() * decayDistX + tau.p4().py() * decayDistY + tau.p4().pz() * decayDistZ) /
              (pAOneMag * decayDistMag);
          if (std::abs(argumentThetaGJmax) <= 1. && std::abs(argumentThetaGJmeasured) <= 1.) {
            double thetaGJmax = std::asin(argumentThetaGJmax);
            double thetaGJmeasured = std::acos(argumentThetaGJmeasured);
            gjAngleDiff = thetaGJmeasured - thetaGJmax;
          }
        }

        if (mvaOpt == kOldDMwoLT || mvaOpt == kNewDMwoLT) {
          mvaInput[0] = std::log(std::max(1.f, (float)tau.pt()));
          mvaInput[1] = std::abs((float)tau.eta());
          mvaInput[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
          mvaInput[3] = std::log(std::max(1.e-2f, neutralIsoPtSum - 0.125f * puCorrPtSum));
          mvaInput[4] = std::log(std::max(1.e-2f, puCorrPtSum));
          mvaInput[5] = tauDecayMode;
        } else if (mvaOpt == kOldDMwLT || mvaOpt == kNewDMwLT) {
          mvaInput[0] = std::log(std::max(1.f, (float)tau.pt()));
          mvaInput[1] = std::abs((float)tau.eta());
          mvaInput[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
          mvaInput[3] = std::log(std::max(1.e-2f, neutralIsoPtSum - 0.125f * puCorrPtSum));
          mvaInput[4] = std::log(std::max(1.e-2f, puCorrPtSum));
          mvaInput[5] = tauDecayMode;
          mvaInput[6] = std::copysign(+1.f, tau.dxy());
          mvaInput[7] = std::sqrt(std::min(1.f, std::abs(tau.dxy())));
          mvaInput[8] = std::min(10.f, std::abs(tau.dxy_Sig()));
          mvaInput[9] = (tau.hasSecondaryVertex()) ? 1. : 0.;
          mvaInput[10] = std::sqrt(decayDistMag);
          mvaInput[11] = std::min(10.f, tau.flightLengthSig());
        } else if (mvaOpt == kDBoldDMwLT || mvaOpt == kDBnewDMwLT) {
          mvaInput[0] = std::log(std::max(1.f, (float)tau.pt()));
          mvaInput[1] = std::abs((float)tau.eta());
          mvaInput[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
          mvaInput[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
          mvaInput[4] = std::log(std::max(1.e-2f, puCorrPtSum));
          mvaInput[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
          mvaInput[6] = tauDecayMode;
          mvaInput[7] = std::min(30.f, nPhoton);
          mvaInput[8] = std::min(0.5f, ptWeightedDetaStrip);
          mvaInput[9] = std::min(0.5f, ptWeightedDphiStrip);
          mvaInput[10] = std::min(0.5f, ptWeightedDrSignal);
          mvaInput[11] = std::min(0.5f, ptWeightedDrIsolation);
          mvaInput[12] = std::min(100.f, leadingTrackChi2);
          mvaInput[13] = std::min(1.f, eRatio);
          mvaInput[14] = std::copysign(+1.f, tau.dxy());
          mvaInput[15] = std::sqrt(std::min(1.f, std::abs(tau.dxy())));
          mvaInput[16] = std::min(10.f, std::abs(tau.dxy_Sig()));
          mvaInput[17] = std::copysign(+1.f, tau.ip3d());
          mvaInput[18] = std::sqrt(std::min(1.f, std::abs(tau.ip3d())));
          mvaInput[19] = std::min(10.f, std::abs(tau.ip3d_Sig()));
          mvaInput[20] = (tau.hasSecondaryVertex()) ? 1. : 0.;
          mvaInput[21] = std::sqrt(decayDistMag);
          mvaInput[22] = std::min(10.f, tau.flightLengthSig());
        } else if (mvaOpt == kPWoldDMwLT || mvaOpt == kPWnewDMwLT) {
          mvaInput[0] = std::log(std::max(1.f, (float)tau.pt()));
          mvaInput[1] = std::abs((float)tau.eta());
          mvaInput[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
          mvaInput[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
          mvaInput[4] = std::log(std::max(1.e-2f, footprintCorrection));
          mvaInput[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
          mvaInput[6] = tauDecayMode;
          mvaInput[7] = std::min(30.f, nPhoton);
          mvaInput[8] = std::min(0.5f, ptWeightedDetaStrip);
          mvaInput[9] = std::min(0.5f, ptWeightedDphiStrip);
          mvaInput[10] = std::min(0.5f, ptWeightedDrSignal);
          mvaInput[11] = std::min(0.5f, ptWeightedDrIsolation);
          mvaInput[12] = std::min(100.f, leadingTrackChi2);
          mvaInput[13] = std::min(1.f, eRatio);
          mvaInput[14] = std::copysign(+1.f, tau.dxy());
          mvaInput[15] = std::sqrt(std::min(1.f, std::abs(tau.dxy())));
          mvaInput[16] = std::min(10.f, std::abs(tau.dxy_Sig()));
          mvaInput[17] = std::copysign(+1.f, tau.ip3d());
          mvaInput[18] = std::sqrt(std::min(1.f, std::abs(tau.ip3d())));
          mvaInput[19] = std::min(10.f, std::abs(tau.ip3d_Sig()));
          mvaInput[20] = (tau.hasSecondaryVertex()) ? 1. : 0.;
          mvaInput[21] = std::sqrt(decayDistMag);
          mvaInput[22] = std::min(10.f, tau.flightLengthSig());
        } else if (mvaOpt == kDBoldDMwLTwGJ || mvaOpt == kDBnewDMwLTwGJ) {
          mvaInput[0] = std::log(std::max(1.f, (float)tau.pt()));
          mvaInput[1] = std::abs((float)tau.eta());
          mvaInput[2] = std::log(std::max(1.e-2f, chargedIsoPtSum));
          mvaInput[3] = std::log(std::max(1.e-2f, neutralIsoPtSum));
          mvaInput[4] = std::log(std::max(1.e-2f, puCorrPtSum));
          mvaInput[5] = std::log(std::max(1.e-2f, photonPtSumOutsideSignalCone));
          mvaInput[6] = tauDecayMode;
          mvaInput[7] = std::min(30.f, nPhoton);
          mvaInput[8] = std::min(0.5f, ptWeightedDetaStrip);
          mvaInput[9] = std::min(0.5f, ptWeightedDphiStrip);
          mvaInput[10] = std::min(0.5f, ptWeightedDrSignal);
          mvaInput[11] = std::min(0.5f, ptWeightedDrIsolation);
          mvaInput[12] = std::min(1.f, eRatio);
          mvaInput[13] = std::copysign(+1.f, tau.dxy());
          mvaInput[14] = std::sqrt(std::min(1.f, std::abs(tau.dxy())));
          mvaInput[15] = std::min(10.f, std::abs(tau.dxy_Sig()));
          mvaInput[16] = std::copysign(+1.f, tau.ip3d());
          mvaInput[17] = std::sqrt(std::min(1.f, std::abs(tau.ip3d())));
          mvaInput[18] = std::min(10.f, std::abs(tau.ip3d_Sig()));
          mvaInput[19] = (tau.hasSecondaryVertex()) ? 1. : 0.;
          mvaInput[20] = std::sqrt(decayDistMag);
          mvaInput[21] = std::min(10.f, tau.flightLengthSig());
          mvaInput[22] = std::max(-1.f, gjAngleDiff);
        }

        return true;
      }
      return false;
    }

  }  // namespace tau
}  // namespace reco