Tau.cc

Go to the documentation of this file.

//
//

#include "DataFormats/PatCandidates/interface/Tau.h"
#include "DataFormats/JetReco/interface/GenJet.h"
#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
#include <algorithm>
#include <boost/bind.hpp>

using namespace pat;


Tau::Tau() :
    Lepton<reco::BaseTau>()
    ,embeddedIsolationTracks_(false)
    ,embeddedLeadTrack_(false)
    ,embeddedSignalTracks_(false)
    ,embeddedLeadPFCand_(false)
    ,embeddedLeadPFChargedHadrCand_(false)
    ,embeddedLeadPFNeutralCand_(false)
    ,embeddedSignalPFCands_(false)
    ,embeddedSignalPFChargedHadrCands_(false)
    ,embeddedSignalPFNeutralHadrCands_(false)
    ,embeddedSignalPFGammaCands_(false)
    ,embeddedIsolationPFCands_(false)
    ,embeddedIsolationPFChargedHadrCands_(false)
    ,embeddedIsolationPFNeutralHadrCands_(false)
    ,embeddedIsolationPFGammaCands_(false)
{
}

Tau::Tau(const reco::BaseTau & aTau) :
    Lepton<reco::BaseTau>(aTau)
    ,embeddedIsolationTracks_(false)
    ,embeddedLeadTrack_(false)
    ,embeddedSignalTracks_(false)
    ,embeddedLeadPFCand_(false)
    ,embeddedLeadPFChargedHadrCand_(false)
    ,embeddedLeadPFNeutralCand_(false)
    ,embeddedSignalPFCands_(false)
    ,embeddedSignalPFChargedHadrCands_(false)
    ,embeddedSignalPFNeutralHadrCands_(false)
    ,embeddedSignalPFGammaCands_(false)
    ,embeddedIsolationPFCands_(false)
    ,embeddedIsolationPFChargedHadrCands_(false)
    ,embeddedIsolationPFNeutralHadrCands_(false)
    ,embeddedIsolationPFGammaCands_(false)
{
    initFromBaseTau(aTau);
}

Tau::Tau(const edm::RefToBase<reco::BaseTau> & aTauRef) :
    Lepton<reco::BaseTau>(aTauRef)
    ,embeddedIsolationTracks_(false)
    ,embeddedLeadTrack_(false)
    ,embeddedSignalTracks_(false)
    ,embeddedLeadPFCand_(false)
    ,embeddedLeadPFChargedHadrCand_(false)
    ,embeddedLeadPFNeutralCand_(false)
    ,embeddedSignalPFCands_(false)
    ,embeddedSignalPFChargedHadrCands_(false)
    ,embeddedSignalPFNeutralHadrCands_(false)
    ,embeddedSignalPFGammaCands_(false)
    ,embeddedIsolationPFCands_(false)
    ,embeddedIsolationPFChargedHadrCands_(false)
    ,embeddedIsolationPFNeutralHadrCands_(false)
    ,embeddedIsolationPFGammaCands_(false)
{
    initFromBaseTau(*aTauRef);
}

Tau::Tau(const edm::Ptr<reco::BaseTau> & aTauRef) :
    Lepton<reco::BaseTau>(aTauRef)
    ,embeddedIsolationTracks_(false)
    ,embeddedLeadTrack_(false)
    ,embeddedSignalTracks_(false)
    ,embeddedLeadPFCand_(false)
    ,embeddedLeadPFChargedHadrCand_(false)
    ,embeddedLeadPFNeutralCand_(false)
    ,embeddedSignalPFCands_(false)
    ,embeddedSignalPFChargedHadrCands_(false)
    ,embeddedSignalPFNeutralHadrCands_(false)
    ,embeddedSignalPFGammaCands_(false)
    ,embeddedIsolationPFCands_(false)
    ,embeddedIsolationPFChargedHadrCands_(false)
    ,embeddedIsolationPFNeutralHadrCands_(false)
    ,embeddedIsolationPFGammaCands_(false)
{
    initFromBaseTau(*aTauRef);
}

void Tau::initFromBaseTau(const reco::BaseTau& aTau) {
    const reco::PFTau * pfTau = dynamic_cast<const reco::PFTau *>(&aTau);
    if (pfTau != nullptr){
      // If PFTau is made from PackedCandidates, directly fill slimmed version
      // without PFSpecific
      const pat::PackedCandidate* pc = dynamic_cast<const pat::PackedCandidate*>(pfTau->leadChargedHadrCand().get());
      if (pc != nullptr) {
        for (const auto& ptr : pfTau->signalChargedHadrCands())
          signalChargedHadrCandPtrs_.push_back(ptr);

        for (const auto& ptr : pfTau->signalNeutrHadrCands())
          signalNeutralHadrCandPtrs_.push_back(ptr);

        for (const auto& ptr : pfTau->signalGammaCands())
          signalGammaCandPtrs_.push_back(ptr);

        for (const auto& ptr : pfTau->isolationChargedHadrCands())
          isolationChargedHadrCandPtrs_.push_back(ptr);

        for (const auto& ptr : pfTau->isolationNeutrHadrCands())
          isolationNeutralHadrCandPtrs_.push_back(ptr);

        for (const auto& ptr : pfTau->isolationGammaCands())
          isolationGammaCandPtrs_.push_back(ptr);
      }
      else {
        pfSpecific_.push_back(pat::tau::TauPFSpecific(*pfTau));
      }
      pfEssential_.push_back(pat::tau::TauPFEssential(*pfTau));
    }
    const reco::CaloTau * caloTau = dynamic_cast<const reco::CaloTau *>(&aTau);
    if (caloTau != nullptr) caloSpecific_.push_back(pat::tau::TauCaloSpecific(*caloTau));
}

