PFRecoTauDiscriminationAgainstMuonSimple.cc

Go to the documentation of this file.

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

#include "FWCore/Utilities/interface/Exception.h"

#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/Math/interface/deltaR.h"

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

#include <vector>
#include <string>
#include <iostream>
#include <atomic>

using reco::tau::format_vint;

namespace {

  class PFRecoTauDiscriminationAgainstMuonSimple final : public PFTauDiscriminationProducerBase {
  public:
    explicit PFRecoTauDiscriminationAgainstMuonSimple(const edm::ParameterSet& cfg)
        : PFTauDiscriminationProducerBase(cfg), moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
      hop_ = cfg.getParameter<double>("HoPMin");
      doCaloMuonVeto_ = cfg.getParameter<bool>("doCaloMuonVeto");
      srcPatMuons_ = cfg.getParameter<edm::InputTag>("srcPatMuons");
      Muons_token = consumes<pat::MuonCollection>(srcPatMuons_);
      dRmuonMatch_ = cfg.getParameter<double>("dRmuonMatch");
      dRmuonMatchLimitedToJetArea_ = cfg.getParameter<bool>("dRmuonMatchLimitedToJetArea");
      minPtMatchedMuon_ = cfg.getParameter<double>("minPtMatchedMuon");
      maxNumberOfMatches_ = cfg.getParameter<int>("maxNumberOfMatches");
      maxNumberOfHitsLast2Stations_ = cfg.getParameter<int>("maxNumberOfHitsLast2Stations");
      typedef std::vector<int> vint;
      maskMatchesDT_ = cfg.getParameter<vint>("maskMatchesDT");
      maskMatchesCSC_ = cfg.getParameter<vint>("maskMatchesCSC");
      maskMatchesRPC_ = cfg.getParameter<vint>("maskMatchesRPC");
      maskHitsDT_ = cfg.getParameter<vint>("maskHitsDT");
      maskHitsCSC_ = cfg.getParameter<vint>("maskHitsCSC");
      maskHitsRPC_ = cfg.getParameter<vint>("maskHitsRPC");
      maxNumberOfSTAMuons_ = cfg.getParameter<int>("maxNumberOfSTAMuons");
      maxNumberOfRPCMuons_ = cfg.getParameter<int>("maxNumberOfRPCMuons");

      verbosity_ = cfg.exists("verbosity") ? cfg.getParameter<int>("verbosity") : 0;
    }
    ~PFRecoTauDiscriminationAgainstMuonSimple() override {}

    void beginEvent(const edm::Event&, const edm::EventSetup&) override;

    double discriminate(const reco::PFTauRef&) const override;

  private:
    std::string moduleLabel_;
    int discriminatorOption_;
    double hop_;
    bool doCaloMuonVeto_;
    edm::InputTag srcPatMuons_;
    edm::Handle<pat::MuonCollection> muons_;
    edm::EDGetTokenT<pat::MuonCollection> Muons_token;
    double dRmuonMatch_;
    bool dRmuonMatchLimitedToJetArea_;
    double minPtMatchedMuon_;
    int maxNumberOfMatches_;
    std::vector<int> maskMatchesDT_;
    std::vector<int> maskMatchesCSC_;
    std::vector<int> maskMatchesRPC_;
    int maxNumberOfHitsLast2Stations_;
    std::vector<int> maskHitsDT_;
    std::vector<int> maskHitsCSC_;
    std::vector<int> maskHitsRPC_;
    int maxNumberOfSTAMuons_, maxNumberOfRPCMuons_;

