7 #include "fastjet/ClusterSequence.hh" 16 minLeptonPt_(pset.getParameter<double>(
"minLeptonPt")),
17 maxLeptonEta_(pset.getParameter<double>(
"maxLeptonEta")),
18 minJetPt_(pset.getParameter<double>(
"minJetPt")),
19 maxJetEta_(pset.getParameter<double>(
"maxJetEta")),
20 wMass_(pset.getParameter<double>(
"wMass")),
21 tMass_(pset.getParameter<double>(
"tMass")),
22 minLeptonPtDilepton_(pset.getParameter<double>(
"minLeptonPtDilepton")),
23 maxLeptonEtaDilepton_(pset.getParameter<double>(
"maxLeptonEtaDilepton")),
24 minDileptonMassDilepton_(pset.getParameter<double>(
"minDileptonMassDilepton")),
25 minLeptonPtSemilepton_(pset.getParameter<double>(
"minLeptonPtSemilepton")),
26 maxLeptonEtaSemilepton_(pset.getParameter<double>(
"maxLeptonEtaSemilepton")),
27 minVetoLeptonPtSemilepton_(pset.getParameter<double>(
"minVetoLeptonPtSemilepton")),
28 maxVetoLeptonEtaSemilepton_(pset.getParameter<double>(
"maxVetoLeptonEtaSemilepton")),
29 minMETSemiLepton_(pset.getParameter<double>(
"minMETSemiLepton")),
30 minMtWSemiLepton_(pset.getParameter<double>(
"minMtWSemiLepton"))
39 produces<reco::GenParticleCollection>(
"neutrinos");
40 produces<reco::GenJetCollection>(
"leptons");
41 produces<reco::GenJetCollection>(
"jets");
43 produces<reco::GenParticleCollection>();
66 std::set<int> bHadronIdxs;
67 for (
size_t i=0,
n=genParticleHandle->size();
i<
n; ++
i ) {
70 if ( status == 1 )
continue;
78 std::vector<size_t> leptonIdxs;
79 for (
size_t i=0, n=finalStateHandle->size();
i<
n; ++
i ) {
82 if ( p.
status() != 1 )
continue;
91 leptonIdxs.push_back(
i);
93 case 12:
case 14:
case 16:
104 std::vector<fastjet::PseudoJet> fjLepInputs;
105 fjLepInputs.reserve(leptonIdxs.size());
106 for (
auto index : leptonIdxs ) {
110 fjLepInputs.push_back(fastjet::PseudoJet(p.
px(), p.
py(), p.
pz(), p.
energy()));
111 fjLepInputs.back().set_user_index(
index);
115 fastjet::ClusterSequence fjLepClusterSeq(fjLepInputs, *
fjLepDef_);
116 std::vector<fastjet::PseudoJet> fjLepJets = fastjet::sorted_by_pt(fjLepClusterSeq.inclusive_jets(
minLeptonPt_));
119 leptons->reserve(fjLepJets.size());
120 std::set<size_t> lepDauIdxs;
121 for (
auto& fjJet : fjLepJets ) {
125 const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
127 std::vector<reco::CandidatePtr> constituents;
129 for (
auto& fjConstituent : fjConstituents ) {
130 const size_t index = fjConstituent.user_index();
132 const int absPdgId =
abs(cand->
pdgId());
133 if ( absPdgId == 11
or absPdgId == 13 ) {
134 if ( lepCand.
isNonnull() and lepCand->
pt() > cand->
pt() )
continue;
137 constituents.push_back(cand);
139 if ( lepCand.
isNull() )
continue;
141 const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
145 lepJet.setPdgId(lepCand->
pdgId());
146 lepJet.setCharge(lepCand->
charge());
148 const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
149 lepJet.setJetArea(jetArea);
151 leptons->push_back(lepJet);
154 for (
auto& fjConstituent : fjConstituents ) {
155 lepDauIdxs.insert(fjConstituent.user_index());
162 std::vector<fastjet::PseudoJet> fjJetInputs;
163 fjJetInputs.reserve(nStables);
164 for (
size_t i=0, n=finalStateHandle->size();
i<
n; ++
i ) {
166 if ( p.
status() != 1 )
continue;
170 if ( absId == 12
or absId == 14
or absId == 16 )
continue;
171 if ( lepDauIdxs.find(
i) != lepDauIdxs.end() )
continue;
173 fjJetInputs.push_back(fastjet::PseudoJet(p.
px(), p.
py(), p.
pz(), p.
energy()));
174 fjJetInputs.back().set_user_index(
i);
177 for (
auto index : bHadronIdxs ) {
181 const double scale = 1
e-20/p.
p();
183 fjJetInputs.back().set_user_index(
index);
187 fastjet::ClusterSequence fjJetClusterSeq(fjJetInputs, *
fjJetDef_);
188 std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(fjJetClusterSeq.inclusive_jets(
minJetPt_));
191 jets->reserve(fjJets.size());
192 std::vector<size_t> bjetIdxs, ljetIdxs;
193 for (
auto& fjJet : fjJets ) {
197 const std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
199 std::vector<reco::CandidatePtr> constituents;
200 bool hasBHadron =
false;
201 for (
size_t j=0,
m=fjConstituents.size(); j<
m; ++j ) {
202 const int index = fjConstituents[j].user_index();
203 if ( bHadronIdxs.find(index) != bHadronIdxs.end() ) {
208 constituents.push_back(cand);
211 const LorentzVector jetP4(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E());
215 const double jetArea = fjJet.has_area() ? fjJet.area() : 0;
216 genJet.setJetArea(jetArea);
219 bjetIdxs.push_back(jets->size());
222 ljetIdxs.push_back(jets->size());
225 jets->push_back(genJet);
231 if ( bjetIdxs.size() < 2 )
break;
234 if ( leptons->size() == 2 and neutrinos->size() >= 2 ) {
236 const int q1 = leptons->at(0).charge();
237 const int q2 = leptons->at(1).charge();
238 if ( q1*q2 > 0 )
break;
240 const auto& lepton1 = q1 > 0 ? leptons->at(0) : leptons->at(1);
241 const auto& lepton2 = q1 > 0 ? leptons->at(1) : leptons->at(0);
248 int selNu1 = -1, selNu2 = -1;
249 for (
int i=0;
i<2; ++
i ) {
251 for (
int j=0; j<2; ++j ) {
252 if (
i == j )
continue;
254 const double newDm = dm1+dm2;
256 if ( newDm < dm ) { dm = newDm; selNu1 =
i; selNu2 = j; }
259 if ( dm >= 1e9 )
break;
261 const auto& nu1 = neutrinos->at(selNu1);
262 const auto& nu2 = neutrinos->at(selNu2);
263 const auto w1LVec = lepton1.p4()+nu1.p4();
264 const auto w2LVec = lepton2.p4()+nu2.p4();
268 int selB1 = -1, selB2 = -1;
269 for (
auto i : bjetIdxs ) {
271 for (
auto j : bjetIdxs ) {
272 if (
i == j )
continue;
274 const double newDm = dm1+dm2;
276 if ( newDm < dm ) { dm = newDm; selB1 =
i; selB2 = j; }
279 if ( dm >= 1e9 )
break;
281 const auto&
bJet1 = jets->at(selB1);
282 const auto&
bJet2 = jets->at(selB2);
283 const auto t1LVec = w1LVec +
bJet1.p4();
284 const auto t2LVec = w2LVec +
bJet2.p4();
287 const int lepQ1 = lepton1.charge();
288 const int lepQ2 = lepton2.charge();
303 pseudoTop->push_back(t1);
304 pseudoTop->push_back(t2);
306 pseudoTop->push_back(w1);
307 pseudoTop->push_back(b1);
308 pseudoTop->push_back(l1);
309 pseudoTop->push_back(n1);
311 pseudoTop->push_back(w2);
312 pseudoTop->push_back(b2);
313 pseudoTop->push_back(l2);
314 pseudoTop->push_back(n2);
316 else if ( !leptons->empty() and !neutrinos->empty() and ljetIdxs.size() >= 2 ) {
318 const auto& lepton = leptons->at(0);
322 bool hasVetoLepton =
false;
323 for (
auto vetoLeptonItr =
next(leptons->begin()); vetoLeptonItr != leptons->end(); ++vetoLeptonItr ) {
326 hasVetoLepton =
true;
329 if ( hasVetoLepton )
break;
332 double metX = 0, metY = 0;
333 for (
auto neutrino : *neutrinos ) {
334 metX += neutrino.
px();
335 metY += neutrino.py();
337 const double metPt = std::hypot(metX, metY);
340 const double mtW =
std::sqrt( 2*lepton.pt()*metPt*
cos(lepton.phi()-atan2(metX, metY)) );
349 const double lpz = lepton.pz(), le = lepton.energy(), lpt = lepton.pt();
350 const double cf = (
wMass_*
wMass_)/2 + lepton.px()*metX + lepton.py()*metY;
351 const double det = cf*cf*lpz*lpz - lpt*lpt*(le*le*metPt*metPt - cf*cf);
352 const double pz = (cf*lpz + (cf < 0 ? -
sqrt(det) :
sqrt(det)))/lpt/lpt;
354 const auto w1LVec = lepton.p4()+nu1P4;
359 for (
auto b1Itr = bjetIdxs.begin(); b1Itr != bjetIdxs.end(); ++b1Itr ) {
360 const double dR =
deltaR(jets->at(*b1Itr).p4(), w1LVec);
361 if ( dR < minDR ) { selB1 = *b1Itr; minDR =
dR; }
363 if ( selB1 == -1 )
break;
364 const auto&
bJet1 = jets->at(selB1);
365 const auto t1LVec = w1LVec +
bJet1.p4();
368 const auto& wJet1 = jets->at(ljetIdxs[0]);
369 const auto& wJet2 = jets->at(ljetIdxs[1]);
370 const auto w2LVec = wJet1.p4() + wJet2.p4();
375 for (
auto i : bjetIdxs ) {
376 if (
int(
i) == selB1 )
continue;
378 if ( newDm < dm ) { dm = newDm; selB2 =
i; }
380 if ( dm >= 1e9 )
break;
382 const auto&
bJet2 = jets->at(selB2);
383 const auto t2LVec = w2LVec +
bJet2.p4();
385 const int q = lepton.charge();
399 pseudoTop->push_back(t1);
400 pseudoTop->push_back(t2);
402 pseudoTop->push_back(w1);
403 pseudoTop->push_back(b1);
404 pseudoTop->push_back(l1);
405 pseudoTop->push_back(n1);
407 pseudoTop->push_back(w2);
408 pseudoTop->push_back(b2);
409 pseudoTop->push_back(
u2);
410 pseudoTop->push_back(d2);
414 if ( pseudoTop->size() == 10 )
441 event.put(
std::move(neutrinos),
"neutrinos");
442 event.put(
std::move(leptons),
"leptons");
463 if ( !mother )
continue;
467 if ( pdgId > 100 )
return true;
475 const unsigned int absPdgId =
abs(p->
pdgId());
476 if ( !
isBHadron(absPdgId) )
return false;
491 if ( absPdgId <= 100 )
return false;
492 if ( absPdgId >= 1000000000 )
return false;
497 const int nq3 = (absPdgId / 10) % 10;
498 const int nq2 = (absPdgId / 100) % 10;
499 const int nq1 = (absPdgId / 1000) % 10;
501 if ( nq3 == 0 )
return false;
502 if ( nq1 == 0 and nq2 == 5 )
return true;
503 if ( nq1 == 5 )
return true;
509 std::auto_ptr<reco::GenParticleCollection>&
outColl)
const 517 outColl->push_back(pOut);
const edm::EDGetTokenT< edm::View< reco::Candidate > > genParticleToken_
bool isBHadron(const reco::Candidate *p) const
PseudoTopProducer(const edm::ParameterSet &pset)
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
T getParameter(std::string const &) const
reco::Particle::Point genVertex_
virtual double pz() const =0
z coordinate of momentum vector
const double maxLeptonEtaSemilepton_
common ppss p3p6s2 common epss epspn46 common const1 w2
const double maxLeptonEta_
const double minVetoLeptonPtSemilepton_
const double maxLeptonEtaDilepton_
std::shared_ptr< JetDef > fjJetDef_
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
double px() const final
x coordinate of momentum vector
edm::Ref< GenParticleCollection > GenParticleRef
persistent reference to a GenParticle
std::vector< GenJet > GenJetCollection
collection of GenJet objects
bool isFromHadron(const reco::Candidate *p) const
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
std::shared_ptr< JetDef > fjLepDef_
const double minMtWSemiLepton_
const double minMETSemiLepton_
reco::Particle::LorentzVector LorentzVector
void clearDaughters()
clear daughter references
void addDaughter(const typename daughters::value_type &)
add a daughter via a reference
void resetDaughters(const edm::ProductID &id)
set daughters product ID
const double minLeptonPtSemilepton_
virtual int status() const =0
status word
virtual double energy() const =0
energy
virtual double py() const =0
y coordinate of momentum vector
bJet1
do the razor with one or two b jets (medium CSV)
auto const T2 &decltype(t1.eta()) t2
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
virtual int pdgId() const =0
PDG identifier.
const reco::Candidate * getLast(const reco::Candidate *p)
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
Cos< T >::type cos(const T &t)
void produce(edm::Event &event, const edm::EventSetup &eventSetup) override
bool isNull() const
Checks for null.
math::XYZPoint Point
point in the space
Abs< T >::type abs(const T &t)
Jets made from MC generator particles.
virtual double p() const =0
magnitude of momentum vector
bool isNonnull() const
Checks for non-null.
double deltaR(double eta1, double eta2, double phi1, double phi2)
ProductID id() const
Accessor for product ID.
reco::GenParticleRef buildGenParticle(const reco::Candidate *p, reco::GenParticleRefProd &refHandle, std::auto_ptr< reco::GenParticleCollection > &outColl) const
virtual double pt() const =0
transverse momentum
math::XYZTLorentzVector LorentzVector
Lorentz vector.
virtual int charge() const =0
electric charge
const double maxVetoLeptonEtaSemilepton_
const double minLeptonPt_
virtual const Point & vertex() const =0
vertex position
const double minDileptonMassDilepton_
virtual double px() const =0
x coordinate of momentum vector
void resetMothers(const edm::ProductID &id)
set mother product ID
const edm::EDGetTokenT< edm::View< reco::Candidate > > finalStateToken_
const double minLeptonPtDilepton_
virtual size_type numberOfDaughters() const =0
number of daughters
void clearMothers()
clear mother references
void writeSpecific(reco::CaloJet &jet, reco::Particle::LorentzVector const &p4, reco::Particle::Point const &point, std::vector< reco::CandidatePtr > const &constituents, edm::EventSetup const &c)