Tau::~Tau() {
}

std::ostream& 
reco::operator<<(std::ostream& out, const pat::Tau& obj) 
{
  if(!out) return out;
  
  out << "\tpat::Tau: ";
  out << std::setiosflags(std::ios::right);
  out << std::setiosflags(std::ios::fixed);
  out << std::setprecision(3);
  out << " E/pT/eta/phi " 
      << obj.energy()<<"/"
      << obj.pt()<<"/"
      << obj.eta()<<"/"
      << obj.phi();
  return out; 
}

const reco::TrackRefVector & Tau::isolationTracks() const {
  if (embeddedIsolationTracks_) {
    if (!isolationTracksTransientRefVector_.isSet()) {
    std::unique_ptr<reco::TrackRefVector> trackRefVec{ new reco::TrackRefVector{}};
    trackRefVec->reserve(isolationTracks_.size());
        for (unsigned int i = 0; i < isolationTracks_.size(); i++) {
          trackRefVec->push_back(reco::TrackRef(&isolationTracks_, i));
        }
    isolationTracksTransientRefVector_.set(std::move(trackRefVec));
    }
    return *isolationTracksTransientRefVector_;
  } else {
    return reco::BaseTau::isolationTracks();
  }
}


reco::TrackRef Tau::leadTrack() const {
  if (embeddedLeadTrack_) {
    return reco::TrackRef(&leadTrack_, 0);
  } else {
    return reco::BaseTau::leadTrack();
  }
}


const reco::TrackRefVector & Tau::signalTracks() const {
  if (embeddedSignalTracks_) {
    if (!signalTracksTransientRefVector_.isSet()) {
        std::unique_ptr<reco::TrackRefVector> trackRefVec{ new reco::TrackRefVector{} };
        trackRefVec->reserve(signalTracks_.size());
        for (unsigned int i = 0; i < signalTracks_.size(); i++) {
          trackRefVec->push_back(reco::TrackRef(&signalTracks_, i));
        }
    signalTracksTransientRefVector_.set(std::move(trackRefVec));
    }
    return *signalTracksTransientRefVector_;
  } else {
    return reco::BaseTau::signalTracks();
  }
}


void Tau::embedIsolationTracks() {
  isolationTracks_.clear();
  reco::TrackRefVector trackRefVec = reco::BaseTau::isolationTracks();
  for (unsigned int i = 0; i < trackRefVec.size(); i++) {
    isolationTracks_.push_back(*trackRefVec.at(i));
  }
  embeddedIsolationTracks_ = true;
}


void Tau::embedLeadTrack() {
  leadTrack_.clear();
  if (reco::BaseTau::leadTrack().isNonnull()) {
      leadTrack_.push_back(*reco::BaseTau::leadTrack());
      embeddedLeadTrack_ = true;
  }
}


void Tau::embedSignalTracks(){
  signalTracks_.clear();
  reco::TrackRefVector trackRefVec = reco::BaseTau::signalTracks();
  for (unsigned int i = 0; i < trackRefVec.size(); i++) {
    signalTracks_.push_back(*trackRefVec.at(i));
  }
  embeddedSignalTracks_ = true;
}


void Tau::setGenJet(const reco::GenJetRef& gj) {
  genJet_.clear();
  genJet_.push_back(*gj);
}

const reco::GenJet * Tau::genJet() const {
  return (!genJet_.empty() ? &genJet_.front() : nullptr);
}


// method to retrieve a tau ID (or throw)
float Tau::tauID(const std::string & name) const {
  for (std::vector<IdPair>::const_iterator it = tauIDs_.begin(), ed = tauIDs_.end(); it != ed; ++it) {
    if (it->first == name) return it->second;
  }
  cms::Exception ex("Key not found");
  ex << "pat::Tau: the ID " << name << " can't be found in this pat::Tau.\n";
  ex << "The available IDs are: ";
  for (std::vector<IdPair>::const_iterator it = tauIDs_.begin(), ed = tauIDs_.end(); it != ed; ++it) {
    ex << "'" << it->first << "' ";
  }
  ex << ".\n";
  throw ex;
}
// check if an ID is there
bool Tau::isTauIDAvailable(const std::string & name) const {
  for (std::vector<IdPair>::const_iterator it = tauIDs_.begin(), ed = tauIDs_.end(); it != ed; ++it) {
    if (it->first == name) return true;
  }
  return false;
}


const pat::tau::TauPFSpecific & Tau::pfSpecific() const {
  if (!isPFTau()) throw cms::Exception("Type Error") << "Requesting a PFTau-specific information from a pat::Tau which wasn't made from a PFTau.\n";
  return pfSpecific_[0]; 
}

const pat::tau::TauPFEssential & Tau::pfEssential() const {
  if (pfEssential_.empty()) throw cms::Exception("Type Error") << "Requesting a PFTau-specific information from a pat::Tau which wasn't made from a PFTau.\n";
  return pfEssential_[0]; 
}


const pat::tau::TauCaloSpecific & Tau::caloSpecific() const {
  if (!isCaloTau()) throw cms::Exception("Type Error") << "Requesting a CaloTau-specific information from a pat::Tau which wasn't made from a CaloTau.\n";
  return caloSpecific_[0]; 
}

