/data/refman/pasoursint/CMSSW_5_3_4/src/L1Trigger/L1ExtraFromDigis/src/L1ExtraParticlesProd.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    L1ExtraParticlesProd
// Class:      L1ExtraParticlesProd
// 
//
// Original Author:  Werner Sun
//         Created:  Mon Oct  2 22:45:32 EDT 2006
// $Id: L1ExtraParticlesProd.cc,v 1.32 2009/12/14 22:23:14 wmtan Exp $
//
//


// system include files
#include <memory>

// user include files
#include "L1Trigger/L1ExtraFromDigis/interface/L1ExtraParticlesProd.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "DataFormats/Common/interface/OrphanHandle.h"

#include "CondFormats/L1TObjects/interface/L1CaloGeometry.h"
#include "CondFormats/DataRecord/interface/L1CaloGeometryRecord.h"

#include "CondFormats/L1TObjects/interface/L1CaloEtScale.h"
#include "CondFormats/DataRecord/interface/L1EmEtScaleRcd.h"
#include "CondFormats/DataRecord/interface/L1JetEtScaleRcd.h"
#include "CondFormats/DataRecord/interface/L1HtMissScaleRcd.h"
#include "CondFormats/DataRecord/interface/L1HfRingEtScaleRcd.h"
#include "CondFormats/L1TObjects/interface/L1GctJetFinderParams.h"
#include "CondFormats/DataRecord/interface/L1GctJetFinderParamsRcd.h"

#include "CondFormats/L1TObjects/interface/L1MuTriggerScales.h"
#include "CondFormats/DataRecord/interface/L1MuTriggerScalesRcd.h"
#include "CondFormats/L1TObjects/interface/L1MuTriggerPtScale.h"
#include "CondFormats/DataRecord/interface/L1MuTriggerPtScaleRcd.h"
#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

// #include "FWCore/Utilities/interface/EDMException.h"

//
// class decleration
//


//
// constants, enums and typedefs
//

//
// static data member definitions
//

double L1ExtraParticlesProd::muonMassGeV_ = 0.105658369 ; // PDG06

//
// constructors and destructor
//
L1ExtraParticlesProd::L1ExtraParticlesProd(const edm::ParameterSet& iConfig)
   : produceMuonParticles_( iConfig.getParameter< bool >(
      "produceMuonParticles" ) ),
     muonSource_( iConfig.getParameter< edm::InputTag >(
        "muonSource" ) ),
     produceCaloParticles_( iConfig.getParameter< bool >(
        "produceCaloParticles" ) ),
     isoEmSource_( iConfig.getParameter< edm::InputTag >(
        "isolatedEmSource" ) ),
     nonIsoEmSource_( iConfig.getParameter< edm::InputTag >(
        "nonIsolatedEmSource" ) ),
     cenJetSource_( iConfig.getParameter< edm::InputTag >(
        "centralJetSource" ) ),
     forJetSource_( iConfig.getParameter< edm::InputTag >(
        "forwardJetSource" ) ),
     tauJetSource_( iConfig.getParameter< edm::InputTag >(
        "tauJetSource" ) ),
     etTotSource_( iConfig.getParameter< edm::InputTag >(
        "etTotalSource" ) ),
     etHadSource_( iConfig.getParameter< edm::InputTag >(
        "etHadSource" ) ),
     etMissSource_( iConfig.getParameter< edm::InputTag >(
        "etMissSource" ) ),
     htMissSource_( iConfig.getParameter< edm::InputTag >(
        "htMissSource" ) ),
     hfRingEtSumsSource_( iConfig.getParameter< edm::InputTag >(
        "hfRingEtSumsSource" ) ),
     hfRingBitCountsSource_( iConfig.getParameter< edm::InputTag >(
        "hfRingBitCountsSource" ) ),
     centralBxOnly_( iConfig.getParameter< bool >( "centralBxOnly" ) ),
     ignoreHtMiss_( iConfig.getParameter< bool >( "ignoreHtMiss" ) )
{
   using namespace l1extra ;

   //register your products
   produces< L1EmParticleCollection >( "Isolated" ) ;
   produces< L1EmParticleCollection >( "NonIsolated" ) ;
   produces< L1JetParticleCollection >( "Central" ) ;
   produces< L1JetParticleCollection >( "Forward" ) ;
   produces< L1JetParticleCollection >( "Tau" ) ;
   produces< L1MuonParticleCollection >() ;
   produces< L1EtMissParticleCollection >( "MET" ) ;
   produces< L1EtMissParticleCollection >( "MHT" ) ;
   produces< L1HFRingsCollection >() ;

   //now do what ever other initialization is needed
}


