Go to the documentation of this file.00001 #include "RecoEgamma/EgammaIsolationAlgos/plugins/EleIsoDetIdCollectionProducer.h"
00002
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 #include "FWCore/Framework/interface/EventSetup.h"
00005 #include "FWCore/Framework/interface/ESHandle.h"
00006
00007 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00008
00009 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00010 #include "RecoCaloTools/Selectors/interface/CaloDualConeSelector.h"
00011 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00012 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00013 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00014 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00015 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00016 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00017 #include "RecoCaloTools/MetaCollections/interface/CaloRecHitMetaCollections.h"
00018 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00019
00020 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00021 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00022
00023 #include "DataFormats/DetId/interface/DetIdCollection.h"
00024
00025 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
00026 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00027
00028 EleIsoDetIdCollectionProducer::EleIsoDetIdCollectionProducer(const edm::ParameterSet& iConfig) :
00029 recHitsLabel_(iConfig.getParameter< edm::InputTag > ("recHitsLabel")),
00030 emObjectLabel_(iConfig.getParameter< edm::InputTag > ("emObjectLabel")),
00031 energyCut_(iConfig.getParameter<double>("energyCut")),
00032 etCut_(iConfig.getParameter<double>("etCut")),
00033 etCandCut_(iConfig.getParameter<double> ("etCandCut")),
00034 outerRadius_(iConfig.getParameter<double>("outerRadius")),
00035 innerRadius_(iConfig.getParameter<double>("innerRadius")),
00036 interestingDetIdCollection_(iConfig.getParameter<std::string>("interestingDetIdCollection")),
00037 severityLevelCut_(iConfig.getParameter<int>("severityLevelCut")),
00038 severityRecHitThreshold_(iConfig.getParameter<double>("severityRecHitThreshold")),
00039 spIdString_(iConfig.getParameter<std::string>("spikeIdString")),
00040 spIdThreshold_(iConfig.getParameter<double>("spikeIdThreshold")),
00041 v_chstatus_(iConfig.getParameter<std::vector<int> >("recHitFlagsToBeExcluded")) {
00042
00043
00044 if ( !spIdString_.compare("kE1OverE9") ) spId_ = EcalSeverityLevelAlgo::kE1OverE9;
00045 else if( !spIdString_.compare("kSwissCross") ) spId_ = EcalSeverityLevelAlgo::kSwissCross;
00046 else if( !spIdString_.compare("kSwissCrossBordersIncluded") ) spId_ = EcalSeverityLevelAlgo::kSwissCrossBordersIncluded;
00047 else spId_ = EcalSeverityLevelAlgo::kSwissCross;
00048
00049
00050 produces< DetIdCollection > (interestingDetIdCollection_) ;
00051 }
00052
00053 EleIsoDetIdCollectionProducer::~EleIsoDetIdCollectionProducer() {}
00054
00055 void EleIsoDetIdCollectionProducer::beginJob ()
00056 {}
00057
00058
00059 void
00060 EleIsoDetIdCollectionProducer::produce (edm::Event& iEvent,
00061 const edm::EventSetup& iSetup)
00062 {
00063 using namespace edm;
00064 using namespace std;
00065
00066
00067 Handle<reco::GsfElectronCollection> emObjectH;
00068 iEvent.getByLabel(emObjectLabel_,emObjectH);
00069
00070
00071 Handle<EcalRecHitCollection> recHitsH;
00072 iEvent.getByLabel(recHitsLabel_,recHitsH);
00073 std::auto_ptr<CaloRecHitMetaCollectionV> recHits_(0);
00074 recHits_ = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*recHitsH));
00075
00076 edm::ESHandle<CaloGeometry> pG;
00077 iSetup.get<CaloGeometryRecord>().get(pG);
00078 const CaloGeometry* caloGeom = pG.product();
00079
00080
00081 edm::ESHandle<EcalChannelStatus> chStatus;
00082 iSetup.get<EcalChannelStatusRcd>().get(chStatus);
00083
00084 CaloDualConeSelector *doubleConeSel_ = 0;
00085 if(recHitsLabel_.instance() == "EcalRecHitsEB")
00086 doubleConeSel_= new CaloDualConeSelector(innerRadius_,outerRadius_, &*pG, DetId::Ecal, EcalBarrel);
00087 else if(recHitsLabel_.instance() == "EcalRecHitsEE")
00088 doubleConeSel_= new CaloDualConeSelector(innerRadius_,outerRadius_, &*pG, DetId::Ecal, EcalEndcap);
00089
00090
00091 std::auto_ptr< DetIdCollection > detIdCollection (new DetIdCollection() ) ;
00092
00093 reco::GsfElectronCollection::const_iterator emObj;
00094 if(doubleConeSel_) {
00095 for (emObj = emObjectH->begin(); emObj != emObjectH->end(); emObj++) {
00096
00097 if(emObj->et() < etCandCut_) continue;
00098
00099 GlobalPoint pclu (emObj->caloPosition().x(),emObj->caloPosition().y(),emObj->caloPosition().z());
00100 std::auto_ptr<CaloRecHitMetaCollectionV> chosen = doubleConeSel_->select(pclu,*recHits_);
00101
00102 CaloRecHitMetaCollectionV::const_iterator recIt;
00103 for (recIt = chosen->begin(); recIt!= chosen->end () ; ++recIt) {
00104
00105 if ( fabs(recIt->energy()) < energyCut_) continue;
00106
00107
00108 double et = recIt->energy() *
00109 caloGeom->getPosition(recIt->detid()).perp() /
00110 caloGeom->getPosition(recIt->detid()).mag();
00111
00112 if ( fabs(et) < etCut_) continue;
00113
00114 if(recHitsLabel_.instance() == "EcalRecHitsEB" &&
00115 EcalSeverityLevelAlgo::severityLevel(
00116 EBDetId(recIt->detid()),
00117 *recHitsH,
00118 *chStatus,
00119 severityRecHitThreshold_,
00120 spId_,
00121 spIdThreshold_
00122 ) >= severityLevelCut_) continue;
00123
00124
00125
00126 std::vector<int>::const_iterator vit = std::find( v_chstatus_.begin(), v_chstatus_.end(), ((EcalRecHit*)(&*recIt))->recoFlag() );
00127 if ( vit != v_chstatus_.end() ) continue;
00128
00129
00130 if(std::find(detIdCollection->begin(),detIdCollection->end(),recIt->detid()) == detIdCollection->end())
00131 detIdCollection->push_back(recIt->detid());
00132 }
00133
00134 }
00135
00136 delete doubleConeSel_;
00137 }
00138
00139 iEvent.put( detIdCollection, interestingDetIdCollection_ );
00140 }