CMS 3D CMS Logo

L1ExtraParticlesProd.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: L1ExtraParticlesProd
4 // Class: L1ExtraParticlesProd
5 //
9 //
10 // Original Author: Werner Sun
11 // Created: Mon Oct 2 22:45:32 EDT 2006
12 //
13 //
14 
15 // system include files
16 #include <memory>
17 
18 // user include files
20 
23 
27 
29 
31 
32 // #include "FWCore/Utilities/interface/EDMException.h"
33 
34 //
35 // class decleration
36 //
37 
38 //
39 // constants, enums and typedefs
40 //
41 
42 //
43 // static data member definitions
44 //
45 
46 double const L1ExtraParticlesProd::muonMassGeV_ = 0.105658369; // PDG06
47 
48 //
49 // constructors and destructor
50 //
52  : produceMuonParticles_(iConfig.getParameter<bool>("produceMuonParticles")),
53  muonSource_(iConfig.getParameter<edm::InputTag>("muonSource")),
54  produceCaloParticles_(iConfig.getParameter<bool>("produceCaloParticles")),
55  isoEmSource_(iConfig.getParameter<edm::InputTag>("isolatedEmSource")),
56  nonIsoEmSource_(iConfig.getParameter<edm::InputTag>("nonIsolatedEmSource")),
57  cenJetSource_(iConfig.getParameter<edm::InputTag>("centralJetSource")),
58  forJetSource_(iConfig.getParameter<edm::InputTag>("forwardJetSource")),
59  tauJetSource_(iConfig.getParameter<edm::InputTag>("tauJetSource")),
60  isoTauJetSource_(iConfig.getParameter<edm::InputTag>("isoTauJetSource")),
61  etTotSource_(iConfig.getParameter<edm::InputTag>("etTotalSource")),
62  etHadSource_(iConfig.getParameter<edm::InputTag>("etHadSource")),
63  etMissSource_(iConfig.getParameter<edm::InputTag>("etMissSource")),
64  htMissSource_(iConfig.getParameter<edm::InputTag>("htMissSource")),
65  hfRingEtSumsSource_(iConfig.getParameter<edm::InputTag>("hfRingEtSumsSource")),
66  hfRingBitCountsSource_(iConfig.getParameter<edm::InputTag>("hfRingBitCountsSource")),
67  centralBxOnly_(iConfig.getParameter<bool>("centralBxOnly")),
68  ignoreHtMiss_(iConfig.getParameter<bool>("ignoreHtMiss")) {
69  using namespace l1extra;
70 
71  // register your products
72  produces<L1EmParticleCollection>("Isolated");
73  produces<L1EmParticleCollection>("NonIsolated");
74  produces<L1JetParticleCollection>("Central");
75  produces<L1JetParticleCollection>("Forward");
76  produces<L1JetParticleCollection>("Tau");
77  produces<L1JetParticleCollection>("IsoTau");
78  produces<L1MuonParticleCollection>();
79  produces<L1EtMissParticleCollection>("MET");
80  produces<L1EtMissParticleCollection>("MHT");
81  produces<L1HFRingsCollection>();
82 
83  // now do what ever other initialization is needed
84  consumes<L1MuGMTReadoutCollection>(muonSource_);
85  consumes<L1GctEmCandCollection>(isoEmSource_);
86  consumes<L1GctEmCandCollection>(nonIsoEmSource_);
87  consumes<L1GctJetCandCollection>(cenJetSource_);
88  consumes<L1GctJetCandCollection>(forJetSource_);
89  consumes<L1GctJetCandCollection>(tauJetSource_);
90  consumes<L1GctJetCandCollection>(isoTauJetSource_);
91  consumes<L1GctEtTotalCollection>(etTotSource_);
92  consumes<L1GctEtMissCollection>(etMissSource_);
93  consumes<L1GctEtHadCollection>(etHadSource_);
94  consumes<L1GctHtMissCollection>(htMissSource_);
95  consumes<L1GctHFRingEtSumsCollection>(hfRingEtSumsSource_);
96  consumes<L1GctHFBitCountsCollection>(hfRingBitCountsSource_);
97 
99  muScalesToken_ = esConsumes<L1MuTriggerScales, L1MuTriggerScalesRcd>();
100  muPtScaleToken_ = esConsumes<L1MuTriggerPtScale, L1MuTriggerPtScaleRcd>();
101  }
102  if (produceCaloParticles_) {
103  caloGeomToken_ = esConsumes<L1CaloGeometry, L1CaloGeometryRecord>();
104  emScaleToken_ = esConsumes<L1CaloEtScale, L1EmEtScaleRcd>();
105  jetScaleToken_ = esConsumes<L1CaloEtScale, L1JetEtScaleRcd>();
106  jetFinderParamsToken_ = esConsumes<L1GctJetFinderParams, L1GctJetFinderParamsRcd>();
107  hfRingEtScaleToken_ = esConsumes<L1CaloEtScale, L1HfRingEtScaleRcd>();
108  if (!ignoreHtMiss_) {
109  htMissScaleToken_ = esConsumes<L1CaloEtScale, L1HtMissScaleRcd>();
110  }
111  }
112 }
113 
115  // do anything here that needs to be done at desctruction time
116  // (e.g. close files, deallocate resources etc.)
117 }
118 
119 //
120 // member functions
121 //
122 
123 // ------------ method called to produce the data ------------
125  using namespace edm;
126  using namespace l1extra;
127  using namespace std;
128 
129  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
130  // ~~~~~~~~~~~~~~~~~~~~ Muons ~~~~~~~~~~~~~~~~~~~~
131  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132 
133  unique_ptr<L1MuonParticleCollection> muColl(new L1MuonParticleCollection);
134 
135  if (produceMuonParticles_) {
137 
139 
140  Handle<L1MuGMTReadoutCollection> hwMuCollection;
141  iEvent.getByLabel(muonSource_, hwMuCollection);
142 
143  vector<L1MuGMTExtendedCand> hwMuCands;
144 
145  if (!hwMuCollection.isValid()) {
146  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1MuGMTReadoutCollection with " << muonSource_
147  << "\nrequested in configuration, but not found in the event." << std::endl;
148  } else {
149  if (centralBxOnly_) {
150  // Get GMT candidates from central bunch crossing only
151  hwMuCands = hwMuCollection->getRecord().getGMTCands();
152  } else {
153  // Get GMT candidates from all bunch crossings
154  vector<L1MuGMTReadoutRecord> records = hwMuCollection->getRecords();
155  vector<L1MuGMTReadoutRecord>::const_iterator rItr = records.begin();
156  vector<L1MuGMTReadoutRecord>::const_iterator rEnd = records.end();
157 
158  for (; rItr != rEnd; ++rItr) {
159  vector<L1MuGMTExtendedCand> tmpCands = rItr->getGMTCands();
160 
161  hwMuCands.insert(hwMuCands.end(), tmpCands.begin(), tmpCands.end());
162  }
163  }
164 
165  // cout << "HW muons" << endl ;
166  vector<L1MuGMTExtendedCand>::const_iterator muItr = hwMuCands.begin();
167  vector<L1MuGMTExtendedCand>::const_iterator muEnd = hwMuCands.end();
168  for (int i = 0; muItr != muEnd; ++muItr, ++i) {
169  // cout << "#" << i
170  // << " name " << muItr->name()
171  // << " empty " << muItr->empty()
172  // << " pt " << muItr->ptIndex()
173  // << " eta " << muItr->etaIndex()
174  // << " phi " << muItr->phiIndex()
175  // << " iso " << muItr->isol()
176  // << " mip " << muItr->mip()
177  // << " bx " << muItr->bx()
178  // << endl ;
179 
180  if (!muItr->empty()) {
181  // keep x and y components non-zero and protect against roundoff.
182  double pt = muPtScale->getPtScale()->getLowEdge(muItr->ptIndex()) + 1.e-6;
183 
184  // cout << "L1Extra pt " << pt << endl ;
185 
186  double eta = muScales->getGMTEtaScale()->getCenter(muItr->etaIndex());
187 
188  double phi = muScales->getPhiScale()->getLowEdge(muItr->phiIndex());
189 
191 
192  muColl->push_back(L1MuonParticle(muItr->charge(), p4, *muItr, muItr->bx()));
193  }
194  }
195  }
196  }
197 
199 
200  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201  // ~~~~~~~~~~~~~~~~~~~~ Calorimeter ~~~~~~~~~~~~~~~~~~~~
202  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
203 
204  unique_ptr<L1EmParticleCollection> isoEmColl(new L1EmParticleCollection);
205 
206  unique_ptr<L1EmParticleCollection> nonIsoEmColl(new L1EmParticleCollection);
207 
208  unique_ptr<L1JetParticleCollection> cenJetColl(new L1JetParticleCollection);
209 
210  unique_ptr<L1JetParticleCollection> forJetColl(new L1JetParticleCollection);
211 
212  unique_ptr<L1JetParticleCollection> tauJetColl(new L1JetParticleCollection);
213 
214  unique_ptr<L1JetParticleCollection> isoTauJetColl(new L1JetParticleCollection);
215 
216  unique_ptr<L1EtMissParticleCollection> etMissColl(new L1EtMissParticleCollection);
217 
218  unique_ptr<L1EtMissParticleCollection> htMissColl(new L1EtMissParticleCollection);
219 
220  unique_ptr<L1HFRingsCollection> hfRingsColl(new L1HFRingsCollection);
221 
222  if (produceCaloParticles_) {
223  // ~~~~~~~~~~~~~~~~~~~~ Geometry ~~~~~~~~~~~~~~~~~~~~
224 
225  ESHandle<L1CaloGeometry> caloGeomESH = iSetup.getHandle(caloGeomToken_);
226  const L1CaloGeometry *caloGeom = &(*caloGeomESH);
227 
228  // ~~~~~~~~~~~~~~~~~~~~ EM ~~~~~~~~~~~~~~~~~~~~
229 
231 
232  // Isolated EM
233  Handle<L1GctEmCandCollection> hwIsoEmCands;
234  iEvent.getByLabel(isoEmSource_, hwIsoEmCands);
235 
236  if (!hwIsoEmCands.isValid()) {
237  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << isoEmSource_
238  << "\nrequested in configuration, but not found in the event." << std::endl;
239  } else {
240  // cout << "HW iso EM" << endl ;
241 
242  L1GctEmCandCollection::const_iterator emItr = hwIsoEmCands->begin();
243  L1GctEmCandCollection::const_iterator emEnd = hwIsoEmCands->end();
244  for (int i = 0; emItr != emEnd; ++emItr, ++i) {
245  // cout << "#" << i
246  // << " name " << emItr->name()
247  // << " empty " << emItr->empty()
248  // << " rank " << emItr->rank()
249  // << " eta " << emItr->etaIndex()
250  // << " sign " << emItr->etaSign()
251  // << " phi " << emItr->phiIndex()
252  // << " iso " << emItr->isolated()
253  // << " bx " << emItr->bx()
254  // << endl ;
255 
256  if (!emItr->empty() && (!centralBxOnly_ || emItr->bx() == 0)) {
257  double et = emScale->et(emItr->rank());
258 
259  // cout << "L1Extra et " << et << endl ;
260 
261  isoEmColl->push_back(L1EmParticle(
262  gctLorentzVector(et, *emItr, caloGeom, true), Ref<L1GctEmCandCollection>(hwIsoEmCands, i), emItr->bx()));
263  }
264  }
265  }
266 
267  // Non-isolated EM
268  Handle<L1GctEmCandCollection> hwNonIsoEmCands;
269  iEvent.getByLabel(nonIsoEmSource_, hwNonIsoEmCands);
270 
271  if (!hwNonIsoEmCands.isValid()) {
272  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << nonIsoEmSource_
273  << "\nrequested in configuration, but not found in the event." << std::endl;
274  } else {
275  // cout << "HW non-iso EM" << endl ;
276  L1GctEmCandCollection::const_iterator emItr = hwNonIsoEmCands->begin();
277  L1GctEmCandCollection::const_iterator emEnd = hwNonIsoEmCands->end();
278  for (int i = 0; emItr != emEnd; ++emItr, ++i) {
279  // cout << "#" << i
280  // << " name " << emItr->name()
281  // << " empty " << emItr->empty()
282  // << " rank " << emItr->rank()
283  // << " eta " << emItr->etaIndex()
284  // << " sign " << emItr->etaSign()
285  // << " phi " << emItr->phiIndex()
286  // << " iso " << emItr->isolated()
287  // << " bx " << emItr->bx()
288  // << endl ;
289 
290  if (!emItr->empty() && (!centralBxOnly_ || emItr->bx() == 0)) {
291  double et = emScale->et(emItr->rank());
292 
293  // cout << "L1Extra et " << et << endl ;
294 
295  nonIsoEmColl->push_back(L1EmParticle(gctLorentzVector(et, *emItr, caloGeom, true),
296  Ref<L1GctEmCandCollection>(hwNonIsoEmCands, i),
297  emItr->bx()));
298  }
299  }
300  }
301 
302  // ~~~~~~~~~~~~~~~~~~~~ Jets ~~~~~~~~~~~~~~~~~~~~
303 
305 
306  // Central jets.
307  Handle<L1GctJetCandCollection> hwCenJetCands;
308  iEvent.getByLabel(cenJetSource_, hwCenJetCands);
309 
310  if (!hwCenJetCands.isValid()) {
311  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << cenJetSource_
312  << "\nrequested in configuration, but not found in the event." << std::endl;
313  } else {
314  // cout << "HW central jets" << endl ;
315  L1GctJetCandCollection::const_iterator jetItr = hwCenJetCands->begin();
316  L1GctJetCandCollection::const_iterator jetEnd = hwCenJetCands->end();
317  for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
318  // cout << "#" << i
319  // << " name " << jetItr->name()
320  // << " empty " << jetItr->empty()
321  // << " rank " << jetItr->rank()
322  // << " eta " << jetItr->etaIndex()
323  // << " sign " << jetItr->etaSign()
324  // << " phi " << jetItr->phiIndex()
325  // << " cen " << jetItr->isCentral()
326  // << " for " << jetItr->isForward()
327  // << " tau " << jetItr->isTau()
328  // << " bx " << jetItr->bx()
329  // << endl ;
330 
331  if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
332  double et = jetScale->et(jetItr->rank());
333 
334  // cout << "L1Extra et " << et << endl ;
335 
336  cenJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
337  Ref<L1GctJetCandCollection>(hwCenJetCands, i),
338  jetItr->bx()));
339  }
340  }
341  }
342 
343  // Forward jets.
344  Handle<L1GctJetCandCollection> hwForJetCands;
345  iEvent.getByLabel(forJetSource_, hwForJetCands);
346 
347  if (!hwForJetCands.isValid()) {
348  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << forJetSource_
349  << "\nrequested in configuration, but not found in the event." << std::endl;
350  } else {
351  // cout << "HW forward jets" << endl ;
352  L1GctJetCandCollection::const_iterator jetItr = hwForJetCands->begin();
353  L1GctJetCandCollection::const_iterator jetEnd = hwForJetCands->end();
354  for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
355  // cout << "#" << i
356  // << " name " << jetItr->name()
357  // << " empty " << jetItr->empty()
358  // << " rank " << jetItr->rank()
359  // << " eta " << jetItr->etaIndex()
360  // << " sign " << jetItr->etaSign()
361  // << " phi " << jetItr->phiIndex()
362  // << " cen " << jetItr->isCentral()
363  // << " for " << jetItr->isForward()
364  // << " tau " << jetItr->isTau()
365  // << " bx " << jetItr->bx()
366  // << endl ;
367 
368  if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
369  double et = jetScale->et(jetItr->rank());
370 
371  // cout << "L1Extra et " << et << endl ;
372 
373  forJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, false),
374  Ref<L1GctJetCandCollection>(hwForJetCands, i),
375  jetItr->bx()));
376  }
377  }
378  }
379 
380  // Tau jets.
381  // cout << "HW tau jets" << endl ;
382  Handle<L1GctJetCandCollection> hwTauJetCands;
383  iEvent.getByLabel(tauJetSource_, hwTauJetCands);
384 
385  if (!hwTauJetCands.isValid()) {
386  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << tauJetSource_
387  << "\nrequested in configuration, but not found in the event." << std::endl;
388  } else {
389  L1GctJetCandCollection::const_iterator jetItr = hwTauJetCands->begin();
390  L1GctJetCandCollection::const_iterator jetEnd = hwTauJetCands->end();
391  for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
392  // cout << "#" << i
393  // << " name " << jetItr->name()
394  // << " empty " << jetItr->empty()
395  // << " rank " << jetItr->rank()
396  // << " eta " << jetItr->etaIndex()
397  // << " sign " << jetItr->etaSign()
398  // << " phi " << jetItr->phiIndex()
399  // << " cen " << jetItr->isCentral()
400  // << " for " << jetItr->isForward()
401  // << " tau " << jetItr->isTau()
402  // << " bx " << jetItr->bx()
403  // << endl ;
404 
405  if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
406  double et = jetScale->et(jetItr->rank());
407 
408  // cout << "L1Extra et " << et << endl ;
409 
410  tauJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
411  Ref<L1GctJetCandCollection>(hwTauJetCands, i),
412  jetItr->bx()));
413  }
414  }
415  }
416 
417  // Isolated Tau jets.
418  // cout << "HW tau jets" << endl ;
419  Handle<L1GctJetCandCollection> hwIsoTauJetCands;
420  iEvent.getByLabel(isoTauJetSource_, hwIsoTauJetCands);
421 
422  if (!hwIsoTauJetCands.isValid()) {
423  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << isoTauJetSource_
424  << "\nrequested in configuration, but not found in the event." << std::endl;
425  } else {
426  L1GctJetCandCollection::const_iterator jetItr = hwIsoTauJetCands->begin();
427  L1GctJetCandCollection::const_iterator jetEnd = hwIsoTauJetCands->end();
428  for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
429  // cout << "#" << i
430  // << " name " << jetItr->name()
431  // << " empty " << jetItr->empty()
432  // << " rank " << jetItr->rank()
433  // << " eta " << jetItr->etaIndex()
434  // << " sign " << jetItr->etaSign()
435  // << " phi " << jetItr->phiIndex()
436  // << " cen " << jetItr->isCentral()
437  // << " for " << jetItr->isForward()
438  // << " tau " << jetItr->isTau()
439  // << " bx " << jetItr->bx()
440  // << endl ;
441 
442  if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
443  double et = jetScale->et(jetItr->rank());
444 
445  // cout << "L1Extra et " << et << endl ;
446 
447  isoTauJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
448  Ref<L1GctJetCandCollection>(hwIsoTauJetCands, i),
449  jetItr->bx()));
450  }
451  }
452  }
453 
454  // ~~~~~~~~~~~~~~~~~~~~ ET Sums ~~~~~~~~~~~~~~~~~~~~
455 
456  double etSumLSB = jetScale->linearLsb();
457 
458  Handle<L1GctEtTotalCollection> hwEtTotColl;
459  iEvent.getByLabel(etTotSource_, hwEtTotColl);
460 
461  Handle<L1GctEtMissCollection> hwEtMissColl;
462  iEvent.getByLabel(etMissSource_, hwEtMissColl);
463 
464  if (!hwEtTotColl.isValid()) {
465  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtTotalCollection with " << etTotSource_
466  << "\nrequested in configuration, but not found in the event." << std::endl;
467  } else if (!hwEtMissColl.isValid()) {
468  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtMissCollection with " << etMissSource_
469  << "\nrequested in configuration, but not found in the event." << std::endl;
470  } else {
471  // Make a L1EtMissParticle even if either L1GctEtTotal or L1GctEtMiss
472  // is missing for a given bx. Keep track of which L1GctEtMiss objects
473  // have a corresponding L1GctEtTotal object.
474  std::vector<bool> etMissMatched;
475 
476  L1GctEtMissCollection::const_iterator hwEtMissItr = hwEtMissColl->begin();
477  L1GctEtMissCollection::const_iterator hwEtMissEnd = hwEtMissColl->end();
478  for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr) {
479  etMissMatched.push_back(false);
480  }
481 
482  // Collate energy sums by bx
483  L1GctEtTotalCollection::const_iterator hwEtTotItr = hwEtTotColl->begin();
484  L1GctEtTotalCollection::const_iterator hwEtTotEnd = hwEtTotColl->end();
485 
486  int iTot = 0;
487  for (; hwEtTotItr != hwEtTotEnd; ++hwEtTotItr, ++iTot) {
488  int bx = hwEtTotItr->bx();
489 
490  if (!centralBxOnly_ || bx == 0) {
491  // ET bin low edge
492  double etTot =
493  (hwEtTotItr->overFlow() ? (double)L1GctEtTotal::kEtTotalMaxValue : (double)hwEtTotItr->et()) * etSumLSB +
494  1.e-6;
495 
496  int iMiss = 0;
497  hwEtMissItr = hwEtMissColl->begin();
498  hwEtMissEnd = hwEtMissColl->end();
499  for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr, ++iMiss) {
500  if (hwEtMissItr->bx() == bx) {
501  etMissMatched[iMiss] = true;
502  break;
503  }
504  }
505 
506  double etMiss = 0.;
507  double phi = 0.;
510 
511  // If a L1GctEtMiss with the right bx is not found, itr == end.
512  if (hwEtMissItr != hwEtMissEnd) {
513  // ET bin low edge
514  etMiss = (hwEtMissItr->overFlow() ? (double)L1GctEtMiss::kEtMissMaxValue : (double)hwEtMissItr->et()) *
515  etSumLSB +
516  1.e-6;
517  // keep x and y components non-zero and
518  // protect against roundoff.
519 
520  phi = caloGeom->etSumPhiBinCenter(hwEtMissItr->phi());
521 
522  p4 = math::PtEtaPhiMLorentzVector(etMiss, 0., phi, 0.);
523 
524  metRef = Ref<L1GctEtMissCollection>(hwEtMissColl, iMiss);
525 
526  // cout << "HW ET Sums " << endl
527  // << "MET: phi " << hwEtMissItr->phi() << " = " << phi
528  // << " et " << hwEtMissItr->et() << " = " << etMiss
529  // << " EtTot " << hwEtTotItr->et() << " = " << etTot
530  // << " bx " << bx
531  // << endl ;
532  }
533  // else
534  // {
535  // cout << "HW ET Sums " << endl
536  // << "MET: phi " << phi
537  // << " et "<< etMiss
538  // << " EtTot " << hwEtTotItr->et() << " = " << etTot
539  // << " bx " << bx
540  // << endl ;
541  // }
542 
543  etMissColl->push_back(L1EtMissParticle(p4,
544  L1EtMissParticle::kMET,
545  etTot,
546  metRef,
547  Ref<L1GctEtTotalCollection>(hwEtTotColl, iTot),
550  bx));
551  }
552  }
553 
554  if (!centralBxOnly_) {
555  // Make L1EtMissParticles for those L1GctEtMiss objects without
556  // a matched L1GctEtTotal object.
557 
558  double etTot = 0.;
559 
560  hwEtMissItr = hwEtMissColl->begin();
561  hwEtMissEnd = hwEtMissColl->end();
562  int iMiss = 0;
563  for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr, ++iMiss) {
564  if (!etMissMatched[iMiss]) {
565  int bx = hwEtMissItr->bx();
566 
567  // ET bin low edge
568  double etMiss =
569  (hwEtMissItr->overFlow() ? (double)L1GctEtMiss::kEtMissMaxValue : (double)hwEtMissItr->et()) *
570  etSumLSB +
571  1.e-6;
572  // keep x and y components non-zero and
573  // protect against roundoff.
574 
575  double phi = caloGeom->etSumPhiBinCenter(hwEtMissItr->phi());
576 
577  math::PtEtaPhiMLorentzVector p4(etMiss, 0., phi, 0.);
578 
579  // cout << "HW ET Sums " << endl
580  // << "MET: phi " << hwEtMissItr->phi() << " = " <<
581  // phi
582  // << " et " << hwEtMissItr->et() << " = " << etMiss
583  // << " EtTot " << etTot
584  // << " bx " << bx
585  // << endl ;
586 
587  etMissColl->push_back(L1EtMissParticle(p4,
588  L1EtMissParticle::kMET,
589  etTot,
590  Ref<L1GctEtMissCollection>(hwEtMissColl, iMiss),
594  bx));
595  }
596  }
597  }
598  }
599 
600  // ~~~~~~~~~~~~~~~~~~~~ HT Sums ~~~~~~~~~~~~~~~~~~~~
601 
602  Handle<L1GctEtHadCollection> hwEtHadColl;
603  iEvent.getByLabel(etHadSource_, hwEtHadColl);
604 
605  Handle<L1GctHtMissCollection> hwHtMissColl;
606  if (!ignoreHtMiss_) {
607  iEvent.getByLabel(htMissSource_, hwHtMissColl);
608  }
609 
611  double htSumLSB = jetFinderParams->getHtLsbGeV();
612 
613  ESHandle<L1CaloEtScale> htMissScale;
614  std::vector<bool> htMissMatched;
615  if (!ignoreHtMiss_) {
616  htMissScale = iSetup.getHandle(htMissScaleToken_);
617 
618  if (!hwEtHadColl.isValid()) {
619  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtHadCollection with " << etHadSource_
620  << "\nrequested in configuration, but not found in the event." << std::endl;
621  } else if (!hwHtMissColl.isValid()) {
622  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHtMissCollection with " << htMissSource_
623  << "\nrequested in configuration, but not found in the event." << std::endl;
624  } else {
625  // Make a L1EtMissParticle even if either L1GctEtHad or L1GctHtMiss
626  // is missing for a given bx. Keep track of which L1GctHtMiss objects
627  // have a corresponding L1GctHtTotal object.
628  L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
629  L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
630  for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr) {
631  htMissMatched.push_back(false);
632  }
633  }
634  }
635 
636  if (!hwEtHadColl.isValid()) {
637  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtHadCollection with " << etHadSource_
638  << "\nrequested in configuration, but not found in the event." << std::endl;
639  } else if (!hwHtMissColl.isValid()) {
640  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHtMissCollection with " << htMissSource_
641  << "\nrequested in configuration, but not found in the event." << std::endl;
642  } else {
643  L1GctEtHadCollection::const_iterator hwEtHadItr = hwEtHadColl->begin();
644  L1GctEtHadCollection::const_iterator hwEtHadEnd = hwEtHadColl->end();
645 
646  int iHad = 0;
647  for (; hwEtHadItr != hwEtHadEnd; ++hwEtHadItr, ++iHad) {
648  int bx = hwEtHadItr->bx();
649 
650  if (!centralBxOnly_ || bx == 0) {
651  // HT bin low edge
652  double htTot =
653  (hwEtHadItr->overFlow() ? (double)L1GctEtHad::kEtHadMaxValue : (double)hwEtHadItr->et()) * htSumLSB +
654  1.e-6;
655 
656  double htMiss = 0.;
657  double phi = 0.;
660 
661  if (!ignoreHtMiss_) {
662  L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
663  L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
664 
665  int iMiss = 0;
666  for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr, ++iMiss) {
667  if (hwHtMissItr->bx() == bx) {
668  htMissMatched[iMiss] = true;
669  break;
670  }
671  }
672 
673  // If a L1GctHtMiss with the right bx is not found, itr == end.
674  if (hwHtMissItr != hwHtMissEnd) {
675  // HT bin low edge
676  htMiss =
677  htMissScale->et(hwHtMissItr->overFlow() ? htMissScale->rankScaleMax() : hwHtMissItr->et()) + 1.e-6;
678  // keep x and y components non-zero and
679  // protect against roundoff.
680 
681  phi = caloGeom->htSumPhiBinCenter(hwHtMissItr->phi());
682 
683  p4 = math::PtEtaPhiMLorentzVector(htMiss, 0., phi, 0.);
684 
685  mhtRef = Ref<L1GctHtMissCollection>(hwHtMissColl, iMiss);
686 
687  // cout << "HW HT Sums " << endl
688  // << "MHT: phi " << hwHtMissItr->phi() <<
689  // " = " << phi
690  // << " ht " << hwHtMissItr->et() << " = "
691  // << htMiss
692  // << " HtTot " << hwEtHadItr->et() << " =
693  // "
694  // << htTot
695  // << " bx " << bx
696  // << endl ;
697  }
698  // else
699  // {
700  // cout << "HW HT Sums " << endl
701  // << "MHT: phi " << phi
702  // << " ht " << htMiss
703  // << " HtTot " << hwEtHadItr->et() << " = " <<
704  // htTot
705  // << " bx " << bx
706  // << endl ;
707  // }
708  }
709 
710  htMissColl->push_back(L1EtMissParticle(p4,
711  L1EtMissParticle::kMHT,
712  htTot,
715  mhtRef,
716  Ref<L1GctEtHadCollection>(hwEtHadColl, iHad),
717  bx));
718  }
719  }
720 
721  if (!centralBxOnly_ && !ignoreHtMiss_) {
722  // Make L1EtMissParticles for those L1GctHtMiss objects without
723  // a matched L1GctHtTotal object.
724  double htTot = 0.;
725 
726  L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
727  L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
728 
729  int iMiss = 0;
730  for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr, ++iMiss) {
731  if (!htMissMatched[iMiss]) {
732  int bx = hwHtMissItr->bx();
733 
734  // HT bin low edge
735  double htMiss =
736  htMissScale->et(hwHtMissItr->overFlow() ? htMissScale->rankScaleMax() : hwHtMissItr->et()) + 1.e-6;
737  // keep x and y components non-zero and
738  // protect against roundoff.
739 
740  double phi = caloGeom->htSumPhiBinCenter(hwHtMissItr->phi());
741 
742  math::PtEtaPhiMLorentzVector p4(htMiss, 0., phi, 0.);
743 
744  // cout << "HW HT Sums " << endl
745  // << "MHT: phi " << hwHtMissItr->phi() << " = " <<
746  // phi
747  // << " ht " << hwHtMissItr->et() << " = " << htMiss
748  // << " HtTot " << htTot
749  // << " bx " << bx
750  // << endl ;
751 
752  htMissColl->push_back(L1EtMissParticle(p4,
753  L1EtMissParticle::kMHT,
754  htTot,
757  Ref<L1GctHtMissCollection>(hwHtMissColl, iMiss),
759  bx));
760  }
761  }
762  }
763  }
764 
765  // ~~~~~~~~~~~~~~~~~~~~ HF Rings ~~~~~~~~~~~~~~~~~~~~
766 
768  iEvent.getByLabel(hfRingEtSumsSource_, hwHFEtSumsColl);
769 
770  Handle<L1GctHFBitCountsCollection> hwHFBitCountsColl;
771  iEvent.getByLabel(hfRingBitCountsSource_, hwHFBitCountsColl);
772 
773  ESHandle<L1CaloEtScale> hfRingEtScale = iSetup.getHandle(hfRingEtScaleToken_);
774 
775  if (!hwHFEtSumsColl.isValid()) {
776  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHFRingEtSumsCollection with " << hfRingEtSumsSource_
777  << "\nrequested in configuration, but not found in the event." << std::endl;
778  } else if (!hwHFBitCountsColl.isValid()) {
779  LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHFBitCountsCollection with " << hfRingBitCountsSource_
780  << "\nrequested in configuration, but not found in the event." << std::endl;
781  } else {
782  L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsItr = hwHFEtSumsColl->begin();
783  L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsEnd = hwHFEtSumsColl->end();
784 
785  int iEtSums = 0;
786  for (; hwHFEtSumsItr != hwHFEtSumsEnd; ++hwHFEtSumsItr, ++iEtSums) {
787  int bx = hwHFEtSumsItr->bx();
788 
789  if (!centralBxOnly_ || bx == 0) {
790  L1GctHFBitCountsCollection::const_iterator hwHFBitCountsItr = hwHFBitCountsColl->begin();
791  L1GctHFBitCountsCollection::const_iterator hwHFBitCountsEnd = hwHFBitCountsColl->end();
792 
793  int iBitCounts = 0;
794  for (; hwHFBitCountsItr != hwHFBitCountsEnd; ++hwHFBitCountsItr, ++iBitCounts) {
795  if (hwHFBitCountsItr->bx() == bx) {
796  break;
797  }
798  }
799 
800  // If a L1GctHFBitCounts with the right bx is not found, itr == end.
801  if (hwHFBitCountsItr != hwHFBitCountsEnd) {
802  // Construct L1HFRings only if both HW objects are present.
803 
804  // cout << "HF Rings " << endl
805 
806  // HF Et sums
807  double etSums[L1HFRings::kNumRings];
808  etSums[L1HFRings::kRing1PosEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(0)) + 1.e-6;
809  etSums[L1HFRings::kRing1NegEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(1)) + 1.e-6;
810  etSums[L1HFRings::kRing2PosEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(2)) + 1.e-6;
811  etSums[L1HFRings::kRing2NegEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(3)) + 1.e-6;
812  // protect against roundoff.
813 
814  // cout << "HF Et Sums "
815  // << hwHFEtSumsItr->etSum( 0 ) << " = "
816  // << etSums[ L1HFRings::kRing1PosEta ] << " "
817  // << hwHFEtSumsItr->etSum( 1 ) << " = "
818  // << etSums[ L1HFRings::kRing1NegEta ] << " "
819  // << hwHFEtSumsItr->etSum( 2 ) << " = "
820  // << etSums[ L1HFRings::kRing2PosEta ] << " "
821  // << hwHFEtSumsItr->etSum( 3 ) << " = "
822  // << etSums[ L1HFRings::kRing2NegEta ]
823  // << endl ;
824 
825  // HF bit counts
826  int bitCounts[L1HFRings::kNumRings];
827  bitCounts[L1HFRings::kRing1PosEta] = hwHFBitCountsItr->bitCount(0);
828  bitCounts[L1HFRings::kRing1NegEta] = hwHFBitCountsItr->bitCount(1);
829  bitCounts[L1HFRings::kRing2PosEta] = hwHFBitCountsItr->bitCount(2);
830  bitCounts[L1HFRings::kRing2NegEta] = hwHFBitCountsItr->bitCount(3);
831 
832  // cout << "HF bit counts "
833  // << hwHFBitCountsItr->bitCount( 0 ) << " "
834  // << hwHFBitCountsItr->bitCount( 1 ) << " "
835  // << hwHFBitCountsItr->bitCount( 2 ) << " "
836  // << hwHFBitCountsItr->bitCount( 3 ) << " "
837  // << endl ;
838 
839  hfRingsColl->push_back(L1HFRings(etSums,
840  bitCounts,
841  Ref<L1GctHFRingEtSumsCollection>(hwHFEtSumsColl, iEtSums),
842  Ref<L1GctHFBitCountsCollection>(hwHFBitCountsColl, iBitCounts),
843  bx));
844  }
845  }
846  }
847  }
848  }
849 
850  OrphanHandle<L1EmParticleCollection> isoEmHandle = iEvent.put(std::move(isoEmColl), "Isolated");
851 
852  OrphanHandle<L1EmParticleCollection> nonIsoEmHandle = iEvent.put(std::move(nonIsoEmColl), "NonIsolated");
853 
854  OrphanHandle<L1JetParticleCollection> cenJetHandle = iEvent.put(std::move(cenJetColl), "Central");
855 
856  OrphanHandle<L1JetParticleCollection> forJetHandle = iEvent.put(std::move(forJetColl), "Forward");
857 
858  OrphanHandle<L1JetParticleCollection> tauJetHandle = iEvent.put(std::move(tauJetColl), "Tau");
859 
860  OrphanHandle<L1JetParticleCollection> IsoTauJetHandle = iEvent.put(std::move(isoTauJetColl), "IsoTau");
861 
862  OrphanHandle<L1EtMissParticleCollection> etMissCollHandle = iEvent.put(std::move(etMissColl), "MET");
863 
864  OrphanHandle<L1EtMissParticleCollection> htMissCollHandle = iEvent.put(std::move(htMissColl), "MHT");
865 
866  OrphanHandle<L1HFRingsCollection> hfRingsCollHandle = iEvent.put(std::move(hfRingsColl));
867 }
868 
869 // math::XYZTLorentzVector
871  const L1GctCand &cand,
872  const L1CaloGeometry *geom,
873  bool central) {
874  // To keep x and y components non-zero.
875  double etCorr = et + 1.e-6; // protect against roundoff, not only for et=0
876 
877  double eta = geom->etaBinCenter(cand.etaIndex(), central);
878 
879  // double tanThOver2 = exp( -eta ) ;
880  // double ez = etCorr * ( 1. - tanThOver2 * tanThOver2 ) / ( 2. *
881  // tanThOver2 ); double e = etCorr * ( 1. + tanThOver2 * tanThOver2 ) /
882  // ( 2. * tanThOver2 );
883 
884  double phi = geom->emJetPhiBinCenter(cand.phiIndex());
885 
886  // return math::XYZTLorentzVector( etCorr * cos( phi ),
887  // etCorr * sin( phi ),
888  // ez,
889  // e ) ;
890  return math::PtEtaPhiMLorentzVector(etCorr, eta, phi, 0.);
891 }
892 
893 // define this as a plug-in
const L1MuScale * getGMTEtaScale() const
get the GMT eta scale
edm::ESGetToken< L1CaloEtScale, L1JetEtScaleRcd > jetScaleToken_
edm::ESGetToken< L1GctJetFinderParams, L1GctJetFinderParamsRcd > jetFinderParamsToken_
math::PtEtaPhiMLorentzVector gctLorentzVector(const double &et, const L1GctCand &cand, const L1CaloGeometry *geom, bool central)
edm::ESGetToken< L1CaloEtScale, L1EmEtScaleRcd > emScaleToken_
edm::ESGetToken< L1CaloEtScale, L1HfRingEtScaleRcd > hfRingEtScaleToken_
double htSumPhiBinCenter(unsigned int phiIndex) const
virtual float getLowEdge(unsigned packed) const =0
get the low edge of bin represented by packed
double etSumPhiBinCenter(unsigned int phiIndex) const
const L1MuScale * getPhiScale() const
get the phi scale
virtual float getCenter(unsigned packed) const =0
get the center of bin represented by packed
static const double muonMassGeV_
const L1MuScale * getPtScale() const
get the Pt scale
edm::ESGetToken< L1MuTriggerPtScale, L1MuTriggerPtScaleRcd > muPtScaleToken_
PtEtaPhiMLorentzVectorD PtEtaPhiMLorentzVector
Lorentz vector with cartesian internal representation.
Definition: LorentzVector.h:25
int iEvent
Definition: GenABIO.cc:224
double et(const uint16_t rank) const
convert from rank to physically meaningful quantity
edm::ESGetToken< L1CaloEtScale, L1HtMissScaleRcd > htMissScaleToken_
edm::InputTag hfRingEtSumsSource_
edm::InputTag hfRingBitCountsSource_
void produce(edm::Event &, const edm::EventSetup &) override
edm::ESGetToken< L1CaloGeometry, L1CaloGeometryRecord > caloGeomToken_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
edm::InputTag isoTauJetSource_
std::vector< L1MuGMTReadoutRecord > const & getRecords() const
edm::ESGetToken< L1MuTriggerScales, L1MuTriggerScalesRcd > muScalesToken_
bool isValid() const
Definition: HandleBase.h:70
std::vector< L1EtMissParticle > L1EtMissParticleCollection
HLT enums.
unsigned rankScaleMax() const
Definition: L1CaloEtScale.h:51
ABC for GCT EM and jet candidates.
Definition: L1GctCand.h:12
L1MuGMTReadoutRecord const & getRecord(int bx=0) const
L1ExtraParticlesProd(const edm::ParameterSet &)
std::vector< L1MuGMTExtendedCand > getGMTCands() const
get GMT candidates vector
def move(src, dest)
Definition: eostools.py:511
#define LogDebug(id)