reco::Candidate::LorentzVector Tau::p4Jet() const
{
  if ( isCaloTau() ) return caloSpecific().p4Jet_;
  if ( isPFTau()   ) return reco::Candidate::LorentzVector(pfEssential().p4Jet_);
  throw cms::Exception("Type Error") << "Requesting a CaloTau/PFTau-specific information from a pat::Tau which wasn't made from either a CaloTau or a PFTau.\n";
}

float Tau::dxy_Sig() const
{
  if ( pfEssential().dxy_error_ != 0 ) return (pfEssential().dxy_/pfEssential().dxy_error_);
  else return 0.;
}

pat::tau::TauPFEssential::CovMatrix Tau::flightLengthCov() const
{
  pat::tau::TauPFEssential::CovMatrix cov;
  const pat::tau::TauPFEssential::CovMatrix& sv = secondaryVertexCov();
  const pat::tau::TauPFEssential::CovMatrix& pv = primaryVertexCov();
  for ( int i = 0; i < dimension; ++i ) {
    for ( int j = 0; j < dimension; ++j ) {
      cov(i,j) = sv(i,j) + pv(i,j);
    }
  }
  return cov;
}

float Tau::ip3d_Sig() const
{
    if( pfEssential().ip3d_error_ != 0 ) return (pfEssential().ip3d_/pfEssential().ip3d_error_);
    else return 0.;
}

float Tau::etaetaMoment() const
{
  if ( isCaloTau() ) return caloSpecific().etaetaMoment_;
  if ( isPFTau()   ) return pfSpecific().etaetaMoment_;
  throw cms::Exception("Type Error") << "Requesting a CaloTau/PFTau-specific information from a pat::Tau which wasn't made from either a CaloTau or a PFTau.\n";
}

float Tau::phiphiMoment() const
{
  if ( isCaloTau() ) return caloSpecific().phiphiMoment_;
  if ( isPFTau()   ) return pfSpecific().phiphiMoment_;
  throw cms::Exception("Type Error") << "Requesting a CaloTau/PFTau-specific information from a pat::Tau which wasn't made from either a CaloTau or a PFTau.\n";
}

float Tau::etaphiMoment() const
{
  if ( isCaloTau() ) return caloSpecific().etaphiMoment_;
  if ( isPFTau()   ) return pfSpecific().etaphiMoment_;
  throw cms::Exception("Type Error") << "Requesting a CaloTau/PFTau-specific information from a pat::Tau which wasn't made from either a CaloTau or a PFTau.\n";
}

void Tau::setDecayMode(int decayMode)
{
  if (!isPFTau()) throw cms::Exception("Type Error") << "Requesting a PFTau-specific information from a pat::Tau which wasn't made from a PFTau.\n";
  pfEssential_[0].decayMode_ = decayMode;
} 

void Tau::embedLeadPFCand() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  leadPFCand_.clear();
  if (pfSpecific_[0].leadPFCand_.isNonnull() ) {
    leadPFCand_.push_back(*static_cast<const reco::PFCandidate*>(&*pfSpecific_[0].leadPFCand_)); //already set in C-tor
    embeddedLeadPFCand_ = true;
  }
}
void Tau::embedLeadPFChargedHadrCand() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  leadPFChargedHadrCand_.clear();
  if (pfSpecific_[0].leadPFChargedHadrCand_.isNonnull() ) {
    leadPFChargedHadrCand_.push_back(*static_cast<const reco::PFCandidate*>(&*pfSpecific_[0].leadPFChargedHadrCand_)); //already set in C-tor
    embeddedLeadPFChargedHadrCand_ = true;
  }
}
void Tau::embedLeadPFNeutralCand() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  leadPFNeutralCand_.clear();
  if (pfSpecific_[0].leadPFNeutralCand_.isNonnull() ) {
    leadPFNeutralCand_.push_back(*static_cast<const reco::PFCandidate*>(&*pfSpecific_[0].leadPFNeutralCand_)); //already set in C-tor
    embeddedLeadPFNeutralCand_ = true;
  }
}

void Tau::embedSignalPFCands() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  std::vector<reco::PFCandidatePtr> candPtrs = pfSpecific_[0].selectedSignalPFCands_;
  for (unsigned int i = 0; i < candPtrs.size(); i++) {
    signalPFCands_.push_back(candPtrs.at(i));
  }
  embeddedSignalPFCands_ = true;
}
void Tau::embedSignalPFChargedHadrCands() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  std::vector<reco::PFCandidatePtr> candPtrs = pfSpecific_[0].selectedSignalPFChargedHadrCands_;
  for (unsigned int i = 0; i < candPtrs.size(); i++) {
    signalPFChargedHadrCands_.push_back(candPtrs.at(i));
  }
  embeddedSignalPFChargedHadrCands_ = true;
}
void Tau::embedSignalPFNeutralHadrCands() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  std::vector<reco::PFCandidatePtr> candPtrs = pfSpecific_[0].selectedSignalPFNeutrHadrCands_;
  for (unsigned int i = 0; i < candPtrs.size(); i++) {
    signalPFNeutralHadrCands_.push_back(candPtrs.at(i));
  }
  embeddedSignalPFNeutralHadrCands_ = true;
}
void Tau::embedSignalPFGammaCands() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  std::vector<reco::PFCandidatePtr> candPtrs = pfSpecific_[0].selectedSignalPFGammaCands_;
  for (unsigned int i = 0; i < candPtrs.size(); i++) {
    signalPFGammaCands_.push_back(candPtrs.at(i));
  }
  embeddedSignalPFGammaCands_ = true;
}

