23 #include "fastjet/PseudoJet.hh"
24 #include <fastjet/JetDefinition.hh>
25 #include <TLorentzVector.h>
37 produces<edm::ValueMap<float>>(
"lsf3");
38 produces<edm::ValueMap<int>>(
"muIdx3SJ");
39 produces<edm::ValueMap<int>>(
"eleIdx3SJ");
48 static bool orderPseudoJet(fastjet::PseudoJet j1, fastjet::PseudoJet j2);
49 std::tuple<float, float>
calculateLSF(std::vector<fastjet::PseudoJet> iCParticles,
50 std::vector<fastjet::PseudoJet> &ljets,
74 unsigned int nJet =
srcJet->size();
75 unsigned int nEle =
srcEle->size();
76 unsigned int nMu =
srcMu->size();
78 std::vector<float> vlsf3;
79 std::vector<int> vmuIdx3SJ;
80 std::vector<int> veleIdx3SJ;
82 int ele_pfmatch_index = -1;
83 int mu_pfmatch_index = -1;
86 for (
unsigned int ij = 0; ij < nJet; ij++) {
87 const pat::Jet &itJet = (*srcJet)[ij];
90 std::vector<fastjet::PseudoJet> lClusterParticles;
91 float lepPt(-1), lepEta(-1), lepPhi(-1);
94 fastjet::PseudoJet
p(
d->px(),
d->py(),
d->pz(),
d->energy());
95 lClusterParticles.emplace_back(
p);
99 double dRmin(0.8), dRele(999), dRmu(999), dRtmp(999);
100 for (
unsigned int il(0); il < nEle; il++) {
101 auto itLep =
srcEle->ptrAt(il);
104 if (dRtmp < dRmin && dRtmp < dRele && itLep->
pt() > lepPt) {
106 lepEta = itLep->eta();
107 lepPhi = itLep->phi();
109 ele_pfmatch_index = il;
115 for (
unsigned int il(0); il < nMu; il++) {
116 auto itLep =
srcMu->ptrAt(il);
119 if (dRtmp < dRmin && dRtmp < dRele && dRtmp < dRmu && itLep->
pt() > lepPt) {
121 lepEta = itLep->eta();
122 lepPhi = itLep->phi();
124 ele_pfmatch_index = -1;
125 mu_pfmatch_index = il;
132 std::vector<fastjet::PseudoJet> psub_3;
133 std::sort(lClusterParticles.begin(), lClusterParticles.end(), orderPseudoJet);
134 auto lsf_3 = calculateLSF(lClusterParticles, psub_3, lepPt, lepEta, lepPhi, lepId, 2.0, 3);
135 vlsf3.push_back(std::get<0>(lsf_3));
136 veleIdx3SJ.push_back(ele_pfmatch_index);
137 vmuIdx3SJ.push_back(mu_pfmatch_index);
149 fillermuIdx3SJ.
insert(
srcJet, vmuIdx3SJ.begin(), vmuIdx3SJ.end());
150 fillermuIdx3SJ.
fill();
155 fillereleIdx3SJ.
insert(
srcJet, veleIdx3SJ.begin(), veleIdx3SJ.end());
156 fillereleIdx3SJ.
fill();
160 template <
typename T>
162 return j1.perp2() > j2.perp2();
165 template <
typename T>
167 std::vector<fastjet::PseudoJet> &lsubjets,
174 float lsf(-1), lmd(-1);
175 if (ilPt > 0 && (ilId == 11 || ilId == 13)) {
178 ilep.SetPtEtaPhiM(ilPt, ilEta, ilPhi, 0.000511);
180 ilep.SetPtEtaPhiM(ilPt, ilEta, ilPhi, 0.105658);
181 fastjet::JetDefinition lCJet_def(fastjet::kt_algorithm,
dr);
182 fastjet::ClusterSequence lCClust_seq(iCParticles, lCJet_def);
184 lsubjets = sorted_by_pt(lCClust_seq.exclusive_jets_up_to(nsj));
186 lsubjets = sorted_by_pt(lCClust_seq.inclusive_jets());
190 for (
unsigned int i0 = 0; i0 < lsubjets.size(); i0++) {
199 pVec.SetPtEtaPhiM(lsubjets[lId].
pt(), lsubjets[lId].
eta(), lsubjets[lId].
phi(), lsubjets[lId].
m());
200 lsf = ilep.Pt() / pVec.Pt();
201 lmd = (ilep - pVec).M() / pVec.M();
204 return std::tuple<float, float>(lsf, lmd);
207 template <
typename T>