CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
BJetEnergyRegressionVarProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: PhysicsTools/NanoAOD
4 // Class: BJetEnergyRegressionVarProducer
5 //
13 //
14 // Original Author: Marco Peruzzi
15 // Created: Tue, 05 Sep 2017 12:24:38 GMT
16 //
17 //
18 
19 // system include files
20 #include <memory>
21 
22 // user include files
25 
28 
31 
33 
37 
40 
44 
45 #include <vector>
46 
47 // user include files
50 
53 
56 
61 
65 
67 
68 //
69 // class declaration
70 //
71 
72 template <typename T>
74 public:
76  : srcJet_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("src"))),
77  srcVtx_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("pvsrc"))),
78  srcSV_(consumes<edm::View<reco::VertexCompositePtrCandidate>>(iConfig.getParameter<edm::InputTag>("svsrc"))),
79  srcGP_(consumes<std::vector<reco::GenParticle>>(iConfig.getParameter<edm::InputTag>("gpsrc"))) {
80  //un prodotto da copiare
81  produces<edm::ValueMap<float>>("leptonPtRel");
82  produces<edm::ValueMap<float>>("leptonPtRatio");
83  produces<edm::ValueMap<float>>("leptonPtRelInv"); //wrong variable?
84  produces<edm::ValueMap<float>>("leptonPtRelv0");
85  produces<edm::ValueMap<float>>("leptonPtRatiov0");
86  produces<edm::ValueMap<float>>("leptonPtRelInvv0"); //v0 ~ heppy?
87  produces<edm::ValueMap<float>>("leptonPt");
88  produces<edm::ValueMap<int>>("leptonPdgId");
89  produces<edm::ValueMap<float>>("leptonDeltaR");
90  produces<edm::ValueMap<float>>("leadTrackPt");
91  produces<edm::ValueMap<float>>("vtxPt");
92  produces<edm::ValueMap<float>>("vtxMass");
93  produces<edm::ValueMap<float>>("vtx3dL");
94  produces<edm::ValueMap<float>>("vtx3deL");
95  produces<edm::ValueMap<int>>("vtxNtrk");
96  produces<edm::ValueMap<float>>("ptD");
97  produces<edm::ValueMap<float>>("genPtwNu");
98  }
100 
101  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
102 
103 private:
104  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
105 
106  std::tuple<float, float, float> calculatePtRatioRel(edm::Ptr<reco::Candidate> lep, edm::Ptr<pat::Jet> jet) const;
107  std::tuple<float, float, float> calculatePtRatioRelSimple(edm::Ptr<reco::Candidate> lep,
108  edm::Ptr<pat::Jet> jet) const; //old version?
109 
110  // ----------member data ---------------------------
111 
116 };
117 
118 //
119 // constants, enums and typedefs
120 //
121 
122 //
123 // static data member definitions
124 //
125 
126 //
127 // member functions
128 //
129 
130 // ------------ method called to produce the data ------------
131 template <typename T>
134  const edm::EventSetup& iSetup) const {
136  iEvent.getByToken(srcJet_, srcJet);
138  iEvent.getByToken(srcVtx_, srcVtx);
140  iEvent.getByToken(srcSV_, srcSV);
142  iEvent.getByToken(srcGP_, srcGP);
143 
144  unsigned int nJet = srcJet->size();
145  // unsigned int nLep = srcLep->size();
146 
147  std::vector<float> leptonPtRel(nJet, 0);
148  std::vector<float> leptonPtRatio(nJet, 0);
149  std::vector<float> leptonPtRelInv(nJet, 0);
150  std::vector<float> leptonPtRel_v0(nJet, 0);
151  std::vector<float> leptonPtRatio_v0(nJet, 0);
152  std::vector<float> leptonPtRelInv_v0(nJet, 0);
153  std::vector<int> leptonPdgId(nJet, 0);
154  std::vector<float> leptonPt(nJet, 0);
155  std::vector<float> leptonDeltaR(nJet, 0);
156  std::vector<float> leadTrackPt(nJet, 0);
157  std::vector<float> vtxPt(nJet, 0);
158  std::vector<float> vtxMass(nJet, 0);
159  std::vector<float> vtx3dL(nJet, 0);
160  std::vector<float> vtx3deL(nJet, 0);
161  std::vector<int> vtxNtrk(nJet, 0);
162  std::vector<float> ptD(nJet, 0);
163  std::vector<float> genPtwNu(nJet, 0);
164 
165  const auto& pv = (*srcVtx)[0];
166  for (unsigned int ij = 0; ij < nJet; ij++) {
167  auto jet = srcJet->ptrAt(ij);
168 
169  if (jet->genJet() != nullptr) {
170  auto genp4 = jet->genJet()->p4();
171  auto gep4wNu = genp4;
172  for (const auto& gp : *srcGP) {
173  if ((abs(gp.pdgId()) == 12 || abs(gp.pdgId()) == 14 || abs(gp.pdgId()) == 16) && gp.status() == 1) {
174  if (reco::deltaR(genp4, gp.p4()) < 0.4) {
175  // std::cout<<" from "<<gep4wNu.pt()<<std::endl;
176  gep4wNu = gep4wNu + gp.p4();
177  // std::cout<<" to "<<gep4wNu.pt()<<std::endl;
178  }
179  }
180  }
181 
182  genPtwNu[ij] = gep4wNu.pt();
183  }
184 
185  float ptMax = 0;
186  float sumWeight = 0;
187  float sumPt = 0;
188 
189  for (const auto& d : jet->daughterPtrVector()) {
190  sumWeight += (d->pt()) * (d->pt());
191  sumPt += d->pt();
192  if (d->pt() > ptMax)
193  ptMax = d->pt();
194  }
195  leadTrackPt[ij] = ptMax;
196  ptD[ij] = (sumWeight > 0 ? sqrt(sumWeight) / sumPt : 0);
197 
198  //lepton properties
199  float maxLepPt = 0;
200  leptonPtRel[ij] = 0;
201 
202  for (const auto& d : jet->daughterPtrVector()) {
203  if (abs(d->pdgId()) == 11 || abs(d->pdgId()) == 13) {
204  if (d->pt() < maxLepPt)
205  continue;
206  auto res = calculatePtRatioRel(d, jet);
207  leptonPtRatio[ij] = std::get<0>(res);
208  leptonPtRel[ij] = std::get<1>(res);
209  leptonPtRelInv[ij] = std::get<2>(res);
210  auto res2 = calculatePtRatioRelSimple(d, jet);
211  leptonPtRatio_v0[ij] = std::get<0>(res2);
212  leptonPtRel_v0[ij] = std::get<1>(res2);
213  leptonPtRelInv_v0[ij] = std::get<2>(res2);
214  leptonPdgId[ij] = d->pdgId();
215  leptonDeltaR[ij] = reco::deltaR(jet->p4(), d->p4());
216  leptonPt[ij] = d->pt();
217  maxLepPt = d->pt();
218  }
219  }
220 
221  //Fill vertex properties
222  VertexDistance3D vdist;
223  float maxFoundSignificance = 0;
224 
225  vtxPt[ij] = 0;
226  vtxMass[ij] = 0;
227  vtx3dL[ij] = 0;
228  vtx3deL[ij] = 0;
229  vtxNtrk[ij] = 0;
230 
231  for (const auto& sv : *srcSV) {
232  GlobalVector flightDir(sv.vertex().x() - pv.x(), sv.vertex().y() - pv.y(), sv.vertex().z() - pv.z());
233  GlobalVector jetDir(jet->px(), jet->py(), jet->pz());
234  if (reco::deltaR2(flightDir, jetDir) < 0.09) {
235  Measurement1D dl = vdist.distance(
236  pv, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
237  if (dl.significance() > maxFoundSignificance) {
238  maxFoundSignificance = dl.significance();
239  vtxPt[ij] = sv.pt();
240  vtxMass[ij] = sv.p4().M();
241  vtx3dL[ij] = dl.value();
242  vtx3deL[ij] = dl.error();
243  vtxNtrk[ij] = sv.numberOfSourceCandidatePtrs();
244  }
245  }
246  }
247  }
248 
249  std::unique_ptr<edm::ValueMap<float>> leptonPtRelV(new edm::ValueMap<float>());
250  edm::ValueMap<float>::Filler fillerRel(*leptonPtRelV);
251  fillerRel.insert(srcJet, leptonPtRel.begin(), leptonPtRel.end());
252  fillerRel.fill();
253  iEvent.put(std::move(leptonPtRelV), "leptonPtRel");
254 
255  std::unique_ptr<edm::ValueMap<float>> leptonPtRatioV(new edm::ValueMap<float>());
256  edm::ValueMap<float>::Filler fillerRatio(*leptonPtRatioV);
257  fillerRatio.insert(srcJet, leptonPtRatio.begin(), leptonPtRatio.end());
258  fillerRatio.fill();
259  iEvent.put(std::move(leptonPtRatioV), "leptonPtRatio");
260 
261  std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV(new edm::ValueMap<float>());
262  edm::ValueMap<float>::Filler fillerRelInv(*leptonPtRelInvV);
263  fillerRelInv.insert(srcJet, leptonPtRelInv.begin(), leptonPtRelInv.end());
264  fillerRelInv.fill();
265  iEvent.put(std::move(leptonPtRelInvV), "leptonPtRelInv");
266 
267  std::unique_ptr<edm::ValueMap<float>> leptonPtRelV_v0(new edm::ValueMap<float>());
268  edm::ValueMap<float>::Filler fillerRel_v0(*leptonPtRelV_v0);
269  fillerRel_v0.insert(srcJet, leptonPtRel_v0.begin(), leptonPtRel_v0.end());
270  fillerRel_v0.fill();
271  iEvent.put(std::move(leptonPtRelV_v0), "leptonPtRelv0");
272 
273  std::unique_ptr<edm::ValueMap<float>> leptonPtRatioV_v0(new edm::ValueMap<float>());
274  edm::ValueMap<float>::Filler fillerRatio_v0(*leptonPtRatioV_v0);
275  fillerRatio_v0.insert(srcJet, leptonPtRatio_v0.begin(), leptonPtRatio_v0.end());
276  fillerRatio_v0.fill();
277  iEvent.put(std::move(leptonPtRatioV_v0), "leptonPtRatiov0");
278 
279  std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV_v0(new edm::ValueMap<float>());
280  edm::ValueMap<float>::Filler fillerRelInv_v0(*leptonPtRelInvV_v0);
281  fillerRelInv_v0.insert(srcJet, leptonPtRelInv_v0.begin(), leptonPtRelInv_v0.end());
282  fillerRelInv_v0.fill();
283  iEvent.put(std::move(leptonPtRelInvV_v0), "leptonPtRelInvv0");
284 
285  std::unique_ptr<edm::ValueMap<float>> leptonPtV(new edm::ValueMap<float>());
286  edm::ValueMap<float>::Filler fillerLpt(*leptonPtV);
287  fillerLpt.insert(srcJet, leptonPt.begin(), leptonPt.end());
288  fillerLpt.fill();
289  iEvent.put(std::move(leptonPtV), "leptonPt");
290 
291  std::unique_ptr<edm::ValueMap<float>> leptonDeltaRV(new edm::ValueMap<float>());
292  edm::ValueMap<float>::Filler fillerLdR(*leptonDeltaRV);
293  fillerLdR.insert(srcJet, leptonDeltaR.begin(), leptonDeltaR.end());
294  fillerLdR.fill();
295  iEvent.put(std::move(leptonDeltaRV), "leptonDeltaR");
296 
297  std::unique_ptr<edm::ValueMap<int>> leptonPtPdgIdV(new edm::ValueMap<int>());
298  edm::ValueMap<int>::Filler fillerId(*leptonPtPdgIdV);
299  fillerId.insert(srcJet, leptonPdgId.begin(), leptonPdgId.end());
300  fillerId.fill();
301  iEvent.put(std::move(leptonPtPdgIdV), "leptonPdgId");
302 
303  std::unique_ptr<edm::ValueMap<float>> leadTrackPtV(new edm::ValueMap<float>());
304  edm::ValueMap<float>::Filler fillerLT(*leadTrackPtV);
305  fillerLT.insert(srcJet, leadTrackPt.begin(), leadTrackPt.end());
306  fillerLT.fill();
307  iEvent.put(std::move(leadTrackPtV), "leadTrackPt");
308 
309  std::unique_ptr<edm::ValueMap<float>> vtxPtV(new edm::ValueMap<float>());
310  edm::ValueMap<float>::Filler fillerVtxPt(*vtxPtV);
311  fillerVtxPt.insert(srcJet, vtxPt.begin(), vtxPt.end());
312  fillerVtxPt.fill();
313  iEvent.put(std::move(vtxPtV), "vtxPt");
314 
315  std::unique_ptr<edm::ValueMap<float>> vtxMassV(new edm::ValueMap<float>());
316  edm::ValueMap<float>::Filler fillerVtxMass(*vtxMassV);
317  fillerVtxMass.insert(srcJet, vtxMass.begin(), vtxMass.end());
318  fillerVtxMass.fill();
319  iEvent.put(std::move(vtxMassV), "vtxMass");
320 
321  std::unique_ptr<edm::ValueMap<float>> vtx3dLV(new edm::ValueMap<float>());
322  edm::ValueMap<float>::Filler fillerVtx3dL(*vtx3dLV);
323  fillerVtx3dL.insert(srcJet, vtx3dL.begin(), vtx3dL.end());
324  fillerVtx3dL.fill();
325  iEvent.put(std::move(vtx3dLV), "vtx3dL");
326 
327  std::unique_ptr<edm::ValueMap<float>> vtx3deLV(new edm::ValueMap<float>());
328  edm::ValueMap<float>::Filler fillerVtx3deL(*vtx3deLV);
329  fillerVtx3deL.insert(srcJet, vtx3deL.begin(), vtx3deL.end());
330  fillerVtx3deL.fill();
331  iEvent.put(std::move(vtx3deLV), "vtx3deL");
332 
333  std::unique_ptr<edm::ValueMap<int>> vtxNtrkV(new edm::ValueMap<int>());
334  edm::ValueMap<int>::Filler fillerVtxNT(*vtxNtrkV);
335  fillerVtxNT.insert(srcJet, vtxNtrk.begin(), vtxNtrk.end());
336  fillerVtxNT.fill();
337  iEvent.put(std::move(vtxNtrkV), "vtxNtrk");
338 
339  std::unique_ptr<edm::ValueMap<float>> ptDV(new edm::ValueMap<float>());
340  edm::ValueMap<float>::Filler fillerPtD(*ptDV);
341  fillerPtD.insert(srcJet, ptD.begin(), ptD.end());
342  fillerPtD.fill();
343  iEvent.put(std::move(ptDV), "ptD");
344 
345  std::unique_ptr<edm::ValueMap<float>> genptV(new edm::ValueMap<float>());
346  edm::ValueMap<float>::Filler fillergenpt(*genptV);
347  fillergenpt.insert(srcJet, genPtwNu.begin(), genPtwNu.end());
348  fillergenpt.fill();
349  iEvent.put(std::move(genptV), "genPtwNu");
350 }
351 
352 template <typename T>
354  edm::Ptr<pat::Jet> jet) const {
355  auto rawp4 = jet->correctedP4("Uncorrected");
356  auto lepp4 = lep->p4();
357 
358  if ((rawp4 - lepp4).R() < 1e-4)
359  return std::tuple<float, float, float>(1.0, 0.0, 0.0);
360 
361  auto jetp4 = (rawp4 - lepp4 * (1.0 / jet->jecFactor("L1FastJet"))) * (jet->pt() / rawp4.pt()) + lepp4;
362  auto ptratio = lepp4.pt() / jetp4.pt();
363  auto ptrel = lepp4.Vect().Cross((jetp4 - lepp4).Vect().Unit()).R();
364  auto ptrelinv = (jetp4 - lepp4).Vect().Cross((lepp4).Vect().Unit()).R();
365 
366  return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
367 }
368 
369 template <typename T>
372  auto lepp4 = lep->p4();
373  auto rawp4 = jet->correctedP4("Uncorrected");
374 
375  if ((rawp4 - lepp4).R() < 1e-4)
376  return std::tuple<float, float, float>(1.0, 0.0, 0.0);
377 
378  auto jetp4 = jet->p4(); //(rawp4 - lepp4*(1.0/jet->jecFactor("L1FastJet")))*(jet->pt()/rawp4.pt())+lepp4;
379  auto ptratio = lepp4.pt() / jetp4.pt();
380  auto ptrel = lepp4.Vect().Cross((jetp4).Vect().Unit()).R();
381  auto ptrelinv = jetp4.Vect().Cross((lepp4).Vect().Unit()).R();
382 
383  return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
384 }
385 
386 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
387 template <typename T>
389  //The following says we do not know what parameters are allowed so do no validation
390  // Please change this to state exactly what you do use, even if it is no parameters
392  desc.add<edm::InputTag>("src")->setComment("jet input collection");
393  // desc.add<edm::InputTag>("musrc")->setComment("muons input collection");
394  // desc.add<edm::InputTag>("elesrc")->setComment("electrons input collection");
395  desc.add<edm::InputTag>("pvsrc")->setComment("primary vertex input collection");
396  desc.add<edm::InputTag>("svsrc")->setComment("secondary vertex input collection");
397  desc.add<edm::InputTag>("gpsrc")->setComment("genparticles for nu recovery");
399  if (typeid(T) == typeid(pat::Jet))
400  modname += "Jet";
401  modname += "RegressionVarProducer";
402  descriptions.add(modname, desc);
403 }
404 
406 
407 //define this as a plug-in
reco::Vertex::Point convertPos(const GlobalPoint &p)
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:133
edm::EDGetTokenT< edm::View< reco::VertexCompositePtrCandidate > > srcSV_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
edm::EDGetTokenT< std::vector< reco::GenParticle > > srcGP_
std::tuple< float, float, float > calculatePtRatioRel(edm::Ptr< reco::Candidate > lep, edm::Ptr< pat::Jet > jet) const
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
reco::Vertex::Error convertError(const GlobalError &ge)
Definition: ConvertError.h:8
void insert(const H &h, I begin, I end)
Definition: ValueMap.h:53
double error() const
Definition: Measurement1D.h:27
BJetEnergyRegressionVarProducer(const edm::ParameterSet &iConfig)
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
tuple d
Definition: ztail.py:151
edm::EDGetTokenT< edm::View< pat::Jet > > srcJet_
int iEvent
Definition: GenABIO.cc:224
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
T sqrt(T t)
Definition: SSEVec.h:19
BJetEnergyRegressionVarProducer< pat::Jet > JetRegressionVarProducer
def move
Definition: eostools.py:511
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ParameterDescriptionBase * add(U const &iLabel, T const &value)
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
Measurement1D distance(const GlobalPoint &vtx1Position, const GlobalError &vtx1PositionError, const GlobalPoint &vtx2Position, const GlobalError &vtx2PositionError) const override
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
double significance() const
Definition: Measurement1D.h:29
edm::EDGetTokenT< std::vector< reco::Vertex > > srcVtx_
Analysis-level calorimeter jet class.
Definition: Jet.h:77
void add(std::string const &label, ParameterSetDescription const &psetDescription)
double value() const
Definition: Measurement1D.h:25
std::tuple< float, float, float > calculatePtRatioRelSimple(edm::Ptr< reco::Candidate > lep, edm::Ptr< pat::Jet > jet) const
constexpr char Jet[]
Definition: modules.cc:9
void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override
long double T