Reference Manual

Go to the documentation of this file.
00001 // C/C++ headers
00002 #include <iostream>
00003 #include <vector>
00004 #include <memory>
00005 
00006 // Framework
00007 #include "FWCore/Framework/interface/Event.h"
00008 #include "FWCore/Framework/interface/EventSetup.h"
00009 #include "DataFormats/Common/interface/Handle.h"
00010 #include "FWCore/Framework/interface/ESHandle.h"
00011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00012 #include "FWCore/Utilities/interface/Exception.h"
00013 
00014 // Reconstruction Classes
00015 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00016 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00017 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00018 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00019 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00020 
00021 // Geometry
00022 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00023 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00024 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00025 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00026 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
00027 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
00028 
00029 // Level 1 Trigger
00030 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
00031 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
00032 #include "CondFormats/L1TObjects/interface/L1CaloGeometry.h"
00033 #include "CondFormats/DataRecord/interface/L1CaloGeometryRecord.h"
00034 
00035 // EgammaCoreTools
00036 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
00037 #include "RecoEcal/EgammaCoreTools/interface/EcalEtaPhiRegion.h"
00038 
00039 // Class header file
00040 #include "RecoEgamma/EgammaHLTProducers/interface/EgammaHLTIslandClusterProducer.h"
00041 
00042 EgammaHLTIslandClusterProducer::EgammaHLTIslandClusterProducer(const edm::ParameterSet& ps)
00043 {
00044   // The verbosity level
00045   std::string verbosityString = ps.getParameter<std::string>("VerbosityLevel");
00046   if      (verbosityString == "DEBUG")   verbosity = IslandClusterAlgo::pDEBUG;
00047   else if (verbosityString == "WARNING") verbosity = IslandClusterAlgo::pWARNING;
00048   else if (verbosityString == "INFO")    verbosity = IslandClusterAlgo::pINFO;
00049   else                                   verbosity = IslandClusterAlgo::pERROR;
00050 
00051   doBarrel_   = ps.getParameter<bool>("doBarrel");
00052   doEndcaps_   = ps.getParameter<bool>("doEndcaps");
00053   doIsolated_   = ps.getParameter<bool>("doIsolated");
00054 
00055   // Parameters to identify the hit collections
00056   barrelHitProducer_   = ps.getParameter<edm::InputTag>("barrelHitProducer");
00057   endcapHitProducer_   = ps.getParameter<edm::InputTag>("endcapHitProducer");
00058   barrelHitCollection_ = ps.getParameter<std::string>("barrelHitCollection");
00059   endcapHitCollection_ = ps.getParameter<std::string>("endcapHitCollection");
00060 
00061   // The names of the produced cluster collections
00062   barrelClusterCollection_  = ps.getParameter<std::string>("barrelClusterCollection");
00063   endcapClusterCollection_  = ps.getParameter<std::string>("endcapClusterCollection");
00064 
00065   // Island algorithm parameters
00066   double barrelSeedThreshold = ps.getParameter<double>("IslandBarrelSeedThr");
00067   double endcapSeedThreshold = ps.getParameter<double>("IslandEndcapSeedThr");
00068 
00069   // L1 matching parameters
00070   l1TagIsolated_ = ps.getParameter< edm::InputTag > ("l1TagIsolated");
00071   l1TagNonIsolated_ = ps.getParameter< edm::InputTag > ("l1TagNonIsolated");
00072   l1LowerThr_ = ps.getParameter<double> ("l1LowerThr");
00073   l1UpperThr_ = ps.getParameter<double> ("l1UpperThr");
00074   l1LowerThrIgnoreIsolation_ = ps.getParameter<double> ("l1LowerThrIgnoreIsolation");
00075 
00076   regionEtaMargin_   = ps.getParameter<double>("regionEtaMargin");
00077   regionPhiMargin_   = ps.getParameter<double>("regionPhiMargin");
00078 
00079    // Parameters for the position calculation:
00080   posCalculator_ = PositionCalc( ps.getParameter<edm::ParameterSet>("posCalcParameters") );
00081 
00082   // Produces a collection of barrel and a collection of endcap clusters
00083 
00084   produces< reco::BasicClusterCollection >(endcapClusterCollection_);
00085   produces< reco::BasicClusterCollection >(barrelClusterCollection_);
00086 
00087   island_p = new IslandClusterAlgo(barrelSeedThreshold, endcapSeedThreshold, posCalculator_,verbosity);
00088 
00089   nEvt_ = 0;
00090 }
00091 
00092 
00093 EgammaHLTIslandClusterProducer::~EgammaHLTIslandClusterProducer()
00094 {
00095   delete island_p;
00096 }
00097 
00098 
00099 void EgammaHLTIslandClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es)
00100 {
00101   //Get the L1 EM Particle Collection
00102   edm::Handle< l1extra::L1EmParticleCollection > emIsolColl ;
00103   if(doIsolated_)
00104     evt.getByLabel(l1TagIsolated_, emIsolColl);
00105   //Get the L1 EM Particle Collection
00106   edm::Handle< l1extra::L1EmParticleCollection > emNonIsolColl ;
00107   evt.getByLabel(l1TagNonIsolated_, emNonIsolColl);
00108   // Get the CaloGeometry
00109   edm::ESHandle<L1CaloGeometry> l1CaloGeom ;
00110   es.get<L1CaloGeometryRecord>().get(l1CaloGeom) ;
00111 
00112   std::vector<EcalEtaPhiRegion> barrelRegions;
00113   std::vector<EcalEtaPhiRegion> endcapRegions;
00114 
00115   if(doIsolated_) {
00116     for( l1extra::L1EmParticleCollection::const_iterator emItr = emIsolColl->begin(); emItr != emIsolColl->end() ;++emItr ){
00117 
00118       if (emItr->et() > l1LowerThr_ && emItr->et() < l1UpperThr_) {
00119         
00120         // Access the GCT hardware object corresponding to the L1Extra EM object.
00121         int etaIndex = emItr->gctEmCand()->etaIndex() ;
00122         
00123         
00124         int phiIndex = emItr->gctEmCand()->phiIndex() ;
00125         // Use the L1CaloGeometry to find the eta, phi bin boundaries.
00126         double etaLow  = l1CaloGeom->etaBinLowEdge( etaIndex ) ;
00127         double etaHigh = l1CaloGeom->etaBinHighEdge( etaIndex ) ;
00128         double phiLow  = l1CaloGeom->emJetPhiBinLowEdge( phiIndex ) ;
00129         double phiHigh = l1CaloGeom->emJetPhiBinHighEdge( phiIndex ) ;
00130 
00131         //Attention isForward does not work
00132         int isforw=0;
00133         int isbarl=0;
00134         if((float)(etaHigh)>1.479 || (float)(etaLow)<-1.479) isforw=1;
00135         if(((float)(etaLow)>-1.479 && (float)(etaLow)<1.479) || 
00136            ((float)(etaHigh)>-1.479 && (float)(etaHigh)<1.479)) isbarl=1;
00137 
00138         //std::cout<<"Island etaindex "<<etaIndex<<" low hig : "<<etaLow<<" "<<etaHigh<<" phi low hig" <<phiLow<<" " << phiHigh<<" isforw "<<emItr->gctEmCand()->regionId().isForward()<<" isforwnew" <<isforw<< std::endl;
00139         
00140         etaLow -= regionEtaMargin_;
00141         etaHigh += regionEtaMargin_;
00142         phiLow -= regionPhiMargin_;
00143         phiHigh += regionPhiMargin_;
00144 
00145         //if (emItr->gctEmCand()->regionId().isForward()) {
00146         if (isforw) {
00147           if (etaHigh>-1.479 && etaHigh<1.479) etaHigh=-1.479;
00148           if ( etaLow>-1.479 &&  etaLow<1.479) etaLow=1.479;
00149           EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00150           endcapRegions.push_back(region);
00151         }
00152         if (isbarl) {
00153           if (etaHigh>1.479) etaHigh=1.479;
00154           if (etaLow<-1.479) etaLow=-1.479;
00155           EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00156           barrelRegions.push_back(region);
00157         }
00158         EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00159         
00160       }
00161     }
00162   }
00163 
00164 
00165   if(!doIsolated_||l1LowerThrIgnoreIsolation_<64) {
00166     for( l1extra::L1EmParticleCollection::const_iterator emItr = emNonIsolColl->begin(); emItr != emNonIsolColl->end() ;++emItr ){
00167 
00168       if(doIsolated_&&emItr->et()<l1LowerThrIgnoreIsolation_) continue;
00169 
00170       if (emItr->et() > l1LowerThr_ && emItr->et() < l1UpperThr_) {
00171         
00172         // Access the GCT hardware object corresponding to the L1Extra EM object.
00173         int etaIndex = emItr->gctEmCand()->etaIndex() ;
00174         
00175         
00176         int phiIndex = emItr->gctEmCand()->phiIndex() ;
00177         // Use the L1CaloGeometry to find the eta, phi bin boundaries.
00178         double etaLow  = l1CaloGeom->etaBinLowEdge( etaIndex ) ;
00179         double etaHigh = l1CaloGeom->etaBinHighEdge( etaIndex ) ;
00180         double phiLow  = l1CaloGeom->emJetPhiBinLowEdge( phiIndex ) ;
00181         double phiHigh = l1CaloGeom->emJetPhiBinHighEdge( phiIndex ) ;
00182 
00183 
00184         int isforw=0;
00185         int isbarl=0;
00186         if((float)(etaHigh)>1.479 || (float)(etaLow)<-1.479) isforw=1;
00187         if(((float)(etaLow)>-1.479 && (float)(etaLow)<1.479) || 
00188            ((float)(etaHigh)>-1.479 && (float)(etaHigh)<1.479)) isbarl=1;
00189 
00190         //std::cout<<"Island etaindex "<<etaIndex<<" low hig : "<<etaLow<<" "<<etaHigh<<" phi low hig" <<phiLow<<" " << phiHigh<<" isforw "<<emItr->gctEmCand()->regionId().isForward()<<" isforwnew" <<isforw<< std::endl;
00191         
00192         etaLow -= regionEtaMargin_;
00193         etaHigh += regionEtaMargin_;
00194         phiLow -= regionPhiMargin_;
00195         phiHigh += regionPhiMargin_;
00196 
00197         //if (emItr->gctEmCand()->regionId().isForward()) {
00198         if (isforw) {
00199           if (etaHigh>-1.479 && etaHigh<1.479) etaHigh=-1.479;
00200           if ( etaLow>-1.479 &&  etaLow<1.479) etaLow=1.479;
00201           EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00202           endcapRegions.push_back(region);
00203         }
00204         if (isbarl) {
00205           if (etaHigh>1.479) etaHigh=1.479;
00206           if (etaLow<-1.479) etaLow=-1.479;
00207           EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00208           barrelRegions.push_back(region);
00209         }
00210         
00211       }
00212     }
00213   }
00214 
00215   if (doEndcaps_ 
00216       //&&endcapRegions.size()!=0
00217       ) {
00218 
00219     clusterizeECALPart(evt, es, endcapHitProducer_.label(), endcapHitCollection_, endcapClusterCollection_, endcapRegions, IslandClusterAlgo::endcap);
00220   }
00221   if (doBarrel_ 
00222       //&& barrelRegions.size()!=0
00223       ) {
00224     clusterizeECALPart(evt, es, barrelHitProducer_.label(), barrelHitCollection_, barrelClusterCollection_, barrelRegions, IslandClusterAlgo::barrel);
00225   }
00226   nEvt_++;
00227 }
00228 
00229 
00230 const EcalRecHitCollection * EgammaHLTIslandClusterProducer::getCollection(edm::Event& evt,
00231                                                                   const std::string& hitProducer_,
00232                                                                   const std::string& hitCollection_)
00233 {
00234   edm::Handle<EcalRecHitCollection> rhcHandle;
00235 
00236   evt.getByLabel(hitProducer_, hitCollection_, rhcHandle);
00237   if (!(rhcHandle.isValid())) 
00238     {
00239       std::cout << "could not get a handle on the EcalRecHitCollection!" << std::endl;
00240       edm::LogError("EgammaHLTIslandClusterProducerError") << "Error! can't get the product " << hitCollection_.c_str() ;
00241       return 0;
00242     } 
00243   return rhcHandle.product();
00244 }
00245 
00246 
00247 void EgammaHLTIslandClusterProducer::clusterizeECALPart(edm::Event &evt, const edm::EventSetup &es,
00248                                                const std::string& hitProducer,
00249                                                const std::string& hitCollection,
00250                                                const std::string& clusterCollection,
00251                                                const std::vector<EcalEtaPhiRegion>& regions,
00252                                                const IslandClusterAlgo::EcalPart& ecalPart)
00253 {
00254   // get the hit collection from the event:
00255   const EcalRecHitCollection *hitCollection_p = getCollection(evt, hitProducer, hitCollection);
00256 
00257   // get the geometry and topology from the event setup:
00258   edm::ESHandle<CaloGeometry> geoHandle;
00259   es.get<CaloGeometryRecord>().get(geoHandle);
00260 
00261   const CaloSubdetectorGeometry *geometry_p;
00262   CaloSubdetectorTopology *topology_p;
00263 
00264   if (ecalPart == IslandClusterAlgo::barrel) 
00265     {
00266       geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
00267       topology_p = new EcalBarrelTopology(geoHandle);
00268     }
00269   else
00270     {
00271       geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
00272       topology_p = new EcalEndcapTopology(geoHandle); 
00273    }
00274 
00275   const CaloSubdetectorGeometry *geometryES_p;
00276   geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00277 
00278   // Run the clusterization algorithm:
00279   reco::BasicClusterCollection clusters;
00280   clusters = island_p->makeClusters(hitCollection_p, geometry_p, topology_p, geometryES_p, ecalPart, true, regions);
00281 
00282   // create an auto_ptr to a BasicClusterCollection, copy the barrel clusters into it and put in the Event:
00283   std::auto_ptr< reco::BasicClusterCollection > clusters_p(new reco::BasicClusterCollection);
00284   clusters_p->assign(clusters.begin(), clusters.end());
00285   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle;
00286   if (ecalPart == IslandClusterAlgo::barrel) 
00287     bccHandle = evt.put(clusters_p, barrelClusterCollection_);
00288   else
00289     bccHandle = evt.put(clusters_p, endcapClusterCollection_);
00290 
00291   delete topology_p;
00292 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoEgamma/EgammaHLTProducers/src/EgammaHLTIslandClusterProducer.cc
