Go to the documentation of this file.00001 #include "Validation/RecoParticleFlow/plugins/PFJetFilter.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "FWCore/ServiceRegistry/interface/Service.h"
00004
00005 #include "FWCore/Framework/interface/Event.h"
00006 #include "DataFormats/Common/interface/Handle.h"
00007 #include "FWCore/Framework/interface/ESHandle.h"
00008 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00009 #include "DataFormats/Candidate/interface/Candidate.h"
00010 #include "DataFormats/JetReco/interface/PFJet.h"
00011
00012 using namespace reco;
00013 using namespace edm;
00014 using namespace std;
00015
00016 PFJetFilter::PFJetFilter(const edm::ParameterSet& iConfig)
00017 {
00018 inputTruthLabel_ = iConfig.getParameter<edm::InputTag>("InputTruthLabel");
00019 inputRecoLabel_ = iConfig.getParameter<edm::InputTag>("InputRecoLabel");
00020
00021 recPt_cut = iConfig.getParameter<double>("recPt");
00022 genPt_cut = iConfig.getParameter<double>("genPt");
00023
00024 eta_min = iConfig.getParameter<double>("minEta");
00025 eta_max = iConfig.getParameter<double>("maxEta");
00026
00027 deltaR_min = iConfig.getParameter<double>("deltaRMin");
00028 deltaR_max = iConfig.getParameter<double>("deltaRMax");
00029
00030 deltaEt_min = iConfig.getParameter<double>("minDeltaEt");
00031 deltaEt_max = iConfig.getParameter<double>("maxDeltaEt");
00032
00033 verbose = iConfig.getParameter<bool>("verbose");
00034
00035 entry = 0;
00036 }
00037
00038 PFJetFilter::~PFJetFilter()
00039 {}
00040
00041
00042 void
00043 PFJetFilter::beginJob() {}
00044
00045 void
00046 PFJetFilter::endJob() {}
00047
00048 bool
00049 PFJetFilter::filter(edm::Event& iEvent, const edm::EventSetup& iESetup)
00050 {
00051
00052
00053 typedef edm::View<reco::Candidate> candidateCollection ;
00054 typedef edm::View<reco::Candidate> candidateCollection ;
00055
00056 const candidateCollection *truth_candidates;
00057 const candidateCollection *reco_candidates;
00058
00059
00060
00061
00062
00063
00064 Handle<candidateCollection> truth_hnd;
00065 bool isGen = iEvent.getByLabel(inputTruthLabel_, truth_hnd);
00066 if ( !isGen ) {
00067 std::cout << "Warning : no Gen jets in input !" << std::endl;
00068 return false;
00069 }
00070
00071 truth_candidates = truth_hnd.product();
00072
00073
00074 Handle<candidateCollection> reco_hnd;
00075 bool isReco = iEvent.getByLabel(inputRecoLabel_, reco_hnd);
00076 if ( !isReco ) {
00077 std::cout << "Warning : no Reco jets in input !" << std::endl;
00078 return false;
00079 }
00080 reco_candidates = reco_hnd.product();
00081 if (!truth_candidates || !reco_candidates) return false;
00082
00083
00084 bool pass = false;
00085
00086 for (unsigned int i = 0; i < reco_candidates->size(); i++) {
00087
00088 const reco::Candidate *particle = &(*reco_candidates)[i];
00089
00090
00091 assert( particle!=NULL );
00092
00093 double rec_pt = particle->pt();
00094 double rec_eta = particle->eta();
00095 double rec_phi = particle->phi();
00096
00097
00098 if (rec_pt < recPt_cut ) continue;
00099
00100
00101 if ( fabs(rec_eta) > eta_max ) continue;
00102 if ( fabs(rec_eta) < eta_min ) continue;
00103
00104 bool Barrel = false;
00105 bool Endcap = false;
00106 bool Forward = false;
00107 if (std::abs(rec_eta) < 1.4 ) Barrel = true;
00108 if (std::abs (rec_eta) > 1.4 && std::abs (rec_eta) < 2.6 ) Endcap = true;
00109 if (std::abs (rec_eta) > 3.5 && std::abs (rec_eta) < 4.5 ) Forward = true;
00110
00111
00112 if ( !Barrel && !Endcap ) continue;
00113
00114
00115 double deltaRmin = 999.;
00116 double ptmin = 0.;
00117 for (unsigned int j = 0; j < reco_candidates->size(); j++) {
00118
00119 if ( i == j ) continue;
00120 const reco::Candidate *other = &(*reco_candidates)[j];
00121 double deltaR = algo_->deltaR(particle, other);
00122 if ( deltaR < deltaRmin &&
00123 other->pt() > 0.25*particle->pt() &&
00124 other->pt() > recPt_cut ) {
00125 deltaRmin = deltaR;
00126 ptmin = other->pt();
00127 }
00128 if ( deltaRmin < deltaR_min ) break;
00129 }
00130 if ( deltaRmin < deltaR_min ) continue;
00131
00132
00133
00134 const reco::Candidate *gen_particle =
00135 algo_->matchByDeltaR(particle,
00136 truth_candidates);
00137
00138
00139 if(gen_particle==NULL) continue;
00140
00141
00142 double deltaR = algo_->deltaR(particle, gen_particle);
00143 if ( deltaR > deltaR_max ) continue;
00144
00145
00146 double true_pt = gen_particle->pt();
00147 double true_eta = gen_particle->eta();
00148 double true_phi = gen_particle->phi();
00149
00150
00151 if (true_pt < genPt_cut ) continue;
00152
00153
00154 double resPt = (rec_pt -true_pt)/true_pt;
00155 double sigma = resolution(true_pt,Barrel);
00156 double avera = response(true_pt,Barrel);
00157 double nSig = (resPt-avera)/sigma;
00158
00159 if ( nSig > deltaEt_max || nSig < deltaEt_min ) {
00160
00161 if ( verbose )
00162 std::cout << "Entry " << entry
00163 << " resPt = " << resPt
00164 << " sigma/avera/nSig = " << sigma << "/" << avera <<"/" << nSig
00165 << " pT (T/R) " << true_pt << "/" << rec_pt
00166 << " Eta (T/R) " << true_eta << "/" << rec_eta
00167 << " Phi (T/R) " << true_phi << "/" << rec_phi
00168 << " deltaRMin/ptmin " << deltaRmin << "/" << ptmin
00169 << std::endl;
00170
00171 pass = true;
00172 }
00173
00174 if ( pass ) break;
00175
00176 }
00177
00178 entry++;
00179 return pass;
00180
00181 }
00182
00183 double
00184 PFJetFilter::resolution(double pt, bool barrel) {
00185
00186 double p0 = barrel ? 1.19200e-02 : 8.45341e-03;
00187 double p1 = barrel ? 1.06138e+00 : 7.96855e-01;
00188 double p2 = barrel ? -2.05929e+00 : -3.12076e-01;
00189
00190 double resp = p0 + p1/sqrt(pt) + p2/pt;
00191 return resp;
00192
00193 }
00194
00195 double
00196 PFJetFilter::response(double pt, bool barrel) {
00197
00198 double p0 = barrel ? 1.09906E-1 : 6.91939E+1;
00199 double p1 = barrel ? -1.61443E-1 : -6.92733E+1;
00200 double p2 = barrel ? 3.45489E+3 : 1.58207E+6;
00201
00202 double resp = p0 + p1 * std::exp(-pt/p2);
00203 return resp;
00204
00205 }