    static std::atomic<unsigned int> numWarnings_;
    static constexpr unsigned int maxWarnings_ = 3;
    int verbosity_;
  };

  std::atomic<unsigned int> PFRecoTauDiscriminationAgainstMuonSimple::numWarnings_{0};

  void PFRecoTauDiscriminationAgainstMuonSimple::beginEvent(const edm::Event& evt, const edm::EventSetup& es) {
    evt.getByToken(Muons_token, muons_);
  }

  double PFRecoTauDiscriminationAgainstMuonSimple::discriminate(const reco::PFTauRef& pfTau) const {
    if (verbosity_) {
      edm::LogPrint("PFTauAgainstMuonSimple") << "<PFRecoTauDiscriminationAgainstMuonSimple::discriminate>:";
      edm::LogPrint("PFTauAgainstMuonSimple") << " moduleLabel = " << moduleLabel_;
      edm::LogPrint("PFTauAgainstMuonSimple")
          << "tau #" << pfTau.key() << ": Pt = " << pfTau->pt() << ", eta = " << pfTau->eta()
          << ", phi = " << pfTau->phi() << ", decay mode = " << pfTau->decayMode();
    }

    //if (pfTau->decayMode() >= 5) return true; //MB: accept all multi-prongs??

    const reco::CandidatePtr& pfLeadChargedHadron = pfTau->leadChargedHadrCand();
    bool passesCaloMuonVeto = true;
    if (pfLeadChargedHadron.isNonnull()) {
      const pat::PackedCandidate* pCand = dynamic_cast<const pat::PackedCandidate*>(pfLeadChargedHadron.get());
      if (pCand != nullptr) {
        double rawCaloEnergyFraction = pCand->rawCaloFraction();
        //if ( !(rawCaloEnergyFraction > 0.) ) rawCaloEnergyFraction = 99; //MB: hack against cases when rawCaloEnergyFraction is not stored; it makes performance of the H/P cut rather poor
        if (verbosity_) {
          edm::LogPrint("PFTauAgainstMuonSimple")
              << "decayMode = " << pfTau->decayMode()
              << ", rawCaloEnergy(ECAL+HCAL)Fraction = " << rawCaloEnergyFraction
              << ", leadPFChargedHadronP = " << pCand->p() << ", leadPFChargedHadron pdgId = " << pCand->pdgId();
        }
        if (pfTau->decayMode() == 0 && rawCaloEnergyFraction < hop_)
          passesCaloMuonVeto = false;
      }
    }
    if (doCaloMuonVeto_ && !passesCaloMuonVeto) {
      if (verbosity_)
        edm::LogPrint("PFTauAgainstMuonSimple") << "--> CaloMuonVeto failed, returning 0";
      return 0.;
    }

    //iterate over muons to find ones to use for discrimination
    std::vector<const pat::Muon*> muonsToCheck;
    size_t numMuons = muons_->size();
    for (size_t idxMuon = 0; idxMuon < numMuons; ++idxMuon) {
      bool matched = false;
      const pat::MuonRef muon(muons_, idxMuon);
      if (verbosity_)
        edm::LogPrint("PFTauAgainstMuonSimple") << "muon #" << muon.key() << ": Pt = " << muon->pt()
                                                << ", eta = " << muon->eta() << ", phi = " << muon->phi();
      //firsty check if the muon corrsponds with the leading tau particle
      for (size_t iCand = 0; iCand < muon->numberOfSourceCandidatePtrs(); ++iCand) {
        const reco::CandidatePtr& srcCand = muon->sourceCandidatePtr(iCand);
        if (srcCand.isNonnull() && pfLeadChargedHadron.isNonnull() && srcCand.id() == pfLeadChargedHadron.id() &&
            srcCand.key() == pfLeadChargedHadron.key()) {
          muonsToCheck.push_back(muon.get());
          matched = true;
          if (verbosity_)
            edm::LogPrint("PFTauAgainstMuonSimple") << " tau has muonRef.";
          break;
        }
      }
      if (matched)
        continue;
      //check dR matching
      if (!(muon->pt() > minPtMatchedMuon_)) {
        if (verbosity_) {
          edm::LogPrint("PFTauAgainstMuonSimple") << " fails Pt cut --> skipping it.";
        }
        continue;
      }
      double dR = deltaR(muon->p4(), pfTau->p4());
      double dRmatch = dRmuonMatch_;
      if (dRmuonMatchLimitedToJetArea_) {
        double jetArea = 0.;
        if (pfTau->jetRef().isNonnull())
          jetArea = pfTau->jetRef()->jetArea();
        if (jetArea > 0.) {
          dRmatch = TMath::Min(dRmatch, TMath::Sqrt(jetArea / TMath::Pi()));
        } else {
          if (numWarnings_ < maxWarnings_) {
            edm::LogInfo("PFRecoTauDiscriminationAgainstMuonSimple::discriminate")
                << "Jet associated to Tau: Pt = " << pfTau->pt() << ", eta = " << pfTau->eta()
                << ", phi = " << pfTau->phi() << " has area = " << jetArea << " !!";
            ++numWarnings_;
          }
          dRmatch = pfTau->signalConeSize();
        }
        dRmatch = TMath::Max(dRmatch, pfTau->signalConeSize());
        if (dR < dRmatch) {
          if (verbosity_)
            edm::LogPrint("PFTauAgainstMuonSimple") << " overlaps with tau, dR = " << dR;
          muonsToCheck.push_back(muon.get());
        }
      }
    }
    //now examine selected muons
    bool pass = true;
    std::vector<int> numMatchesDT(4);
    std::vector<int> numMatchesCSC(4);
    std::vector<int> numMatchesRPC(4);
    std::vector<int> numHitsDT(4);
    std::vector<int> numHitsCSC(4);
    std::vector<int> numHitsRPC(4);
    //MB: clear counters of matched segments and hits globally as in the AgainstMuon2 discriminant, but note that they will sum for all matched muons. However it is not likely that there is more than one matched muon.
    for (int iStation = 0; iStation < 4; ++iStation) {
      numMatchesDT[iStation] = 0;
      numMatchesCSC[iStation] = 0;
      numMatchesRPC[iStation] = 0;
      numHitsDT[iStation] = 0;
      numHitsCSC[iStation] = 0;
      numHitsRPC[iStation] = 0;
    }
    int numSTAMuons = 0, numRPCMuons = 0;
    if (verbosity_ && !muonsToCheck.empty())
      edm::LogPrint("PFTauAgainstMuonSimple") << "Muons to check (" << muonsToCheck.size() << "):";
    size_t iMu = 0;
    for (const auto& mu : muonsToCheck) {
      if (mu->isStandAloneMuon())
        numSTAMuons++;
      if (mu->muonID("RPCMuLoose"))
        numRPCMuons++;
      reco::tau::countMatches(*mu, numMatchesDT, numMatchesCSC, numMatchesRPC);
      int numStationsWithMatches = 0;
      for (int iStation = 0; iStation < 4; ++iStation) {
        if (numMatchesDT[iStation] > 0 && !maskMatchesDT_[iStation])
          ++numStationsWithMatches;
        if (numMatchesCSC[iStation] > 0 && !maskMatchesCSC_[iStation])
          ++numStationsWithMatches;
        if (numMatchesRPC[iStation] > 0 && !maskMatchesRPC_[iStation])
          ++numStationsWithMatches;
      }
      reco::tau::countHits(*mu, numHitsDT, numHitsCSC, numHitsRPC);
      int numLast2StationsWithHits = 0;
      for (int iStation = 2; iStation < 4; ++iStation) {
        if (numHitsDT[iStation] > 0 && !maskHitsDT_[iStation])
          ++numLast2StationsWithHits;
        if (numHitsCSC[iStation] > 0 && !maskHitsCSC_[iStation])
          ++numLast2StationsWithHits;
        if (numHitsRPC[iStation] > 0 && !maskHitsRPC_[iStation])
          ++numLast2StationsWithHits;
      }
      if (verbosity_)
        edm::LogPrint("PFTauAgainstMuonSimple")
            << "\t" << iMu << ": Pt = " << mu->pt() << ", eta = " << mu->eta() << ", phi = " << mu->phi() << "\n\t"
            << "   isSTA: " << mu->isStandAloneMuon() << ", isRPCLoose: " << mu->muonID("RPCMuLoose")
            << "\n\t   numMatchesDT  = " << format_vint(numMatchesDT)
            << "\n\t   numMatchesCSC = " << format_vint(numMatchesCSC)
            << "\n\t   numMatchesRPC = " << format_vint(numMatchesRPC)
            << "\n\t   --> numStationsWithMatches = " << numStationsWithMatches
            << "\n\t   numHitsDT  = " << format_vint(numHitsDT) << "\n\t   numHitsCSC = " << format_vint(numHitsCSC)
            << "\n\t   numHitsRPC = " << format_vint(numHitsRPC)
            << "\n\t   --> numLast2StationsWithHits = " << numLast2StationsWithHits;
      if (maxNumberOfMatches_ >= 0 && numStationsWithMatches > maxNumberOfMatches_) {
        pass = false;
        break;
      }
      if (maxNumberOfHitsLast2Stations_ >= 0 && numLast2StationsWithHits > maxNumberOfHitsLast2Stations_) {
        pass = false;
        break;
      }
      if (maxNumberOfSTAMuons_ >= 0 && numSTAMuons > maxNumberOfSTAMuons_) {
        pass = false;
        break;
      }
      if (maxNumberOfRPCMuons_ >= 0 && numRPCMuons > maxNumberOfRPCMuons_) {
        pass = false;
        break;
      }
      iMu++;
    }

    double discriminatorValue = pass ? 1. : 0.;
    if (verbosity_)
      edm::LogPrint("PFTauAgainstMuonSimple") << "--> returning discriminatorValue = " << discriminatorValue;

    return discriminatorValue;
  }

}  // namespace

DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstMuonSimple);