|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
#include "DQMOffline/RecoB/interface/AcceptJet.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include <iostream>
using namespace std;
AcceptJet::AcceptJet(const double& etaMin_,
const double& etaMax_,
const double& ptMin_,
const double& ptMax_,
const double& pMin_,
const double& pMax_,
const double& ratioMin_,
const double& ratioMax_,
const bool& doJetID_)
: etaMin(etaMin_),
etaMax(etaMax_),
ptRecJetMin(ptMin_),
ptRecJetMax(ptMax_),
pRecJetMin(pMin_),
pRecJetMax(pMax_),
ratioMin(ratioMin_),
ratioMax(ratioMax_),
doJetID(doJetID_) {}
bool AcceptJet::operator()(const reco::Jet& jet,
const int& jetFlavour,
const edm::Handle<reco::SoftLeptonTagInfoCollection>& infos,
const double jec) const {
// temporary fudge to correct for double loop error
// jetPartonMomentum /= 2.0;
if (fabs(jet.eta()) < etaMin || fabs(jet.eta()) > etaMax)
return false;
// if ( jetFlavour.underlyingParton4Vec().P() < pPartonMin ||
// jetFlavour.underlyingParton4Vec().P() > pPartonMax ) accept = false;
//
// if ( jetFlavour.underlyingParton4Vec().Pt() < ptPartonMin ||
// jetFlavour.underlyingParton4Vec().Pt() > ptPartonMax ) accept = false;
if (jet.pt() * jec < ptRecJetMin || jet.pt() * jec > ptRecJetMax)
return false;
if (jet.p() * jec < pRecJetMin || jet.p() * jec > pRecJetMax)
return false;
if (!infos.isValid()) {
LogDebug("infos not valid") << "A valid SoftLeptonTagInfoCollection was not found!"
<< " Skipping ratio check.";
} else {
double pToEratio = ratio(jet, infos);
if (pToEratio / jec < ratioMin || pToEratio / jec > ratioMax)
return false;
}
if (doJetID) {
const reco::PFJet* pfjet = dynamic_cast<const reco::PFJet*>(&jet);
if (pfjet) {
double neutralHadronEnergyFraction = pfjet->neutralHadronEnergy() / (jet.energy() * jec);
double neutralEmEnergyFraction = pfjet->neutralEmEnergy() / (jet.energy() * jec);
int nConstituents = pfjet->getPFConstituents().size();
double chargedHadronEnergyFraction = pfjet->chargedHadronEnergy() / (jet.energy() * jec);
int chargedMultiplicity = pfjet->chargedMultiplicity();
double chargedEmEnergyFraction = pfjet->chargedEmEnergy() / (jet.energy() * jec);
if (!(neutralHadronEnergyFraction < 0.99 && neutralEmEnergyFraction < 0.99 && nConstituents > 1 &&
chargedHadronEnergyFraction > 0.0 && chargedMultiplicity > 0.0 && chargedEmEnergyFraction < 0.99))
return false; //2012 values
}
//Looks to not work -> to be fixed
//else std::cout<<"Jets are not PF jets, put 'doJetID' to False."<<std::endl;
}
return true;
}
double AcceptJet::ratio(const reco::Jet& jet, const edm::Handle<reco::SoftLeptonTagInfoCollection>& infos) const {
double jetRatio = 0.0;
reco::SoftLeptonTagInfoCollection::const_iterator infoiter = infos->begin();
for (; infoiter != infos->end(); ++infoiter) {
if (reco::deltaR(jet.eta(), jet.phi(), infoiter->jet()->eta(), infoiter->jet()->phi()) > 1e-4)
continue;
if (infoiter->leptons() == 0)
break;
for (unsigned int k = 0; k != infoiter->leptons(); ++k) {
double tempRatio = infoiter->properties(k).ratio;
if (tempRatio > jetRatio)
jetRatio = tempRatio;
}
}
return jetRatio;
}
|