CMSSW Reference Manual

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