Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-05 03:16:36 

0001 // -*- C++ -*-
0002 //
0003 // Package:    L1ExtraParticlesProd
0004 // Class:      L1ExtraParticlesProd
0005 //
0006 /**\class L1ExtraParticlesProd \file L1ExtraParticlesProd.cc
0007  * src/L1ExtraParticlesProd/src/L1ExtraParticlesProd.cc
0008  */
0009 //
0010 // Original Author:  Werner Sun
0011 //         Created:  Mon Oct  2 22:45:32 EDT 2006
0012 //
0013 //
0014 
0015 // system include files
0016 #include <memory>
0017 
0018 // user include files
0019 #include "L1Trigger/L1ExtraFromDigis/interface/L1ExtraParticlesProd.h"
0020 
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 
0024 #include "DataFormats/Common/interface/Handle.h"
0025 #include "DataFormats/Common/interface/OrphanHandle.h"
0026 #include "FWCore/Framework/interface/ESHandle.h"
0027 
0028 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"
0029 
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 
0032 // #include "FWCore/Utilities/interface/EDMException.h"
0033 
0034 //
0035 // class decleration
0036 //
0037 
0038 //
0039 // constants, enums and typedefs
0040 //
0041 
0042 //
0043 // static data member definitions
0044 //
0045 
0046 double const L1ExtraParticlesProd::muonMassGeV_ = 0.105658369;  // PDG06
0047 
0048 //
0049 // constructors and destructor
0050 //
0051 L1ExtraParticlesProd::L1ExtraParticlesProd(const edm::ParameterSet &iConfig)
0052     : produceMuonParticles_(iConfig.getParameter<bool>("produceMuonParticles")),
0053       muonSource_(iConfig.getParameter<edm::InputTag>("muonSource")),
0054       produceCaloParticles_(iConfig.getParameter<bool>("produceCaloParticles")),
0055       isoEmSource_(iConfig.getParameter<edm::InputTag>("isolatedEmSource")),
0056       nonIsoEmSource_(iConfig.getParameter<edm::InputTag>("nonIsolatedEmSource")),
0057       cenJetSource_(iConfig.getParameter<edm::InputTag>("centralJetSource")),
0058       forJetSource_(iConfig.getParameter<edm::InputTag>("forwardJetSource")),
0059       tauJetSource_(iConfig.getParameter<edm::InputTag>("tauJetSource")),
0060       isoTauJetSource_(iConfig.getParameter<edm::InputTag>("isoTauJetSource")),
0061       etTotSource_(iConfig.getParameter<edm::InputTag>("etTotalSource")),
0062       etHadSource_(iConfig.getParameter<edm::InputTag>("etHadSource")),
0063       etMissSource_(iConfig.getParameter<edm::InputTag>("etMissSource")),
0064       htMissSource_(iConfig.getParameter<edm::InputTag>("htMissSource")),
0065       hfRingEtSumsSource_(iConfig.getParameter<edm::InputTag>("hfRingEtSumsSource")),
0066       hfRingBitCountsSource_(iConfig.getParameter<edm::InputTag>("hfRingBitCountsSource")),
0067       centralBxOnly_(iConfig.getParameter<bool>("centralBxOnly")),
0068       ignoreHtMiss_(iConfig.getParameter<bool>("ignoreHtMiss")) {
0069   using namespace l1extra;
0070 
0071   // register your products
0072   produces<L1EmParticleCollection>("Isolated");
0073   produces<L1EmParticleCollection>("NonIsolated");
0074   produces<L1JetParticleCollection>("Central");
0075   produces<L1JetParticleCollection>("Forward");
0076   produces<L1JetParticleCollection>("Tau");
0077   produces<L1JetParticleCollection>("IsoTau");
0078   produces<L1MuonParticleCollection>();
0079   produces<L1EtMissParticleCollection>("MET");
0080   produces<L1EtMissParticleCollection>("MHT");
0081   produces<L1HFRingsCollection>();
0082 
0083   // now do what ever other initialization is needed
0084   consumes<L1MuGMTReadoutCollection>(muonSource_);
0085   consumes<L1GctEmCandCollection>(isoEmSource_);
0086   consumes<L1GctEmCandCollection>(nonIsoEmSource_);
0087   consumes<L1GctJetCandCollection>(cenJetSource_);
0088   consumes<L1GctJetCandCollection>(forJetSource_);
0089   consumes<L1GctJetCandCollection>(tauJetSource_);
0090   consumes<L1GctJetCandCollection>(isoTauJetSource_);
0091   consumes<L1GctEtTotalCollection>(etTotSource_);
0092   consumes<L1GctEtMissCollection>(etMissSource_);
0093   consumes<L1GctEtHadCollection>(etHadSource_);
0094   consumes<L1GctHtMissCollection>(htMissSource_);
0095   consumes<L1GctHFRingEtSumsCollection>(hfRingEtSumsSource_);
0096   consumes<L1GctHFBitCountsCollection>(hfRingBitCountsSource_);
0097 
0098   if (produceMuonParticles_) {
0099     muScalesToken_ = esConsumes<L1MuTriggerScales, L1MuTriggerScalesRcd>();
0100     muPtScaleToken_ = esConsumes<L1MuTriggerPtScale, L1MuTriggerPtScaleRcd>();
0101   }
0102   if (produceCaloParticles_) {
0103     caloGeomToken_ = esConsumes<L1CaloGeometry, L1CaloGeometryRecord>();
0104     emScaleToken_ = esConsumes<L1CaloEtScale, L1EmEtScaleRcd>();
0105     jetScaleToken_ = esConsumes<L1CaloEtScale, L1JetEtScaleRcd>();
0106     jetFinderParamsToken_ = esConsumes<L1GctJetFinderParams, L1GctJetFinderParamsRcd>();
0107     hfRingEtScaleToken_ = esConsumes<L1CaloEtScale, L1HfRingEtScaleRcd>();
0108     if (!ignoreHtMiss_) {
0109       htMissScaleToken_ = esConsumes<L1CaloEtScale, L1HtMissScaleRcd>();
0110     }
0111   }
0112 }
0113 
0114 L1ExtraParticlesProd::~L1ExtraParticlesProd() {
0115   // do anything here that needs to be done at desctruction time
0116   // (e.g. close files, deallocate resources etc.)
0117 }
0118 
0119 //
0120 // member functions
0121 //
0122 
0123 // ------------ method called to produce the data  ------------
0124 void L1ExtraParticlesProd::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0125   using namespace edm;
0126   using namespace l1extra;
0127   using namespace std;
0128 
0129   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0130   // ~~~~~~~~~~~~~~~~~~~~ Muons ~~~~~~~~~~~~~~~~~~~~
0131   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0132 
0133   unique_ptr<L1MuonParticleCollection> muColl(new L1MuonParticleCollection);
0134 
0135   if (produceMuonParticles_) {
0136     ESHandle<L1MuTriggerScales> muScales = iSetup.getHandle(muScalesToken_);
0137 
0138     ESHandle<L1MuTriggerPtScale> muPtScale = iSetup.getHandle(muPtScaleToken_);
0139 
0140     Handle<L1MuGMTReadoutCollection> hwMuCollection;
0141     iEvent.getByLabel(muonSource_, hwMuCollection);
0142 
0143     vector<L1MuGMTExtendedCand> hwMuCands;
0144 
0145     if (!hwMuCollection.isValid()) {
0146       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1MuGMTReadoutCollection with " << muonSource_
0147                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0148     } else {
0149       if (centralBxOnly_) {
0150         // Get GMT candidates from central bunch crossing only
0151         hwMuCands = hwMuCollection->getRecord().getGMTCands();
0152       } else {
0153         // Get GMT candidates from all bunch crossings
0154         vector<L1MuGMTReadoutRecord> records = hwMuCollection->getRecords();
0155         vector<L1MuGMTReadoutRecord>::const_iterator rItr = records.begin();
0156         vector<L1MuGMTReadoutRecord>::const_iterator rEnd = records.end();
0157 
0158         for (; rItr != rEnd; ++rItr) {
0159           vector<L1MuGMTExtendedCand> tmpCands = rItr->getGMTCands();
0160 
0161           hwMuCands.insert(hwMuCands.end(), tmpCands.begin(), tmpCands.end());
0162         }
0163       }
0164 
0165       //       cout << "HW muons" << endl ;
0166       vector<L1MuGMTExtendedCand>::const_iterator muItr = hwMuCands.begin();
0167       vector<L1MuGMTExtendedCand>::const_iterator muEnd = hwMuCands.end();
0168       for (int i = 0; muItr != muEnd; ++muItr, ++i) {
0169         //   cout << "#" << i
0170         //        << " name " << muItr->name()
0171         //        << " empty " << muItr->empty()
0172         //        << " pt " << muItr->ptIndex()
0173         //        << " eta " << muItr->etaIndex()
0174         //        << " phi " << muItr->phiIndex()
0175         //        << " iso " << muItr->isol()
0176         //        << " mip " << muItr->mip()
0177         //        << " bx " << muItr->bx()
0178         //        << endl ;
0179 
0180         if (!muItr->empty()) {
0181           // keep x and y components non-zero and protect against roundoff.
0182           double pt = muPtScale->getPtScale()->getLowEdge(muItr->ptIndex()) + 1.e-6;
0183 
0184           //        cout << "L1Extra pt " << pt << endl ;
0185 
0186           double eta = muScales->getGMTEtaScale()->getCenter(muItr->etaIndex());
0187 
0188           double phi = muScales->getPhiScale()->getLowEdge(muItr->phiIndex());
0189 
0190           math::PtEtaPhiMLorentzVector p4(pt, eta, phi, muonMassGeV_);
0191 
0192           muColl->push_back(L1MuonParticle(muItr->charge(), p4, *muItr, muItr->bx()));
0193         }
0194       }
0195     }
0196   }
0197 
0198   iEvent.put(std::move(muColl));
0199 
0200   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0201   // ~~~~~~~~~~~~~~~~~~~~ Calorimeter ~~~~~~~~~~~~~~~~~~~~
0202   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0203 
0204   unique_ptr<L1EmParticleCollection> isoEmColl(new L1EmParticleCollection);
0205   unique_ptr<L1EmParticleCollection> nonIsoEmColl(new L1EmParticleCollection);
0206   unique_ptr<L1JetParticleCollection> cenJetColl(new L1JetParticleCollection);
0207   unique_ptr<L1JetParticleCollection> forJetColl(new L1JetParticleCollection);
0208   unique_ptr<L1JetParticleCollection> tauJetColl(new L1JetParticleCollection);
0209   unique_ptr<L1JetParticleCollection> isoTauJetColl(new L1JetParticleCollection);
0210   unique_ptr<L1EtMissParticleCollection> etMissColl(new L1EtMissParticleCollection);
0211   unique_ptr<L1EtMissParticleCollection> htMissColl(new L1EtMissParticleCollection);
0212   unique_ptr<L1HFRingsCollection> hfRingsColl(new L1HFRingsCollection);
0213 
0214   if (produceCaloParticles_) {
0215     // ~~~~~~~~~~~~~~~~~~~~ Geometry ~~~~~~~~~~~~~~~~~~~~
0216 
0217     ESHandle<L1CaloGeometry> caloGeomESH = iSetup.getHandle(caloGeomToken_);
0218     const L1CaloGeometry *caloGeom = &(*caloGeomESH);
0219 
0220     // ~~~~~~~~~~~~~~~~~~~~ EM ~~~~~~~~~~~~~~~~~~~~
0221 
0222     ESHandle<L1CaloEtScale> emScale = iSetup.getHandle(emScaleToken_);
0223 
0224     // Isolated EM
0225     Handle<L1GctEmCandCollection> hwIsoEmCands;
0226     iEvent.getByLabel(isoEmSource_, hwIsoEmCands);
0227 
0228     if (!hwIsoEmCands.isValid()) {
0229       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << isoEmSource_
0230                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0231     } else {
0232       //       cout << "HW iso EM" << endl ;
0233 
0234       L1GctEmCandCollection::const_iterator emItr = hwIsoEmCands->begin();
0235       L1GctEmCandCollection::const_iterator emEnd = hwIsoEmCands->end();
0236       for (int i = 0; emItr != emEnd; ++emItr, ++i) {
0237         //   cout << "#" << i
0238         //        << " name " << emItr->name()
0239         //        << " empty " << emItr->empty()
0240         //        << " rank " << emItr->rank()
0241         //        << " eta " << emItr->etaIndex()
0242         //        << " sign " << emItr->etaSign()
0243         //        << " phi " << emItr->phiIndex()
0244         //        << " iso " << emItr->isolated()
0245         //        << " bx " << emItr->bx()
0246         //        << endl ;
0247 
0248         if (!emItr->empty() && (!centralBxOnly_ || emItr->bx() == 0)) {
0249           double et = emScale->et(emItr->rank());
0250 
0251           //        cout << "L1Extra et " << et << endl ;
0252 
0253           isoEmColl->push_back(L1EmParticle(
0254               gctLorentzVector(et, *emItr, caloGeom, true), Ref<L1GctEmCandCollection>(hwIsoEmCands, i), emItr->bx()));
0255         }
0256       }
0257     }
0258 
0259     // Non-isolated EM
0260     Handle<L1GctEmCandCollection> hwNonIsoEmCands;
0261     iEvent.getByLabel(nonIsoEmSource_, hwNonIsoEmCands);
0262 
0263     if (!hwNonIsoEmCands.isValid()) {
0264       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << nonIsoEmSource_
0265                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0266     } else {
0267       //       cout << "HW non-iso EM" << endl ;
0268       L1GctEmCandCollection::const_iterator emItr = hwNonIsoEmCands->begin();
0269       L1GctEmCandCollection::const_iterator emEnd = hwNonIsoEmCands->end();
0270       for (int i = 0; emItr != emEnd; ++emItr, ++i) {
0271         //   cout << "#" << i
0272         //        << " name " << emItr->name()
0273         //        << " empty " << emItr->empty()
0274         //        << " rank " << emItr->rank()
0275         //        << " eta " << emItr->etaIndex()
0276         //        << " sign " << emItr->etaSign()
0277         //        << " phi " << emItr->phiIndex()
0278         //        << " iso " << emItr->isolated()
0279         //        << " bx " << emItr->bx()
0280         //        << endl ;
0281 
0282         if (!emItr->empty() && (!centralBxOnly_ || emItr->bx() == 0)) {
0283           double et = emScale->et(emItr->rank());
0284 
0285           //        cout << "L1Extra et " << et << endl ;
0286 
0287           nonIsoEmColl->push_back(L1EmParticle(gctLorentzVector(et, *emItr, caloGeom, true),
0288                                                Ref<L1GctEmCandCollection>(hwNonIsoEmCands, i),
0289                                                emItr->bx()));
0290         }
0291       }
0292     }
0293 
0294     // ~~~~~~~~~~~~~~~~~~~~ Jets ~~~~~~~~~~~~~~~~~~~~
0295 
0296     ESHandle<L1CaloEtScale> jetScale = iSetup.getHandle(jetScaleToken_);
0297 
0298     // Central jets.
0299     Handle<L1GctJetCandCollection> hwCenJetCands;
0300     iEvent.getByLabel(cenJetSource_, hwCenJetCands);
0301 
0302     if (!hwCenJetCands.isValid()) {
0303       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << cenJetSource_
0304                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0305     } else {
0306       //       cout << "HW central jets" << endl ;
0307       L1GctJetCandCollection::const_iterator jetItr = hwCenJetCands->begin();
0308       L1GctJetCandCollection::const_iterator jetEnd = hwCenJetCands->end();
0309       for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
0310         //   cout << "#" << i
0311         //        << " name " << jetItr->name()
0312         //        << " empty " << jetItr->empty()
0313         //        << " rank " << jetItr->rank()
0314         //        << " eta " << jetItr->etaIndex()
0315         //        << " sign " << jetItr->etaSign()
0316         //        << " phi " << jetItr->phiIndex()
0317         //        << " cen " << jetItr->isCentral()
0318         //        << " for " << jetItr->isForward()
0319         //        << " tau " << jetItr->isTau()
0320         //        << " bx " << jetItr->bx()
0321         //        << endl ;
0322 
0323         if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
0324           double et = jetScale->et(jetItr->rank());
0325 
0326           //        cout << "L1Extra et " << et << endl ;
0327 
0328           cenJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
0329                                               Ref<L1GctJetCandCollection>(hwCenJetCands, i),
0330                                               jetItr->bx()));
0331         }
0332       }
0333     }
0334 
0335     // Forward jets.
0336     Handle<L1GctJetCandCollection> hwForJetCands;
0337     iEvent.getByLabel(forJetSource_, hwForJetCands);
0338 
0339     if (!hwForJetCands.isValid()) {
0340       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << forJetSource_
0341                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0342     } else {
0343       //       cout << "HW forward jets" << endl ;
0344       L1GctJetCandCollection::const_iterator jetItr = hwForJetCands->begin();
0345       L1GctJetCandCollection::const_iterator jetEnd = hwForJetCands->end();
0346       for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
0347         //   cout << "#" << i
0348         //        << " name " << jetItr->name()
0349         //        << " empty " << jetItr->empty()
0350         //        << " rank " << jetItr->rank()
0351         //        << " eta " << jetItr->etaIndex()
0352         //        << " sign " << jetItr->etaSign()
0353         //        << " phi " << jetItr->phiIndex()
0354         //        << " cen " << jetItr->isCentral()
0355         //        << " for " << jetItr->isForward()
0356         //        << " tau " << jetItr->isTau()
0357         //        << " bx " << jetItr->bx()
0358         //        << endl ;
0359 
0360         if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
0361           double et = jetScale->et(jetItr->rank());
0362 
0363           //        cout << "L1Extra et " << et << endl ;
0364 
0365           forJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, false),
0366                                               Ref<L1GctJetCandCollection>(hwForJetCands, i),
0367                                               jetItr->bx()));
0368         }
0369       }
0370     }
0371 
0372     // Tau jets.
0373     //       cout << "HW tau jets" << endl ;
0374     Handle<L1GctJetCandCollection> hwTauJetCands;
0375     iEvent.getByLabel(tauJetSource_, hwTauJetCands);
0376 
0377     if (!hwTauJetCands.isValid()) {
0378       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << tauJetSource_
0379                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0380     } else {
0381       L1GctJetCandCollection::const_iterator jetItr = hwTauJetCands->begin();
0382       L1GctJetCandCollection::const_iterator jetEnd = hwTauJetCands->end();
0383       for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
0384         //   cout << "#" << i
0385         //        << " name " << jetItr->name()
0386         //        << " empty " << jetItr->empty()
0387         //        << " rank " << jetItr->rank()
0388         //        << " eta " << jetItr->etaIndex()
0389         //        << " sign " << jetItr->etaSign()
0390         //        << " phi " << jetItr->phiIndex()
0391         //        << " cen " << jetItr->isCentral()
0392         //        << " for " << jetItr->isForward()
0393         //        << " tau " << jetItr->isTau()
0394         //        << " bx " << jetItr->bx()
0395         //        << endl ;
0396 
0397         if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
0398           double et = jetScale->et(jetItr->rank());
0399 
0400           //        cout << "L1Extra et " << et << endl ;
0401 
0402           tauJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
0403                                               Ref<L1GctJetCandCollection>(hwTauJetCands, i),
0404                                               jetItr->bx()));
0405         }
0406       }
0407     }
0408 
0409     // Isolated Tau jets.
0410     //       cout << "HW tau jets" << endl ;
0411     Handle<L1GctJetCandCollection> hwIsoTauJetCands;
0412     iEvent.getByLabel(isoTauJetSource_, hwIsoTauJetCands);
0413 
0414     if (!hwIsoTauJetCands.isValid()) {
0415       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << isoTauJetSource_
0416                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0417     } else {
0418       L1GctJetCandCollection::const_iterator jetItr = hwIsoTauJetCands->begin();
0419       L1GctJetCandCollection::const_iterator jetEnd = hwIsoTauJetCands->end();
0420       for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
0421         //   cout << "#" << i
0422         //        << " name " << jetItr->name()
0423         //        << " empty " << jetItr->empty()
0424         //        << " rank " << jetItr->rank()
0425         //        << " eta " << jetItr->etaIndex()
0426         //        << " sign " << jetItr->etaSign()
0427         //        << " phi " << jetItr->phiIndex()
0428         //        << " cen " << jetItr->isCentral()
0429         //        << " for " << jetItr->isForward()
0430         //        << " tau " << jetItr->isTau()
0431         //        << " bx " << jetItr->bx()
0432         //        << endl ;
0433 
0434         if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
0435           double et = jetScale->et(jetItr->rank());
0436 
0437           //        cout << "L1Extra et " << et << endl ;
0438 
0439           isoTauJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
0440                                                  Ref<L1GctJetCandCollection>(hwIsoTauJetCands, i),
0441                                                  jetItr->bx()));
0442         }
0443       }
0444     }
0445 
0446     // ~~~~~~~~~~~~~~~~~~~~ ET Sums ~~~~~~~~~~~~~~~~~~~~
0447 
0448     double etSumLSB = jetScale->linearLsb();
0449 
0450     Handle<L1GctEtTotalCollection> hwEtTotColl;
0451     iEvent.getByLabel(etTotSource_, hwEtTotColl);
0452 
0453     Handle<L1GctEtMissCollection> hwEtMissColl;
0454     iEvent.getByLabel(etMissSource_, hwEtMissColl);
0455 
0456     if (!hwEtTotColl.isValid()) {
0457       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtTotalCollection with " << etTotSource_
0458                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0459     } else if (!hwEtMissColl.isValid()) {
0460       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtMissCollection with " << etMissSource_
0461                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0462     } else {
0463       // Make a L1EtMissParticle even if either L1GctEtTotal or L1GctEtMiss
0464       // is missing for a given bx.  Keep track of which L1GctEtMiss objects
0465       // have a corresponding L1GctEtTotal object.
0466       std::vector<bool> etMissMatched;
0467 
0468       L1GctEtMissCollection::const_iterator hwEtMissItr = hwEtMissColl->begin();
0469       L1GctEtMissCollection::const_iterator hwEtMissEnd = hwEtMissColl->end();
0470       for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr) {
0471         etMissMatched.push_back(false);
0472       }
0473 
0474       // Collate energy sums by bx
0475       L1GctEtTotalCollection::const_iterator hwEtTotItr = hwEtTotColl->begin();
0476       L1GctEtTotalCollection::const_iterator hwEtTotEnd = hwEtTotColl->end();
0477 
0478       int iTot = 0;
0479       for (; hwEtTotItr != hwEtTotEnd; ++hwEtTotItr, ++iTot) {
0480         int bx = hwEtTotItr->bx();
0481 
0482         if (!centralBxOnly_ || bx == 0) {
0483           // ET bin low edge
0484           double etTot =
0485               (hwEtTotItr->overFlow() ? (double)L1GctEtTotal::kEtTotalMaxValue : (double)hwEtTotItr->et()) * etSumLSB +
0486               1.e-6;
0487 
0488           int iMiss = 0;
0489           hwEtMissItr = hwEtMissColl->begin();
0490           hwEtMissEnd = hwEtMissColl->end();
0491           for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr, ++iMiss) {
0492             if (hwEtMissItr->bx() == bx) {
0493               etMissMatched[iMiss] = true;
0494               break;
0495             }
0496           }
0497 
0498           double etMiss = 0.;
0499           double phi = 0.;
0500           math::PtEtaPhiMLorentzVector p4;
0501           Ref<L1GctEtMissCollection> metRef;
0502 
0503           // If a L1GctEtMiss with the right bx is not found, itr == end.
0504           if (hwEtMissItr != hwEtMissEnd) {
0505             // ET bin low edge
0506             etMiss = (hwEtMissItr->overFlow() ? (double)L1GctEtMiss::kEtMissMaxValue : (double)hwEtMissItr->et()) *
0507                          etSumLSB +
0508                      1.e-6;
0509             // keep x and y components non-zero and
0510             // protect against roundoff.
0511 
0512             phi = caloGeom->etSumPhiBinCenter(hwEtMissItr->phi());
0513 
0514             p4 = math::PtEtaPhiMLorentzVector(etMiss, 0., phi, 0.);
0515 
0516             metRef = Ref<L1GctEtMissCollection>(hwEtMissColl, iMiss);
0517 
0518             //         cout << "HW ET Sums " << endl
0519             //          << "MET: phi " << hwEtMissItr->phi() << " = " << phi
0520             //          << " et " << hwEtMissItr->et() << " = " << etMiss
0521             //          << " EtTot " << hwEtTotItr->et() << " = " << etTot
0522             //          << " bx " << bx
0523             //          << endl ;
0524           }
0525           //       else
0526           //         {
0527           //           cout << "HW ET Sums " << endl
0528           //            << "MET: phi " << phi
0529           //            << " et "<< etMiss
0530           //            << " EtTot " << hwEtTotItr->et() << " = " << etTot
0531           //            << " bx " << bx
0532           //            << endl ;
0533           //         }
0534 
0535           etMissColl->push_back(L1EtMissParticle(p4,
0536                                                  L1EtMissParticle::kMET,
0537                                                  etTot,
0538                                                  metRef,
0539                                                  Ref<L1GctEtTotalCollection>(hwEtTotColl, iTot),
0540                                                  Ref<L1GctHtMissCollection>(),
0541                                                  Ref<L1GctEtHadCollection>(),
0542                                                  bx));
0543         }
0544       }
0545 
0546       if (!centralBxOnly_) {
0547         // Make L1EtMissParticles for those L1GctEtMiss objects without
0548         // a matched L1GctEtTotal object.
0549 
0550         double etTot = 0.;
0551 
0552         hwEtMissItr = hwEtMissColl->begin();
0553         hwEtMissEnd = hwEtMissColl->end();
0554         int iMiss = 0;
0555         for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr, ++iMiss) {
0556           if (!etMissMatched[iMiss]) {
0557             int bx = hwEtMissItr->bx();
0558 
0559             // ET bin low edge
0560             double etMiss =
0561                 (hwEtMissItr->overFlow() ? (double)L1GctEtMiss::kEtMissMaxValue : (double)hwEtMissItr->et()) *
0562                     etSumLSB +
0563                 1.e-6;
0564             // keep x and y components non-zero and
0565             // protect against roundoff.
0566 
0567             double phi = caloGeom->etSumPhiBinCenter(hwEtMissItr->phi());
0568 
0569             math::PtEtaPhiMLorentzVector p4(etMiss, 0., phi, 0.);
0570 
0571             //        cout << "HW ET Sums " << endl
0572             //             << "MET: phi " << hwEtMissItr->phi() << " = " <<
0573             // phi
0574             //             << " et " << hwEtMissItr->et() << " = " << etMiss
0575             //             << " EtTot " << etTot
0576             //             << " bx " << bx
0577             //             << endl ;
0578 
0579             etMissColl->push_back(L1EtMissParticle(p4,
0580                                                    L1EtMissParticle::kMET,
0581                                                    etTot,
0582                                                    Ref<L1GctEtMissCollection>(hwEtMissColl, iMiss),
0583                                                    Ref<L1GctEtTotalCollection>(),
0584                                                    Ref<L1GctHtMissCollection>(),
0585                                                    Ref<L1GctEtHadCollection>(),
0586                                                    bx));
0587           }
0588         }
0589       }
0590     }
0591 
0592     // ~~~~~~~~~~~~~~~~~~~~ HT Sums ~~~~~~~~~~~~~~~~~~~~
0593 
0594     Handle<L1GctEtHadCollection> hwEtHadColl;
0595     iEvent.getByLabel(etHadSource_, hwEtHadColl);
0596 
0597     Handle<L1GctHtMissCollection> hwHtMissColl;
0598     if (!ignoreHtMiss_) {
0599       iEvent.getByLabel(htMissSource_, hwHtMissColl);
0600     }
0601 
0602     ESHandle<L1GctJetFinderParams> jetFinderParams = iSetup.getHandle(jetFinderParamsToken_);
0603     double htSumLSB = jetFinderParams->getHtLsbGeV();
0604 
0605     ESHandle<L1CaloEtScale> htMissScale;
0606     std::vector<bool> htMissMatched;
0607     if (!ignoreHtMiss_) {
0608       htMissScale = iSetup.getHandle(htMissScaleToken_);
0609 
0610       if (!hwEtHadColl.isValid()) {
0611         LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtHadCollection with " << etHadSource_
0612                                          << "\nrequested in configuration, but not found in the event." << std::endl;
0613       } else if (!hwHtMissColl.isValid()) {
0614         LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHtMissCollection with " << htMissSource_
0615                                          << "\nrequested in configuration, but not found in the event." << std::endl;
0616       } else {
0617         // Make a L1EtMissParticle even if either L1GctEtHad or L1GctHtMiss
0618         // is missing for a given bx. Keep track of which L1GctHtMiss objects
0619         // have a corresponding L1GctHtTotal object.
0620         L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
0621         L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
0622         for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr) {
0623           htMissMatched.push_back(false);
0624         }
0625       }
0626     }
0627 
0628     if (!hwEtHadColl.isValid()) {
0629       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtHadCollection with " << etHadSource_
0630                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0631     } else if (!hwHtMissColl.isValid()) {
0632       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHtMissCollection with " << htMissSource_
0633                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0634     } else {
0635       L1GctEtHadCollection::const_iterator hwEtHadItr = hwEtHadColl->begin();
0636       L1GctEtHadCollection::const_iterator hwEtHadEnd = hwEtHadColl->end();
0637 
0638       int iHad = 0;
0639       for (; hwEtHadItr != hwEtHadEnd; ++hwEtHadItr, ++iHad) {
0640         int bx = hwEtHadItr->bx();
0641 
0642         if (!centralBxOnly_ || bx == 0) {
0643           // HT bin low edge
0644           double htTot =
0645               (hwEtHadItr->overFlow() ? (double)L1GctEtHad::kEtHadMaxValue : (double)hwEtHadItr->et()) * htSumLSB +
0646               1.e-6;
0647 
0648           double htMiss = 0.;
0649           double phi = 0.;
0650           math::PtEtaPhiMLorentzVector p4;
0651           Ref<L1GctHtMissCollection> mhtRef;
0652 
0653           if (!ignoreHtMiss_) {
0654             L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
0655             L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
0656 
0657             int iMiss = 0;
0658             for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr, ++iMiss) {
0659               if (hwHtMissItr->bx() == bx) {
0660                 htMissMatched[iMiss] = true;
0661                 break;
0662               }
0663             }
0664 
0665             // If a L1GctHtMiss with the right bx is not found, itr == end.
0666             if (hwHtMissItr != hwHtMissEnd) {
0667               // HT bin low edge
0668               htMiss =
0669                   htMissScale->et(hwHtMissItr->overFlow() ? htMissScale->rankScaleMax() : hwHtMissItr->et()) + 1.e-6;
0670               // keep x and y components non-zero and
0671               // protect against roundoff.
0672 
0673               phi = caloGeom->htSumPhiBinCenter(hwHtMissItr->phi());
0674 
0675               p4 = math::PtEtaPhiMLorentzVector(htMiss, 0., phi, 0.);
0676 
0677               mhtRef = Ref<L1GctHtMissCollection>(hwHtMissColl, iMiss);
0678 
0679               //           cout << "HW HT Sums " << endl
0680               //            << "MHT: phi " << hwHtMissItr->phi() <<
0681               // " = " << phi
0682               //            << " ht " << hwHtMissItr->et() << " = "
0683               // << htMiss
0684               //            << " HtTot " << hwEtHadItr->et() << " =
0685               // "
0686               // << htTot
0687               //            << " bx " << bx
0688               //            << endl ;
0689             }
0690             //         else
0691             //       {
0692             //         cout << "HW HT Sums " << endl
0693             //          << "MHT: phi " << phi
0694             //          << " ht " << htMiss
0695             //          << " HtTot " << hwEtHadItr->et() << " = " <<
0696             // htTot
0697             //          << " bx " << bx
0698             //          << endl ;
0699             //       }
0700           }
0701 
0702           htMissColl->push_back(L1EtMissParticle(p4,
0703                                                  L1EtMissParticle::kMHT,
0704                                                  htTot,
0705                                                  Ref<L1GctEtMissCollection>(),
0706                                                  Ref<L1GctEtTotalCollection>(),
0707                                                  mhtRef,
0708                                                  Ref<L1GctEtHadCollection>(hwEtHadColl, iHad),
0709                                                  bx));
0710         }
0711       }
0712 
0713       if (!centralBxOnly_ && !ignoreHtMiss_) {
0714         // Make L1EtMissParticles for those L1GctHtMiss objects without
0715         // a matched L1GctHtTotal object.
0716         double htTot = 0.;
0717 
0718         L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
0719         L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
0720 
0721         int iMiss = 0;
0722         for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr, ++iMiss) {
0723           if (!htMissMatched[iMiss]) {
0724             int bx = hwHtMissItr->bx();
0725 
0726             // HT bin low edge
0727             double htMiss =
0728                 htMissScale->et(hwHtMissItr->overFlow() ? htMissScale->rankScaleMax() : hwHtMissItr->et()) + 1.e-6;
0729             // keep x and y components non-zero and
0730             // protect against roundoff.
0731 
0732             double phi = caloGeom->htSumPhiBinCenter(hwHtMissItr->phi());
0733 
0734             math::PtEtaPhiMLorentzVector p4(htMiss, 0., phi, 0.);
0735 
0736             //        cout << "HW HT Sums " << endl
0737             //             << "MHT: phi " << hwHtMissItr->phi() << " = " <<
0738             // phi
0739             //             << " ht " << hwHtMissItr->et() << " = " << htMiss
0740             //             << " HtTot " << htTot
0741             //             << " bx " << bx
0742             //             << endl ;
0743 
0744             htMissColl->push_back(L1EtMissParticle(p4,
0745                                                    L1EtMissParticle::kMHT,
0746                                                    htTot,
0747                                                    Ref<L1GctEtMissCollection>(),
0748                                                    Ref<L1GctEtTotalCollection>(),
0749                                                    Ref<L1GctHtMissCollection>(hwHtMissColl, iMiss),
0750                                                    Ref<L1GctEtHadCollection>(),
0751                                                    bx));
0752           }
0753         }
0754       }
0755     }
0756 
0757     // ~~~~~~~~~~~~~~~~~~~~ HF Rings ~~~~~~~~~~~~~~~~~~~~
0758 
0759     Handle<L1GctHFRingEtSumsCollection> hwHFEtSumsColl;
0760     iEvent.getByLabel(hfRingEtSumsSource_, hwHFEtSumsColl);
0761 
0762     Handle<L1GctHFBitCountsCollection> hwHFBitCountsColl;
0763     iEvent.getByLabel(hfRingBitCountsSource_, hwHFBitCountsColl);
0764 
0765     ESHandle<L1CaloEtScale> hfRingEtScale = iSetup.getHandle(hfRingEtScaleToken_);
0766 
0767     if (!hwHFEtSumsColl.isValid()) {
0768       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHFRingEtSumsCollection with " << hfRingEtSumsSource_
0769                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0770     } else if (!hwHFBitCountsColl.isValid()) {
0771       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHFBitCountsCollection with " << hfRingBitCountsSource_
0772                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0773     } else {
0774       L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsItr = hwHFEtSumsColl->begin();
0775       L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsEnd = hwHFEtSumsColl->end();
0776 
0777       int iEtSums = 0;
0778       for (; hwHFEtSumsItr != hwHFEtSumsEnd; ++hwHFEtSumsItr, ++iEtSums) {
0779         int bx = hwHFEtSumsItr->bx();
0780 
0781         if (!centralBxOnly_ || bx == 0) {
0782           L1GctHFBitCountsCollection::const_iterator hwHFBitCountsItr = hwHFBitCountsColl->begin();
0783           L1GctHFBitCountsCollection::const_iterator hwHFBitCountsEnd = hwHFBitCountsColl->end();
0784 
0785           int iBitCounts = 0;
0786           for (; hwHFBitCountsItr != hwHFBitCountsEnd; ++hwHFBitCountsItr, ++iBitCounts) {
0787             if (hwHFBitCountsItr->bx() == bx) {
0788               break;
0789             }
0790           }
0791 
0792           // If a L1GctHFBitCounts with the right bx is not found, itr == end.
0793           if (hwHFBitCountsItr != hwHFBitCountsEnd) {
0794             // Construct L1HFRings only if both HW objects are present.
0795 
0796             // cout << "HF Rings " << endl
0797 
0798             // HF Et sums
0799             double etSums[L1HFRings::kNumRings];
0800             etSums[L1HFRings::kRing1PosEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(0)) + 1.e-6;
0801             etSums[L1HFRings::kRing1NegEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(1)) + 1.e-6;
0802             etSums[L1HFRings::kRing2PosEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(2)) + 1.e-6;
0803             etSums[L1HFRings::kRing2NegEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(3)) + 1.e-6;
0804             // protect against roundoff.
0805 
0806             //         cout << "HF Et Sums "
0807             //          << hwHFEtSumsItr->etSum( 0 ) << " = "
0808             //          << etSums[ L1HFRings::kRing1PosEta ] << " "
0809             //          << hwHFEtSumsItr->etSum( 1 ) << " = "
0810             //          << etSums[ L1HFRings::kRing1NegEta ] << " "
0811             //          << hwHFEtSumsItr->etSum( 2 ) << " = "
0812             //          << etSums[ L1HFRings::kRing2PosEta ] << " "
0813             //          << hwHFEtSumsItr->etSum( 3 ) << " = "
0814             //          << etSums[ L1HFRings::kRing2NegEta ]
0815             //          << endl ;
0816 
0817             // HF bit counts
0818             int bitCounts[L1HFRings::kNumRings];
0819             bitCounts[L1HFRings::kRing1PosEta] = hwHFBitCountsItr->bitCount(0);
0820             bitCounts[L1HFRings::kRing1NegEta] = hwHFBitCountsItr->bitCount(1);
0821             bitCounts[L1HFRings::kRing2PosEta] = hwHFBitCountsItr->bitCount(2);
0822             bitCounts[L1HFRings::kRing2NegEta] = hwHFBitCountsItr->bitCount(3);
0823 
0824             //         cout << "HF bit counts "
0825             //          << hwHFBitCountsItr->bitCount( 0 ) << " "
0826             //          << hwHFBitCountsItr->bitCount( 1 ) << " "
0827             //          << hwHFBitCountsItr->bitCount( 2 ) << " "
0828             //          << hwHFBitCountsItr->bitCount( 3 ) << " "
0829             //          << endl ;
0830 
0831             hfRingsColl->push_back(L1HFRings(etSums,
0832                                              bitCounts,
0833                                              Ref<L1GctHFRingEtSumsCollection>(hwHFEtSumsColl, iEtSums),
0834                                              Ref<L1GctHFBitCountsCollection>(hwHFBitCountsColl, iBitCounts),
0835                                              bx));
0836           }
0837         }
0838       }
0839     }
0840   }
0841 
0842   iEvent.put(std::move(isoEmColl), "Isolated");
0843   iEvent.put(std::move(nonIsoEmColl), "NonIsolated");
0844   iEvent.put(std::move(cenJetColl), "Central");
0845   iEvent.put(std::move(forJetColl), "Forward");
0846   iEvent.put(std::move(tauJetColl), "Tau");
0847   iEvent.put(std::move(isoTauJetColl), "IsoTau");
0848   iEvent.put(std::move(etMissColl), "MET");
0849   iEvent.put(std::move(htMissColl), "MHT");
0850   iEvent.put(std::move(hfRingsColl));
0851 }
0852 
0853 // math::XYZTLorentzVector
0854 math::PtEtaPhiMLorentzVector L1ExtraParticlesProd::gctLorentzVector(const double &et,
0855                                                                     const L1GctCand &cand,
0856                                                                     const L1CaloGeometry *geom,
0857                                                                     bool central) {
0858   // To keep x and y components non-zero.
0859   double etCorr = et + 1.e-6;  // protect against roundoff, not only for et=0
0860 
0861   double eta = geom->etaBinCenter(cand.etaIndex(), central);
0862 
0863   //    double tanThOver2 = exp( -eta ) ;
0864   //    double ez = etCorr * ( 1. - tanThOver2 * tanThOver2 ) / ( 2. *
0865   //    tanThOver2 ); double e  = etCorr * ( 1. + tanThOver2 * tanThOver2 ) /
0866   //    ( 2. * tanThOver2 );
0867 
0868   double phi = geom->emJetPhiBinCenter(cand.phiIndex());
0869 
0870   //    return math::XYZTLorentzVector( etCorr * cos( phi ),
0871   //                   etCorr * sin( phi ),
0872   //                   ez,
0873   //                   e ) ;
0874   return math::PtEtaPhiMLorentzVector(etCorr, eta, phi, 0.);
0875 }
0876 
0877 // define this as a plug-in
0878 DEFINE_FWK_MODULE(L1ExtraParticlesProd);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!