Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
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;
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 }}
00195
00196 #include "FWCore/Framework/interface/MakerMacros.h"
00197 DEFINE_EDM_PLUGIN(RecoTauModifierPluginFactory,
00198 reco::tau::RecoTauEnergyRecoveryPlugin,
00199 "RecoTauEnergyRecoveryPlugin");