CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/Validation/RecoParticleFlow/plugins/PFJetFilter.cc

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 // Typedefs to use views
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   // Retrieve!
00061   // ==========================================================
00062 
00063  // Get Truth Candidates (GenCandidates, GenJets, etc.)
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    // Get Reco Candidates (PFlow, CaloJet, etc.)
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       // This protection is meant at not being used !
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       // skip PFjets with pt < recPt_cut GeV
00098       if (rec_pt < recPt_cut ) continue;
00099 
00100       // skip PFjets with eta > maxEta_cut or eta < minEta_cut
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       // Keep only jets in the barrel or the endcaps, within the tracker acceptance
00112       if ( !Barrel && !Endcap ) continue;
00113 
00114       // Find the closets recJet
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       // Find the closest genJet.
00134       const reco::Candidate *gen_particle = 
00135         algo_->matchByDeltaR(particle,
00136                              truth_candidates);
00137 
00138       // Check there is a genJet associated to the recoJet
00139       if(gen_particle==NULL) continue; 
00140 
00141       // check deltaR is small enough
00142       double deltaR = algo_->deltaR(particle, gen_particle);
00143       if ( deltaR > deltaR_max ) continue;
00144 
00145       // double true_E = gen_particle->p();
00146       double true_pt = gen_particle->pt();
00147       double true_eta = gen_particle->eta();
00148       double true_phi = gen_particle->phi();
00149 
00150       // skip PFjets with pt < genPt_cut GeV
00151       if (true_pt < genPt_cut ) continue;
00152       
00153       // Find the pT residual
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 }