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/MessageLogger/interface/MessageLogger.h"
00011 #include "FWCore/Framework/interface/ESHandle.h"
00012 
00013 // Reconstruction Classes
00014 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00015 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00016 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00017 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00018 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00019 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00020 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00021 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
00022 
00023 // Geometry
00024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00025 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00026 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00027 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00028 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
00029 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
00030 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
00031 
00032 
00033 // Level 1 Trigger
00034 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
00035 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
00036 #include "CondFormats/L1TObjects/interface/L1CaloGeometry.h"
00037 #include "CondFormats/DataRecord/interface/L1CaloGeometryRecord.h"
00038 
00039 // EgammaCoreTools
00040 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
00041 #include "RecoEcal/EgammaCoreTools/interface/EcalEtaPhiRegion.h"
00042 
00043 // Class header file
00044 #include "RecoEgamma/EgammaHLTProducers/interface/EgammaHLTHybridClusterProducer.h"
00045 
00046 
00047 EgammaHLTHybridClusterProducer::EgammaHLTHybridClusterProducer(const edm::ParameterSet& ps)
00048 {
00049 
00050     // The debug level
00051   std::string debugString = ps.getParameter<std::string>("debugLevel");
00052   if      (debugString == "DEBUG")   debugL = HybridClusterAlgo::pDEBUG;
00053   else if (debugString == "INFO")    debugL = HybridClusterAlgo::pINFO;
00054   else                               debugL = HybridClusterAlgo::pERROR;
00055 
00056   basicclusterCollection_ = ps.getParameter<std::string>("basicclusterCollection");
00057   superclusterCollection_ = ps.getParameter<std::string>("superclusterCollection");
00058   hitproducer_ = ps.getParameter<edm::InputTag>("ecalhitproducer");
00059   hitcollection_ =ps.getParameter<std::string>("ecalhitcollection");
00060 
00061 
00062 
00063   // L1 matching parameters
00064   l1TagIsolated_ = ps.getParameter< edm::InputTag > ("l1TagIsolated");
00065   l1TagNonIsolated_ = ps.getParameter< edm::InputTag > ("l1TagNonIsolated");
00066 
00067   doIsolated_   = ps.getParameter<bool>("doIsolated");
00068 
00069   l1LowerThr_ = ps.getParameter<double> ("l1LowerThr");
00070   l1UpperThr_ = ps.getParameter<double> ("l1UpperThr");
00071   l1LowerThrIgnoreIsolation_ = ps.getParameter<double> ("l1LowerThrIgnoreIsolation");
00072 
00073   regionEtaMargin_   = ps.getParameter<double>("regionEtaMargin");
00074   regionPhiMargin_   = ps.getParameter<double>("regionPhiMargin");
00075 
00076   // Parameters for the position calculation:
00077   posCalculator_ = PositionCalc( ps.getParameter<edm::ParameterSet>("posCalcParameters") );
00078   
00079 
00080   hybrid_p = new HybridClusterAlgo(ps.getParameter<double>("HybridBarrelSeedThr"), 
00081                                    ps.getParameter<int>("step"),
00082                                    ps.getParameter<double>("ethresh"),
00083                                    ps.getParameter<double>("eseed"),
00084                                    ps.getParameter<double>("ewing"),
00085                                    ps.getParameter<std::vector<int> >("RecHitFlagToBeExcluded"),
00086                                    posCalculator_,
00087                                    debugL,
00088                                    ps.getParameter<bool>("dynamicEThresh"),
00089                                    ps.getParameter<double>("eThreshA"),
00090                                    ps.getParameter<double>("eThreshB"),
00091                                    ps.getParameter<std::vector<int> >("RecHitSeverityToBeExcluded"),
00092                                    ps.getParameter<double>("severityRecHitThreshold"),
00093                                    ps.getParameter<int>("severitySpikeId"),
00094                                    ps.getParameter<double>("severitySpikeThreshold"),
00095                                    ps.getParameter<bool>("excludeFlagged")
00096                                    );
00097 
00098   bool dynamicPhiRoad = ps.getParameter<bool>("dynamicPhiRoad");
00099     if (dynamicPhiRoad) {
00100      edm::ParameterSet bremRecoveryPset = ps.getParameter<edm::ParameterSet>("bremRecoveryPset");
00101      hybrid_p->setDynamicPhiRoad(bremRecoveryPset);
00102   }
00103 
00104 
00105   produces< reco::BasicClusterCollection >(basicclusterCollection_);
00106   produces< reco::SuperClusterCollection >(superclusterCollection_);
00107   nEvt_ = 0;
00108 }
00109 
00110 
00111 EgammaHLTHybridClusterProducer::~EgammaHLTHybridClusterProducer()
00112 {
00113   delete hybrid_p;
00114 }
00115 
00116 
00117 void EgammaHLTHybridClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es)
00118 {
00119   // get the hit collection from the event:
00120   edm::Handle<EcalRecHitCollection> rhcHandle;
00121   //  evt.getByType(rhcHandle);
00122   evt.getByLabel(hitproducer_.label(), hitcollection_, rhcHandle);
00123   if (!(rhcHandle.isValid())) 
00124     {
00125       if (debugL <= HybridClusterAlgo::pINFO)
00126         std::cout << "could not get a handle on the EcalRecHitCollection!" << std::endl;
00127       return;
00128     }
00129   const EcalRecHitCollection *hit_collection = rhcHandle.product();
00130 
00131   // get the collection geometry:
00132   edm::ESHandle<CaloGeometry> geoHandle;
00133   es.get<CaloGeometryRecord>().get(geoHandle);
00134   const CaloGeometry& geometry = *geoHandle;
00135   const CaloSubdetectorGeometry *geometry_p;
00136   std::auto_ptr<const CaloSubdetectorTopology> topology;
00137 
00138   edm::ESHandle<EcalChannelStatus> chStatus;
00139   es.get<EcalChannelStatusRcd>().get(chStatus);
00140   const EcalChannelStatus* theEcalChStatus = chStatus.product();
00141 
00142   //if (debugL == HybridClusterAlgo::pDEBUG)
00143   //std::cout << "\n\n\n" << hitcollection_ << "\n\n" << std::endl;
00144 
00145   if(hitcollection_ == "EcalRecHitsEB") {
00146     geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
00147     topology.reset(new EcalBarrelTopology(geoHandle));
00148   } else if(hitcollection_ == "EcalRecHitsEE") {
00149     geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
00150     topology.reset(new EcalEndcapTopology(geoHandle));
00151   } else if(hitcollection_ == "EcalRecHitsPS") {
00152     geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00153     topology.reset(new EcalPreshowerTopology (geoHandle));
00154   } else throw(std::runtime_error("\n\nHybrid Cluster Producer encountered invalied ecalhitcollection type.\n\n"));
00155     
00156   //Get the L1 EM Particle Collection
00157   //Get the L1 EM Particle Collection
00158   edm::Handle< l1extra::L1EmParticleCollection > emIsolColl ;
00159   if(doIsolated_)
00160     evt.getByLabel(l1TagIsolated_, emIsolColl);
00161   //Get the L1 EM Particle Collection
00162   edm::Handle< l1extra::L1EmParticleCollection > emNonIsolColl ;
00163   evt.getByLabel(l1TagNonIsolated_, emNonIsolColl);
00164 
00165   // Get the CaloGeometry
00166   edm::ESHandle<L1CaloGeometry> l1CaloGeom ;
00167   es.get<L1CaloGeometryRecord>().get(l1CaloGeom) ;
00168 
00169   std::vector<EcalEtaPhiRegion> regions;
00170 
00171   if(doIsolated_) {
00172     for( l1extra::L1EmParticleCollection::const_iterator emItr = emIsolColl->begin(); emItr != emIsolColl->end() ;++emItr ){
00173 
00174     if (emItr->et() > l1LowerThr_ && emItr->et() < l1UpperThr_
00175         //&&
00177 ) {
00178 
00179       //bool isolated = emItr->gctEmCand()->isolated();
00180       //if ((l1Isolated_ &&isolated) || (!l1Isolated_ &&!isolated)) {
00181 
00182       // Access the GCT hardware object corresponding to the L1Extra EM object.
00183       int etaIndex = emItr->gctEmCand()->etaIndex() ;
00184       int phiIndex = emItr->gctEmCand()->phiIndex() ;
00185       // Use the L1CaloGeometry to find the eta, phi bin boundaries.
00186       double etaLow  = l1CaloGeom->etaBinLowEdge( etaIndex ) ;
00187       double etaHigh = l1CaloGeom->etaBinHighEdge( etaIndex ) ;
00188       double phiLow  = l1CaloGeom->emJetPhiBinLowEdge( phiIndex ) ;
00189       double phiHigh = l1CaloGeom->emJetPhiBinHighEdge( phiIndex ) ;
00190 
00191        int isbarl=0;
00192        //Part of the region is in the barel if either the upper or lower
00193        //edge of the region is within the barrel
00194        if(((float)(etaLow)>-1.479 && (float)(etaLow)<1.479) || 
00195           ((float)(etaHigh)>-1.479 && (float)(etaHigh)<1.479)) isbarl=1;
00196 
00197         //std::cout<<"Hybrid etaindex "<<etaIndex<<" low hig : "<<etaLow<<" "<<etaHigh<<" phi low hig" <<phiLow<<" " << phiHigh<<" isforw "<<emItr->gctEmCand()->regionId().isForward()<<" isforwnew" <<isforw<< std::endl;
00198 
00199       etaLow -= regionEtaMargin_;
00200       etaHigh += regionEtaMargin_;
00201       phiLow -= regionPhiMargin_;
00202       phiHigh += regionPhiMargin_;
00203 
00204       if (etaHigh>1.479) etaHigh=1.479;
00205       if (etaLow<-1.479) etaLow=-1.479;
00206 
00207       if(isbarl) regions.push_back(EcalEtaPhiRegion(etaLow,etaHigh,phiLow,phiHigh));
00208 
00209     }
00210   }
00211   }
00212 
00213   if(!doIsolated_||l1LowerThrIgnoreIsolation_<64) {
00214     for( l1extra::L1EmParticleCollection::const_iterator emItr = emNonIsolColl->begin(); emItr != emNonIsolColl->end() ;++emItr ){
00215 
00216       if(doIsolated_&&emItr->et()<l1LowerThrIgnoreIsolation_) continue;
00217 
00218     if (emItr->et() > l1LowerThr_ && emItr->et() < l1UpperThr_
00219         //&&
00221 ) {
00222 
00223       //bool isolated = emItr->gctEmCand()->isolated();
00224       //if ((l1Isolated_ &&isolated) || (!l1Isolated_ &&!isolated)) {
00225 
00226       // Access the GCT hardware object corresponding to the L1Extra EM object.
00227       int etaIndex = emItr->gctEmCand()->etaIndex() ;
00228       int phiIndex = emItr->gctEmCand()->phiIndex() ;
00229       // Use the L1CaloGeometry to find the eta, phi bin boundaries.
00230       double etaLow  = l1CaloGeom->etaBinLowEdge( etaIndex ) ;
00231       double etaHigh = l1CaloGeom->etaBinHighEdge( etaIndex ) ;
00232       double phiLow  = l1CaloGeom->emJetPhiBinLowEdge( phiIndex ) ;
00233       double phiHigh = l1CaloGeom->emJetPhiBinHighEdge( phiIndex ) ;
00234 
00235        int isbarl=0;
00236        //Part of the region is in the barel if either the upper or lower
00237        //edge of the region is within the barrel
00238        if(((float)(etaLow)>-1.479 && (float)(etaLow)<1.479) || 
00239           ((float)(etaHigh)>-1.479 && (float)(etaHigh)<1.479)) isbarl=1;
00240 
00241         //std::cout<<"Hybrid etaindex "<<etaIndex<<" low hig : "<<etaLow<<" "<<etaHigh<<" phi low hig" <<phiLow<<" " << phiHigh<<" isforw "<<emItr->gctEmCand()->regionId().isForward()<<" isforwnew" <<isforw<< std::endl;
00242 
00243       etaLow -= regionEtaMargin_;
00244       etaHigh += regionEtaMargin_;
00245       phiLow -= regionPhiMargin_;
00246       phiHigh += regionPhiMargin_;
00247 
00248       if (etaHigh>1.479) etaHigh=1.479;
00249       if (etaLow<-1.479) etaLow=-1.479;
00250 
00251       if(isbarl) regions.push_back(EcalEtaPhiRegion(etaLow,etaHigh,phiLow,phiHigh));
00252 
00253     }
00254   }
00255   }
00256 
00257   // make the Basic clusters!
00258   reco::BasicClusterCollection basicClusters;
00259   hybrid_p->makeClusters(hit_collection, geometry_p, basicClusters, true, regions,theEcalChStatus);
00260   //if (debugL == HybridClusterAlgo::pDEBUG)
00261   //std::cout << "Hybrid Finished clustering - BasicClusterCollection returned to producer..." << std::endl;
00262 
00263   // create an auto_ptr to a BasicClusterCollection, copy the clusters into it and put in the Event:
00264   std::auto_ptr< reco::BasicClusterCollection > basicclusters_p(new reco::BasicClusterCollection);
00265   basicclusters_p->assign(basicClusters.begin(), basicClusters.end());
00266   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle =  evt.put(basicclusters_p, 
00267                                                                        basicclusterCollection_);
00268   //Basic clusters now in the event.
00269   //if (debugL == HybridClusterAlgo::pDEBUG)
00270   //std::cout << "Basic Clusters now put into event." << std::endl;
00271   
00272   //Weird though it is, get the BasicClusters back out of the event.  We need the
00273   //edm::Ref to these guys to make our superclusters for Hybrid.
00274 //  edm::Handle<reco::BasicClusterCollection> bccHandle;
00275  // evt.getByLabel("clusterproducer",basicclusterCollection_, bccHandle);
00276   if (!(bccHandle.isValid())) {
00277     //if (debugL <= HybridClusterAlgo::pINFO)
00278     //std::cout << "could not get a handle on the BasicClusterCollection!" << std::endl;
00279     return;
00280   }
00281   reco::BasicClusterCollection clusterCollection = *bccHandle;
00282   //if (debugL == HybridClusterAlgo::pDEBUG)
00283   //std::cout << "Got the BasicClusterCollection" << std::endl;
00284 
00285   reco::CaloClusterPtrVector clusterRefVector;
00286   for (unsigned int i = 0; i < clusterCollection.size(); i++){
00287     clusterRefVector.push_back(reco::CaloClusterPtr(bccHandle, i));
00288   }
00289 
00290   reco::SuperClusterCollection superClusters = hybrid_p->makeSuperClusters(clusterRefVector);
00291   //if (debugL == HybridClusterAlgo::pDEBUG)
00292   //std::cout << "Found: " << superClusters.size() << " superclusters." << std::endl;
00293 
00294   std::auto_ptr< reco::SuperClusterCollection > superclusters_p(new reco::SuperClusterCollection);
00295   superclusters_p->assign(superClusters.begin(), superClusters.end());
00296   evt.put(superclusters_p, superclusterCollection_);
00297 
00298   //if (debugL == HybridClusterAlgo::pDEBUG)
00299   //std::cout << "Hybrid Clusters (Basic/Super) added to the Event! :-)" << std::endl;
00300 
00301   nEvt_++;
00302 }
00303
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/RecoEgamma/EgammaHLTProducers/src/EgammaHLTHybridClusterProducer.cc