void Tau::embedIsolationPFCands() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  std::vector<reco::PFCandidatePtr> candPtrs = pfSpecific_[0].selectedIsolationPFCands_;
  for (unsigned int i = 0; i < candPtrs.size(); i++) {
    isolationPFCands_.push_back(candPtrs.at(i));
  }
  embeddedIsolationPFCands_ = true;
}

void Tau::embedIsolationPFChargedHadrCands() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  std::vector<reco::PFCandidatePtr> candPtrs = pfSpecific_[0].selectedIsolationPFChargedHadrCands_;
  for (unsigned int i = 0; i < candPtrs.size(); i++) {
    isolationPFChargedHadrCands_.push_back(candPtrs.at(i));
  }
  embeddedIsolationPFChargedHadrCands_ = true;
}
void Tau::embedIsolationPFNeutralHadrCands() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  std::vector<reco::PFCandidatePtr> candPtrs = pfSpecific_[0].selectedIsolationPFNeutrHadrCands_;
  for (unsigned int i = 0; i < candPtrs.size(); i++) {
    isolationPFNeutralHadrCands_.push_back(candPtrs.at(i));
  }
  embeddedIsolationPFNeutralHadrCands_ = true;
}
void Tau::embedIsolationPFGammaCands() {
  if (!isPFTau() ) {//additional check with warning in pat::tau producer
    return;
  }
  std::vector<reco::PFCandidatePtr> candPtrs = pfSpecific_[0].selectedIsolationPFGammaCands_;
  for (unsigned int i = 0; i < candPtrs.size(); i++) {
    isolationPFGammaCands_.push_back(candPtrs.at(i));
  }
  embeddedIsolationPFGammaCands_ = true;
}

reco::PFRecoTauChargedHadronRef Tau::leadTauChargedHadronCandidate() const {
  if(!isPFTau() ) throw cms::Exception("Type Error") << "Requesting content that is not stored in miniAOD.\n";
  if ( !pfSpecific().signalTauChargedHadronCandidates_.empty() ) {
    return reco::PFRecoTauChargedHadronRef(&pfSpecific().signalTauChargedHadronCandidates_,0);
  } else {
    return reco::PFRecoTauChargedHadronRef();
  }
}

const reco::PFCandidatePtr convertToPFCandidatePtr(const reco::CandidatePtr& ptr) {
  const reco::PFCandidate* pf_cand = dynamic_cast<const reco::PFCandidate*>(&*ptr);
  if (pf_cand)
    return edm::Ptr<reco::PFCandidate>(ptr);
  return reco::PFCandidatePtr();
}

const reco::PFCandidatePtr Tau::leadPFChargedHadrCand() const { 
  if(!embeddedLeadPFChargedHadrCand_){
    if(pfSpecific_.empty()) return reco::PFCandidatePtr(); 
    else return convertToPFCandidatePtr(pfSpecific().leadPFChargedHadrCand_); 
  }else
    return reco::PFCandidatePtr(&leadPFChargedHadrCand_,0);
}

const reco::PFCandidatePtr Tau::leadPFNeutralCand() const { 
  if(!embeddedLeadPFNeutralCand_){
    if(pfSpecific_.empty()) return reco::PFCandidatePtr();
    else return convertToPFCandidatePtr(pfSpecific().leadPFNeutralCand_);
  }else
    return reco::PFCandidatePtr(&leadPFNeutralCand_,0);
}

const reco::PFCandidatePtr Tau::leadPFCand() const { 
  if(!embeddedLeadPFCand_){
    if(pfSpecific_.empty()) return reco::PFCandidatePtr();
    return convertToPFCandidatePtr(pfSpecific().leadPFCand_);
  }else
    return reco::PFCandidatePtr(&leadPFCand_,0);
}

const std::vector<reco::PFCandidatePtr>& Tau::signalPFCands() const { 
  if (embeddedSignalPFCands_) {
   if (!signalPFCandsTransientPtrs_.isSet()) {
      std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
      aPtrs->reserve(signalPFCands_.size());
      for (unsigned int i = 0; i < signalPFCands_.size(); i++) {
    aPtrs->push_back(reco::PFCandidatePtr(&signalPFCands_, i) );
      }
      signalPFCandsTransientPtrs_.set(std::move(aPtrs));
    }
    return *signalPFCandsTransientPtrs_;
  } else {
    if(pfSpecific_.empty() || pfSpecific().selectedSignalPFCands_.empty() || !pfSpecific().selectedSignalPFCands_.front().isAvailable()){
      // this part of code is called when reading from patTuple or miniAOD
      // it returns empty collection in correct format so it can be substituted by reco::Candidates if available
      if (!signalPFCandsTransientPtrs_.isSet()) {
        std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
        signalPFCandsTransientPtrs_.set(std::move(aPtrs));
      }
      return *signalPFCandsTransientPtrs_;
    } else return pfSpecific().selectedSignalPFCands_;
  }
}

