CMS 3D CMS Logo

InterestingDetIdCollectionProducer.cc

Go to the documentation of this file.
00001 #include "RecoEcal/EgammaClusterProducers/interface/InterestingDetIdCollectionProducer.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/EcalDetId/interface/EBDetId.h"
00010 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00011 
00012 #include "DataFormats/DetId/interface/DetIdCollection.h"
00013 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00014 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00015 
00016 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00017 #include "Geometry/CaloTopology/interface/CaloTopology.h"
00018 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
00019 
00020 
00021 
00022 InterestingDetIdCollectionProducer::InterestingDetIdCollectionProducer(const edm::ParameterSet& iConfig) 
00023 {
00024 
00025   recHitsLabel_ = iConfig.getParameter< edm::InputTag > ("recHitsLabel");
00026   basicClusters_ = iConfig.getParameter< edm::InputTag > ("basicClustersLabel");
00027 
00028   interestingDetIdCollection_ = iConfig.getParameter<std::string>("interestingDetIdCollection");
00029   
00030   minimalEtaSize_ = iConfig.getParameter<int> ("etaSize");
00031   minimalPhiSize_ = iConfig.getParameter<int> ("phiSize");
00032   if ( minimalPhiSize_ % 2 == 0 ||  minimalEtaSize_ % 2 == 0)
00033     edm::LogError("InterestingDetIdCollectionProducerError") << "Size of eta/phi should be odd numbers";
00034  
00035    //register your products
00036   produces< DetIdCollection > (interestingDetIdCollection_) ;
00037   
00038 }
00039 
00040 
00041 InterestingDetIdCollectionProducer::~InterestingDetIdCollectionProducer()
00042 {}
00043 
00044 void InterestingDetIdCollectionProducer::beginJob (const edm::EventSetup& iSetup)  
00045 {
00046   edm::ESHandle<CaloTopology> theCaloTopology;
00047   iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
00048   caloTopology_ = &(*theCaloTopology); 
00049 }
00050 
00051 // ------------ method called to produce the data  ------------
00052 void
00053 InterestingDetIdCollectionProducer::produce (edm::Event& iEvent, 
00054                                 const edm::EventSetup& iSetup)
00055 {
00056    using namespace edm;
00057    using namespace std;
00058 
00059    // take BasicClusters
00060   Handle<reco::BasicClusterCollection> pClusters;
00061   iEvent.getByLabel(basicClusters_, pClusters);
00062   
00063   // take EcalRecHits
00064   Handle<EcalRecHitCollection> recHitsHandle;
00065   iEvent.getByLabel(recHitsLabel_,recHitsHandle);
00066 
00067   //Create empty output collections
00068   std::auto_ptr< DetIdCollection > detIdCollection (new DetIdCollection() ) ;
00069   
00070   reco::BasicClusterCollection::const_iterator clusIt;
00071   
00072   for (clusIt=pClusters->begin(); clusIt!=pClusters->end(); clusIt++) {
00073     //PG barrel
00074     
00075     float eMax=0.;
00076     DetId eMaxId(0);
00077 
00078     std::vector<DetId> clusterDetIds = (*clusIt).getHitsByDetId();
00079     std::vector<DetId>::iterator posCurrent;
00080 
00081     EcalRecHit testEcalRecHit;
00082     
00083     for(posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); posCurrent++)
00084       {
00085         EcalRecHitCollection::const_iterator itt = recHitsHandle->find(*posCurrent);
00086         if ((!((*posCurrent).null())) && (itt != recHitsHandle->end()) && ((*itt).energy() > eMax) )
00087           {
00088             eMax = (*itt).energy();
00089             eMaxId = (*itt).id();
00090           }
00091       }
00092     
00093     if (eMaxId.null())
00094     continue;
00095     
00096     const CaloSubdetectorTopology* topology  = caloTopology_->getSubdetectorTopology(eMaxId.det(),eMaxId.subdetId());
00097     std::vector<DetId> xtalsToStore=topology->getWindow(eMaxId,minimalEtaSize_,minimalPhiSize_);
00098     std::vector<DetId> xtalsInClus=(*clusIt).getHitsByDetId();
00099     
00100     for (unsigned int ii=0;ii<xtalsInClus.size();ii++)
00101       {
00102         if (std::find(xtalsToStore.begin(),xtalsToStore.end(),xtalsInClus[ii]) == xtalsToStore.end())
00103           xtalsToStore.push_back(xtalsInClus[ii]);
00104       }
00105     
00106     for (unsigned int iCry=0;iCry<xtalsToStore.size();iCry++)
00107       {
00108         if ( 
00109             std::find(detIdCollection->begin(),detIdCollection->end(),xtalsToStore[iCry]) == detIdCollection->end()
00110             )
00111           detIdCollection->push_back(xtalsToStore[iCry]);
00112       }     
00113   }
00114   
00115   //  std::cout << "Interesting DetId Collection size is " << detIdCollection->size() <<  " BCs are " << pClusters->size() << std::endl;
00116   iEvent.put( detIdCollection, interestingDetIdCollection_ );
00117 
00118 }

Generated on Tue Jun 9 17:43:15 2009 for CMSSW by  doxygen 1.5.4