CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_10/src/RecoTauTag/RecoTau/plugins/RecoTauEnergyRecoveryPlugin.cc

Go to the documentation of this file.
00001 /*
00002  * =============================================================================
00003  *       Filename:  RecoTauEnergyRecoveryPlugin.cc
00004  *
00005  *    Description:  Recovery energy of visible decay products
00006  *                  lost due to photon conversions/nuclear interactions
00007  *                  in tracker material.
00008  *        Created:  09/02/2011 10:28:00
00009  *
00010  *         Authors:  Christian Veelken (LLR)
00011  *
00012  * =============================================================================
00013  */
00014 
00015 #include "RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h"
00016 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
00017 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
00018 #include "DataFormats/JetReco/interface/PFJet.h"
00019 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00020 #include "DataFormats/TrackReco/interface/Track.h"
00021 #include "DataFormats/VertexReco/interface/Vertex.h"
00022 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00023 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
00024 
00025 #include "FWCore/Framework/interface/ESHandle.h"
00026 
00027 #include <algorithm>
00028 
00029 namespace reco { namespace tau {
00030 
00031 class RecoTauEnergyRecoveryPlugin : public RecoTauModifierPlugin
00032 {
00033  public:
00034 
00035   explicit RecoTauEnergyRecoveryPlugin(const edm::ParameterSet&);
00036   virtual ~RecoTauEnergyRecoveryPlugin();
00037   void operator()(PFTau&) const;
00038   virtual void beginEvent();
00039 
00040  private:
00041 
00042   edm::InputTag srcVertices_;
00043   const reco::Vertex* theEventVertex_;
00044 
00045   RecoTauQualityCuts* qcuts_;
00046 
00047   unsigned corrLevel_;
00048 
00049   double lev1PhiWindow_;
00050   double lev1EtaWindow_;
00051 };
00052 
00053 RecoTauEnergyRecoveryPlugin::RecoTauEnergyRecoveryPlugin(const edm::ParameterSet& cfg)
00054   : RecoTauModifierPlugin(cfg),
00055     theEventVertex_(0),
00056     corrLevel_(cfg.getParameter<unsigned>("corrLevel")),
00057     lev1PhiWindow_(cfg.getParameter<double>("lev1PhiWindow")),
00058     lev1EtaWindow_(cfg.getParameter<double>("lev1EtaWindow"))
00059 {
00060   edm::ParameterSet cfgQualityCuts = cfg.getParameter<edm::ParameterSet>("qualityCuts");
00061   srcVertices_ = cfgQualityCuts.getParameter<edm::InputTag>("primaryVertexSrc");
00062   qcuts_ = new RecoTauQualityCuts(cfgQualityCuts.getParameter<edm::ParameterSet>("signalQualityCuts"));
00063 }
00064 
00065 RecoTauEnergyRecoveryPlugin::~RecoTauEnergyRecoveryPlugin()
00066 {
00067   delete qcuts_;
00068 }
00069 
00070 void RecoTauEnergyRecoveryPlugin::beginEvent()
00071 {
00072   edm::Handle<reco::VertexCollection> vertices;
00073   evt()->getByLabel(srcVertices_, vertices);
00074 
00075   if ( vertices->size() >= 1 ) {
00076     qcuts_->setPV(reco::VertexRef(vertices, 0));
00077     theEventVertex_ = &vertices->at(0);
00078   } else {
00079     theEventVertex_ = 0;
00080   }
00081 }
00082 
00083 const reco::TrackBaseRef getTrack(const reco::PFCandidate& cand)
00084 {
00085   if      ( cand.trackRef().isNonnull()    ) return reco::TrackBaseRef(cand.trackRef());
00086   else if ( cand.gsfTrackRef().isNonnull() ) return reco::TrackBaseRef(cand.gsfTrackRef());
00087   else return reco::TrackBaseRef();
00088 }
00089 
00090 bool isTauSignalPFCandidate(const reco::PFTau& tau, const reco::PFCandidatePtr& pfJetConstituent)
00091 {
00092   bool retVal = false;
00093 
00094   const reco::PFCandidateRefVector& signalPFCandidates = tau.signalPFCands();
00095   for ( reco::PFCandidateRefVector::const_iterator signalPFCandidate = signalPFCandidates.begin();
00096         signalPFCandidate != signalPFCandidates.end(); ++signalPFCandidate ) {
00097     if ( pfJetConstituent.key() == signalPFCandidate->key() ) retVal = true;
00098   }
00099 
00100   return retVal;
00101 }
00102 
00103 double square(double x)
00104 {
00105   return x*x;
00106 }
00107 
00108 void RecoTauEnergyRecoveryPlugin::operator()(PFTau& tau) const
00109 {
00110   double tauEnergyCorr = 0.;
00111 
00112   if ( corrLevel_ & 1 && theEventVertex_ ) {
00113 
00114     bool needsCorrLevel1 = false;
00115 
00116     std::vector<reco::PFCandidatePtr> pfJetConstituents = tau.jetRef()->getPFConstituents();
00117     for ( std::vector<reco::PFCandidatePtr>::const_iterator pfJetConstituent = pfJetConstituents.begin();
00118           pfJetConstituent != pfJetConstituents.end(); ++pfJetConstituent ) {
00119       const reco::TrackBaseRef track = getTrack(**pfJetConstituent);
00120       if ( track.isNonnull() ) {
00121         double dIP = fabs(track->dxy(theEventVertex_->position()));
00122         double dZ = fabs(track->dz(theEventVertex_->position()));
00123       if ( track->pt() > 2.0 && fabs(tau.eta() - (*pfJetConstituent)->eta()) < lev1EtaWindow_ &&
00124            !isTauSignalPFCandidate(tau, *pfJetConstituent) && (dZ < 0.2 || dIP > 0.10) ) needsCorrLevel1 = true;
00125       }
00126     }
00127 
00128     if ( needsCorrLevel1 ) {
00129       std::vector<reco::PFCandidatePtr> pfJetConstituents = tau.jetRef()->getPFConstituents();
00130       for ( std::vector<reco::PFCandidatePtr>::const_iterator pfJetConstituent = pfJetConstituents.begin();
00131             pfJetConstituent != pfJetConstituents.end(); ++pfJetConstituent ) {
00132         if ( fabs(tau.eta() - (*pfJetConstituent)->eta()) < lev1EtaWindow_ &&
00133              fabs(tau.phi() - (*pfJetConstituent)->phi()) < lev1PhiWindow_ ) {
00134           if ( (*pfJetConstituent)->particleId() == reco::PFCandidate::h0 ) {
00135             tauEnergyCorr += (*pfJetConstituent)->energy();
00136           } else {
00137             if ( !isTauSignalPFCandidate(tau, *pfJetConstituent) ) {
00138               double caloEn = (*pfJetConstituent)->ecalEnergy() + (*pfJetConstituent)->hcalEnergy();
00139               tauEnergyCorr += caloEn;
00140             }
00141           }
00142         }
00143       }
00144     }
00145   }
00146 
00147   if ( corrLevel_ & 2 && theEventVertex_ ) {
00148 
00149     double leadTrackMom = 0.;
00150     double leadTrackMomErr = 0.;
00151     double jetCaloEn = 0.;
00152 
00153     std::vector<reco::PFCandidatePtr> pfJetConstituents = tau.jetRef()->getPFConstituents();
00154     for ( std::vector<reco::PFCandidatePtr>::const_iterator pfJetConstituent = pfJetConstituents.begin();
00155           pfJetConstituent != pfJetConstituents.end(); ++pfJetConstituent ) {
00156       const reco::TrackBaseRef track = getTrack(**pfJetConstituent);
00157       if ( track.isNonnull() ) {
00158         double trackPt = track->pt();
00159         double trackPtErr = track->ptError();
00160         if ( qcuts_->filter(**pfJetConstituent) &&
00161              trackPtErr < (0.20*trackPt) && track->normalizedChi2() < 5.0 && track->hitPattern().numberOfValidPixelHits() >= 1 &&
00162              (trackPt - 3.*trackPtErr) > (*pfJetConstituent)->pt() && trackPt < (3.*tau.jetRef()->pt()) ) {
00163           if ( track->p() > leadTrackMom ) {
00164             leadTrackMom = track->p();
00165             leadTrackMomErr = leadTrackMom*(trackPtErr/trackPt);
00166           }
00167         }
00168       }
00169 
00170       double caloEn = (*pfJetConstituent)->ecalEnergy() + (*pfJetConstituent)->hcalEnergy();
00171       jetCaloEn += caloEn;
00172     }
00173 
00174     if ( leadTrackMom > tau.p() ) {
00175       const double chargedPionMass = 0.13957; // GeV
00176       double leadTrackEn = sqrt(square(leadTrackMom) + square(chargedPionMass));
00177       double jetCaloEnErr = 1.00*sqrt(std::max(jetCaloEn, leadTrackEn));
00178       double combEn = ((1./square(jetCaloEnErr))*jetCaloEn + (1./square(leadTrackMomErr))*leadTrackEn)/
00179                       ((1./square(jetCaloEnErr)) + (1./square(leadTrackMomErr)));
00180       tauEnergyCorr = std::max(tauEnergyCorr, combEn - tau.energy());
00181     }
00182   }
00183 
00184   if ( tau.energy() > 0. ) {
00185     double tauEnergy_corrected = tau.energy() + tauEnergyCorr;
00186     double scaleFactor = tauEnergy_corrected/tau.energy();
00187     double tauPx_corrected = scaleFactor*tau.px();
00188     double tauPy_corrected = scaleFactor*tau.py();
00189     double tauPz_corrected = scaleFactor*tau.pz();
00190     tau.setalternatLorentzVect(reco::Candidate::LorentzVector(tauPx_corrected, tauPy_corrected, tauPz_corrected, tauEnergy_corrected));
00191   }
00192 }
00193 
00194 }} // end namespace reco::tau
00195 
00196 #include "FWCore/Framework/interface/MakerMacros.h"
00197 DEFINE_EDM_PLUGIN(RecoTauModifierPluginFactory,
00198     reco::tau::RecoTauEnergyRecoveryPlugin,
00199     "RecoTauEnergyRecoveryPlugin");