const std::vector<reco::PFCandidatePtr>& Tau::signalPFChargedHadrCands() const {
  if (embeddedSignalPFChargedHadrCands_) {
   if (!signalPFChargedHadrCandsTransientPtrs_.isSet()) {
      std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{ new std::vector<reco::PFCandidatePtr>{}};
      aPtrs->reserve(signalPFChargedHadrCands_.size());
      for (unsigned int i = 0; i < signalPFChargedHadrCands_.size(); i++) {
    aPtrs->push_back(reco::PFCandidatePtr(&signalPFChargedHadrCands_, i) );
      }
      signalPFChargedHadrCandsTransientPtrs_.set(std::move(aPtrs));
    }
    return *signalPFChargedHadrCandsTransientPtrs_;
  } else {
    if(pfSpecific_.empty() || pfSpecific().selectedSignalPFChargedHadrCands_.empty() || !pfSpecific().selectedSignalPFChargedHadrCands_.front().isAvailable()){
      // this part of code is called when reading from patTuple or miniAOD
      // it returns empty collection in correct format so it can be substituted by reco::Candidates if available
      if (!signalPFChargedHadrCandsTransientPtrs_.isSet()) {
        std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
        signalPFChargedHadrCandsTransientPtrs_.set(std::move(aPtrs));
      }
      return *signalPFChargedHadrCandsTransientPtrs_;
    } else return pfSpecific().selectedSignalPFChargedHadrCands_;
  }
} 

const std::vector<reco::PFCandidatePtr>& Tau::signalPFNeutrHadrCands() const {
  if (embeddedSignalPFNeutralHadrCands_) {
   if (!signalPFNeutralHadrCandsTransientPtrs_.isSet()) {
      std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
      aPtrs->reserve(signalPFNeutralHadrCands_.size());
      for (unsigned int i = 0; i < signalPFNeutralHadrCands_.size(); i++) {
    aPtrs->push_back(reco::PFCandidatePtr(&signalPFNeutralHadrCands_, i) );
      }
      signalPFNeutralHadrCandsTransientPtrs_.set(std::move(aPtrs));
    }
    return *signalPFNeutralHadrCandsTransientPtrs_;
  } else {
    if(pfSpecific_.empty() || pfSpecific().selectedSignalPFNeutrHadrCands_.empty() || !pfSpecific().selectedSignalPFNeutrHadrCands_.front().isAvailable()){
      // this part of code is called when reading from patTuple or miniAOD
      // it returns empty collection in correct format so it can be substituted by reco::Candidates if available
      if (!signalPFNeutralHadrCandsTransientPtrs_.isSet()) {
        std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
        signalPFNeutralHadrCandsTransientPtrs_.set(std::move(aPtrs));
      }
      return *signalPFNeutralHadrCandsTransientPtrs_;
    } else return pfSpecific().selectedSignalPFNeutrHadrCands_;
  }
} 

const std::vector<reco::PFCandidatePtr>& Tau::signalPFGammaCands() const {
  if (embeddedSignalPFGammaCands_) {
   if (!signalPFGammaCandsTransientPtrs_.isSet()) {
      std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{ new std::vector<reco::PFCandidatePtr>{}};
      aPtrs->reserve(signalPFGammaCands_.size());
      for (unsigned int i = 0; i < signalPFGammaCands_.size(); i++) {
    aPtrs->push_back(reco::PFCandidatePtr(&signalPFGammaCands_, i) );
      }
      signalPFGammaCandsTransientPtrs_.set(std::move(aPtrs));
    }
    return *signalPFGammaCandsTransientPtrs_;
  } else {
    if(pfSpecific_.empty() || pfSpecific().selectedSignalPFGammaCands_.empty() || !pfSpecific().selectedSignalPFGammaCands_.front().isAvailable()){
      // this part of code is called when reading from patTuple or miniAOD
      // it returns empty collection in correct format so it can be substituted by reco::Candidates if available
      if (!signalPFGammaCandsTransientPtrs_.isSet()) {
        std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
        signalPFGammaCandsTransientPtrs_.set(std::move(aPtrs));
      }
      return *signalPFGammaCandsTransientPtrs_;
    } else return pfSpecific().selectedSignalPFGammaCands_;
  }
}

const std::vector<reco::PFRecoTauChargedHadron> & Tau::signalTauChargedHadronCandidates() const {
  if(pfSpecific_.empty()) throw cms::Exception("Type Error") << "Requesting content that is not stored in miniAOD.\n";
  return pfSpecific().signalTauChargedHadronCandidates_;
}

const std::vector<reco::RecoTauPiZero> & Tau::signalPiZeroCandidates() const {
  if(pfSpecific_.empty()) throw cms::Exception("Type Error") << "Requesting content that is not stored in miniAOD.\n";
  return pfSpecific().signalPiZeroCandidates_;
}

const std::vector<reco::PFCandidatePtr>& Tau::isolationPFCands() const {
  if (embeddedIsolationPFCands_) {
  if (!isolationPFCandsTransientPtrs_.isSet()) {
      std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{ new std::vector<reco::PFCandidatePtr>{}};
      aPtrs->reserve(isolationPFCands_.size());
      for (unsigned int i = 0; i < isolationPFCands_.size(); i++) {
    aPtrs->push_back(reco::PFCandidatePtr(&isolationPFCands_, i) );
      }
      isolationPFCandsTransientPtrs_.set(std::move(aPtrs));
    }
    return *isolationPFCandsTransientPtrs_;
  } else {
    if(pfSpecific_.empty() || pfSpecific().selectedIsolationPFCands_.empty() || !pfSpecific().selectedIsolationPFCands_.front().isAvailable()){
      // this part of code is called when reading from patTuple or miniAOD
      // it returns empty collection in correct format so it can be substituted by reco::Candidates if available
      if (!isolationPFCandsTransientPtrs_.isSet()) {
        std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
        isolationPFCandsTransientPtrs_.set(std::move(aPtrs));
      }
      return *isolationPFCandsTransientPtrs_;
    } else return pfSpecific().selectedIsolationPFCands_;
  }
} 

