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