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/EgammaHLTMulti5x5ClusterProducer.h"
00041 
00042 EgammaHLTMulti5x5ClusterProducer::EgammaHLTMulti5x5ClusterProducer(const edm::ParameterSet& ps)
00043 {
00044   // The verbosity level
00045   std::string verbosityString = ps.getParameter<std::string>("VerbosityLevel");
00046   if      (verbosityString == "DEBUG")   verbosity = Multi5x5ClusterAlgo::pDEBUG;
00047   else if (verbosityString == "WARNING") verbosity = Multi5x5ClusterAlgo::pWARNING;
00048   else if (verbosityString == "INFO")    verbosity = Multi5x5ClusterAlgo::pINFO;
00049   else                                   verbosity = Multi5x5ClusterAlgo::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   // Multi5x5 algorithm parameters
00066   double barrelSeedThreshold = ps.getParameter<double>("Multi5x5BarrelSeedThr");
00067   double endcapSeedThreshold = ps.getParameter<double>("Multi5x5EndcapSeedThr");
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   // exclude recHit flags from seeding
00083   std::vector<int> v_chstatus = ps.getParameter<std::vector<int> >("RecHitFlagToBeExcluded");
00084 
00085   // Produces a collection of barrel and a collection of endcap clusters
00086 
00087   produces< reco::BasicClusterCollection >(endcapClusterCollection_);
00088   produces< reco::BasicClusterCollection >(barrelClusterCollection_);
00089 
00090   Multi5x5_p = new Multi5x5ClusterAlgo(barrelSeedThreshold, endcapSeedThreshold, v_chstatus, posCalculator_,verbosity);
00091   
00092   /*
00093   shapeAlgo_ = ClusterShapeAlgo(providedParameters);//new
00094 
00095 
00096   clustershapecollectionEB_ = ps.getParameter<std::string>("clustershapecollectionEB");
00097   clustershapecollectionEE_ = ps.getParameter<std::string>("clustershapecollectionEE");
00098 
00099   //AssociationMap
00100   barrelClusterShapeAssociation_ = ps.getParameter<std::string>("barrelShapeAssociation");
00101   endcapClusterShapeAssociation_ = ps.getParameter<std::string>("endcapShapeAssociation");
00102   
00103   // Produces a collection of barrel and a collection of endcap clusters
00104 
00105   produces< reco::ClusterShapeCollection>(clustershapecollectionEE_);//new
00106   //produces< reco::BasicClusterCollection >(endcapClusterCollection_);
00107   produces< reco::ClusterShapeCollection>(clustershapecollectionEB_);//new
00108  // produces< reco::BasicClusterCollection >(barrelClusterCollection_);
00109   produces< reco::BasicClusterShapeAssociationCollection >(barrelClusterShapeAssociation_);//new
00110   produces< reco::BasicClusterShapeAssociationCollection >(endcapClusterShapeAssociation_);//new
00111 
00112   */
00113 
00114   nEvt_ = 0;
00115 }
00116 
00117 
00118 EgammaHLTMulti5x5ClusterProducer::~EgammaHLTMulti5x5ClusterProducer()
00119 {
00120   delete Multi5x5_p;
00121 }
00122 
00123 
00124 void EgammaHLTMulti5x5ClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es)
00125 {
00126   //Get the L1 EM Particle Collection
00127   edm::Handle< l1extra::L1EmParticleCollection > emIsolColl ;
00128   if(doIsolated_)
00129     evt.getByLabel(l1TagIsolated_, emIsolColl);
00130   //Get the L1 EM Particle Collection
00131   edm::Handle< l1extra::L1EmParticleCollection > emNonIsolColl ;
00132   evt.getByLabel(l1TagNonIsolated_, emNonIsolColl);
00133   // Get the CaloGeometry
00134   edm::ESHandle<L1CaloGeometry> l1CaloGeom ;
00135   es.get<L1CaloGeometryRecord>().get(l1CaloGeom) ;
00136 
00137   std::vector<EcalEtaPhiRegion> barrelRegions;
00138   std::vector<EcalEtaPhiRegion> endcapRegions;
00139 
00140   if(doIsolated_) {
00141     for( l1extra::L1EmParticleCollection::const_iterator emItr = emIsolColl->begin(); emItr != emIsolColl->end() ;++emItr ){
00142 
00143       if (emItr->et() > l1LowerThr_ && emItr->et() < l1UpperThr_) {
00144         
00145         // Access the GCT hardware object corresponding to the L1Extra EM object.
00146         int etaIndex = emItr->gctEmCand()->etaIndex() ;
00147         
00148         
00149         int phiIndex = emItr->gctEmCand()->phiIndex() ;
00150         // Use the L1CaloGeometry to find the eta, phi bin boundaries.
00151         double etaLow  = l1CaloGeom->etaBinLowEdge( etaIndex ) ;
00152         double etaHigh = l1CaloGeom->etaBinHighEdge( etaIndex ) ;
00153         double phiLow  = l1CaloGeom->emJetPhiBinLowEdge( phiIndex ) ;
00154         double phiHigh = l1CaloGeom->emJetPhiBinHighEdge( phiIndex ) ;
00155 
00156         //Attention isForward does not work
00157         int isforw=0;
00158         int isbarl=0;
00159         if((float)(etaHigh)>1.479 || (float)(etaLow)<-1.479) isforw=1;
00160         if(((float)(etaLow)>-1.479 && (float)(etaLow)<1.479) || 
00161            ((float)(etaHigh)>-1.479 && (float)(etaHigh)<1.479)) isbarl=1;
00162 
00163         //std::cout<<"Multi5x5 etaindex "<<etaIndex<<" low hig : "<<etaLow<<" "<<etaHigh<<" phi low hig" <<phiLow<<" " << phiHigh<<" isforw "<<emItr->gctEmCand()->regionId().isForward()<<" isforwnew" <<isforw<< std::endl;
00164         
00165         etaLow -= regionEtaMargin_;
00166         etaHigh += regionEtaMargin_;
00167         phiLow -= regionPhiMargin_;
00168         phiHigh += regionPhiMargin_;
00169 
00170         //if (emItr->gctEmCand()->regionId().isForward()) {
00171         if (isforw) {
00172           if (etaHigh>-1.479 && etaHigh<1.479) etaHigh=-1.479;
00173           if ( etaLow>-1.479 &&  etaLow<1.479) etaLow=1.479;
00174           EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00175           endcapRegions.push_back(region);
00176         }
00177         if (isbarl) {
00178           if (etaHigh>1.479) etaHigh=1.479;
00179           if (etaLow<-1.479) etaLow=-1.479;
00180           EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00181           barrelRegions.push_back(region);
00182         }
00183         EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00184         
00185       }
00186     }
00187   }
00188 
00189 
00190   if(!doIsolated_||l1LowerThrIgnoreIsolation_<64) {
00191     for( l1extra::L1EmParticleCollection::const_iterator emItr = emNonIsolColl->begin(); emItr != emNonIsolColl->end() ;++emItr ){
00192 
00193       if(doIsolated_&&emItr->et()<l1LowerThrIgnoreIsolation_) continue;
00194 
00195       if (emItr->et() > l1LowerThr_ && emItr->et() < l1UpperThr_) {
00196         
00197         // Access the GCT hardware object corresponding to the L1Extra EM object.
00198         int etaIndex = emItr->gctEmCand()->etaIndex() ;
00199         
00200         
00201         int phiIndex = emItr->gctEmCand()->phiIndex() ;
00202         // Use the L1CaloGeometry to find the eta, phi bin boundaries.
00203         double etaLow  = l1CaloGeom->etaBinLowEdge( etaIndex ) ;
00204         double etaHigh = l1CaloGeom->etaBinHighEdge( etaIndex ) ;
00205         double phiLow  = l1CaloGeom->emJetPhiBinLowEdge( phiIndex ) ;
00206         double phiHigh = l1CaloGeom->emJetPhiBinHighEdge( phiIndex ) ;
00207 
00208 
00209         int isforw=0;
00210         int isbarl=0;
00211         if((float)(etaHigh)>1.479 || (float)(etaLow)<-1.479) isforw=1;
00212         if(((float)(etaLow)>-1.479 && (float)(etaLow)<1.479) || 
00213            ((float)(etaHigh)>-1.479 && (float)(etaHigh)<1.479)) isbarl=1;
00214 
00215         //std::cout<<"Multi5x5 etaindex "<<etaIndex<<" low hig : "<<etaLow<<" "<<etaHigh<<" phi low hig" <<phiLow<<" " << phiHigh<<" isforw "<<emItr->gctEmCand()->regionId().isForward()<<" isforwnew" <<isforw<< std::endl;
00216         
00217         etaLow -= regionEtaMargin_;
00218         etaHigh += regionEtaMargin_;
00219         phiLow -= regionPhiMargin_;
00220         phiHigh += regionPhiMargin_;
00221 
00222         //if (emItr->gctEmCand()->regionId().isForward()) {
00223         if (isforw) {
00224           if (etaHigh>-1.479 && etaHigh<1.479) etaHigh=-1.479;
00225           if ( etaLow>-1.479 &&  etaLow<1.479) etaLow=1.479;
00226           EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00227           endcapRegions.push_back(region);
00228         }
00229         if (isbarl) {
00230           if (etaHigh>1.479) etaHigh=1.479;
00231           if (etaLow<-1.479) etaLow=-1.479;
00232           EcalEtaPhiRegion region(etaLow,etaHigh,phiLow,phiHigh);
00233           barrelRegions.push_back(region);
00234         }
00235         
00236       }
00237     }
00238   }
00239 
00240 
00241   if (doEndcaps_ 
00242       //&&endcapRegions.size()!=0
00243       ) {
00244 
00245     clusterizeECALPart(evt, es, endcapHitProducer_.label(), endcapHitCollection_, endcapClusterCollection_, endcapRegions, reco::CaloID::DET_ECAL_ENDCAP);//old
00246    }
00247   if (doBarrel_ 
00248       //&& barrelRegions.size()!=0
00249       ) {
00250     clusterizeECALPart(evt, es, barrelHitProducer_.label(), barrelHitCollection_, barrelClusterCollection_, barrelRegions, reco::CaloID::DET_ECAL_BARREL);//old
00251  
00252   }
00253   nEvt_++;
00254 }
00255 
00256 
00257 const EcalRecHitCollection * EgammaHLTMulti5x5ClusterProducer::getCollection(edm::Event& evt,
00258                                                                   const std::string& hitProducer_,
00259                                                                   const std::string& hitCollection_)
00260 {
00261   edm::Handle<EcalRecHitCollection> rhcHandle;
00262 
00263   evt.getByLabel(hitProducer_, hitCollection_, rhcHandle);
00264   if (!(rhcHandle.isValid())) 
00265     {
00266       std::cout << "could not get a handle on the EcalRecHitCollection!" << std::endl;
00267       edm::LogError("EgammaHLTMulti5x5ClusterProducerError") << "Error! can't get the product " << hitCollection_.c_str() ;
00268       return 0;
00269     } 
00270   return rhcHandle.product();
00271 }
00272 
00273 
00274 void EgammaHLTMulti5x5ClusterProducer::clusterizeECALPart(edm::Event &evt, const edm::EventSetup &es,
00275                                                const std::string& hitProducer,
00276                                                const std::string& hitCollection,
00277                                                const std::string& clusterCollection,
00278                                                const std::vector<EcalEtaPhiRegion>& regions,
00279                                                 const reco::CaloID::Detectors detector)
00280 {
00281 
00282   // get the hit collection from the event:
00283   const EcalRecHitCollection *hitCollection_p = getCollection(evt, hitProducer, hitCollection);
00284 
00285   // get the geometry and topology from the event setup:
00286   edm::ESHandle<CaloGeometry> geoHandle;
00287   es.get<CaloGeometryRecord>().get(geoHandle);
00288 
00289   const CaloSubdetectorGeometry *geometry_p;
00290   CaloSubdetectorTopology *topology_p;
00291 
00292   if (detector == reco::CaloID::DET_ECAL_BARREL) 
00293     {
00294       geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
00295       topology_p = new EcalBarrelTopology(geoHandle);
00296     }
00297   else
00298     {
00299       geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
00300       topology_p = new EcalEndcapTopology(geoHandle); 
00301    }
00302 
00303 
00304   const CaloSubdetectorGeometry *geometryES_p;
00305   geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00306 
00307   // Run the clusterization algorithm:
00308   reco::BasicClusterCollection clusters;
00309   clusters = Multi5x5_p->makeClusters(hitCollection_p, geometry_p, topology_p, geometryES_p, detector, true, regions);
00310 
00311   // create an auto_ptr to a BasicClusterCollection, copy the barrel clusters into it and put in the Event:
00312   std::auto_ptr< reco::BasicClusterCollection > clusters_p(new reco::BasicClusterCollection);
00313   clusters_p->assign(clusters.begin(), clusters.end());
00314   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle;
00315   if (detector == reco::CaloID::DET_ECAL_BARREL) 
00316     bccHandle = evt.put(clusters_p, barrelClusterCollection_);
00317   else
00318     bccHandle = evt.put(clusters_p, endcapClusterCollection_);
00319 
00320   delete topology_p;
00321 }
/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoEgamma/EgammaHLTProducers/src/EgammaHLTMulti5x5ClusterProducer.cc