const std::vector<reco::PFCandidatePtr>& Tau::isolationPFChargedHadrCands() const {
  if (embeddedIsolationPFChargedHadrCands_) {
   if (!isolationPFChargedHadrCandsTransientPtrs_.isSet()) {
      std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{ new std::vector<reco::PFCandidatePtr>{}};
      aPtrs->reserve(isolationPFChargedHadrCands_.size());
      for (unsigned int i = 0; i < isolationPFChargedHadrCands_.size(); i++) {
    aPtrs->push_back(reco::PFCandidatePtr(&isolationPFChargedHadrCands_, i) );
      }
      isolationPFChargedHadrCandsTransientPtrs_.set(std::move(aPtrs));
    }
    return *isolationPFChargedHadrCandsTransientPtrs_;
  } else {
    if(pfSpecific_.empty() || pfSpecific().selectedIsolationPFChargedHadrCands_.empty() || !pfSpecific().selectedIsolationPFChargedHadrCands_.front().isAvailable()){
      // this part of code is called when reading from patTuple or miniAOD
      // it returns empty collection in correct format so it can be substituted by reco::Candidates if available
      if (!isolationPFChargedHadrCandsTransientPtrs_.isSet()) {
        std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
        isolationPFChargedHadrCandsTransientPtrs_.set(std::move(aPtrs));
      }
      return *isolationPFChargedHadrCandsTransientPtrs_;
    } else return pfSpecific().selectedIsolationPFChargedHadrCands_;
  }
}

const std::vector<reco::PFCandidatePtr>& Tau::isolationPFNeutrHadrCands() const {
  if (embeddedIsolationPFNeutralHadrCands_) {
    if (!isolationPFNeutralHadrCandsTransientPtrs_.isSet()) {
      std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{ new std::vector<reco::PFCandidatePtr>{}};
      aPtrs->reserve(isolationPFNeutralHadrCands_.size());
      for (unsigned int i = 0; i < isolationPFNeutralHadrCands_.size(); i++) {
    aPtrs->push_back(reco::PFCandidatePtr(&isolationPFNeutralHadrCands_, i) );
      }
      isolationPFNeutralHadrCandsTransientPtrs_.set(std::move(aPtrs));
    }
    return *isolationPFNeutralHadrCandsTransientPtrs_;
  } else {
    if(pfSpecific_.empty() || pfSpecific().selectedIsolationPFNeutrHadrCands_.empty() || !pfSpecific().selectedIsolationPFNeutrHadrCands_.front().isAvailable()){
      // this part of code is called when reading from patTuple or miniAOD
      // it returns empty collection in correct format so it can be substituted by reco::Candidates if available
      if (!isolationPFNeutralHadrCandsTransientPtrs_.isSet()) {
        std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
        isolationPFNeutralHadrCandsTransientPtrs_.set(std::move(aPtrs));
      }
      return *isolationPFNeutralHadrCandsTransientPtrs_;
    } else return pfSpecific().selectedIsolationPFNeutrHadrCands_;
  }
} 

const std::vector<reco::PFCandidatePtr>& Tau::isolationPFGammaCands() const {
  if (embeddedIsolationPFGammaCands_) {
    if (!isolationPFGammaCandsTransientPtrs_.isSet()) {
      std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
      aPtrs->reserve(isolationPFGammaCands_.size());
      for (unsigned int i = 0; i < isolationPFGammaCands_.size(); i++) {
    aPtrs->push_back(reco::PFCandidatePtr(&isolationPFGammaCands_, i) );
      }
      isolationPFGammaCandsTransientPtrs_.set(std::move(aPtrs));
    }
    return *isolationPFGammaCandsTransientPtrs_;
  } else {
    if(pfSpecific_.empty() || pfSpecific().selectedIsolationPFGammaCands_.empty() || !pfSpecific().selectedIsolationPFGammaCands_.front().isAvailable()){
      // this part of code is called when reading from patTuple or miniAOD
      // it returns empty collection in correct format so it can be substituted by reco::Candidates if available
      if (!isolationPFGammaCandsTransientPtrs_.isSet()) {
        std::unique_ptr<std::vector<reco::PFCandidatePtr> > aPtrs{new std::vector<reco::PFCandidatePtr>{}};
        isolationPFGammaCandsTransientPtrs_.set(std::move(aPtrs));
      }
      return *isolationPFGammaCandsTransientPtrs_;
    } else return pfSpecific().selectedIsolationPFGammaCands_;
  }
}

const std::vector<reco::PFRecoTauChargedHadron> & Tau::isolationTauChargedHadronCandidates() const {
  if(pfSpecific_.empty()) throw cms::Exception("Type Error") << "Requesting content that is not stored in miniAOD.\n";
  return pfSpecific().isolationTauChargedHadronCandidates_;
}

const std::vector<reco::RecoTauPiZero> & Tau::isolationPiZeroCandidates() const {
  if(pfSpecific_.empty()) throw cms::Exception("Type Error") << "Requesting content that is not stored in miniAOD.\n";
  return pfSpecific().isolationPiZeroCandidates_;
}


// initialize the jet to a given JEC level during creation starting from Uncorrected
void Tau::initializeJEC(unsigned int level, unsigned int set)
{
  currentJECSet(set);
  currentJECLevel(level);
  setP4(jec_[set].correction(level)*p4());
}

