CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoEgamma/EgammaIsolationAlgos/plugins/EleIsoDetIdCollectionProducer.cc

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     //register your products
00050     produces< DetIdCollection > (interestingDetIdCollection_) ;
00051 }
00052 
00053 EleIsoDetIdCollectionProducer::~EleIsoDetIdCollectionProducer() {}
00054 
00055 void EleIsoDetIdCollectionProducer::beginJob () 
00056 {}
00057 
00058 // ------------ method called to produce the data  ------------
00059     void
00060 EleIsoDetIdCollectionProducer::produce (edm::Event& iEvent, 
00061         const edm::EventSetup& iSetup)
00062 {
00063     using namespace edm;
00064     using namespace std;
00065 
00066     //Get EM Object
00067     Handle<reco::GsfElectronCollection> emObjectH;
00068     iEvent.getByLabel(emObjectLabel_,emObjectH);
00069 
00070     // take EcalRecHits
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     //Get the channel status from the db
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     //Create empty output collections
00091     std::auto_ptr< DetIdCollection > detIdCollection (new DetIdCollection() ) ;
00092 
00093     reco::GsfElectronCollection::const_iterator emObj;
00094     if(doubleConeSel_) { //if cone selector was created
00095         for (emObj = emObjectH->begin(); emObj != emObjectH->end();  emObj++) { //Loop over candidates
00096 
00097             if(emObj->et() < etCandCut_) continue; //don't calculate if object hasn't enough energy
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) { // Select RecHits 
00104 
00105                 if ( fabs(recIt->energy()) < energyCut_) continue;  //dont fill if below E noise value
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;  //dont fill if below ET noise value
00113 
00114                 if(recHitsLabel_.instance() == "EcalRecHitsEB" && //make sure we have a barrel rechit
00115                    EcalSeverityLevelAlgo::severityLevel(          //call the severity level method
00116                        EBDetId(recIt->detid()),                   //passing the EBDetId
00117                        *recHitsH,                                 //the rechit collection in order to calculate the swiss crss
00118                        *chStatus,                                 //and the EcalChannelRecHitRcd
00119                         severityRecHitThreshold_,                 //only consider rechits with ET >
00120                         spId_,                                    //the SpikeId method (currently kE1OverE9 or kSwissCross)
00121                         spIdThreshold_                            //cut value for above
00122                    ) >= severityLevelCut_) continue;              //then if the severity level is too high, we continue to the next rechit
00123 
00124                 //Check based on flags to protect from recovered channels from non-read towers
00125                 //Assumption is that v_chstatus_ is empty unless doFlagChecks() has been called
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; // the recHit has to be excluded from the iso sum
00128 
00129 
00130                 if(std::find(detIdCollection->begin(),detIdCollection->end(),recIt->detid()) == detIdCollection->end()) 
00131                             detIdCollection->push_back(recIt->detid()); 
00132             } //end rechits
00133 
00134         } //end candidates
00135 
00136         delete doubleConeSel_;
00137     } //end if cone selector was created
00138     
00139     iEvent.put( detIdCollection, interestingDetIdCollection_ );
00140 }