PatUtils/src/CaloJetSelector.cc

0001
0002 #include "PhysicsTools/PatUtils/interface/CaloJetSelector.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004
0005 using pat::CaloJetSelector;
0006
0007 //______________________________________________________________________________
0008 const pat::ParticleStatus CaloJetSelector::filter(  //const unsigned int&        index,
0009                                                     // const edm::View<reco::CaloJet>& Jets
0010     const reco::CaloJet& Jet) const {
0011   ParticleStatus result = GOOD;
0012
0013   ///Retrieve information
0014   ///Pt Jet
0015   if (Jet.p4().Pt() < config_.Ptmin)
0016     result = BAD;
0017
0018   ///eta region
0019   double eta = fabs(Jet.p4().Eta());
0020   if (eta > config_.Etamax)
0021     result = BAD;
0022
0023   ///electromagnetic fraction
0024   double EMF = Jet.emEnergyFraction();
0025   if (EMF < config_.EMFmin || EMF > config_.EMFmax)
0026     result = BAD;
0027
0028   ///(EMCalEnergyFraction + HadCalEnergyFraction) / (EMCalEnergyFraction - HadCalEnergyFraction)
0029   double HadF = Jet.emEnergyFraction();
0030   double EMvsHadF = 0.;
0031   if (EMF - HadF != 0.)
0032     EMvsHadF = (EMF + HadF) / (EMF - HadF);
0033   if (EMvsHadF < config_.EMvsHadFmin || EMvsHadF > config_.EMvsHadFmax)
0034     result = BAD;
0035
0036   //ratio Energy over Momentum (both from calorimeter)
0037   //Useful? Both come from a lorentz-vector
0038   //double EoverP = 0.;
0039   //if (Jet.p4().P()!=0.) EoverP = Jet.energy() / Jet.p4().P();
0040   //if (EoverP > config_.EoverPmax) result = BAD;
0041
0042   ///n90: number of towers containing 90% of the jet's energy
0043   double n90 = Jet.n90();
0044   if (n90 < config_.N90min || n90 > config_.N90max)
0045     result = BAD;
0046
0047   ///Tower Number
0048   std::vector<CaloTowerPtr> jetTowers = Jet.getCaloConstituents();
0049   if (jetTowers.size() < config_.NCaloTowersmin || jetTowers.size() > config_.NCaloTowersmax)
0050     result = BAD;
0051
0052   //calculate tower related variables:
0053   double MaxEnergyTower = 0., SumTowPt = 0., SumTowPtR = 0.;
0054   for (std::vector<CaloTowerPtr>::const_iterator tow = jetTowers.begin(), towend = jetTowers.end(); tow != towend;
0055        ++tow) {
0056     SumTowPt += (*tow)->et();
0057     SumTowPtR += (*tow)->et() * deltaR(Jet.p4().Eta(), Jet.p4().Phi(), (*tow)->eta(), (*tow)->phi());
0058     if ((*tow)->et() > MaxEnergyTower)
0059       MaxEnergyTower = (*tow)->et();
0060   }
0061
0062   ///Highest Et Tower / Et Jet
0063   double EtTowerMaxOverEtJet = 0.;
0064   if (Jet.p4().Et() != 0.)
0065     EtTowerMaxOverEtJet = MaxEnergyTower / Jet.p4().Et();
0066   if (EtTowerMaxOverEtJet < config_.HighestTowerOverJetmin || EtTowerMaxOverEtJet > config_.HighestTowerOverJetmax)
0067     result = BAD;
0068
0069   ///Sum(E Twr * DeltaR(Twr-Jet)) / Sum(E Twr)
0070   double RWidth = 0.;
0071   if (SumTowPt != 0.)
0072     RWidth = SumTowPtR / SumTowPt;
0073   if (RWidth < config_.RWidthmin || RWidth > config_.RWidthmax)
0074     result = BAD;
0075
0076   ///Pt Jet / Towers Area
0077   double PtJetoverArea = 0.;
0078   if (Jet.towersArea() != 0.)
0079     PtJetoverArea = Jet.p4().Pt() / Jet.towersArea();
0080   if (PtJetoverArea < config_.PtJetOverArea_min || PtJetoverArea > config_.PtJetOverArea_max)
0081     result = BAD;
0082
0083   ///Highest Et Tower / Towers Area
0084   double PtToweroverArea = 0.;
0085   if (Jet.towersArea() != 0.)
0086     PtToweroverArea = MaxEnergyTower / Jet.towersArea();
0087   if (PtToweroverArea < config_.PtTowerOverArea_min || PtToweroverArea > config_.PtTowerOverArea_max)
0088     result = BAD;
0089
0090   return result;
0091 }