int Tau::jecSet(const std::string& set) const
{
  for ( std::vector<pat::TauJetCorrFactors>::const_iterator corrFactor = jec_.begin(); 
    corrFactor != jec_.end(); ++corrFactor ) {
    if ( corrFactor->jecSet() == set ) return corrFactor-jec_.begin(); 
  }
  return -1;
}

const std::vector<std::string> Tau::availableJECSets() const
{
  std::vector<std::string> sets;
  for ( std::vector<pat::TauJetCorrFactors>::const_iterator corrFactor = jec_.begin(); 
    corrFactor != jec_.end(); ++corrFactor ) {
    sets.push_back(corrFactor->jecSet());
  }
  return sets;
}

const std::vector<std::string> Tau::availableJECLevels(const int& set) const
{
  return set>=0 ? jec_.at(set).correctionLabels() : std::vector<std::string>();
}

float Tau::jecFactor(const std::string& level, const std::string& set) const
{
  for ( unsigned int idx = 0; idx < jec_.size(); ++idx ) {
    if ( set.empty() || jec_.at(idx).jecSet() == set ){
      if ( jec_[idx].jecLevel(level) >= 0 ) 
    return jecFactor(jec_[idx].jecLevel(level), idx);
      else
    throw cms::Exception("InvalidRequest") 
      << "This JEC level " << level << " does not exist. \n";
    }
  }
  throw cms::Exception("InvalidRequest") 
    << "This jet does not carry any jet energy correction factor information \n"
    << "for a jet energy correction set with label " << set << "\n";
}

float Tau::jecFactor(const unsigned int& level, const unsigned int& set) const
{
  if ( !jecSetsAvailable() )
    throw cms::Exception("InvalidRequest") 
      << "This jet does not carry any jet energy correction factor information \n";
  if ( !jecSetAvailable(set) )
    throw cms::Exception("InvalidRequest") 
      << "This jet does not carry any jet energy correction factor information \n"
      << "for a jet energy correction set with index " << set << "\n";
  return jec_.at(set).correction(level)/jec_.at(currentJECSet_).correction(currentJECLevel_);
}

Tau Tau::correctedTauJet(const std::string& level, const std::string& set) const
{
  // rescale p4 of the jet; the update of current values is
  // done within the called jecFactor function
  for ( unsigned int idx = 0; idx < jec_.size(); ++idx ) {
    if ( set.empty() || jec_.at(idx).jecSet() == set ) {
      if ( jec_[idx].jecLevel(level) >= 0 ) 
    return correctedTauJet(jec_[idx].jecLevel(level), idx);
      else
    throw cms::Exception("InvalidRequest") 
      << "This JEC level " << level << " does not exist. \n";
    }
  }
  throw cms::Exception("InvalidRequest") 
    << "This JEC set " << set << " does not exist. \n";
}

Tau Tau::correctedTauJet(const unsigned int& level, const unsigned int& set) const
{
  Tau correctedTauJet(*this);
  //rescale p4 of the jet
  correctedTauJet.setP4(jecFactor(level, set)*p4());
  // update current level and set
  correctedTauJet.currentJECSet(set); 
  correctedTauJet.currentJECLevel(level); 
  return correctedTauJet;
}



// return the leading candidate from signal(PF)ChargedHadrCandPtrs_ collection
const reco::CandidatePtr Tau::leadChargedHadrCand() const {
    const reco::PFCandidatePtr leadPF = leadPFChargedHadrCand();
    if (leadPF.isAvailable() || signalChargedHadrCandPtrs_.isNull()) return leadPF;
    reco::CandidatePtr ret;
    for (const reco::CandidatePtr & p : signalChargedHadrCandPtrs_) {
        if (ret.isNull() || (p->pt() > ret->pt())) ret = p;
    }
    return ret;

}

const reco::CandidatePtr Tau::leadNeutralCand() const {
    const reco::PFCandidatePtr leadPF = leadPFNeutralCand();
    if (leadPF.isAvailable() || signalNeutralHadrCandPtrs_.isNull()) return leadPF;
    reco::CandidatePtr ret;
    for (const reco::CandidatePtr & p : signalNeutralHadrCandPtrs_) {
        if (ret.isNull() || (p->pt() > ret->pt())) ret = p;
    }
    return ret;

}

const reco::CandidatePtr Tau::leadCand() const {
    const reco::PFCandidatePtr leadPF = leadPFCand();
    if (leadPF.isAvailable() || !Tau::ExistSignalCands()) return leadPF;
    else return Tau::signalCands()[0];
}


bool Tau::ExistSignalCands() const {
  return !(signalChargedHadrCandPtrs_.isNull() && signalNeutralHadrCandPtrs_.isNull() && signalGammaCandPtrs_.isNull());
}

bool Tau::ExistIsolationCands() const {
  return !(isolationChargedHadrCandPtrs_.isNull() && isolationNeutralHadrCandPtrs_.isNull() && isolationGammaCandPtrs_.isNull());
}


