Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:20:12 

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   OrphanHandle<L1MuonParticleCollection> muHandle = iEvent.put(std::move(muColl));
0199 
0200   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0201   // ~~~~~~~~~~~~~~~~~~~~ Calorimeter ~~~~~~~~~~~~~~~~~~~~
0202   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0203 
0204   unique_ptr<L1EmParticleCollection> isoEmColl(new L1EmParticleCollection);
0205 
0206   unique_ptr<L1EmParticleCollection> nonIsoEmColl(new L1EmParticleCollection);
0207 
0208   unique_ptr<L1JetParticleCollection> cenJetColl(new L1JetParticleCollection);
0209 
0210   unique_ptr<L1JetParticleCollection> forJetColl(new L1JetParticleCollection);
0211 
0212   unique_ptr<L1JetParticleCollection> tauJetColl(new L1JetParticleCollection);
0213 
0214   unique_ptr<L1JetParticleCollection> isoTauJetColl(new L1JetParticleCollection);
0215 
0216   unique_ptr<L1EtMissParticleCollection> etMissColl(new L1EtMissParticleCollection);
0217 
0218   unique_ptr<L1EtMissParticleCollection> htMissColl(new L1EtMissParticleCollection);
0219 
0220   unique_ptr<L1HFRingsCollection> hfRingsColl(new L1HFRingsCollection);
0221 
0222   if (produceCaloParticles_) {
0223     // ~~~~~~~~~~~~~~~~~~~~ Geometry ~~~~~~~~~~~~~~~~~~~~
0224 
0225     ESHandle<L1CaloGeometry> caloGeomESH = iSetup.getHandle(caloGeomToken_);
0226     const L1CaloGeometry *caloGeom = &(*caloGeomESH);
0227 
0228     // ~~~~~~~~~~~~~~~~~~~~ EM ~~~~~~~~~~~~~~~~~~~~
0229 
0230     ESHandle<L1CaloEtScale> emScale = iSetup.getHandle(emScaleToken_);
0231 
0232     // Isolated EM
0233     Handle<L1GctEmCandCollection> hwIsoEmCands;
0234     iEvent.getByLabel(isoEmSource_, hwIsoEmCands);
0235 
0236     if (!hwIsoEmCands.isValid()) {
0237       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << isoEmSource_
0238                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0239     } else {
0240       //       cout << "HW iso EM" << endl ;
0241 
0242       L1GctEmCandCollection::const_iterator emItr = hwIsoEmCands->begin();
0243       L1GctEmCandCollection::const_iterator emEnd = hwIsoEmCands->end();
0244       for (int i = 0; emItr != emEnd; ++emItr, ++i) {
0245         //   cout << "#" << i
0246         //        << " name " << emItr->name()
0247         //        << " empty " << emItr->empty()
0248         //        << " rank " << emItr->rank()
0249         //        << " eta " << emItr->etaIndex()
0250         //        << " sign " << emItr->etaSign()
0251         //        << " phi " << emItr->phiIndex()
0252         //        << " iso " << emItr->isolated()
0253         //        << " bx " << emItr->bx()
0254         //        << endl ;
0255 
0256         if (!emItr->empty() && (!centralBxOnly_ || emItr->bx() == 0)) {
0257           double et = emScale->et(emItr->rank());
0258 
0259           //        cout << "L1Extra et " << et << endl ;
0260 
0261           isoEmColl->push_back(L1EmParticle(
0262               gctLorentzVector(et, *emItr, caloGeom, true), Ref<L1GctEmCandCollection>(hwIsoEmCands, i), emItr->bx()));
0263         }
0264       }
0265     }
0266 
0267     // Non-isolated EM
0268     Handle<L1GctEmCandCollection> hwNonIsoEmCands;
0269     iEvent.getByLabel(nonIsoEmSource_, hwNonIsoEmCands);
0270 
0271     if (!hwNonIsoEmCands.isValid()) {
0272       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << nonIsoEmSource_
0273                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0274     } else {
0275       //       cout << "HW non-iso EM" << endl ;
0276       L1GctEmCandCollection::const_iterator emItr = hwNonIsoEmCands->begin();
0277       L1GctEmCandCollection::const_iterator emEnd = hwNonIsoEmCands->end();
0278       for (int i = 0; emItr != emEnd; ++emItr, ++i) {
0279         //   cout << "#" << i
0280         //        << " name " << emItr->name()
0281         //        << " empty " << emItr->empty()
0282         //        << " rank " << emItr->rank()
0283         //        << " eta " << emItr->etaIndex()
0284         //        << " sign " << emItr->etaSign()
0285         //        << " phi " << emItr->phiIndex()
0286         //        << " iso " << emItr->isolated()
0287         //        << " bx " << emItr->bx()
0288         //        << endl ;
0289 
0290         if (!emItr->empty() && (!centralBxOnly_ || emItr->bx() == 0)) {
0291           double et = emScale->et(emItr->rank());
0292 
0293           //        cout << "L1Extra et " << et << endl ;
0294 
0295           nonIsoEmColl->push_back(L1EmParticle(gctLorentzVector(et, *emItr, caloGeom, true),
0296                                                Ref<L1GctEmCandCollection>(hwNonIsoEmCands, i),
0297                                                emItr->bx()));
0298         }
0299       }
0300     }
0301 
0302     // ~~~~~~~~~~~~~~~~~~~~ Jets ~~~~~~~~~~~~~~~~~~~~
0303 
0304     ESHandle<L1CaloEtScale> jetScale = iSetup.getHandle(jetScaleToken_);
0305 
0306     // Central jets.
0307     Handle<L1GctJetCandCollection> hwCenJetCands;
0308     iEvent.getByLabel(cenJetSource_, hwCenJetCands);
0309 
0310     if (!hwCenJetCands.isValid()) {
0311       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << cenJetSource_
0312                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0313     } else {
0314       //       cout << "HW central jets" << endl ;
0315       L1GctJetCandCollection::const_iterator jetItr = hwCenJetCands->begin();
0316       L1GctJetCandCollection::const_iterator jetEnd = hwCenJetCands->end();
0317       for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
0318         //   cout << "#" << i
0319         //        << " name " << jetItr->name()
0320         //        << " empty " << jetItr->empty()
0321         //        << " rank " << jetItr->rank()
0322         //        << " eta " << jetItr->etaIndex()
0323         //        << " sign " << jetItr->etaSign()
0324         //        << " phi " << jetItr->phiIndex()
0325         //        << " cen " << jetItr->isCentral()
0326         //        << " for " << jetItr->isForward()
0327         //        << " tau " << jetItr->isTau()
0328         //        << " bx " << jetItr->bx()
0329         //        << endl ;
0330 
0331         if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
0332           double et = jetScale->et(jetItr->rank());
0333 
0334           //        cout << "L1Extra et " << et << endl ;
0335 
0336           cenJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
0337                                               Ref<L1GctJetCandCollection>(hwCenJetCands, i),
0338                                               jetItr->bx()));
0339         }
0340       }
0341     }
0342 
0343     // Forward jets.
0344     Handle<L1GctJetCandCollection> hwForJetCands;
0345     iEvent.getByLabel(forJetSource_, hwForJetCands);
0346 
0347     if (!hwForJetCands.isValid()) {
0348       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEmCandCollection with " << forJetSource_
0349                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0350     } else {
0351       //       cout << "HW forward jets" << endl ;
0352       L1GctJetCandCollection::const_iterator jetItr = hwForJetCands->begin();
0353       L1GctJetCandCollection::const_iterator jetEnd = hwForJetCands->end();
0354       for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
0355         //   cout << "#" << i
0356         //        << " name " << jetItr->name()
0357         //        << " empty " << jetItr->empty()
0358         //        << " rank " << jetItr->rank()
0359         //        << " eta " << jetItr->etaIndex()
0360         //        << " sign " << jetItr->etaSign()
0361         //        << " phi " << jetItr->phiIndex()
0362         //        << " cen " << jetItr->isCentral()
0363         //        << " for " << jetItr->isForward()
0364         //        << " tau " << jetItr->isTau()
0365         //        << " bx " << jetItr->bx()
0366         //        << endl ;
0367 
0368         if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
0369           double et = jetScale->et(jetItr->rank());
0370 
0371           //        cout << "L1Extra et " << et << endl ;
0372 
0373           forJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, false),
0374                                               Ref<L1GctJetCandCollection>(hwForJetCands, i),
0375                                               jetItr->bx()));
0376         }
0377       }
0378     }
0379 
0380     // Tau jets.
0381     //       cout << "HW tau jets" << endl ;
0382     Handle<L1GctJetCandCollection> hwTauJetCands;
0383     iEvent.getByLabel(tauJetSource_, hwTauJetCands);
0384 
0385     if (!hwTauJetCands.isValid()) {
0386       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << tauJetSource_
0387                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0388     } else {
0389       L1GctJetCandCollection::const_iterator jetItr = hwTauJetCands->begin();
0390       L1GctJetCandCollection::const_iterator jetEnd = hwTauJetCands->end();
0391       for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
0392         //   cout << "#" << i
0393         //        << " name " << jetItr->name()
0394         //        << " empty " << jetItr->empty()
0395         //        << " rank " << jetItr->rank()
0396         //        << " eta " << jetItr->etaIndex()
0397         //        << " sign " << jetItr->etaSign()
0398         //        << " phi " << jetItr->phiIndex()
0399         //        << " cen " << jetItr->isCentral()
0400         //        << " for " << jetItr->isForward()
0401         //        << " tau " << jetItr->isTau()
0402         //        << " bx " << jetItr->bx()
0403         //        << endl ;
0404 
0405         if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
0406           double et = jetScale->et(jetItr->rank());
0407 
0408           //        cout << "L1Extra et " << et << endl ;
0409 
0410           tauJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
0411                                               Ref<L1GctJetCandCollection>(hwTauJetCands, i),
0412                                               jetItr->bx()));
0413         }
0414       }
0415     }
0416 
0417     // Isolated Tau jets.
0418     //       cout << "HW tau jets" << endl ;
0419     Handle<L1GctJetCandCollection> hwIsoTauJetCands;
0420     iEvent.getByLabel(isoTauJetSource_, hwIsoTauJetCands);
0421 
0422     if (!hwIsoTauJetCands.isValid()) {
0423       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctJetCandCollection with " << isoTauJetSource_
0424                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0425     } else {
0426       L1GctJetCandCollection::const_iterator jetItr = hwIsoTauJetCands->begin();
0427       L1GctJetCandCollection::const_iterator jetEnd = hwIsoTauJetCands->end();
0428       for (int i = 0; jetItr != jetEnd; ++jetItr, ++i) {
0429         //   cout << "#" << i
0430         //        << " name " << jetItr->name()
0431         //        << " empty " << jetItr->empty()
0432         //        << " rank " << jetItr->rank()
0433         //        << " eta " << jetItr->etaIndex()
0434         //        << " sign " << jetItr->etaSign()
0435         //        << " phi " << jetItr->phiIndex()
0436         //        << " cen " << jetItr->isCentral()
0437         //        << " for " << jetItr->isForward()
0438         //        << " tau " << jetItr->isTau()
0439         //        << " bx " << jetItr->bx()
0440         //        << endl ;
0441 
0442         if (!jetItr->empty() && (!centralBxOnly_ || jetItr->bx() == 0)) {
0443           double et = jetScale->et(jetItr->rank());
0444 
0445           //        cout << "L1Extra et " << et << endl ;
0446 
0447           isoTauJetColl->push_back(L1JetParticle(gctLorentzVector(et, *jetItr, caloGeom, true),
0448                                                  Ref<L1GctJetCandCollection>(hwIsoTauJetCands, i),
0449                                                  jetItr->bx()));
0450         }
0451       }
0452     }
0453 
0454     // ~~~~~~~~~~~~~~~~~~~~ ET Sums ~~~~~~~~~~~~~~~~~~~~
0455 
0456     double etSumLSB = jetScale->linearLsb();
0457 
0458     Handle<L1GctEtTotalCollection> hwEtTotColl;
0459     iEvent.getByLabel(etTotSource_, hwEtTotColl);
0460 
0461     Handle<L1GctEtMissCollection> hwEtMissColl;
0462     iEvent.getByLabel(etMissSource_, hwEtMissColl);
0463 
0464     if (!hwEtTotColl.isValid()) {
0465       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtTotalCollection with " << etTotSource_
0466                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0467     } else if (!hwEtMissColl.isValid()) {
0468       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtMissCollection with " << etMissSource_
0469                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0470     } else {
0471       // Make a L1EtMissParticle even if either L1GctEtTotal or L1GctEtMiss
0472       // is missing for a given bx.  Keep track of which L1GctEtMiss objects
0473       // have a corresponding L1GctEtTotal object.
0474       std::vector<bool> etMissMatched;
0475 
0476       L1GctEtMissCollection::const_iterator hwEtMissItr = hwEtMissColl->begin();
0477       L1GctEtMissCollection::const_iterator hwEtMissEnd = hwEtMissColl->end();
0478       for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr) {
0479         etMissMatched.push_back(false);
0480       }
0481 
0482       // Collate energy sums by bx
0483       L1GctEtTotalCollection::const_iterator hwEtTotItr = hwEtTotColl->begin();
0484       L1GctEtTotalCollection::const_iterator hwEtTotEnd = hwEtTotColl->end();
0485 
0486       int iTot = 0;
0487       for (; hwEtTotItr != hwEtTotEnd; ++hwEtTotItr, ++iTot) {
0488         int bx = hwEtTotItr->bx();
0489 
0490         if (!centralBxOnly_ || bx == 0) {
0491           // ET bin low edge
0492           double etTot =
0493               (hwEtTotItr->overFlow() ? (double)L1GctEtTotal::kEtTotalMaxValue : (double)hwEtTotItr->et()) * etSumLSB +
0494               1.e-6;
0495 
0496           int iMiss = 0;
0497           hwEtMissItr = hwEtMissColl->begin();
0498           hwEtMissEnd = hwEtMissColl->end();
0499           for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr, ++iMiss) {
0500             if (hwEtMissItr->bx() == bx) {
0501               etMissMatched[iMiss] = true;
0502               break;
0503             }
0504           }
0505 
0506           double etMiss = 0.;
0507           double phi = 0.;
0508           math::PtEtaPhiMLorentzVector p4;
0509           Ref<L1GctEtMissCollection> metRef;
0510 
0511           // If a L1GctEtMiss with the right bx is not found, itr == end.
0512           if (hwEtMissItr != hwEtMissEnd) {
0513             // ET bin low edge
0514             etMiss = (hwEtMissItr->overFlow() ? (double)L1GctEtMiss::kEtMissMaxValue : (double)hwEtMissItr->et()) *
0515                          etSumLSB +
0516                      1.e-6;
0517             // keep x and y components non-zero and
0518             // protect against roundoff.
0519 
0520             phi = caloGeom->etSumPhiBinCenter(hwEtMissItr->phi());
0521 
0522             p4 = math::PtEtaPhiMLorentzVector(etMiss, 0., phi, 0.);
0523 
0524             metRef = Ref<L1GctEtMissCollection>(hwEtMissColl, iMiss);
0525 
0526             //         cout << "HW ET Sums " << endl
0527             //          << "MET: phi " << hwEtMissItr->phi() << " = " << phi
0528             //          << " et " << hwEtMissItr->et() << " = " << etMiss
0529             //          << " EtTot " << hwEtTotItr->et() << " = " << etTot
0530             //          << " bx " << bx
0531             //          << endl ;
0532           }
0533           //       else
0534           //         {
0535           //           cout << "HW ET Sums " << endl
0536           //            << "MET: phi " << phi
0537           //            << " et "<< etMiss
0538           //            << " EtTot " << hwEtTotItr->et() << " = " << etTot
0539           //            << " bx " << bx
0540           //            << endl ;
0541           //         }
0542 
0543           etMissColl->push_back(L1EtMissParticle(p4,
0544                                                  L1EtMissParticle::kMET,
0545                                                  etTot,
0546                                                  metRef,
0547                                                  Ref<L1GctEtTotalCollection>(hwEtTotColl, iTot),
0548                                                  Ref<L1GctHtMissCollection>(),
0549                                                  Ref<L1GctEtHadCollection>(),
0550                                                  bx));
0551         }
0552       }
0553 
0554       if (!centralBxOnly_) {
0555         // Make L1EtMissParticles for those L1GctEtMiss objects without
0556         // a matched L1GctEtTotal object.
0557 
0558         double etTot = 0.;
0559 
0560         hwEtMissItr = hwEtMissColl->begin();
0561         hwEtMissEnd = hwEtMissColl->end();
0562         int iMiss = 0;
0563         for (; hwEtMissItr != hwEtMissEnd; ++hwEtMissItr, ++iMiss) {
0564           if (!etMissMatched[iMiss]) {
0565             int bx = hwEtMissItr->bx();
0566 
0567             // ET bin low edge
0568             double etMiss =
0569                 (hwEtMissItr->overFlow() ? (double)L1GctEtMiss::kEtMissMaxValue : (double)hwEtMissItr->et()) *
0570                     etSumLSB +
0571                 1.e-6;
0572             // keep x and y components non-zero and
0573             // protect against roundoff.
0574 
0575             double phi = caloGeom->etSumPhiBinCenter(hwEtMissItr->phi());
0576 
0577             math::PtEtaPhiMLorentzVector p4(etMiss, 0., phi, 0.);
0578 
0579             //        cout << "HW ET Sums " << endl
0580             //             << "MET: phi " << hwEtMissItr->phi() << " = " <<
0581             // phi
0582             //             << " et " << hwEtMissItr->et() << " = " << etMiss
0583             //             << " EtTot " << etTot
0584             //             << " bx " << bx
0585             //             << endl ;
0586 
0587             etMissColl->push_back(L1EtMissParticle(p4,
0588                                                    L1EtMissParticle::kMET,
0589                                                    etTot,
0590                                                    Ref<L1GctEtMissCollection>(hwEtMissColl, iMiss),
0591                                                    Ref<L1GctEtTotalCollection>(),
0592                                                    Ref<L1GctHtMissCollection>(),
0593                                                    Ref<L1GctEtHadCollection>(),
0594                                                    bx));
0595           }
0596         }
0597       }
0598     }
0599 
0600     // ~~~~~~~~~~~~~~~~~~~~ HT Sums ~~~~~~~~~~~~~~~~~~~~
0601 
0602     Handle<L1GctEtHadCollection> hwEtHadColl;
0603     iEvent.getByLabel(etHadSource_, hwEtHadColl);
0604 
0605     Handle<L1GctHtMissCollection> hwHtMissColl;
0606     if (!ignoreHtMiss_) {
0607       iEvent.getByLabel(htMissSource_, hwHtMissColl);
0608     }
0609 
0610     ESHandle<L1GctJetFinderParams> jetFinderParams = iSetup.getHandle(jetFinderParamsToken_);
0611     double htSumLSB = jetFinderParams->getHtLsbGeV();
0612 
0613     ESHandle<L1CaloEtScale> htMissScale;
0614     std::vector<bool> htMissMatched;
0615     if (!ignoreHtMiss_) {
0616       htMissScale = iSetup.getHandle(htMissScaleToken_);
0617 
0618       if (!hwEtHadColl.isValid()) {
0619         LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtHadCollection with " << etHadSource_
0620                                          << "\nrequested in configuration, but not found in the event." << std::endl;
0621       } else if (!hwHtMissColl.isValid()) {
0622         LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHtMissCollection with " << htMissSource_
0623                                          << "\nrequested in configuration, but not found in the event." << std::endl;
0624       } else {
0625         // Make a L1EtMissParticle even if either L1GctEtHad or L1GctHtMiss
0626         // is missing for a given bx. Keep track of which L1GctHtMiss objects
0627         // have a corresponding L1GctHtTotal object.
0628         L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
0629         L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
0630         for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr) {
0631           htMissMatched.push_back(false);
0632         }
0633       }
0634     }
0635 
0636     if (!hwEtHadColl.isValid()) {
0637       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctEtHadCollection with " << etHadSource_
0638                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0639     } else if (!hwHtMissColl.isValid()) {
0640       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHtMissCollection with " << htMissSource_
0641                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0642     } else {
0643       L1GctEtHadCollection::const_iterator hwEtHadItr = hwEtHadColl->begin();
0644       L1GctEtHadCollection::const_iterator hwEtHadEnd = hwEtHadColl->end();
0645 
0646       int iHad = 0;
0647       for (; hwEtHadItr != hwEtHadEnd; ++hwEtHadItr, ++iHad) {
0648         int bx = hwEtHadItr->bx();
0649 
0650         if (!centralBxOnly_ || bx == 0) {
0651           // HT bin low edge
0652           double htTot =
0653               (hwEtHadItr->overFlow() ? (double)L1GctEtHad::kEtHadMaxValue : (double)hwEtHadItr->et()) * htSumLSB +
0654               1.e-6;
0655 
0656           double htMiss = 0.;
0657           double phi = 0.;
0658           math::PtEtaPhiMLorentzVector p4;
0659           Ref<L1GctHtMissCollection> mhtRef;
0660 
0661           if (!ignoreHtMiss_) {
0662             L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
0663             L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
0664 
0665             int iMiss = 0;
0666             for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr, ++iMiss) {
0667               if (hwHtMissItr->bx() == bx) {
0668                 htMissMatched[iMiss] = true;
0669                 break;
0670               }
0671             }
0672 
0673             // If a L1GctHtMiss with the right bx is not found, itr == end.
0674             if (hwHtMissItr != hwHtMissEnd) {
0675               // HT bin low edge
0676               htMiss =
0677                   htMissScale->et(hwHtMissItr->overFlow() ? htMissScale->rankScaleMax() : hwHtMissItr->et()) + 1.e-6;
0678               // keep x and y components non-zero and
0679               // protect against roundoff.
0680 
0681               phi = caloGeom->htSumPhiBinCenter(hwHtMissItr->phi());
0682 
0683               p4 = math::PtEtaPhiMLorentzVector(htMiss, 0., phi, 0.);
0684 
0685               mhtRef = Ref<L1GctHtMissCollection>(hwHtMissColl, iMiss);
0686 
0687               //           cout << "HW HT Sums " << endl
0688               //            << "MHT: phi " << hwHtMissItr->phi() <<
0689               // " = " << phi
0690               //            << " ht " << hwHtMissItr->et() << " = "
0691               // << htMiss
0692               //            << " HtTot " << hwEtHadItr->et() << " =
0693               // "
0694               // << htTot
0695               //            << " bx " << bx
0696               //            << endl ;
0697             }
0698             //         else
0699             //       {
0700             //         cout << "HW HT Sums " << endl
0701             //          << "MHT: phi " << phi
0702             //          << " ht " << htMiss
0703             //          << " HtTot " << hwEtHadItr->et() << " = " <<
0704             // htTot
0705             //          << " bx " << bx
0706             //          << endl ;
0707             //       }
0708           }
0709 
0710           htMissColl->push_back(L1EtMissParticle(p4,
0711                                                  L1EtMissParticle::kMHT,
0712                                                  htTot,
0713                                                  Ref<L1GctEtMissCollection>(),
0714                                                  Ref<L1GctEtTotalCollection>(),
0715                                                  mhtRef,
0716                                                  Ref<L1GctEtHadCollection>(hwEtHadColl, iHad),
0717                                                  bx));
0718         }
0719       }
0720 
0721       if (!centralBxOnly_ && !ignoreHtMiss_) {
0722         // Make L1EtMissParticles for those L1GctHtMiss objects without
0723         // a matched L1GctHtTotal object.
0724         double htTot = 0.;
0725 
0726         L1GctHtMissCollection::const_iterator hwHtMissItr = hwHtMissColl->begin();
0727         L1GctHtMissCollection::const_iterator hwHtMissEnd = hwHtMissColl->end();
0728 
0729         int iMiss = 0;
0730         for (; hwHtMissItr != hwHtMissEnd; ++hwHtMissItr, ++iMiss) {
0731           if (!htMissMatched[iMiss]) {
0732             int bx = hwHtMissItr->bx();
0733 
0734             // HT bin low edge
0735             double htMiss =
0736                 htMissScale->et(hwHtMissItr->overFlow() ? htMissScale->rankScaleMax() : hwHtMissItr->et()) + 1.e-6;
0737             // keep x and y components non-zero and
0738             // protect against roundoff.
0739 
0740             double phi = caloGeom->htSumPhiBinCenter(hwHtMissItr->phi());
0741 
0742             math::PtEtaPhiMLorentzVector p4(htMiss, 0., phi, 0.);
0743 
0744             //        cout << "HW HT Sums " << endl
0745             //             << "MHT: phi " << hwHtMissItr->phi() << " = " <<
0746             // phi
0747             //             << " ht " << hwHtMissItr->et() << " = " << htMiss
0748             //             << " HtTot " << htTot
0749             //             << " bx " << bx
0750             //             << endl ;
0751 
0752             htMissColl->push_back(L1EtMissParticle(p4,
0753                                                    L1EtMissParticle::kMHT,
0754                                                    htTot,
0755                                                    Ref<L1GctEtMissCollection>(),
0756                                                    Ref<L1GctEtTotalCollection>(),
0757                                                    Ref<L1GctHtMissCollection>(hwHtMissColl, iMiss),
0758                                                    Ref<L1GctEtHadCollection>(),
0759                                                    bx));
0760           }
0761         }
0762       }
0763     }
0764 
0765     // ~~~~~~~~~~~~~~~~~~~~ HF Rings ~~~~~~~~~~~~~~~~~~~~
0766 
0767     Handle<L1GctHFRingEtSumsCollection> hwHFEtSumsColl;
0768     iEvent.getByLabel(hfRingEtSumsSource_, hwHFEtSumsColl);
0769 
0770     Handle<L1GctHFBitCountsCollection> hwHFBitCountsColl;
0771     iEvent.getByLabel(hfRingBitCountsSource_, hwHFBitCountsColl);
0772 
0773     ESHandle<L1CaloEtScale> hfRingEtScale = iSetup.getHandle(hfRingEtScaleToken_);
0774 
0775     if (!hwHFEtSumsColl.isValid()) {
0776       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHFRingEtSumsCollection with " << hfRingEtSumsSource_
0777                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0778     } else if (!hwHFBitCountsColl.isValid()) {
0779       LogDebug("L1ExtraParticlesProd") << "\nWarning: L1GctHFBitCountsCollection with " << hfRingBitCountsSource_
0780                                        << "\nrequested in configuration, but not found in the event." << std::endl;
0781     } else {
0782       L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsItr = hwHFEtSumsColl->begin();
0783       L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsEnd = hwHFEtSumsColl->end();
0784 
0785       int iEtSums = 0;
0786       for (; hwHFEtSumsItr != hwHFEtSumsEnd; ++hwHFEtSumsItr, ++iEtSums) {
0787         int bx = hwHFEtSumsItr->bx();
0788 
0789         if (!centralBxOnly_ || bx == 0) {
0790           L1GctHFBitCountsCollection::const_iterator hwHFBitCountsItr = hwHFBitCountsColl->begin();
0791           L1GctHFBitCountsCollection::const_iterator hwHFBitCountsEnd = hwHFBitCountsColl->end();
0792 
0793           int iBitCounts = 0;
0794           for (; hwHFBitCountsItr != hwHFBitCountsEnd; ++hwHFBitCountsItr, ++iBitCounts) {
0795             if (hwHFBitCountsItr->bx() == bx) {
0796               break;
0797             }
0798           }
0799 
0800           // If a L1GctHFBitCounts with the right bx is not found, itr == end.
0801           if (hwHFBitCountsItr != hwHFBitCountsEnd) {
0802             // Construct L1HFRings only if both HW objects are present.
0803 
0804             // cout << "HF Rings " << endl
0805 
0806             // HF Et sums
0807             double etSums[L1HFRings::kNumRings];
0808             etSums[L1HFRings::kRing1PosEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(0)) + 1.e-6;
0809             etSums[L1HFRings::kRing1NegEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(1)) + 1.e-6;
0810             etSums[L1HFRings::kRing2PosEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(2)) + 1.e-6;
0811             etSums[L1HFRings::kRing2NegEta] = hfRingEtScale->et(hwHFEtSumsItr->etSum(3)) + 1.e-6;
0812             // protect against roundoff.
0813 
0814             //         cout << "HF Et Sums "
0815             //          << hwHFEtSumsItr->etSum( 0 ) << " = "
0816             //          << etSums[ L1HFRings::kRing1PosEta ] << " "
0817             //          << hwHFEtSumsItr->etSum( 1 ) << " = "
0818             //          << etSums[ L1HFRings::kRing1NegEta ] << " "
0819             //          << hwHFEtSumsItr->etSum( 2 ) << " = "
0820             //          << etSums[ L1HFRings::kRing2PosEta ] << " "
0821             //          << hwHFEtSumsItr->etSum( 3 ) << " = "
0822             //          << etSums[ L1HFRings::kRing2NegEta ]
0823             //          << endl ;
0824 
0825             // HF bit counts
0826             int bitCounts[L1HFRings::kNumRings];
0827             bitCounts[L1HFRings::kRing1PosEta] = hwHFBitCountsItr->bitCount(0);
0828             bitCounts[L1HFRings::kRing1NegEta] = hwHFBitCountsItr->bitCount(1);
0829             bitCounts[L1HFRings::kRing2PosEta] = hwHFBitCountsItr->bitCount(2);
0830             bitCounts[L1HFRings::kRing2NegEta] = hwHFBitCountsItr->bitCount(3);
0831 
0832             //         cout << "HF bit counts "
0833             //          << hwHFBitCountsItr->bitCount( 0 ) << " "
0834             //          << hwHFBitCountsItr->bitCount( 1 ) << " "
0835             //          << hwHFBitCountsItr->bitCount( 2 ) << " "
0836             //          << hwHFBitCountsItr->bitCount( 3 ) << " "
0837             //          << endl ;
0838 
0839             hfRingsColl->push_back(L1HFRings(etSums,
0840                                              bitCounts,
0841                                              Ref<L1GctHFRingEtSumsCollection>(hwHFEtSumsColl, iEtSums),
0842                                              Ref<L1GctHFBitCountsCollection>(hwHFBitCountsColl, iBitCounts),
0843                                              bx));
0844           }
0845         }
0846       }
0847     }
0848   }
0849 
0850   OrphanHandle<L1EmParticleCollection> isoEmHandle = iEvent.put(std::move(isoEmColl), "Isolated");
0851 
0852   OrphanHandle<L1EmParticleCollection> nonIsoEmHandle = iEvent.put(std::move(nonIsoEmColl), "NonIsolated");
0853 
0854   OrphanHandle<L1JetParticleCollection> cenJetHandle = iEvent.put(std::move(cenJetColl), "Central");
0855 
0856   OrphanHandle<L1JetParticleCollection> forJetHandle = iEvent.put(std::move(forJetColl), "Forward");
0857 
0858   OrphanHandle<L1JetParticleCollection> tauJetHandle = iEvent.put(std::move(tauJetColl), "Tau");
0859 
0860   OrphanHandle<L1JetParticleCollection> IsoTauJetHandle = iEvent.put(std::move(isoTauJetColl), "IsoTau");
0861 
0862   OrphanHandle<L1EtMissParticleCollection> etMissCollHandle = iEvent.put(std::move(etMissColl), "MET");
0863 
0864   OrphanHandle<L1EtMissParticleCollection> htMissCollHandle = iEvent.put(std::move(htMissColl), "MHT");
0865 
0866   OrphanHandle<L1HFRingsCollection> hfRingsCollHandle = iEvent.put(std::move(hfRingsColl));
0867 }
0868 
0869 // math::XYZTLorentzVector
0870 math::PtEtaPhiMLorentzVector L1ExtraParticlesProd::gctLorentzVector(const double &et,
0871                                                                     const L1GctCand &cand,
0872                                                                     const L1CaloGeometry *geom,
0873                                                                     bool central) {
0874   // To keep x and y components non-zero.
0875   double etCorr = et + 1.e-6;  // protect against roundoff, not only for et=0
0876 
0877   double eta = geom->etaBinCenter(cand.etaIndex(), central);
0878 
0879   //    double tanThOver2 = exp( -eta ) ;
0880   //    double ez = etCorr * ( 1. - tanThOver2 * tanThOver2 ) / ( 2. *
0881   //    tanThOver2 ); double e  = etCorr * ( 1. + tanThOver2 * tanThOver2 ) /
0882   //    ( 2. * tanThOver2 );
0883 
0884   double phi = geom->emJetPhiBinCenter(cand.phiIndex());
0885 
0886   //    return math::XYZTLorentzVector( etCorr * cos( phi ),
0887   //                   etCorr * sin( phi ),
0888   //                   ez,
0889   //                   e ) ;
0890   return math::PtEtaPhiMLorentzVector(etCorr, eta, phi, 0.);
0891 }
0892 
0893 // define this as a plug-in
0894 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!