RecoB/src/AcceptJet.cc

0001 #include "DQMOffline/RecoB/interface/AcceptJet.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/JetReco/interface/PFJet.h"
0005
0006 #include <iostream>
0007
0008 using namespace std;
0009
0010 AcceptJet::AcceptJet(const double& etaMin_,
0011                      const double& etaMax_,
0012                      const double& ptMin_,
0013                      const double& ptMax_,
0014                      const double& pMin_,
0015                      const double& pMax_,
0016                      const double& ratioMin_,
0017                      const double& ratioMax_,
0018                      const bool& doJetID_)
0019     : etaMin(etaMin_),
0020       etaMax(etaMax_),
0021       ptRecJetMin(ptMin_),
0022       ptRecJetMax(ptMax_),
0023       pRecJetMin(pMin_),
0024       pRecJetMax(pMax_),
0025       ratioMin(ratioMin_),
0026       ratioMax(ratioMax_),
0027       doJetID(doJetID_) {}
0028
0029 bool AcceptJet::operator()(const reco::Jet& jet,
0030                            const int& jetFlavour,
0031                            const edm::Handle<reco::SoftLeptonTagInfoCollection>& infos,
0032                            const double jec) const {
0033   // temporary fudge to correct for double loop error
0034   //  jetPartonMomentum /= 2.0;
0035
0036   if (fabs(jet.eta()) < etaMin || fabs(jet.eta()) > etaMax)
0037     return false;
0038
0039   //   if ( jetFlavour.underlyingParton4Vec().P() < pPartonMin  ||
0040   //        jetFlavour.underlyingParton4Vec().P() > pPartonMax  ) accept = false;
0041   //
0042   //   if ( jetFlavour.underlyingParton4Vec().Pt() < ptPartonMin  ||
0043   //        jetFlavour.underlyingParton4Vec().Pt() > ptPartonMax  ) accept = false;
0044
0045   if (jet.pt() * jec < ptRecJetMin || jet.pt() * jec > ptRecJetMax)
0046     return false;
0047
0048   if (jet.p() * jec < pRecJetMin || jet.p() * jec > pRecJetMax)
0049     return false;
0050
0051   if (!infos.isValid()) {
0052     LogDebug("infos not valid") << "A valid SoftLeptonTagInfoCollection was not found!"
0053                                 << " Skipping ratio check.";
0054   } else {
0055     double pToEratio = ratio(jet, infos);
0056     if (pToEratio / jec < ratioMin || pToEratio / jec > ratioMax)
0057       return false;
0058   }
0059
0060   if (doJetID) {
0061     const reco::PFJet* pfjet = dynamic_cast<const reco::PFJet*>(&jet);
0062     if (pfjet) {
0063       double neutralHadronEnergyFraction = pfjet->neutralHadronEnergy() / (jet.energy() * jec);
0064       double neutralEmEnergyFraction = pfjet->neutralEmEnergy() / (jet.energy() * jec);
0065       int nConstituents = pfjet->getPFConstituents().size();
0066       double chargedHadronEnergyFraction = pfjet->chargedHadronEnergy() / (jet.energy() * jec);
0067       int chargedMultiplicity = pfjet->chargedMultiplicity();
0068       double chargedEmEnergyFraction = pfjet->chargedEmEnergy() / (jet.energy() * jec);
0069       if (!(neutralHadronEnergyFraction < 0.99 && neutralEmEnergyFraction < 0.99 && nConstituents > 1 &&
0070             chargedHadronEnergyFraction > 0.0 && chargedMultiplicity > 0.0 && chargedEmEnergyFraction < 0.99))
0071         return false;  //2012 values
0072     }
0073     //Looks to not work -> to be fixed
0074     //else std::cout<<"Jets are not PF jets, put 'doJetID' to False."<<std::endl;
0075   }
0076
0077   return true;
0078 }
0079
0080 double AcceptJet::ratio(const reco::Jet& jet, const edm::Handle<reco::SoftLeptonTagInfoCollection>& infos) const {
0081   double jetRatio = 0.0;
0082   reco::SoftLeptonTagInfoCollection::const_iterator infoiter = infos->begin();
0083   for (; infoiter != infos->end(); ++infoiter) {
0084     if (reco::deltaR(jet.eta(), jet.phi(), infoiter->jet()->eta(), infoiter->jet()->phi()) > 1e-4)
0085       continue;
0086
0087     if (infoiter->leptons() == 0)
0088       break;
0089
0090     for (unsigned int k = 0; k != infoiter->leptons(); ++k) {
0091       double tempRatio = infoiter->properties(k).ratio;
0092       if (tempRatio > jetRatio)
0093         jetRatio = tempRatio;
0094     }
0095   }
0096
0097   return jetRatio;
0098 }