pfRecoTauChargedHadronAuxFunctions.cc

Go to the documentation of this file.

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

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"

#include <TMath.h>

namespace reco { namespace tau {

const reco::Track* getTrackFromChargedHadron(const reco::PFRecoTauChargedHadron& chargedHadron) {
  // Charged hadron made from track (reco::Track) - RECO/AOD only
  if ( chargedHadron.getTrack().isNonnull()) {
    return chargedHadron.getTrack().get();
  }
  // In MiniAOD, even isolated tracks are saved as candidates, so the track Ptr doesn't exist
  const pat::PackedCandidate* chargedPFPCand = dynamic_cast<const pat::PackedCandidate*> (chargedHadron.getChargedPFCandidate().get());
  if (chargedPFPCand != nullptr) {
    return chargedPFPCand->bestTrack();
  }
  const pat::PackedCandidate* lostTrackCand = dynamic_cast<const pat::PackedCandidate*> (chargedHadron.getLostTrackCandidate().get());
  if (lostTrackCand != nullptr) {
    return lostTrackCand->bestTrack();
  }
  return nullptr;
}
  
void setChargedHadronP4(reco::PFRecoTauChargedHadron& chargedHadron, double scaleFactor_neutralPFCands)
{
  double chargedHadronP  = 0.;
  double chargedHadronPx = 0.;
  double chargedHadronPy = 0.;
  double chargedHadronPz = 0.;
  double SumNeutrals = 0.;
  if ( chargedHadron.algoIs(reco::PFRecoTauChargedHadron::kChargedPFCandidate) ||
       chargedHadron.algoIs(reco::PFRecoTauChargedHadron::kPFNeutralHadron)   ) {
    const reco::CandidatePtr& chargedPFCand = chargedHadron.getChargedPFCandidate();
    assert(chargedPFCand.isNonnull());
    chargedHadronP     += chargedPFCand->p();
    chargedHadronPx     = chargedPFCand->px();
    chargedHadronPy     = chargedPFCand->py();
    chargedHadronPz     = chargedPFCand->pz();
  } else if ( chargedHadron.algoIs(reco::PFRecoTauChargedHadron::kTrack) ) {
    const reco::Track* track = getTrackFromChargedHadron(chargedHadron);
    if (track != nullptr) {
      chargedHadronP     += track->p();
      chargedHadronPx     = track->px();
      chargedHadronPy     = track->py();
      chargedHadronPz     = track->pz();
    } else { // lost tracks from MiniAOD that don't have track information saved
      const reco::CandidatePtr& lostTrack = chargedHadron.getLostTrackCandidate();
      assert(lostTrack.isNonnull());
      chargedHadronP     += lostTrack->p();
      chargedHadronPx     = lostTrack->px();
      chargedHadronPy     = lostTrack->py();
      chargedHadronPz     = lostTrack->pz();
    }
  } else assert(0);
  const std::vector<reco::CandidatePtr>& neutralPFCands = chargedHadron.getNeutralPFCandidates();
  for ( std::vector<reco::CandidatePtr>::const_iterator neutralPFCand = neutralPFCands.begin();
    neutralPFCand != neutralPFCands.end(); ++neutralPFCand ) {
    SumNeutrals += (*neutralPFCand)->p();
  }
  double noNeutrals=chargedHadronP;
  chargedHadronP+=scaleFactor_neutralPFCands*SumNeutrals;
  double ptRatio=chargedHadronP/noNeutrals;
  chargedHadronPx*=ptRatio;
  chargedHadronPy*=ptRatio;
  chargedHadronPz*=ptRatio;
  
      
  reco::Candidate::LorentzVector chargedHadronP4 = compChargedHadronP4fromPxPyPz(chargedHadronPx, chargedHadronPy, chargedHadronPz);
  chargedHadron.setP4(chargedHadronP4);
}

reco::Candidate::LorentzVector compChargedHadronP4fromPxPyPz(double chargedHadronPx, double chargedHadronPy, double chargedHadronPz)
{
  const double chargedPionMass = 0.13957; // GeV  
  double chargedHadronEn = sqrt(chargedHadronPx*chargedHadronPx + chargedHadronPy*chargedHadronPy + chargedHadronPz*chargedHadronPz + chargedPionMass*chargedPionMass);  
  reco::Candidate::LorentzVector chargedHadronP4(chargedHadronPx, chargedHadronPy, chargedHadronPz, chargedHadronEn);
  return chargedHadronP4;
}

reco::Candidate::LorentzVector compChargedHadronP4fromPThetaPhi(double chargedHadronP, double chargedHadronTheta, double chargedHadronPhi)
{
  double chargedHadronPx = chargedHadronP*TMath::Cos(chargedHadronPhi)*TMath::Sin(chargedHadronTheta);
  double chargedHadronPy = chargedHadronP*TMath::Sin(chargedHadronPhi)*TMath::Sin(chargedHadronTheta);
  double chargedHadronPz = chargedHadronP*TMath::Cos(chargedHadronTheta);    
  return compChargedHadronP4fromPxPyPz(chargedHadronPx, chargedHadronPy, chargedHadronPz);
}

}} // end namespace reco::tau