L1ExtraParticlesProd::~L1ExtraParticlesProd()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to produce the data  ------------
void
L1ExtraParticlesProd::produce( edm::Event& iEvent,
                               const edm::EventSetup& iSetup)
{
   using namespace edm ;
   using namespace l1extra ;
   using namespace std ;

   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   // ~~~~~~~~~~~~~~~~~~~~ Muons ~~~~~~~~~~~~~~~~~~~~
   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   auto_ptr< L1MuonParticleCollection > muColl(
     new L1MuonParticleCollection );

   if( produceMuonParticles_ )
   {
      ESHandle< L1MuTriggerScales > muScales ;
      iSetup.get< L1MuTriggerScalesRcd >().get( muScales ) ;

      ESHandle< L1MuTriggerPtScale > muPtScale ;
      iSetup.get< L1MuTriggerPtScaleRcd >().get( muPtScale ) ;

      Handle< L1MuGMTReadoutCollection > hwMuCollection ;
      iEvent.getByLabel( muonSource_, hwMuCollection ) ;

      vector< L1MuGMTExtendedCand > hwMuCands ;

      if( !hwMuCollection.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1MuGMTReadoutCollection with " << muonSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
          if( centralBxOnly_ )
            {
              // Get GMT candidates from central bunch crossing only
              hwMuCands = hwMuCollection->getRecord().getGMTCands() ;
            }
          else
            {
              // Get GMT candidates from all bunch crossings
              vector< L1MuGMTReadoutRecord > records = hwMuCollection->getRecords();
              vector< L1MuGMTReadoutRecord >::const_iterator rItr = records.begin();
              vector< L1MuGMTReadoutRecord >::const_iterator rEnd = records.end();

              for( ; rItr != rEnd ; ++rItr )
                {
                  vector< L1MuGMTExtendedCand > tmpCands = rItr->getGMTCands() ;

                  hwMuCands.insert( hwMuCands.end(),
                                    tmpCands.begin(),
                                    tmpCands.end() ) ;
                }
            }

//       cout << "HW muons" << endl ;
          vector< L1MuGMTExtendedCand >::const_iterator muItr = hwMuCands.begin() ;
          vector< L1MuGMTExtendedCand >::const_iterator muEnd = hwMuCands.end() ;
          for( int i = 0 ; muItr != muEnd ; ++muItr, ++i )
            {
//       cout << "#" << i
//            << " name " << muItr->name()
//            << " empty " << muItr->empty()
//            << " pt " << muItr->ptIndex()
//            << " eta " << muItr->etaIndex()
//            << " phi " << muItr->phiIndex()
//            << " iso " << muItr->isol()
//            << " mip " << muItr->mip()
//            << " bx " << muItr->bx()
//            << endl ;

              if( !muItr->empty() )
                {
                  // keep x and y components non-zero and protect against roundoff.
                  double pt =
                    muPtScale->getPtScale()->getLowEdge( muItr->ptIndex() ) + 1.e-6 ;

//          cout << "L1Extra pt " << pt << endl ;

                  double eta =
                    muScales->getGMTEtaScale()->getCenter( muItr->etaIndex() ) ;

                  double phi =
                    muScales->getPhiScale()->getLowEdge( muItr->phiIndex() ) ;

                  math::PtEtaPhiMLorentzVector p4( pt,
                                                   eta,
                                                   phi,
                                                   muonMassGeV_ ) ;

                  muColl->push_back(
                                    L1MuonParticle( muItr->charge(),
                                                    p4,
                                                    *muItr,
                                                    muItr->bx() )
                                    ) ;
                }
            }
        }
   }
   
   OrphanHandle< L1MuonParticleCollection > muHandle =
     iEvent.put( muColl );


   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   // ~~~~~~~~~~~~~~~~~~~~ Calorimeter ~~~~~~~~~~~~~~~~~~~~
   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   auto_ptr< L1EmParticleCollection > isoEmColl(
      new L1EmParticleCollection );

   auto_ptr< L1EmParticleCollection > nonIsoEmColl(
      new L1EmParticleCollection );

   auto_ptr< L1JetParticleCollection > cenJetColl(
      new L1JetParticleCollection );

   auto_ptr< L1JetParticleCollection > forJetColl(
      new L1JetParticleCollection );

   auto_ptr< L1JetParticleCollection > tauJetColl(
      new L1JetParticleCollection );

   auto_ptr< L1EtMissParticleCollection > etMissColl(
      new L1EtMissParticleCollection );

   auto_ptr< L1EtMissParticleCollection > htMissColl(
      new L1EtMissParticleCollection );

   auto_ptr< L1HFRingsCollection > hfRingsColl( new L1HFRingsCollection );

   if( produceCaloParticles_ )
   {
      // ~~~~~~~~~~~~~~~~~~~~ Geometry ~~~~~~~~~~~~~~~~~~~~

      ESHandle< L1CaloGeometry > caloGeomESH ;
      iSetup.get< L1CaloGeometryRecord >().get( caloGeomESH ) ;
      const L1CaloGeometry* caloGeom = &( *caloGeomESH ) ;

      // ~~~~~~~~~~~~~~~~~~~~ EM ~~~~~~~~~~~~~~~~~~~~

      ESHandle< L1CaloEtScale > emScale ;
      iSetup.get< L1EmEtScaleRcd >().get( emScale ) ;

      // Isolated EM
      Handle< L1GctEmCandCollection > hwIsoEmCands ;
      iEvent.getByLabel( isoEmSource_, hwIsoEmCands ) ;

      if( !hwIsoEmCands.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctEmCandCollection with " << isoEmSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
//       cout << "HW iso EM" << endl ;

          L1GctEmCandCollection::const_iterator emItr = hwIsoEmCands->begin() ;
          L1GctEmCandCollection::const_iterator emEnd = hwIsoEmCands->end() ;
          for( int i = 0 ; emItr != emEnd ; ++emItr, ++i )
            {
//       cout << "#" << i
//            << " name " << emItr->name()
//            << " empty " << emItr->empty()
//            << " rank " << emItr->rank()
//            << " eta " << emItr->etaIndex()
//            << " sign " << emItr->etaSign()
//            << " phi " << emItr->phiIndex()
//            << " iso " << emItr->isolated()
//            << " bx " << emItr->bx()
//            << endl ;

              if( !emItr->empty() &&
                  ( !centralBxOnly_ || emItr->bx() == 0 ) )
                {
                  double et = emScale->et( emItr->rank() ) ;

//          cout << "L1Extra et " << et << endl ;

                  isoEmColl->push_back(
                                       L1EmParticle( gctLorentzVector( et, *emItr, caloGeom, true ),
                                                     Ref< L1GctEmCandCollection >( hwIsoEmCands,
                                                                                   i ),
                                                     emItr->bx() ) ) ;
                }
            }
        }


      // Non-isolated EM
      Handle< L1GctEmCandCollection > hwNonIsoEmCands ;
      iEvent.getByLabel( nonIsoEmSource_, hwNonIsoEmCands ) ;

      if( !hwNonIsoEmCands.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctEmCandCollection with " << nonIsoEmSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
//       cout << "HW non-iso EM" << endl ;
          L1GctEmCandCollection::const_iterator emItr = hwNonIsoEmCands->begin() ;
          L1GctEmCandCollection::const_iterator emEnd = hwNonIsoEmCands->end() ;
          for( int i = 0 ; emItr != emEnd ; ++emItr, ++i )
            {
//       cout << "#" << i
//            << " name " << emItr->name()
//            << " empty " << emItr->empty()
//            << " rank " << emItr->rank()
//            << " eta " << emItr->etaIndex()
//            << " sign " << emItr->etaSign()
//            << " phi " << emItr->phiIndex()
//            << " iso " << emItr->isolated()
//            << " bx " << emItr->bx()
//            << endl ;

              if( !emItr->empty() &&
                  ( !centralBxOnly_ || emItr->bx() == 0 ) )
                {
                  double et = emScale->et( emItr->rank() ) ;

//          cout << "L1Extra et " << et << endl ;

                  nonIsoEmColl->push_back(
                                          L1EmParticle( gctLorentzVector( et, *emItr, caloGeom, true ),
                                                        Ref< L1GctEmCandCollection >( hwNonIsoEmCands,
                                                                                      i ),
                                                        emItr->bx() ) );
                }
            }
        }


      // ~~~~~~~~~~~~~~~~~~~~ Jets ~~~~~~~~~~~~~~~~~~~~

      ESHandle< L1CaloEtScale > jetScale ;
      iSetup.get< L1JetEtScaleRcd >().get( jetScale ) ;

      // Central jets.
      Handle< L1GctJetCandCollection > hwCenJetCands ;
      iEvent.getByLabel( cenJetSource_, hwCenJetCands ) ;

      if( !hwCenJetCands.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctJetCandCollection with " << cenJetSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
//       cout << "HW central jets" << endl ;
          L1GctJetCandCollection::const_iterator jetItr = hwCenJetCands->begin() ;
          L1GctJetCandCollection::const_iterator jetEnd = hwCenJetCands->end() ;
          for( int i = 0 ; jetItr != jetEnd ; ++jetItr, ++i )
            {
//       cout << "#" << i
//            << " name " << jetItr->name()
//            << " empty " << jetItr->empty()
//            << " rank " << jetItr->rank()
//            << " eta " << jetItr->etaIndex()
//            << " sign " << jetItr->etaSign()
//            << " phi " << jetItr->phiIndex()
//            << " cen " << jetItr->isCentral()
//            << " for " << jetItr->isForward()
//            << " tau " << jetItr->isTau()
//            << " bx " << jetItr->bx()
//            << endl ;

              if( !jetItr->empty() &&
                  ( !centralBxOnly_ || jetItr->bx() == 0 ) )
                {
                  double et = jetScale->et( jetItr->rank() ) ;

//          cout << "L1Extra et " << et << endl ;

                  cenJetColl->push_back(
                                        L1JetParticle( gctLorentzVector( et, *jetItr, caloGeom, true ),
                                                       Ref< L1GctJetCandCollection >( hwCenJetCands,
                                                                                      i ),
                                                       jetItr->bx() ) ) ;
                }
            }
        }


      // Forward jets.
      Handle< L1GctJetCandCollection > hwForJetCands ;
      iEvent.getByLabel( forJetSource_, hwForJetCands ) ;

      if( !hwForJetCands.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctEmCandCollection with " << forJetSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
//       cout << "HW forward jets" << endl ;
          L1GctJetCandCollection::const_iterator jetItr = hwForJetCands->begin() ;
          L1GctJetCandCollection::const_iterator jetEnd = hwForJetCands->end() ;
          for( int i = 0 ; jetItr != jetEnd ; ++jetItr, ++i )
            {
//       cout << "#" << i
//            << " name " << jetItr->name()
//            << " empty " << jetItr->empty()
//            << " rank " << jetItr->rank()
//            << " eta " << jetItr->etaIndex()
//            << " sign " << jetItr->etaSign()
//            << " phi " << jetItr->phiIndex()
//            << " cen " << jetItr->isCentral()
//            << " for " << jetItr->isForward()
//            << " tau " << jetItr->isTau()
//            << " bx " << jetItr->bx()
//            << endl ;

              if( !jetItr->empty() &&
                  ( !centralBxOnly_ || jetItr->bx() == 0 ) )
                {
                  double et = jetScale->et( jetItr->rank() ) ;

//          cout << "L1Extra et " << et << endl ;

                  forJetColl->push_back(
                                        L1JetParticle( gctLorentzVector( et, *jetItr, caloGeom, false ),
                                                       Ref< L1GctJetCandCollection >( hwForJetCands,
                                                                                      i ),
                                                       jetItr->bx() ) ) ;
                }
            }
        }


      // Tau jets.
//       cout << "HW tau jets" << endl ;
      Handle< L1GctJetCandCollection > hwTauJetCands ;
      iEvent.getByLabel( tauJetSource_, hwTauJetCands ) ;

      if( !hwTauJetCands.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctJetCandCollection with " << tauJetSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
          L1GctJetCandCollection::const_iterator jetItr = hwTauJetCands->begin() ;
          L1GctJetCandCollection::const_iterator jetEnd = hwTauJetCands->end() ;
          for( int i = 0 ; jetItr != jetEnd ; ++jetItr, ++i )
            {
//       cout << "#" << i
//            << " name " << jetItr->name()
//            << " empty " << jetItr->empty()
//            << " rank " << jetItr->rank()
//            << " eta " << jetItr->etaIndex()
//            << " sign " << jetItr->etaSign()
//            << " phi " << jetItr->phiIndex()
//            << " cen " << jetItr->isCentral()
//            << " for " << jetItr->isForward()
//            << " tau " << jetItr->isTau()
//            << " bx " << jetItr->bx()
//            << endl ;

              if( !jetItr->empty() &&
                  ( !centralBxOnly_ || jetItr->bx() == 0 ) )
                {
                  double et = jetScale->et( jetItr->rank() ) ;

//          cout << "L1Extra et " << et << endl ;

                  tauJetColl->push_back(
                                        L1JetParticle( gctLorentzVector( et, *jetItr, caloGeom, true ),
                                                       Ref< L1GctJetCandCollection >( hwTauJetCands,
                                                                                      i ),
                                                       jetItr->bx() ) ) ;
                }
            }
        }

      // ~~~~~~~~~~~~~~~~~~~~ ET Sums ~~~~~~~~~~~~~~~~~~~~

      double etSumLSB = jetScale->linearLsb() ;

      Handle< L1GctEtTotalCollection > hwEtTotColl ;
      iEvent.getByLabel( etTotSource_, hwEtTotColl ) ;

      Handle< L1GctEtMissCollection > hwEtMissColl ;
      iEvent.getByLabel( etMissSource_, hwEtMissColl ) ;

      if( !hwEtTotColl.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctEtTotalCollection with " << etTotSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else if( !hwEtMissColl.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctEtMissCollection with " << etMissSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
          // Make a L1EtMissParticle even if either L1GctEtTotal or L1GctEtMiss
          // is missing for a given bx.  Keep track of which L1GctEtMiss objects
          // have a corresponding L1GctEtTotal object.
          std::vector< bool > etMissMatched ;

          L1GctEtMissCollection::const_iterator hwEtMissItr =
            hwEtMissColl->begin() ;
          L1GctEtMissCollection::const_iterator hwEtMissEnd =
            hwEtMissColl->end() ;
          for( ; hwEtMissItr != hwEtMissEnd ; ++hwEtMissItr )
            {
              etMissMatched.push_back( false ) ;
            }
      
          // Collate energy sums by bx
          L1GctEtTotalCollection::const_iterator hwEtTotItr =
            hwEtTotColl->begin() ;
          L1GctEtTotalCollection::const_iterator hwEtTotEnd =
            hwEtTotColl->end() ;

          int iTot = 0 ;
          for( ; hwEtTotItr != hwEtTotEnd ; ++hwEtTotItr, ++iTot )
            {
              int bx = hwEtTotItr->bx() ;

              if( !centralBxOnly_ || bx == 0 )
                {
                  // ET bin low edge
                  double etTot =
                    ( hwEtTotItr->overFlow() ?
                      ( double ) L1GctEtTotal::kEtTotalMaxValue :
                      ( double ) hwEtTotItr->et() ) * etSumLSB + 1.e-6 ;

                  int iMiss = 0 ;
                  hwEtMissItr = hwEtMissColl->begin() ;
                  hwEtMissEnd = hwEtMissColl->end() ;
                  for( ; hwEtMissItr != hwEtMissEnd ; ++hwEtMissItr, ++iMiss )
                    {
                      if( hwEtMissItr->bx() == bx )
                        {
                          etMissMatched[ iMiss ] = true ;
                          break ;
                        }
                    }

                  double etMiss = 0. ;
                  double phi = 0. ;
                  math::PtEtaPhiMLorentzVector p4 ;
                  Ref< L1GctEtMissCollection > metRef ;

                  // If a L1GctEtMiss with the right bx is not found, itr == end.
                  if( hwEtMissItr != hwEtMissEnd )
                    {
                      // ET bin low edge
                      etMiss =
                        ( hwEtMissItr->overFlow() ?
                          ( double ) L1GctEtMiss::kEtMissMaxValue :
                          ( double ) hwEtMissItr->et() ) * etSumLSB + 1.e-6 ;
                      // keep x and y components non-zero and
                      // protect against roundoff.

                      phi = caloGeom->etSumPhiBinCenter( hwEtMissItr->phi() ) ;

                      p4 = math::PtEtaPhiMLorentzVector( etMiss,
                                                         0.,
                                                         phi,
                                                         0. ) ;

                      metRef = Ref< L1GctEtMissCollection >( hwEtMissColl, iMiss ) ;

//             cout << "HW ET Sums " << endl
//                  << "MET: phi " << hwEtMissItr->phi() << " = " << phi
//                  << " et " << hwEtMissItr->et() << " = " << etMiss
//                  << " EtTot " << hwEtTotItr->et() << " = " << etTot
//                  << " bx " << bx
//                  << endl ;
                    }
//         else
//           {
//             cout << "HW ET Sums " << endl
//                  << "MET: phi " << phi
//                  << " et "<< etMiss
//                  << " EtTot " << hwEtTotItr->et() << " = " << etTot
//                  << " bx " << bx
//                  << endl ;
//           }

                  etMissColl->push_back(
                                        L1EtMissParticle(
                                                         p4,
                                                         L1EtMissParticle::kMET,
                                                         etTot,
                                                         metRef,
                                                         Ref< L1GctEtTotalCollection >( hwEtTotColl, iTot ),
                                                         Ref< L1GctHtMissCollection >(),
                                                         Ref< L1GctEtHadCollection >(),
                                                         bx
                                                         ) ) ;
                }
            }

          if( !centralBxOnly_ )
            {
              // Make L1EtMissParticles for those L1GctEtMiss objects without
              // a matched L1GctEtTotal object.

              double etTot = 0. ;

              hwEtMissItr = hwEtMissColl->begin() ;
              hwEtMissEnd = hwEtMissColl->end() ;
              int iMiss = 0 ;
              for( ; hwEtMissItr != hwEtMissEnd ; ++hwEtMissItr, ++iMiss )
                {
                  if( !etMissMatched[ iMiss ] )
                    {
                      int bx = hwEtMissItr->bx() ;

                      // ET bin low edge
                      double etMiss =
                        ( hwEtMissItr->overFlow() ?
                          ( double ) L1GctEtMiss::kEtMissMaxValue :
                          ( double ) hwEtMissItr->et() ) * etSumLSB + 1.e-6 ;
                      // keep x and y components non-zero and
                      // protect against roundoff.

                      double phi =
                        caloGeom->etSumPhiBinCenter( hwEtMissItr->phi() ) ;

                      math::PtEtaPhiMLorentzVector p4( etMiss,
                                                       0.,
                                                       phi,
                                                       0. ) ;

//                cout << "HW ET Sums " << endl
//                     << "MET: phi " << hwEtMissItr->phi() << " = " << phi
//                     << " et " << hwEtMissItr->et() << " = " << etMiss
//                     << " EtTot " << etTot
//                     << " bx " << bx
//                     << endl ;

                      etMissColl->push_back(
                                            L1EtMissParticle(
                                                             p4,
                                                             L1EtMissParticle::kMET,
                                                             etTot,
                                                             Ref< L1GctEtMissCollection >( hwEtMissColl, iMiss ),
                                                             Ref< L1GctEtTotalCollection >(),
                                                             Ref< L1GctHtMissCollection >(),
                                                             Ref< L1GctEtHadCollection >(),
                                                             bx
                                                             ) ) ;
                    }
                }
            }
        }

      // ~~~~~~~~~~~~~~~~~~~~ HT Sums ~~~~~~~~~~~~~~~~~~~~

      Handle< L1GctEtHadCollection > hwEtHadColl ;
      iEvent.getByLabel( etHadSource_, hwEtHadColl ) ;

      Handle< L1GctHtMissCollection > hwHtMissColl ;
      if( !ignoreHtMiss_ )
        {
          iEvent.getByLabel( htMissSource_, hwHtMissColl ) ;
        }

      ESHandle< L1GctJetFinderParams > jetFinderParams ;
      iSetup.get< L1GctJetFinderParamsRcd >().get( jetFinderParams ) ;
      double htSumLSB = jetFinderParams->getHtLsbGeV();

      ESHandle< L1CaloEtScale > htMissScale ;
      std::vector< bool > htMissMatched ;
      if( !ignoreHtMiss_ )
        {
          iSetup.get< L1HtMissScaleRcd >().get( htMissScale ) ;

          if( !hwEtHadColl.isValid() )
            {
              LogDebug("L1ExtraParticlesProd")
                << "\nWarning: L1GctEtHadCollection with " << etHadSource_
                << "\nrequested in configuration, but not found in the event."
                << std::endl;
            }
          else if( !hwHtMissColl.isValid() )
            {
              LogDebug("L1ExtraParticlesProd")
                << "\nWarning: L1GctHtMissCollection with " << htMissSource_
                << "\nrequested in configuration, but not found in the event."
                << std::endl;
            }
          else
            {
              // Make a L1EtMissParticle even if either L1GctEtHad or L1GctHtMiss
              // is missing for a given bx. Keep track of which L1GctHtMiss objects
              // have a corresponding L1GctHtTotal object.
              L1GctHtMissCollection::const_iterator hwHtMissItr =
                hwHtMissColl->begin() ;
              L1GctHtMissCollection::const_iterator hwHtMissEnd =
                hwHtMissColl->end() ;
              for( ; hwHtMissItr != hwHtMissEnd ; ++hwHtMissItr )
                {
                  htMissMatched.push_back( false ) ;
                }
            }
        }

      if( !hwEtHadColl.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctEtHadCollection with " << etHadSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else if( !hwHtMissColl.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctHtMissCollection with " << htMissSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
          L1GctEtHadCollection::const_iterator hwEtHadItr =
            hwEtHadColl->begin() ;
          L1GctEtHadCollection::const_iterator hwEtHadEnd =
            hwEtHadColl->end() ;

          int iHad = 0 ;
          for( ; hwEtHadItr != hwEtHadEnd ; ++hwEtHadItr, ++iHad )
            {
              int bx = hwEtHadItr->bx() ;

              if( !centralBxOnly_ || bx == 0 )
                {
                  // HT bin low edge
                  double htTot =
                    ( hwEtHadItr->overFlow() ?
                      ( double ) L1GctEtHad::kEtHadMaxValue :
                      ( double ) hwEtHadItr->et() ) * htSumLSB + 1.e-6 ;

                  double htMiss = 0. ;
                  double phi = 0. ;
                  math::PtEtaPhiMLorentzVector p4 ;
                  Ref< L1GctHtMissCollection > mhtRef ;

                  if( !ignoreHtMiss_ )
                    {
                      L1GctHtMissCollection::const_iterator hwHtMissItr =
                        hwHtMissColl->begin() ;
                      L1GctHtMissCollection::const_iterator hwHtMissEnd =
                        hwHtMissColl->end() ;

                      int iMiss = 0 ;
                      for( ; hwHtMissItr != hwHtMissEnd ; ++hwHtMissItr, ++iMiss )
                        {
                          if( hwHtMissItr->bx() == bx )
                            {
                              htMissMatched[ iMiss ] = true ;
                              break ;
                            }
                        }

                      // If a L1GctHtMiss with the right bx is not found, itr == end.
                      if( hwHtMissItr != hwHtMissEnd )
                        {
                          // HT bin low edge
                          htMiss =
                            htMissScale->et( hwHtMissItr->overFlow() ?
                                             htMissScale->rankScaleMax() :
                                             hwHtMissItr->et() ) + 1.e-6 ;
                          // keep x and y components non-zero and
                          // protect against roundoff.

                          phi =
                            caloGeom->htSumPhiBinCenter( hwHtMissItr->phi() ) ;

                          p4 = math::PtEtaPhiMLorentzVector( htMiss,
                                                             0.,
                                                             phi,
                                                             0. ) ;

                          mhtRef=Ref< L1GctHtMissCollection >( hwHtMissColl, iMiss );

//                 cout << "HW HT Sums " << endl
//                      << "MHT: phi " << hwHtMissItr->phi() << " = " << phi
//                      << " ht " << hwHtMissItr->et() << " = " << htMiss
//                      << " HtTot " << hwEtHadItr->et() << " = " << htTot
//                      << " bx " << bx
//                      << endl ;
                        }
//             else
//               {
//                 cout << "HW HT Sums " << endl
//                      << "MHT: phi " << phi
//                      << " ht " << htMiss
//                      << " HtTot " << hwEtHadItr->et() << " = " << htTot
//                      << " bx " << bx
//                      << endl ;
//               }
                    }

                  htMissColl->push_back(
                                        L1EtMissParticle(
                                                         p4,
                                                         L1EtMissParticle::kMHT,
                                                         htTot,
                                                         Ref< L1GctEtMissCollection >(),
                                                         Ref< L1GctEtTotalCollection >(),
                                                         mhtRef,
                                                         Ref< L1GctEtHadCollection >( hwEtHadColl, iHad ),
                                                         bx
                                                         ) ) ;
                }
            }

          if( !centralBxOnly_  && !ignoreHtMiss_ )
            {
              // Make L1EtMissParticles for those L1GctHtMiss objects without
              // a matched L1GctHtTotal object.
              double htTot = 0. ;

              L1GctHtMissCollection::const_iterator hwHtMissItr =
                hwHtMissColl->begin() ;
              L1GctHtMissCollection::const_iterator hwHtMissEnd =
                hwHtMissColl->end() ;

              int iMiss = 0 ;
              for( ; hwHtMissItr != hwHtMissEnd ; ++hwHtMissItr, ++iMiss )
                {
                  if( !htMissMatched[ iMiss ] )
                    {
                      int bx = hwHtMissItr->bx() ;

                      // HT bin low edge
                      double htMiss =
                        htMissScale->et( hwHtMissItr->overFlow() ?
                                         htMissScale->rankScaleMax() :
                                         hwHtMissItr->et() ) + 1.e-6 ;
                      // keep x and y components non-zero and
                      // protect against roundoff.

                      double phi =
                        caloGeom->htSumPhiBinCenter( hwHtMissItr->phi() ) ;

                      math::PtEtaPhiMLorentzVector p4( htMiss,
                                                       0.,
                                                       phi,
                                                       0. ) ;

//                cout << "HW HT Sums " << endl
//                     << "MHT: phi " << hwHtMissItr->phi() << " = " << phi
//                     << " ht " << hwHtMissItr->et() << " = " << htMiss
//                     << " HtTot " << htTot
//                     << " bx " << bx
//                     << endl ;

                      htMissColl->push_back(
                                            L1EtMissParticle(
                                                             p4,
                                                             L1EtMissParticle::kMHT,
                                                             htTot,
                                                             Ref< L1GctEtMissCollection >(),
                                                             Ref< L1GctEtTotalCollection >(),
                                                             Ref< L1GctHtMissCollection >( hwHtMissColl, iMiss ),
                                                             Ref< L1GctEtHadCollection >(),
                                                             bx
                                                             ) ) ;
                    }
                }
            }
        }

      // ~~~~~~~~~~~~~~~~~~~~ HF Rings ~~~~~~~~~~~~~~~~~~~~

      Handle< L1GctHFRingEtSumsCollection > hwHFEtSumsColl ;
      iEvent.getByLabel( hfRingEtSumsSource_, hwHFEtSumsColl ) ;

      Handle< L1GctHFBitCountsCollection > hwHFBitCountsColl ;
      iEvent.getByLabel( hfRingBitCountsSource_, hwHFBitCountsColl ) ;

      ESHandle< L1CaloEtScale > hfRingEtScale ;
      iSetup.get< L1HfRingEtScaleRcd >().get( hfRingEtScale ) ;

      if( !hwHFEtSumsColl.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctHFRingEtSumsCollection with " << hfRingEtSumsSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else if( !hwHFBitCountsColl.isValid() )
        {
          LogDebug("L1ExtraParticlesProd")
            << "\nWarning: L1GctHFBitCountsCollection with " << hfRingBitCountsSource_
            << "\nrequested in configuration, but not found in the event."
            << std::endl;
        }
      else
        {
          L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsItr =
            hwHFEtSumsColl->begin() ;
          L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsEnd =
            hwHFEtSumsColl->end() ;

          int iEtSums = 0 ;
          for( ; hwHFEtSumsItr != hwHFEtSumsEnd ; ++hwHFEtSumsItr, ++iEtSums )
            {
              int bx = hwHFEtSumsItr->bx() ;

              if( !centralBxOnly_ || bx == 0 )
                {
                  L1GctHFBitCountsCollection::const_iterator hwHFBitCountsItr =
                    hwHFBitCountsColl->begin() ;
                  L1GctHFBitCountsCollection::const_iterator hwHFBitCountsEnd =
                    hwHFBitCountsColl->end() ;

                  int iBitCounts = 0 ;
                  for( ; hwHFBitCountsItr != hwHFBitCountsEnd ;
                       ++hwHFBitCountsItr, ++iBitCounts )
                    {
                      if( hwHFBitCountsItr->bx() == bx )
                        {
                          break ;
                        }
                    }

                  // If a L1GctHFBitCounts with the right bx is not found, itr == end.
                  if( hwHFBitCountsItr != hwHFBitCountsEnd )
                    {
                      // Construct L1HFRings only if both HW objects are present.

                      // cout << "HF Rings " << endl

                      // HF Et sums
                      double etSums[ L1HFRings::kNumRings ] ;
                      etSums[ L1HFRings::kRing1PosEta ] =
                        hfRingEtScale->et( hwHFEtSumsItr->etSum( 0 ) ) + 1.e-6 ;
                      etSums[ L1HFRings::kRing1NegEta ] =
                        hfRingEtScale->et( hwHFEtSumsItr->etSum( 1 ) ) + 1.e-6 ;
                      etSums[ L1HFRings::kRing2PosEta ] =
                        hfRingEtScale->et( hwHFEtSumsItr->etSum( 2 ) ) + 1.e-6 ;
                      etSums[ L1HFRings::kRing2NegEta ] =
                        hfRingEtScale->et( hwHFEtSumsItr->etSum( 3 ) ) + 1.e-6 ;
                      // protect against roundoff.

//             cout << "HF Et Sums "
//                  << hwHFEtSumsItr->etSum( 0 ) << " = "
//                  << etSums[ L1HFRings::kRing1PosEta ] << " "
//                  << hwHFEtSumsItr->etSum( 1 ) << " = "
//                  << etSums[ L1HFRings::kRing1NegEta ] << " "
//                  << hwHFEtSumsItr->etSum( 2 ) << " = "
//                  << etSums[ L1HFRings::kRing2PosEta ] << " "
//                  << hwHFEtSumsItr->etSum( 3 ) << " = "
//                  << etSums[ L1HFRings::kRing2NegEta ]
//                  << endl ;

                      // HF bit counts
                      int bitCounts[ L1HFRings::kNumRings ] ;
                      bitCounts[ L1HFRings::kRing1PosEta ] =
                        hwHFBitCountsItr->bitCount( 0 ) ;
                      bitCounts[ L1HFRings::kRing1NegEta ] =
                        hwHFBitCountsItr->bitCount( 1 ) ;
                      bitCounts[ L1HFRings::kRing2PosEta ] =
                        hwHFBitCountsItr->bitCount( 2 ) ;
                      bitCounts[ L1HFRings::kRing2NegEta ] =
                        hwHFBitCountsItr->bitCount( 3 ) ;

//             cout << "HF bit counts "
//                  << hwHFBitCountsItr->bitCount( 0 ) << " "
//                  << hwHFBitCountsItr->bitCount( 1 ) << " "
//                  << hwHFBitCountsItr->bitCount( 2 ) << " "
//                  << hwHFBitCountsItr->bitCount( 3 ) << " "
//                  << endl ;

                      hfRingsColl->push_back(
                                             L1HFRings( etSums,
                                                        bitCounts,
                                                        Ref< L1GctHFRingEtSumsCollection >( hwHFEtSumsColl,
                                                                                            iEtSums ),
                                                        Ref< L1GctHFBitCountsCollection >( hwHFBitCountsColl,
                                                                                           iBitCounts ),
                                                        bx
                                                        ) ) ;
                    }
                }
            }
        }
   }

   OrphanHandle< L1EmParticleCollection > isoEmHandle =
     iEvent.put( isoEmColl, "Isolated" ) ;

   OrphanHandle< L1EmParticleCollection > nonIsoEmHandle =
     iEvent.put( nonIsoEmColl, "NonIsolated" ) ;

   OrphanHandle< L1JetParticleCollection > cenJetHandle =
     iEvent.put( cenJetColl, "Central" ) ;

   OrphanHandle< L1JetParticleCollection > forJetHandle =
     iEvent.put( forJetColl, "Forward" ) ;

   OrphanHandle< L1JetParticleCollection > tauJetHandle =
     iEvent.put( tauJetColl, "Tau" ) ;

   OrphanHandle< L1EtMissParticleCollection > etMissCollHandle =
     iEvent.put( etMissColl, "MET" ) ;

   OrphanHandle< L1EtMissParticleCollection > htMissCollHandle =
     iEvent.put( htMissColl, "MHT" ) ;

   OrphanHandle< L1HFRingsCollection > hfRingsCollHandle =
     iEvent.put( hfRingsColl ) ;
}

//math::XYZTLorentzVector
math::PtEtaPhiMLorentzVector
L1ExtraParticlesProd::gctLorentzVector( const double& et,
                                        const L1GctCand& cand,
                                        const L1CaloGeometry* geom,
                                        bool central )
{
   // To keep x and y components non-zero.
   double etCorr = et + 1.e-6 ; // protect against roundoff, not only for et=0

   double eta = geom->etaBinCenter( cand.etaIndex(), central ) ;

//    double tanThOver2 = exp( -eta ) ;
//    double ez = etCorr * ( 1. - tanThOver2 * tanThOver2 ) / ( 2. * tanThOver2 );
//    double e  = etCorr * ( 1. + tanThOver2 * tanThOver2 ) / ( 2. * tanThOver2 );

   double phi = geom->emJetPhiBinCenter( cand.phiIndex() ) ;

//    return math::XYZTLorentzVector( etCorr * cos( phi ),
//                                 etCorr * sin( phi ),
//                                 ez,
//                                 e ) ;
   return math::PtEtaPhiMLorentzVector( etCorr,
                                        eta,
                                        phi,
                                        0. ) ;
}     

// ------------ method called once each job just before starting event loop  ------------
void 
L1ExtraParticlesProd::beginJob()
{
}

// ------------ method called once each job just after ending the event loop  ------------
void 
L1ExtraParticlesProd::endJob() {
}

//define this as a plug-in
DEFINE_FWK_MODULE(L1ExtraParticlesProd);