Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    L1ExtraParticlesProd
00004 // Class:      L1ExtraParticlesProd
00005 // 
00008 //
00009 // Original Author:  Werner Sun
00010 //         Created:  Mon Oct  2 22:45:32 EDT 2006
00011 // $Id: L1ExtraParticlesProd.cc,v 1.32 2009/12/14 22:23:14 wmtan Exp $
00012 //
00013 //
00014 
00015 
00016 // system include files
00017 #include <memory>
00018 
00019 // user include files
00020 #include "L1Trigger/L1ExtraFromDigis/interface/L1ExtraParticlesProd.h"
00021 
00022 #include "FWCore/Framework/interface/Frameworkfwd.h"
00023 #include "FWCore/Framework/interface/MakerMacros.h"
00024 
00025 #include "DataFormats/Common/interface/Handle.h"
00026 #include "FWCore/Framework/interface/ESHandle.h"
00027 #include "DataFormats/Common/interface/OrphanHandle.h"
00028 
00029 #include "CondFormats/L1TObjects/interface/L1CaloGeometry.h"
00030 #include "CondFormats/DataRecord/interface/L1CaloGeometryRecord.h"
00031 
00032 #include "CondFormats/L1TObjects/interface/L1CaloEtScale.h"
00033 #include "CondFormats/DataRecord/interface/L1EmEtScaleRcd.h"
00034 #include "CondFormats/DataRecord/interface/L1JetEtScaleRcd.h"
00035 #include "CondFormats/DataRecord/interface/L1HtMissScaleRcd.h"
00036 #include "CondFormats/DataRecord/interface/L1HfRingEtScaleRcd.h"
00037 #include "CondFormats/L1TObjects/interface/L1GctJetFinderParams.h"
00038 #include "CondFormats/DataRecord/interface/L1GctJetFinderParamsRcd.h"
00039 
00040 #include "CondFormats/L1TObjects/interface/L1MuTriggerScales.h"
00041 #include "CondFormats/DataRecord/interface/L1MuTriggerScalesRcd.h"
00042 #include "CondFormats/L1TObjects/interface/L1MuTriggerPtScale.h"
00043 #include "CondFormats/DataRecord/interface/L1MuTriggerPtScaleRcd.h"
00044 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"
00045 
00046 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00047 
00048 // #include "FWCore/Utilities/interface/EDMException.h"
00049 
00050 //
00051 // class decleration
00052 //
00053 
00054 
00055 //
00056 // constants, enums and typedefs
00057 //
00058 
00059 //
00060 // static data member definitions
00061 //
00062 
00063 double L1ExtraParticlesProd::muonMassGeV_ = 0.105658369 ; // PDG06
00064 
00065 //
00066 // constructors and destructor
00067 //
00068 L1ExtraParticlesProd::L1ExtraParticlesProd(const edm::ParameterSet& iConfig)
00069    : produceMuonParticles_( iConfig.getParameter< bool >(
00070       "produceMuonParticles" ) ),
00071      muonSource_( iConfig.getParameter< edm::InputTag >(
00072         "muonSource" ) ),
00073      produceCaloParticles_( iConfig.getParameter< bool >(
00074         "produceCaloParticles" ) ),
00075      isoEmSource_( iConfig.getParameter< edm::InputTag >(
00076         "isolatedEmSource" ) ),
00077      nonIsoEmSource_( iConfig.getParameter< edm::InputTag >(
00078         "nonIsolatedEmSource" ) ),
00079      cenJetSource_( iConfig.getParameter< edm::InputTag >(
00080         "centralJetSource" ) ),
00081      forJetSource_( iConfig.getParameter< edm::InputTag >(
00082         "forwardJetSource" ) ),
00083      tauJetSource_( iConfig.getParameter< edm::InputTag >(
00084         "tauJetSource" ) ),
00085      etTotSource_( iConfig.getParameter< edm::InputTag >(
00086         "etTotalSource" ) ),
00087      etHadSource_( iConfig.getParameter< edm::InputTag >(
00088         "etHadSource" ) ),
00089      etMissSource_( iConfig.getParameter< edm::InputTag >(
00090         "etMissSource" ) ),
00091      htMissSource_( iConfig.getParameter< edm::InputTag >(
00092         "htMissSource" ) ),
00093      hfRingEtSumsSource_( iConfig.getParameter< edm::InputTag >(
00094         "hfRingEtSumsSource" ) ),
00095      hfRingBitCountsSource_( iConfig.getParameter< edm::InputTag >(
00096         "hfRingBitCountsSource" ) ),
00097      centralBxOnly_( iConfig.getParameter< bool >( "centralBxOnly" ) ),
00098      ignoreHtMiss_( iConfig.getParameter< bool >( "ignoreHtMiss" ) )
00099 {
00100    using namespace l1extra ;
00101 
00102    //register your products
00103    produces< L1EmParticleCollection >( "Isolated" ) ;
00104    produces< L1EmParticleCollection >( "NonIsolated" ) ;
00105    produces< L1JetParticleCollection >( "Central" ) ;
00106    produces< L1JetParticleCollection >( "Forward" ) ;
00107    produces< L1JetParticleCollection >( "Tau" ) ;
00108    produces< L1MuonParticleCollection >() ;
00109    produces< L1EtMissParticleCollection >( "MET" ) ;
00110    produces< L1EtMissParticleCollection >( "MHT" ) ;
00111    produces< L1HFRingsCollection >() ;
00112 
00113    //now do what ever other initialization is needed
00114 }
00115 
00116 
00117 L1ExtraParticlesProd::~L1ExtraParticlesProd()
00118 {
00119  
00120    // do anything here that needs to be done at desctruction time
00121    // (e.g. close files, deallocate resources etc.)
00122 
00123 }
00124 
00125 
00126 //
00127 // member functions
00128 //
00129 
00130 // ------------ method called to produce the data  ------------
00131 void
00132 L1ExtraParticlesProd::produce( edm::Event& iEvent,
00133                                const edm::EventSetup& iSetup)
00134 {
00135    using namespace edm ;
00136    using namespace l1extra ;
00137    using namespace std ;
00138 
00139    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00140    // ~~~~~~~~~~~~~~~~~~~~ Muons ~~~~~~~~~~~~~~~~~~~~
00141    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00142 
00143    auto_ptr< L1MuonParticleCollection > muColl(
00144      new L1MuonParticleCollection );
00145 
00146    if( produceMuonParticles_ )
00147    {
00148       ESHandle< L1MuTriggerScales > muScales ;
00149       iSetup.get< L1MuTriggerScalesRcd >().get( muScales ) ;
00150 
00151       ESHandle< L1MuTriggerPtScale > muPtScale ;
00152       iSetup.get< L1MuTriggerPtScaleRcd >().get( muPtScale ) ;
00153 
00154       Handle< L1MuGMTReadoutCollection > hwMuCollection ;
00155       iEvent.getByLabel( muonSource_, hwMuCollection ) ;
00156 
00157       vector< L1MuGMTExtendedCand > hwMuCands ;
00158 
00159       if( !hwMuCollection.isValid() )
00160         {
00161           LogDebug("L1ExtraParticlesProd")
00162             << "\nWarning: L1MuGMTReadoutCollection with " << muonSource_
00163             << "\nrequested in configuration, but not found in the event."
00164             << std::endl;
00165         }
00166       else
00167         {
00168           if( centralBxOnly_ )
00169             {
00170               // Get GMT candidates from central bunch crossing only
00171               hwMuCands = hwMuCollection->getRecord().getGMTCands() ;
00172             }
00173           else
00174             {
00175               // Get GMT candidates from all bunch crossings
00176               vector< L1MuGMTReadoutRecord > records = hwMuCollection->getRecords();
00177               vector< L1MuGMTReadoutRecord >::const_iterator rItr = records.begin();
00178               vector< L1MuGMTReadoutRecord >::const_iterator rEnd = records.end();
00179 
00180               for( ; rItr != rEnd ; ++rItr )
00181                 {
00182                   vector< L1MuGMTExtendedCand > tmpCands = rItr->getGMTCands() ;
00183 
00184                   hwMuCands.insert( hwMuCands.end(),
00185                                     tmpCands.begin(),
00186                                     tmpCands.end() ) ;
00187                 }
00188             }
00189 
00190 //       cout << "HW muons" << endl ;
00191           vector< L1MuGMTExtendedCand >::const_iterator muItr = hwMuCands.begin() ;
00192           vector< L1MuGMTExtendedCand >::const_iterator muEnd = hwMuCands.end() ;
00193           for( int i = 0 ; muItr != muEnd ; ++muItr, ++i )
00194             {
00195 //       cout << "#" << i
00196 //            << " name " << muItr->name()
00197 //            << " empty " << muItr->empty()
00198 //            << " pt " << muItr->ptIndex()
00199 //            << " eta " << muItr->etaIndex()
00200 //            << " phi " << muItr->phiIndex()
00201 //            << " iso " << muItr->isol()
00202 //            << " mip " << muItr->mip()
00203 //            << " bx " << muItr->bx()
00204 //            << endl ;
00205 
00206               if( !muItr->empty() )
00207                 {
00208                   // keep x and y components non-zero and protect against roundoff.
00209                   double pt =
00210                     muPtScale->getPtScale()->getLowEdge( muItr->ptIndex() ) + 1.e-6 ;
00211 
00212 //          cout << "L1Extra pt " << pt << endl ;
00213 
00214                   double eta =
00215                     muScales->getGMTEtaScale()->getCenter( muItr->etaIndex() ) ;
00216 
00217                   double phi =
00218                     muScales->getPhiScale()->getLowEdge( muItr->phiIndex() ) ;
00219 
00220                   math::PtEtaPhiMLorentzVector p4( pt,
00221                                                    eta,
00222                                                    phi,
00223                                                    muonMassGeV_ ) ;
00224 
00225                   muColl->push_back(
00226                                     L1MuonParticle( muItr->charge(),
00227                                                     p4,
00228                                                     *muItr,
00229                                                     muItr->bx() )
00230                                     ) ;
00231                 }
00232             }
00233         }
00234    }
00235    
00236    OrphanHandle< L1MuonParticleCollection > muHandle =
00237      iEvent.put( muColl );
00238 
00239 
00240    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00241    // ~~~~~~~~~~~~~~~~~~~~ Calorimeter ~~~~~~~~~~~~~~~~~~~~
00242    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00243 
00244    auto_ptr< L1EmParticleCollection > isoEmColl(
00245       new L1EmParticleCollection );
00246 
00247    auto_ptr< L1EmParticleCollection > nonIsoEmColl(
00248       new L1EmParticleCollection );
00249 
00250    auto_ptr< L1JetParticleCollection > cenJetColl(
00251       new L1JetParticleCollection );
00252 
00253    auto_ptr< L1JetParticleCollection > forJetColl(
00254       new L1JetParticleCollection );
00255 
00256    auto_ptr< L1JetParticleCollection > tauJetColl(
00257       new L1JetParticleCollection );
00258 
00259    auto_ptr< L1EtMissParticleCollection > etMissColl(
00260       new L1EtMissParticleCollection );
00261 
00262    auto_ptr< L1EtMissParticleCollection > htMissColl(
00263       new L1EtMissParticleCollection );
00264 
00265    auto_ptr< L1HFRingsCollection > hfRingsColl( new L1HFRingsCollection );
00266 
00267    if( produceCaloParticles_ )
00268    {
00269       // ~~~~~~~~~~~~~~~~~~~~ Geometry ~~~~~~~~~~~~~~~~~~~~
00270 
00271       ESHandle< L1CaloGeometry > caloGeomESH ;
00272       iSetup.get< L1CaloGeometryRecord >().get( caloGeomESH ) ;
00273       const L1CaloGeometry* caloGeom = &( *caloGeomESH ) ;
00274 
00275       // ~~~~~~~~~~~~~~~~~~~~ EM ~~~~~~~~~~~~~~~~~~~~
00276 
00277       ESHandle< L1CaloEtScale > emScale ;
00278       iSetup.get< L1EmEtScaleRcd >().get( emScale ) ;
00279 
00280       // Isolated EM
00281       Handle< L1GctEmCandCollection > hwIsoEmCands ;
00282       iEvent.getByLabel( isoEmSource_, hwIsoEmCands ) ;
00283 
00284       if( !hwIsoEmCands.isValid() )
00285         {
00286           LogDebug("L1ExtraParticlesProd")
00287             << "\nWarning: L1GctEmCandCollection with " << isoEmSource_
00288             << "\nrequested in configuration, but not found in the event."
00289             << std::endl;
00290         }
00291       else
00292         {
00293 //       cout << "HW iso EM" << endl ;
00294 
00295           L1GctEmCandCollection::const_iterator emItr = hwIsoEmCands->begin() ;
00296           L1GctEmCandCollection::const_iterator emEnd = hwIsoEmCands->end() ;
00297           for( int i = 0 ; emItr != emEnd ; ++emItr, ++i )
00298             {
00299 //       cout << "#" << i
00300 //            << " name " << emItr->name()
00301 //            << " empty " << emItr->empty()
00302 //            << " rank " << emItr->rank()
00303 //            << " eta " << emItr->etaIndex()
00304 //            << " sign " << emItr->etaSign()
00305 //            << " phi " << emItr->phiIndex()
00306 //            << " iso " << emItr->isolated()
00307 //            << " bx " << emItr->bx()
00308 //            << endl ;
00309 
00310               if( !emItr->empty() &&
00311                   ( !centralBxOnly_ || emItr->bx() == 0 ) )
00312                 {
00313                   double et = emScale->et( emItr->rank() ) ;
00314 
00315 //          cout << "L1Extra et " << et << endl ;
00316 
00317                   isoEmColl->push_back(
00318                                        L1EmParticle( gctLorentzVector( et, *emItr, caloGeom, true ),
00319                                                      Ref< L1GctEmCandCollection >( hwIsoEmCands,
00320                                                                                    i ),
00321                                                      emItr->bx() ) ) ;
00322                 }
00323             }
00324         }
00325 
00326 
00327       // Non-isolated EM
00328       Handle< L1GctEmCandCollection > hwNonIsoEmCands ;
00329       iEvent.getByLabel( nonIsoEmSource_, hwNonIsoEmCands ) ;
00330 
00331       if( !hwNonIsoEmCands.isValid() )
00332         {
00333           LogDebug("L1ExtraParticlesProd")
00334             << "\nWarning: L1GctEmCandCollection with " << nonIsoEmSource_
00335             << "\nrequested in configuration, but not found in the event."
00336             << std::endl;
00337         }
00338       else
00339         {
00340 //       cout << "HW non-iso EM" << endl ;
00341           L1GctEmCandCollection::const_iterator emItr = hwNonIsoEmCands->begin() ;
00342           L1GctEmCandCollection::const_iterator emEnd = hwNonIsoEmCands->end() ;
00343           for( int i = 0 ; emItr != emEnd ; ++emItr, ++i )
00344             {
00345 //       cout << "#" << i
00346 //            << " name " << emItr->name()
00347 //            << " empty " << emItr->empty()
00348 //            << " rank " << emItr->rank()
00349 //            << " eta " << emItr->etaIndex()
00350 //            << " sign " << emItr->etaSign()
00351 //            << " phi " << emItr->phiIndex()
00352 //            << " iso " << emItr->isolated()
00353 //            << " bx " << emItr->bx()
00354 //            << endl ;
00355 
00356               if( !emItr->empty() &&
00357                   ( !centralBxOnly_ || emItr->bx() == 0 ) )
00358                 {
00359                   double et = emScale->et( emItr->rank() ) ;
00360 
00361 //          cout << "L1Extra et " << et << endl ;
00362 
00363                   nonIsoEmColl->push_back(
00364                                           L1EmParticle( gctLorentzVector( et, *emItr, caloGeom, true ),
00365                                                         Ref< L1GctEmCandCollection >( hwNonIsoEmCands,
00366                                                                                       i ),
00367                                                         emItr->bx() ) );
00368                 }
00369             }
00370         }
00371 
00372 
00373       // ~~~~~~~~~~~~~~~~~~~~ Jets ~~~~~~~~~~~~~~~~~~~~
00374 
00375       ESHandle< L1CaloEtScale > jetScale ;
00376       iSetup.get< L1JetEtScaleRcd >().get( jetScale ) ;
00377 
00378       // Central jets.
00379       Handle< L1GctJetCandCollection > hwCenJetCands ;
00380       iEvent.getByLabel( cenJetSource_, hwCenJetCands ) ;
00381 
00382       if( !hwCenJetCands.isValid() )
00383         {
00384           LogDebug("L1ExtraParticlesProd")
00385             << "\nWarning: L1GctJetCandCollection with " << cenJetSource_
00386             << "\nrequested in configuration, but not found in the event."
00387             << std::endl;
00388         }
00389       else
00390         {
00391 //       cout << "HW central jets" << endl ;
00392           L1GctJetCandCollection::const_iterator jetItr = hwCenJetCands->begin() ;
00393           L1GctJetCandCollection::const_iterator jetEnd = hwCenJetCands->end() ;
00394           for( int i = 0 ; jetItr != jetEnd ; ++jetItr, ++i )
00395             {
00396 //       cout << "#" << i
00397 //            << " name " << jetItr->name()
00398 //            << " empty " << jetItr->empty()
00399 //            << " rank " << jetItr->rank()
00400 //            << " eta " << jetItr->etaIndex()
00401 //            << " sign " << jetItr->etaSign()
00402 //            << " phi " << jetItr->phiIndex()
00403 //            << " cen " << jetItr->isCentral()
00404 //            << " for " << jetItr->isForward()
00405 //            << " tau " << jetItr->isTau()
00406 //            << " bx " << jetItr->bx()
00407 //            << endl ;
00408 
00409               if( !jetItr->empty() &&
00410                   ( !centralBxOnly_ || jetItr->bx() == 0 ) )
00411                 {
00412                   double et = jetScale->et( jetItr->rank() ) ;
00413 
00414 //          cout << "L1Extra et " << et << endl ;
00415 
00416                   cenJetColl->push_back(
00417                                         L1JetParticle( gctLorentzVector( et, *jetItr, caloGeom, true ),
00418                                                        Ref< L1GctJetCandCollection >( hwCenJetCands,
00419                                                                                       i ),
00420                                                        jetItr->bx() ) ) ;
00421                 }
00422             }
00423         }
00424 
00425 
00426       // Forward jets.
00427       Handle< L1GctJetCandCollection > hwForJetCands ;
00428       iEvent.getByLabel( forJetSource_, hwForJetCands ) ;
00429 
00430       if( !hwForJetCands.isValid() )
00431         {
00432           LogDebug("L1ExtraParticlesProd")
00433             << "\nWarning: L1GctEmCandCollection with " << forJetSource_
00434             << "\nrequested in configuration, but not found in the event."
00435             << std::endl;
00436         }
00437       else
00438         {
00439 //       cout << "HW forward jets" << endl ;
00440           L1GctJetCandCollection::const_iterator jetItr = hwForJetCands->begin() ;
00441           L1GctJetCandCollection::const_iterator jetEnd = hwForJetCands->end() ;
00442           for( int i = 0 ; jetItr != jetEnd ; ++jetItr, ++i )
00443             {
00444 //       cout << "#" << i
00445 //            << " name " << jetItr->name()
00446 //            << " empty " << jetItr->empty()
00447 //            << " rank " << jetItr->rank()
00448 //            << " eta " << jetItr->etaIndex()
00449 //            << " sign " << jetItr->etaSign()
00450 //            << " phi " << jetItr->phiIndex()
00451 //            << " cen " << jetItr->isCentral()
00452 //            << " for " << jetItr->isForward()
00453 //            << " tau " << jetItr->isTau()
00454 //            << " bx " << jetItr->bx()
00455 //            << endl ;
00456 
00457               if( !jetItr->empty() &&
00458                   ( !centralBxOnly_ || jetItr->bx() == 0 ) )
00459                 {
00460                   double et = jetScale->et( jetItr->rank() ) ;
00461 
00462 //          cout << "L1Extra et " << et << endl ;
00463 
00464                   forJetColl->push_back(
00465                                         L1JetParticle( gctLorentzVector( et, *jetItr, caloGeom, false ),
00466                                                        Ref< L1GctJetCandCollection >( hwForJetCands,
00467                                                                                       i ),
00468                                                        jetItr->bx() ) ) ;
00469                 }
00470             }
00471         }
00472 
00473 
00474       // Tau jets.
00475 //       cout << "HW tau jets" << endl ;
00476       Handle< L1GctJetCandCollection > hwTauJetCands ;
00477       iEvent.getByLabel( tauJetSource_, hwTauJetCands ) ;
00478 
00479       if( !hwTauJetCands.isValid() )
00480         {
00481           LogDebug("L1ExtraParticlesProd")
00482             << "\nWarning: L1GctJetCandCollection with " << tauJetSource_
00483             << "\nrequested in configuration, but not found in the event."
00484             << std::endl;
00485         }
00486       else
00487         {
00488           L1GctJetCandCollection::const_iterator jetItr = hwTauJetCands->begin() ;
00489           L1GctJetCandCollection::const_iterator jetEnd = hwTauJetCands->end() ;
00490           for( int i = 0 ; jetItr != jetEnd ; ++jetItr, ++i )
00491             {
00492 //       cout << "#" << i
00493 //            << " name " << jetItr->name()
00494 //            << " empty " << jetItr->empty()
00495 //            << " rank " << jetItr->rank()
00496 //            << " eta " << jetItr->etaIndex()
00497 //            << " sign " << jetItr->etaSign()
00498 //            << " phi " << jetItr->phiIndex()
00499 //            << " cen " << jetItr->isCentral()
00500 //            << " for " << jetItr->isForward()
00501 //            << " tau " << jetItr->isTau()
00502 //            << " bx " << jetItr->bx()
00503 //            << endl ;
00504 
00505               if( !jetItr->empty() &&
00506                   ( !centralBxOnly_ || jetItr->bx() == 0 ) )
00507                 {
00508                   double et = jetScale->et( jetItr->rank() ) ;
00509 
00510 //          cout << "L1Extra et " << et << endl ;
00511 
00512                   tauJetColl->push_back(
00513                                         L1JetParticle( gctLorentzVector( et, *jetItr, caloGeom, true ),
00514                                                        Ref< L1GctJetCandCollection >( hwTauJetCands,
00515                                                                                       i ),
00516                                                        jetItr->bx() ) ) ;
00517                 }
00518             }
00519         }
00520 
00521       // ~~~~~~~~~~~~~~~~~~~~ ET Sums ~~~~~~~~~~~~~~~~~~~~
00522 
00523       double etSumLSB = jetScale->linearLsb() ;
00524 
00525       Handle< L1GctEtTotalCollection > hwEtTotColl ;
00526       iEvent.getByLabel( etTotSource_, hwEtTotColl ) ;
00527 
00528       Handle< L1GctEtMissCollection > hwEtMissColl ;
00529       iEvent.getByLabel( etMissSource_, hwEtMissColl ) ;
00530 
00531       if( !hwEtTotColl.isValid() )
00532         {
00533           LogDebug("L1ExtraParticlesProd")
00534             << "\nWarning: L1GctEtTotalCollection with " << etTotSource_
00535             << "\nrequested in configuration, but not found in the event."
00536             << std::endl;
00537         }
00538       else if( !hwEtMissColl.isValid() )
00539         {
00540           LogDebug("L1ExtraParticlesProd")
00541             << "\nWarning: L1GctEtMissCollection with " << etMissSource_
00542             << "\nrequested in configuration, but not found in the event."
00543             << std::endl;
00544         }
00545       else
00546         {
00547           // Make a L1EtMissParticle even if either L1GctEtTotal or L1GctEtMiss
00548           // is missing for a given bx.  Keep track of which L1GctEtMiss objects
00549           // have a corresponding L1GctEtTotal object.
00550           std::vector< bool > etMissMatched ;
00551 
00552           L1GctEtMissCollection::const_iterator hwEtMissItr =
00553             hwEtMissColl->begin() ;
00554           L1GctEtMissCollection::const_iterator hwEtMissEnd =
00555             hwEtMissColl->end() ;
00556           for( ; hwEtMissItr != hwEtMissEnd ; ++hwEtMissItr )
00557             {
00558               etMissMatched.push_back( false ) ;
00559             }
00560       
00561           // Collate energy sums by bx
00562           L1GctEtTotalCollection::const_iterator hwEtTotItr =
00563             hwEtTotColl->begin() ;
00564           L1GctEtTotalCollection::const_iterator hwEtTotEnd =
00565             hwEtTotColl->end() ;
00566 
00567           int iTot = 0 ;
00568           for( ; hwEtTotItr != hwEtTotEnd ; ++hwEtTotItr, ++iTot )
00569             {
00570               int bx = hwEtTotItr->bx() ;
00571 
00572               if( !centralBxOnly_ || bx == 0 )
00573                 {
00574                   // ET bin low edge
00575                   double etTot =
00576                     ( hwEtTotItr->overFlow() ?
00577                       ( double ) L1GctEtTotal::kEtTotalMaxValue :
00578                       ( double ) hwEtTotItr->et() ) * etSumLSB + 1.e-6 ;
00579 
00580                   int iMiss = 0 ;
00581                   hwEtMissItr = hwEtMissColl->begin() ;
00582                   hwEtMissEnd = hwEtMissColl->end() ;
00583                   for( ; hwEtMissItr != hwEtMissEnd ; ++hwEtMissItr, ++iMiss )
00584                     {
00585                       if( hwEtMissItr->bx() == bx )
00586                         {
00587                           etMissMatched[ iMiss ] = true ;
00588                           break ;
00589                         }
00590                     }
00591 
00592                   double etMiss = 0. ;
00593                   double phi = 0. ;
00594                   math::PtEtaPhiMLorentzVector p4 ;
00595                   Ref< L1GctEtMissCollection > metRef ;
00596 
00597                   // If a L1GctEtMiss with the right bx is not found, itr == end.
00598                   if( hwEtMissItr != hwEtMissEnd )
00599                     {
00600                       // ET bin low edge
00601                       etMiss =
00602                         ( hwEtMissItr->overFlow() ?
00603                           ( double ) L1GctEtMiss::kEtMissMaxValue :
00604                           ( double ) hwEtMissItr->et() ) * etSumLSB + 1.e-6 ;
00605                       // keep x and y components non-zero and
00606                       // protect against roundoff.
00607 
00608                       phi = caloGeom->etSumPhiBinCenter( hwEtMissItr->phi() ) ;
00609 
00610                       p4 = math::PtEtaPhiMLorentzVector( etMiss,
00611                                                          0.,
00612                                                          phi,
00613                                                          0. ) ;
00614 
00615                       metRef = Ref< L1GctEtMissCollection >( hwEtMissColl, iMiss ) ;
00616 
00617 //             cout << "HW ET Sums " << endl
00618 //                  << "MET: phi " << hwEtMissItr->phi() << " = " << phi
00619 //                  << " et " << hwEtMissItr->et() << " = " << etMiss
00620 //                  << " EtTot " << hwEtTotItr->et() << " = " << etTot
00621 //                  << " bx " << bx
00622 //                  << endl ;
00623                     }
00624 //         else
00625 //           {
00626 //             cout << "HW ET Sums " << endl
00627 //                  << "MET: phi " << phi
00628 //                  << " et "<< etMiss
00629 //                  << " EtTot " << hwEtTotItr->et() << " = " << etTot
00630 //                  << " bx " << bx
00631 //                  << endl ;
00632 //           }
00633 
00634                   etMissColl->push_back(
00635                                         L1EtMissParticle(
00636                                                          p4,
00637                                                          L1EtMissParticle::kMET,
00638                                                          etTot,
00639                                                          metRef,
00640                                                          Ref< L1GctEtTotalCollection >( hwEtTotColl, iTot ),
00641                                                          Ref< L1GctHtMissCollection >(),
00642                                                          Ref< L1GctEtHadCollection >(),
00643                                                          bx
00644                                                          ) ) ;
00645                 }
00646             }
00647 
00648           if( !centralBxOnly_ )
00649             {
00650               // Make L1EtMissParticles for those L1GctEtMiss objects without
00651               // a matched L1GctEtTotal object.
00652 
00653               double etTot = 0. ;
00654 
00655               hwEtMissItr = hwEtMissColl->begin() ;
00656               hwEtMissEnd = hwEtMissColl->end() ;
00657               int iMiss = 0 ;
00658               for( ; hwEtMissItr != hwEtMissEnd ; ++hwEtMissItr, ++iMiss )
00659                 {
00660                   if( !etMissMatched[ iMiss ] )
00661                     {
00662                       int bx = hwEtMissItr->bx() ;
00663 
00664                       // ET bin low edge
00665                       double etMiss =
00666                         ( hwEtMissItr->overFlow() ?
00667                           ( double ) L1GctEtMiss::kEtMissMaxValue :
00668                           ( double ) hwEtMissItr->et() ) * etSumLSB + 1.e-6 ;
00669                       // keep x and y components non-zero and
00670                       // protect against roundoff.
00671 
00672                       double phi =
00673                         caloGeom->etSumPhiBinCenter( hwEtMissItr->phi() ) ;
00674 
00675                       math::PtEtaPhiMLorentzVector p4( etMiss,
00676                                                        0.,
00677                                                        phi,
00678                                                        0. ) ;
00679 
00680 //                cout << "HW ET Sums " << endl
00681 //                     << "MET: phi " << hwEtMissItr->phi() << " = " << phi
00682 //                     << " et " << hwEtMissItr->et() << " = " << etMiss
00683 //                     << " EtTot " << etTot
00684 //                     << " bx " << bx
00685 //                     << endl ;
00686 
00687                       etMissColl->push_back(
00688                                             L1EtMissParticle(
00689                                                              p4,
00690                                                              L1EtMissParticle::kMET,
00691                                                              etTot,
00692                                                              Ref< L1GctEtMissCollection >( hwEtMissColl, iMiss ),
00693                                                              Ref< L1GctEtTotalCollection >(),
00694                                                              Ref< L1GctHtMissCollection >(),
00695                                                              Ref< L1GctEtHadCollection >(),
00696                                                              bx
00697                                                              ) ) ;
00698                     }
00699                 }
00700             }
00701         }
00702 
00703       // ~~~~~~~~~~~~~~~~~~~~ HT Sums ~~~~~~~~~~~~~~~~~~~~
00704 
00705       Handle< L1GctEtHadCollection > hwEtHadColl ;
00706       iEvent.getByLabel( etHadSource_, hwEtHadColl ) ;
00707 
00708       Handle< L1GctHtMissCollection > hwHtMissColl ;
00709       if( !ignoreHtMiss_ )
00710         {
00711           iEvent.getByLabel( htMissSource_, hwHtMissColl ) ;
00712         }
00713 
00714       ESHandle< L1GctJetFinderParams > jetFinderParams ;
00715       iSetup.get< L1GctJetFinderParamsRcd >().get( jetFinderParams ) ;
00716       double htSumLSB = jetFinderParams->getHtLsbGeV();
00717 
00718       ESHandle< L1CaloEtScale > htMissScale ;
00719       std::vector< bool > htMissMatched ;
00720       if( !ignoreHtMiss_ )
00721         {
00722           iSetup.get< L1HtMissScaleRcd >().get( htMissScale ) ;
00723 
00724           if( !hwEtHadColl.isValid() )
00725             {
00726               LogDebug("L1ExtraParticlesProd")
00727                 << "\nWarning: L1GctEtHadCollection with " << etHadSource_
00728                 << "\nrequested in configuration, but not found in the event."
00729                 << std::endl;
00730             }
00731           else if( !hwHtMissColl.isValid() )
00732             {
00733               LogDebug("L1ExtraParticlesProd")
00734                 << "\nWarning: L1GctHtMissCollection with " << htMissSource_
00735                 << "\nrequested in configuration, but not found in the event."
00736                 << std::endl;
00737             }
00738           else
00739             {
00740               // Make a L1EtMissParticle even if either L1GctEtHad or L1GctHtMiss
00741               // is missing for a given bx. Keep track of which L1GctHtMiss objects
00742               // have a corresponding L1GctHtTotal object.
00743               L1GctHtMissCollection::const_iterator hwHtMissItr =
00744                 hwHtMissColl->begin() ;
00745               L1GctHtMissCollection::const_iterator hwHtMissEnd =
00746                 hwHtMissColl->end() ;
00747               for( ; hwHtMissItr != hwHtMissEnd ; ++hwHtMissItr )
00748                 {
00749                   htMissMatched.push_back( false ) ;
00750                 }
00751             }
00752         }
00753 
00754       if( !hwEtHadColl.isValid() )
00755         {
00756           LogDebug("L1ExtraParticlesProd")
00757             << "\nWarning: L1GctEtHadCollection with " << etHadSource_
00758             << "\nrequested in configuration, but not found in the event."
00759             << std::endl;
00760         }
00761       else if( !hwHtMissColl.isValid() )
00762         {
00763           LogDebug("L1ExtraParticlesProd")
00764             << "\nWarning: L1GctHtMissCollection with " << htMissSource_
00765             << "\nrequested in configuration, but not found in the event."
00766             << std::endl;
00767         }
00768       else
00769         {
00770           L1GctEtHadCollection::const_iterator hwEtHadItr =
00771             hwEtHadColl->begin() ;
00772           L1GctEtHadCollection::const_iterator hwEtHadEnd =
00773             hwEtHadColl->end() ;
00774 
00775           int iHad = 0 ;
00776           for( ; hwEtHadItr != hwEtHadEnd ; ++hwEtHadItr, ++iHad )
00777             {
00778               int bx = hwEtHadItr->bx() ;
00779 
00780               if( !centralBxOnly_ || bx == 0 )
00781                 {
00782                   // HT bin low edge
00783                   double htTot =
00784                     ( hwEtHadItr->overFlow() ?
00785                       ( double ) L1GctEtHad::kEtHadMaxValue :
00786                       ( double ) hwEtHadItr->et() ) * htSumLSB + 1.e-6 ;
00787 
00788                   double htMiss = 0. ;
00789                   double phi = 0. ;
00790                   math::PtEtaPhiMLorentzVector p4 ;
00791                   Ref< L1GctHtMissCollection > mhtRef ;
00792 
00793                   if( !ignoreHtMiss_ )
00794                     {
00795                       L1GctHtMissCollection::const_iterator hwHtMissItr =
00796                         hwHtMissColl->begin() ;
00797                       L1GctHtMissCollection::const_iterator hwHtMissEnd =
00798                         hwHtMissColl->end() ;
00799 
00800                       int iMiss = 0 ;
00801                       for( ; hwHtMissItr != hwHtMissEnd ; ++hwHtMissItr, ++iMiss )
00802                         {
00803                           if( hwHtMissItr->bx() == bx )
00804                             {
00805                               htMissMatched[ iMiss ] = true ;
00806                               break ;
00807                             }
00808                         }
00809 
00810                       // If a L1GctHtMiss with the right bx is not found, itr == end.
00811                       if( hwHtMissItr != hwHtMissEnd )
00812                         {
00813                           // HT bin low edge
00814                           htMiss =
00815                             htMissScale->et( hwHtMissItr->overFlow() ?
00816                                              htMissScale->rankScaleMax() :
00817                                              hwHtMissItr->et() ) + 1.e-6 ;
00818                           // keep x and y components non-zero and
00819                           // protect against roundoff.
00820 
00821                           phi =
00822                             caloGeom->htSumPhiBinCenter( hwHtMissItr->phi() ) ;
00823 
00824                           p4 = math::PtEtaPhiMLorentzVector( htMiss,
00825                                                              0.,
00826                                                              phi,
00827                                                              0. ) ;
00828 
00829                           mhtRef=Ref< L1GctHtMissCollection >( hwHtMissColl, iMiss );
00830 
00831 //                 cout << "HW HT Sums " << endl
00832 //                      << "MHT: phi " << hwHtMissItr->phi() << " = " << phi
00833 //                      << " ht " << hwHtMissItr->et() << " = " << htMiss
00834 //                      << " HtTot " << hwEtHadItr->et() << " = " << htTot
00835 //                      << " bx " << bx
00836 //                      << endl ;
00837                         }
00838 //             else
00839 //               {
00840 //                 cout << "HW HT Sums " << endl
00841 //                      << "MHT: phi " << phi
00842 //                      << " ht " << htMiss
00843 //                      << " HtTot " << hwEtHadItr->et() << " = " << htTot
00844 //                      << " bx " << bx
00845 //                      << endl ;
00846 //               }
00847                     }
00848 
00849                   htMissColl->push_back(
00850                                         L1EtMissParticle(
00851                                                          p4,
00852                                                          L1EtMissParticle::kMHT,
00853                                                          htTot,
00854                                                          Ref< L1GctEtMissCollection >(),
00855                                                          Ref< L1GctEtTotalCollection >(),
00856                                                          mhtRef,
00857                                                          Ref< L1GctEtHadCollection >( hwEtHadColl, iHad ),
00858                                                          bx
00859                                                          ) ) ;
00860                 }
00861             }
00862 
00863           if( !centralBxOnly_  && !ignoreHtMiss_ )
00864             {
00865               // Make L1EtMissParticles for those L1GctHtMiss objects without
00866               // a matched L1GctHtTotal object.
00867               double htTot = 0. ;
00868 
00869               L1GctHtMissCollection::const_iterator hwHtMissItr =
00870                 hwHtMissColl->begin() ;
00871               L1GctHtMissCollection::const_iterator hwHtMissEnd =
00872                 hwHtMissColl->end() ;
00873 
00874               int iMiss = 0 ;
00875               for( ; hwHtMissItr != hwHtMissEnd ; ++hwHtMissItr, ++iMiss )
00876                 {
00877                   if( !htMissMatched[ iMiss ] )
00878                     {
00879                       int bx = hwHtMissItr->bx() ;
00880 
00881                       // HT bin low edge
00882                       double htMiss =
00883                         htMissScale->et( hwHtMissItr->overFlow() ?
00884                                          htMissScale->rankScaleMax() :
00885                                          hwHtMissItr->et() ) + 1.e-6 ;
00886                       // keep x and y components non-zero and
00887                       // protect against roundoff.
00888 
00889                       double phi =
00890                         caloGeom->htSumPhiBinCenter( hwHtMissItr->phi() ) ;
00891 
00892                       math::PtEtaPhiMLorentzVector p4( htMiss,
00893                                                        0.,
00894                                                        phi,
00895                                                        0. ) ;
00896 
00897 //                cout << "HW HT Sums " << endl
00898 //                     << "MHT: phi " << hwHtMissItr->phi() << " = " << phi
00899 //                     << " ht " << hwHtMissItr->et() << " = " << htMiss
00900 //                     << " HtTot " << htTot
00901 //                     << " bx " << bx
00902 //                     << endl ;
00903 
00904                       htMissColl->push_back(
00905                                             L1EtMissParticle(
00906                                                              p4,
00907                                                              L1EtMissParticle::kMHT,
00908                                                              htTot,
00909                                                              Ref< L1GctEtMissCollection >(),
00910                                                              Ref< L1GctEtTotalCollection >(),
00911                                                              Ref< L1GctHtMissCollection >( hwHtMissColl, iMiss ),
00912                                                              Ref< L1GctEtHadCollection >(),
00913                                                              bx
00914                                                              ) ) ;
00915                     }
00916                 }
00917             }
00918         }
00919 
00920       // ~~~~~~~~~~~~~~~~~~~~ HF Rings ~~~~~~~~~~~~~~~~~~~~
00921 
00922       Handle< L1GctHFRingEtSumsCollection > hwHFEtSumsColl ;
00923       iEvent.getByLabel( hfRingEtSumsSource_, hwHFEtSumsColl ) ;
00924 
00925       Handle< L1GctHFBitCountsCollection > hwHFBitCountsColl ;
00926       iEvent.getByLabel( hfRingBitCountsSource_, hwHFBitCountsColl ) ;
00927 
00928       ESHandle< L1CaloEtScale > hfRingEtScale ;
00929       iSetup.get< L1HfRingEtScaleRcd >().get( hfRingEtScale ) ;
00930 
00931       if( !hwHFEtSumsColl.isValid() )
00932         {
00933           LogDebug("L1ExtraParticlesProd")
00934             << "\nWarning: L1GctHFRingEtSumsCollection with " << hfRingEtSumsSource_
00935             << "\nrequested in configuration, but not found in the event."
00936             << std::endl;
00937         }
00938       else if( !hwHFBitCountsColl.isValid() )
00939         {
00940           LogDebug("L1ExtraParticlesProd")
00941             << "\nWarning: L1GctHFBitCountsCollection with " << hfRingBitCountsSource_
00942             << "\nrequested in configuration, but not found in the event."
00943             << std::endl;
00944         }
00945       else
00946         {
00947           L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsItr =
00948             hwHFEtSumsColl->begin() ;
00949           L1GctHFRingEtSumsCollection::const_iterator hwHFEtSumsEnd =
00950             hwHFEtSumsColl->end() ;
00951 
00952           int iEtSums = 0 ;
00953           for( ; hwHFEtSumsItr != hwHFEtSumsEnd ; ++hwHFEtSumsItr, ++iEtSums )
00954             {
00955               int bx = hwHFEtSumsItr->bx() ;
00956 
00957               if( !centralBxOnly_ || bx == 0 )
00958                 {
00959                   L1GctHFBitCountsCollection::const_iterator hwHFBitCountsItr =
00960                     hwHFBitCountsColl->begin() ;
00961                   L1GctHFBitCountsCollection::const_iterator hwHFBitCountsEnd =
00962                     hwHFBitCountsColl->end() ;
00963 
00964                   int iBitCounts = 0 ;
00965                   for( ; hwHFBitCountsItr != hwHFBitCountsEnd ;
00966                        ++hwHFBitCountsItr, ++iBitCounts )
00967                     {
00968                       if( hwHFBitCountsItr->bx() == bx )
00969                         {
00970                           break ;
00971                         }
00972                     }
00973 
00974                   // If a L1GctHFBitCounts with the right bx is not found, itr == end.
00975                   if( hwHFBitCountsItr != hwHFBitCountsEnd )
00976                     {
00977                       // Construct L1HFRings only if both HW objects are present.
00978 
00979                       // cout << "HF Rings " << endl
00980 
00981                       // HF Et sums
00982                       double etSums[ L1HFRings::kNumRings ] ;
00983                       etSums[ L1HFRings::kRing1PosEta ] =
00984                         hfRingEtScale->et( hwHFEtSumsItr->etSum( 0 ) ) + 1.e-6 ;
00985                       etSums[ L1HFRings::kRing1NegEta ] =
00986                         hfRingEtScale->et( hwHFEtSumsItr->etSum( 1 ) ) + 1.e-6 ;
00987                       etSums[ L1HFRings::kRing2PosEta ] =
00988                         hfRingEtScale->et( hwHFEtSumsItr->etSum( 2 ) ) + 1.e-6 ;
00989                       etSums[ L1HFRings::kRing2NegEta ] =
00990                         hfRingEtScale->et( hwHFEtSumsItr->etSum( 3 ) ) + 1.e-6 ;
00991                       // protect against roundoff.
00992 
00993 //             cout << "HF Et Sums "
00994 //                  << hwHFEtSumsItr->etSum( 0 ) << " = "
00995 //                  << etSums[ L1HFRings::kRing1PosEta ] << " "
00996 //                  << hwHFEtSumsItr->etSum( 1 ) << " = "
00997 //                  << etSums[ L1HFRings::kRing1NegEta ] << " "
00998 //                  << hwHFEtSumsItr->etSum( 2 ) << " = "
00999 //                  << etSums[ L1HFRings::kRing2PosEta ] << " "
01000 //                  << hwHFEtSumsItr->etSum( 3 ) << " = "
01001 //                  << etSums[ L1HFRings::kRing2NegEta ]
01002 //                  << endl ;
01003 
01004                       // HF bit counts
01005                       int bitCounts[ L1HFRings::kNumRings ] ;
01006                       bitCounts[ L1HFRings::kRing1PosEta ] =
01007                         hwHFBitCountsItr->bitCount( 0 ) ;
01008                       bitCounts[ L1HFRings::kRing1NegEta ] =
01009                         hwHFBitCountsItr->bitCount( 1 ) ;
01010                       bitCounts[ L1HFRings::kRing2PosEta ] =
01011                         hwHFBitCountsItr->bitCount( 2 ) ;
01012                       bitCounts[ L1HFRings::kRing2NegEta ] =
01013                         hwHFBitCountsItr->bitCount( 3 ) ;
01014 
01015 //             cout << "HF bit counts "
01016 //                  << hwHFBitCountsItr->bitCount( 0 ) << " "
01017 //                  << hwHFBitCountsItr->bitCount( 1 ) << " "
01018 //                  << hwHFBitCountsItr->bitCount( 2 ) << " "
01019 //                  << hwHFBitCountsItr->bitCount( 3 ) << " "
01020 //                  << endl ;
01021 
01022                       hfRingsColl->push_back(
01023                                              L1HFRings( etSums,
01024                                                         bitCounts,
01025                                                         Ref< L1GctHFRingEtSumsCollection >( hwHFEtSumsColl,
01026                                                                                             iEtSums ),
01027                                                         Ref< L1GctHFBitCountsCollection >( hwHFBitCountsColl,
01028                                                                                            iBitCounts ),
01029                                                         bx
01030                                                         ) ) ;
01031                     }
01032                 }
01033             }
01034         }
01035    }
01036 
01037    OrphanHandle< L1EmParticleCollection > isoEmHandle =
01038      iEvent.put( isoEmColl, "Isolated" ) ;
01039 
01040    OrphanHandle< L1EmParticleCollection > nonIsoEmHandle =
01041      iEvent.put( nonIsoEmColl, "NonIsolated" ) ;
01042 
01043    OrphanHandle< L1JetParticleCollection > cenJetHandle =
01044      iEvent.put( cenJetColl, "Central" ) ;
01045 
01046    OrphanHandle< L1JetParticleCollection > forJetHandle =
01047      iEvent.put( forJetColl, "Forward" ) ;
01048 
01049    OrphanHandle< L1JetParticleCollection > tauJetHandle =
01050      iEvent.put( tauJetColl, "Tau" ) ;
01051 
01052    OrphanHandle< L1EtMissParticleCollection > etMissCollHandle =
01053      iEvent.put( etMissColl, "MET" ) ;
01054 
01055    OrphanHandle< L1EtMissParticleCollection > htMissCollHandle =
01056      iEvent.put( htMissColl, "MHT" ) ;
01057 
01058    OrphanHandle< L1HFRingsCollection > hfRingsCollHandle =
01059      iEvent.put( hfRingsColl ) ;
01060 }
01061 
01062 //math::XYZTLorentzVector
01063 math::PtEtaPhiMLorentzVector
01064 L1ExtraParticlesProd::gctLorentzVector( const double& et,
01065                                         const L1GctCand& cand,
01066                                         const L1CaloGeometry* geom,
01067                                         bool central )
01068 {
01069    // To keep x and y components non-zero.
01070    double etCorr = et + 1.e-6 ; // protect against roundoff, not only for et=0
01071 
01072    double eta = geom->etaBinCenter( cand.etaIndex(), central ) ;
01073 
01074 //    double tanThOver2 = exp( -eta ) ;
01075 //    double ez = etCorr * ( 1. - tanThOver2 * tanThOver2 ) / ( 2. * tanThOver2 );
01076 //    double e  = etCorr * ( 1. + tanThOver2 * tanThOver2 ) / ( 2. * tanThOver2 );
01077 
01078    double phi = geom->emJetPhiBinCenter( cand.phiIndex() ) ;
01079 
01080 //    return math::XYZTLorentzVector( etCorr * cos( phi ),
01081 //                                 etCorr * sin( phi ),
01082 //                                 ez,
01083 //                                 e ) ;
01084    return math::PtEtaPhiMLorentzVector( etCorr,
01085                                         eta,
01086                                         phi,
01087                                         0. ) ;
01088 }     
01089 
01090 // ------------ method called once each job just before starting event loop  ------------
01091 void 
01092 L1ExtraParticlesProd::beginJob()
01093 {
01094 }
01095 
01096 // ------------ method called once each job just after ending the event loop  ------------
01097 void 
01098 L1ExtraParticlesProd::endJob() {
01099 }
01100 
01101 //define this as a plug-in
01102 DEFINE_FWK_MODULE(L1ExtraParticlesProd);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/L1Trigger/L1ExtraFromDigis/src/L1ExtraParticlesProd.cc