reco::CandidatePtrVector Tau::signalCands() const {
    std::vector<reco::PFCandidatePtr> r0 = signalPFCands();
    reco::CandidatePtrVector ret;
    if(!Tau::ExistSignalCands() || (!r0.empty() && r0.front().isAvailable())) {
      for (const auto & p : r0) ret.push_back(p);
      return ret;
    } else {
      reco::CandidatePtrVector ret2;
      std::vector<std::pair<float,size_t> > pt_index;
      size_t index=0;
      for (const auto & p : signalChargedHadrCandPtrs_){ ret2.push_back(p); pt_index.push_back(std::make_pair(p->pt(),index)); index++;}
      for (const auto & p : signalNeutralHadrCandPtrs_){ ret2.push_back(p); pt_index.push_back(std::make_pair(p->pt(),index)); index++;}
      for (const auto & p : signalGammaCandPtrs_){ ret2.push_back(p); pt_index.push_back(std::make_pair(p->pt(),index)); index++;}
      std::sort(pt_index.begin(),pt_index.end(),
            boost::bind(&std::pair<float,size_t>::first,_1) >
        boost::bind(&std::pair<float,size_t>::first,_2));
      for( const auto & p : pt_index){
    ret.push_back(ret2[p.second]);
      }
      return ret;
    }
}

reco::CandidatePtrVector Tau::signalChargedHadrCands() const {
   std::vector<reco::PFCandidatePtr> r0 = signalPFChargedHadrCands();
   if(signalChargedHadrCandPtrs_.isNull() || (!r0.empty() && r0.front().isAvailable())) {
        reco::CandidatePtrVector ret;
        for (const auto & p : r0) ret.push_back(p);
    return ret;
   } else {
        return signalChargedHadrCandPtrs_;
   }
}


reco::CandidatePtrVector Tau::signalNeutrHadrCands() const {
   std::vector<reco::PFCandidatePtr> r0 = signalPFNeutrHadrCands();
   if(signalNeutralHadrCandPtrs_.isNull() || (!r0.empty() && r0.front().isAvailable())) {
        reco::CandidatePtrVector ret;
        for (const auto & p : r0) ret.push_back(p);
    return ret;
   } else {
        return signalNeutralHadrCandPtrs_;
   }
}

reco::CandidatePtrVector Tau::signalGammaCands() const {
   std::vector<reco::PFCandidatePtr> r0 = signalPFGammaCands();
   if(signalGammaCandPtrs_.isNull() || (!r0.empty() && r0.front().isAvailable())) {
        reco::CandidatePtrVector ret;
        for (const auto & p : r0) ret.push_back(p);
    return ret;
   } else {
        return signalGammaCandPtrs_;
   }
}

reco::CandidatePtrVector Tau::isolationCands() const {
  std::vector<reco::PFCandidatePtr> r0 = isolationPFCands();
  reco::CandidatePtrVector ret;
  if(!Tau::ExistIsolationCands() || (!r0.empty() && r0.front().isAvailable())) {
    for (const auto & p : r0) ret.push_back(p);
    return ret;
  } else {
    reco::CandidatePtrVector ret2;
    std::vector<std::pair<float,size_t> > pt_index;
    size_t index=0;
    for (const auto & p : isolationChargedHadrCandPtrs_){ ret2.push_back(p); pt_index.push_back(std::make_pair(p->pt(),index)); index++;}
    for (const auto & p : isolationNeutralHadrCandPtrs_){ ret2.push_back(p); pt_index.push_back(std::make_pair(p->pt(),index)); index++;}
    for (const auto & p : isolationGammaCandPtrs_){ ret2.push_back(p); pt_index.push_back(std::make_pair(p->pt(),index)); index++;}
    std::sort(pt_index.begin(),pt_index.end(),
          boost::bind(&std::pair<float,size_t>::first,_1) >
          boost::bind(&std::pair<float,size_t>::first,_2));
    for( const auto & p : pt_index){
      ret.push_back(ret2[p.second]);
    }
    return ret;
  }
}

reco::CandidatePtrVector Tau::isolationChargedHadrCands() const {
   std::vector<reco::PFCandidatePtr> r0 = isolationPFChargedHadrCands();
   if(isolationChargedHadrCandPtrs_.isNull() || (!r0.empty() && r0.front().isAvailable())) {
        reco::CandidatePtrVector ret;
        for (const auto & p : r0) ret.push_back(p);
    return ret;
   } else {
        return isolationChargedHadrCandPtrs_;
   }
}

reco::CandidatePtrVector Tau::isolationNeutrHadrCands() const {
   reco::CandidatePtrVector ret;
   std::vector<reco::PFCandidatePtr> r0 = isolationPFNeutrHadrCands();
   if(isolationNeutralHadrCandPtrs_.isNull() || (!r0.empty() && r0.front().isAvailable())) {
        reco::CandidatePtrVector ret;
        for (const auto & p : r0) ret.push_back(p);
    return ret;
   } else {
        return isolationNeutralHadrCandPtrs_;
   }
}

reco::CandidatePtrVector Tau::isolationGammaCands() const {
   std::vector<reco::PFCandidatePtr> r0 = isolationPFGammaCands();
   if(isolationGammaCandPtrs_.isNull() || (!r0.empty() && r0.front().isAvailable())) {
        reco::CandidatePtrVector ret;
        for (const auto & p : r0) ret.push_back(p);
    return ret;
   } else {
        return isolationGammaCandPtrs_;
   }
}



size_t Tau::numberOfSourceCandidatePtrs() const {
  if (Tau::ExistSignalCands()) return Tau::signalCands().size();
  else if(pfSpecific_.empty()) return 0; 
  else return pfSpecific().selectedSignalPFCands_.size();
}
reco::CandidatePtr Tau::sourceCandidatePtr( size_type i ) const {
  if (Tau::ExistSignalCands()) return Tau::signalCands()[i];
  else if(pfSpecific_.empty()) return reco::CandidatePtr();
  else return pfSpecific().selectedSignalPFCands_[i];
}