Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    CentralityProducer
00004 // Class:      CentralityProducer
00005 // 
00013 //
00014 // Original Author:  Yetkin Yilmaz, Young Soo Park
00015 //         Created:  Wed Jun 11 15:31:41 CEST 2008
00016 // $Id: CentralityProducer.cc,v 1.40 2012/12/26 23:38:37 wmtan Exp $
00017 //
00018 //
00019 
00020 
00021 // system include files
00022 #include <memory>
00023 #include <iostream>
00024 
00025 // user include files
00026 #include "FWCore/Framework/interface/Frameworkfwd.h"
00027 #include "FWCore/Framework/interface/EDFilter.h"
00028 #include "FWCore/Framework/interface/Event.h"
00029 #include "FWCore/Framework/interface/EventSetup.h"
00030 
00031 #include "FWCore/Framework/interface/MakerMacros.h"
00032 #include "FWCore/Framework/interface/ESHandle.h"
00033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00034 #include "FWCore/Utilities/interface/InputTag.h"
00035 
00036 #include "DataFormats/Candidate/interface/Candidate.h"
00037 
00038 #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
00039 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00040 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00041 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00042 #include "DataFormats/Common/interface/Ref.h"
00043 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
00044 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
00045 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00046 #include "DataFormats/TrackReco/interface/Track.h"
00047 
00048 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00049 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00050 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00051 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00052 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
00053 
00054 using namespace std;
00055 
00056 //
00057 // class declaration
00058 //
00059 namespace reco{
00060 
00061 class CentralityProducer : public edm::EDFilter {
00062    public:
00063       explicit CentralityProducer(const edm::ParameterSet&);
00064       ~CentralityProducer();
00065 
00066    private:
00067       virtual void beginJob() ;
00068       virtual bool filter(edm::Event&, const edm::EventSetup&);
00069       virtual void endJob() ;
00070       
00071       // ----------member data ---------------------------
00072 
00073    bool recoLevel_;
00074 
00075    bool doFilter_;
00076    bool produceHFhits_;
00077    bool produceHFtowers_;
00078    bool produceEcalhits_;
00079    bool produceBasicClusters_;
00080    bool produceZDChits_;
00081    bool produceETmidRap_;
00082    bool producePixelhits_;
00083    bool produceTracks_;
00084    bool reuseAny_;
00085    bool producePixelTracks_;
00086 
00087   bool doPixelCut_;
00088 
00089   double midRapidityRange_;
00090   double trackPtCut_;
00091   double trackEtaCut_;
00092 
00093    edm::InputTag  srcHFhits_;   
00094    edm::InputTag  srcTowers_;
00095    edm::InputTag srcEEhits_;
00096    edm::InputTag srcEBhits_;
00097    edm::InputTag srcBasicClustersEE_;
00098    edm::InputTag srcBasicClustersEB_;
00099    edm::InputTag srcZDChits_;
00100    edm::InputTag srcPixelhits_;
00101    edm::InputTag srcTracks_;
00102    edm::InputTag srcPixelTracks_;
00103 
00104    edm::InputTag reuseTag_;
00105 
00106   bool useQuality_;
00107   reco::TrackBase::TrackQuality trackQuality_;
00108 
00109   const TrackerGeometry* trackGeo_;
00110   const CaloGeometry* caloGeo_;
00111 
00112 };
00113 
00114 //
00115 // constants, enums and typedefs
00116 //
00117 
00118 
00119 //
00120 // static data member definitions
00121 //
00122 
00123 //
00124 // constructors and destructor
00125 //
00126   CentralityProducer::CentralityProducer(const edm::ParameterSet& iConfig) :
00127     trackGeo_(0),
00128     caloGeo_(0)
00129 {
00130    //register your products
00131    doFilter_ = iConfig.getParameter<bool>("doFilter");
00132    produceHFhits_ = iConfig.getParameter<bool>("produceHFhits");
00133    produceHFtowers_ = iConfig.getParameter<bool>("produceHFtowers");
00134    produceBasicClusters_ = iConfig.getParameter<bool>("produceBasicClusters");
00135    produceEcalhits_ = iConfig.getParameter<bool>("produceEcalhits");
00136    produceZDChits_ = iConfig.getParameter<bool>("produceZDChits");
00137    produceETmidRap_ = iConfig.getParameter<bool>("produceETmidRapidity");
00138    producePixelhits_ = iConfig.getParameter<bool>("producePixelhits");
00139    produceTracks_ = iConfig.getParameter<bool>("produceTracks");
00140    producePixelTracks_ = iConfig.getParameter<bool>("producePixelTracks");
00141 
00142    midRapidityRange_ = iConfig.getParameter<double>("midRapidityRange");
00143    trackPtCut_ = iConfig.getParameter<double>("trackPtCut");
00144    trackEtaCut_ = iConfig.getParameter<double>("trackEtaCut");
00145 
00146    if(produceHFhits_)  srcHFhits_ = iConfig.getParameter<edm::InputTag>("srcHFhits");
00147    if(produceHFtowers_ || produceETmidRap_) srcTowers_ = iConfig.getParameter<edm::InputTag>("srcTowers");
00148 
00149    if(produceEcalhits_){
00150       srcEBhits_ = iConfig.getParameter<edm::InputTag>("srcEBhits");
00151       srcEEhits_ = iConfig.getParameter<edm::InputTag>("srcEEhits");
00152    }
00153    if(produceBasicClusters_){
00154       srcBasicClustersEE_ = iConfig.getParameter<edm::InputTag>("srcBasicClustersEE");
00155       srcBasicClustersEB_ = iConfig.getParameter<edm::InputTag>("srcBasicClustersEB");
00156    }
00157    if(produceZDChits_) srcZDChits_ = iConfig.getParameter<edm::InputTag>("srcZDChits");
00158    if(producePixelhits_){
00159      srcPixelhits_ = iConfig.getParameter<edm::InputTag>("srcPixelhits");
00160      doPixelCut_ = iConfig.getParameter<bool>("doPixelCut");
00161    }
00162    if(produceTracks_) srcTracks_ = iConfig.getParameter<edm::InputTag>("srcTracks");
00163    if(producePixelTracks_) srcPixelTracks_ = iConfig.getParameter<edm::InputTag>("srcPixelTracks");
00164    
00165    reuseAny_ = !produceHFhits_ || !produceHFtowers_ || !produceBasicClusters_ || !produceEcalhits_ || !produceZDChits_;
00166    if(reuseAny_) reuseTag_ = iConfig.getParameter<edm::InputTag>("srcReUse");
00167 
00168    useQuality_   = iConfig.getParameter<bool>("UseQuality");
00169    trackQuality_ = TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));
00170 
00171    produces<reco::Centrality>();
00172    
00173 }
00174 
00175 
00176 CentralityProducer::~CentralityProducer()
00177 {
00178   
00179   // do anything here that needs to be done at desctruction time
00180   // (e.g. close files, deallocate resources etc.)
00181 
00182 }
00183 
00184 
00185 //
00186 // member functions
00187 //
00188 
00189 // ------------ method called to produce the data  ------------
00190 bool
00191 CentralityProducer::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
00192 {
00193   using namespace edm;
00194   using namespace reco;
00195 
00196   if(!trackGeo_ && doPixelCut_){
00197     edm::ESHandle<TrackerGeometry> tGeo;
00198     iSetup.get<TrackerDigiGeometryRecord>().get(tGeo);
00199     trackGeo_ = tGeo.product();
00200   }
00201   
00202   if(!caloGeo_ && 0){
00203     edm::ESHandle<CaloGeometry> cGeo;
00204     iSetup.get<CaloGeometryRecord>().get(cGeo);
00205     caloGeo_ = cGeo.product();
00206   }
00207 
00208   std::auto_ptr<Centrality> creco(new Centrality());
00209   Handle<Centrality> inputCentrality;
00210 
00211   if(reuseAny_) iEvent.getByLabel(reuseTag_,inputCentrality);
00212   
00213   if(produceHFhits_){
00214      creco->etHFhitSumPlus_ = 0;
00215      creco->etHFhitSumMinus_ = 0;
00216 
00217      Handle<HFRecHitCollection> hits;
00218      iEvent.getByLabel(srcHFhits_,hits);
00219      for( size_t ihit = 0; ihit<hits->size(); ++ ihit){
00220         const HFRecHit & rechit = (*hits)[ ihit ];
00221         if(rechit.id().ieta() > 0 )
00222            creco->etHFhitSumPlus_ += rechit.energy();
00223         if(rechit.id().ieta() < 0)
00224            creco->etHFhitSumMinus_ += rechit.energy();
00225      }       
00226   }else{
00227      creco->etHFhitSumMinus_ = inputCentrality->EtHFhitSumMinus();
00228      creco->etHFhitSumPlus_ = inputCentrality->EtHFhitSumPlus();
00229   }
00230   
00231   if(produceHFtowers_ || produceETmidRap_){
00232      creco->etHFtowerSumPlus_ = 0;
00233      creco->etHFtowerSumMinus_ = 0;
00234      creco->etMidRapiditySum_ = 0;
00235      
00236      Handle<CaloTowerCollection> towers;
00237 
00238      iEvent.getByLabel(srcTowers_,towers);
00239 
00240         for( size_t i = 0; i<towers->size(); ++ i){
00241            const CaloTower & tower = (*towers)[ i ];
00242            double eta = tower.eta();
00243            bool isHF = tower.ietaAbs() > 29;
00244            if(produceHFtowers_){
00245               if(isHF && eta > 0){
00246                  creco->etHFtowerSumPlus_ += tower.pt();
00247               }
00248               if(isHF && eta < 0){
00249                  creco->etHFtowerSumMinus_ += tower.pt();
00250               }
00251            }else{
00252               creco->etHFtowerSumMinus_ = inputCentrality->EtHFtowerSumMinus();
00253               creco->etHFtowerSumPlus_ = inputCentrality->EtHFtowerSumPlus();
00254            }
00255            if(produceETmidRap_){
00256               if(fabs(eta) < midRapidityRange_) creco->etMidRapiditySum_ += tower.pt()/(midRapidityRange_*2.);
00257            }else creco->etMidRapiditySum_ = inputCentrality->EtMidRapiditySum();
00258         }
00259   }else{
00260      creco->etHFtowerSumMinus_ = inputCentrality->EtHFtowerSumMinus();
00261      creco->etHFtowerSumPlus_ = inputCentrality->EtHFtowerSumPlus();
00262      creco->etMidRapiditySum_ = inputCentrality->EtMidRapiditySum();
00263   }
00264   
00265   if(produceBasicClusters_){
00266      creco->etEESumPlus_ = 0;
00267      creco->etEESumMinus_ = 0;
00268      creco->etEBSum_ = 0;
00269      
00270      Handle<BasicClusterCollection> clusters;
00271      iEvent.getByLabel(srcBasicClustersEE_, clusters);
00272      for( size_t i = 0; i<clusters->size(); ++ i){
00273         const BasicCluster & cluster = (*clusters)[ i ];
00274         double eta = cluster.eta();
00275         double tg = cluster.position().rho()/cluster.position().r();
00276         double et = cluster.energy()*tg;
00277         if(eta > 0)
00278            creco->etEESumPlus_ += et;
00279         if(eta < 0)
00280            creco->etEESumMinus_ += et;
00281      }
00282      
00283      iEvent.getByLabel(srcBasicClustersEB_, clusters);
00284      for( size_t i = 0; i<clusters->size(); ++ i){
00285         const BasicCluster & cluster = (*clusters)[ i ];
00286         double tg = cluster.position().rho()/cluster.position().r();
00287         double et = cluster.energy()*tg;
00288         creco->etEBSum_ += et;
00289      }
00290   }else{
00291      creco->etEESumMinus_ = inputCentrality->EtEESumMinus();
00292      creco->etEESumPlus_ = inputCentrality->EtEESumPlus();
00293      creco->etEBSum_ = inputCentrality->EtEBSum();
00294   }
00295   
00296   if(producePixelhits_){
00297      creco->pixelMultiplicity_ = 0;
00298      const SiPixelRecHitCollection* rechits;
00299      Handle<SiPixelRecHitCollection> rchts;
00300      iEvent.getByLabel(srcPixelhits_,rchts);
00301      rechits = rchts.product();
00302      int nPixel =0 ;
00303      for (SiPixelRecHitCollection::const_iterator it = rechits->begin(); it!=rechits->end();it++)
00304      {
00305         SiPixelRecHitCollection::DetSet hits = *it;
00306         DetId detId = DetId(hits.detId());
00307         SiPixelRecHitCollection::const_iterator recHitMatch = rechits->find(detId);
00308         const SiPixelRecHitCollection::DetSet recHitRange = *recHitMatch;
00309         for ( SiPixelRecHitCollection::DetSet::const_iterator recHitIterator = recHitRange.begin(); 
00310               recHitIterator != recHitRange.end(); ++recHitIterator) {
00311           // add selection if needed, now all hits.
00312           if(doPixelCut_){
00313             const SiPixelRecHit * recHit = &(*recHitIterator);
00314             const PixelGeomDetUnit* pixelLayer = dynamic_cast<const PixelGeomDetUnit*> (trackGeo_->idToDet(recHit->geographicalId()));
00315             GlobalPoint gpos = pixelLayer->toGlobal(recHit->localPosition());
00316             math::XYZVector rechitPos(gpos.x(),gpos.y(),gpos.z());
00317             double abeta = fabs(rechitPos.eta());
00318             int clusterSize = recHit->cluster()->size();
00319             if (                abeta < 0.5 && clusterSize < 1) continue;
00320             if ( abeta > 0.5 && abeta < 1   && clusterSize < 2) continue;
00321             if ( abeta > 1.  && abeta < 1.5 && clusterSize < 3) continue;
00322             if ( abeta > 1.5 && abeta < 2.  && clusterSize < 4) continue;
00323             if ( abeta > 2.  && abeta < 2.5 && clusterSize < 6) continue;
00324             if ( abeta > 2.5 && abeta < 5   && clusterSize < 9) continue;
00325           }
00326           nPixel++;
00327         } 
00328      }
00329      creco->pixelMultiplicity_ = nPixel;
00330   }else{
00331      creco->pixelMultiplicity_ = inputCentrality->multiplicityPixel();
00332   }
00333 
00334   if(produceTracks_){
00335      edm::Handle<reco::TrackCollection> tracks;
00336      iEvent.getByLabel(srcTracks_,tracks);
00337      int nTracks = 0;
00338 
00339      double trackCounter = 0;
00340      double trackCounterEta = 0;
00341      double trackCounterEtaPt = 0;
00342 
00343      for(unsigned int i = 0 ; i < tracks->size(); ++i){
00344        const Track& track = (*tracks)[i];
00345        if(useQuality_ && !track.quality(trackQuality_)) continue;
00346        nTracks++;
00347 
00348        if( track.pt() > trackPtCut_)  trackCounter++;
00349        if(fabs(track.eta())<trackEtaCut_) {
00350          trackCounterEta++;
00351          if (track.pt() > trackPtCut_) trackCounterEtaPt++;
00352        }
00353      }
00354 
00355      creco->trackMultiplicity_ = nTracks;
00356      creco->ntracksPtCut_ = trackCounter; 
00357      creco->ntracksEtaCut_ = trackCounterEta;
00358      creco->ntracksEtaPtCut_ = trackCounterEtaPt;
00359 
00360   }else{
00361      creco->trackMultiplicity_ = inputCentrality->Ntracks();
00362      creco->ntracksPtCut_= inputCentrality->NtracksPtCut();
00363      creco->ntracksEtaCut_= inputCentrality->NtracksEtaCut();
00364      creco->ntracksEtaPtCut_= inputCentrality->NtracksEtaPtCut();
00365   }
00366 
00367   if(producePixelTracks_){
00368     edm::Handle<reco::TrackCollection> pixeltracks;
00369     iEvent.getByLabel(srcPixelTracks_,pixeltracks);
00370     int nPixelTracks = pixeltracks->size();
00371     creco->nPixelTracks_ = nPixelTracks;
00372   }
00373   else{
00374     creco->nPixelTracks_ = inputCentrality->NpixelTracks();
00375   }
00376 
00377   if(produceZDChits_){
00378      creco->zdcSumPlus_ = 0;
00379      creco->zdcSumMinus_ = 0;
00380      
00381      Handle<ZDCRecHitCollection> hits;
00382      bool zdcAvailable = iEvent.getByLabel(srcZDChits_,hits);
00383      if(zdcAvailable){
00384         for( size_t ihit = 0; ihit<hits->size(); ++ ihit){
00385            const ZDCRecHit & rechit = (*hits)[ ihit ];
00386            if(rechit.id().zside() > 0 )
00387               creco->zdcSumPlus_ += rechit.energy();
00388            if(rechit.id().zside() < 0)
00389               creco->zdcSumMinus_ += rechit.energy();
00390         }
00391      }else{
00392         creco->zdcSumPlus_ = -9;
00393         creco->zdcSumMinus_ = -9;
00394      }
00395   }else{
00396      creco->zdcSumMinus_ = inputCentrality->zdcSumMinus();
00397      creco->zdcSumPlus_ = inputCentrality->zdcSumPlus();
00398   }
00399   
00400   iEvent.put(creco);
00401   return true;
00402 }
00403 
00404 // ------------ method called once each job just before starting event loop  ------------
00405 void 
00406 CentralityProducer::beginJob()
00407 {
00408 }
00409 
00410 // ------------ method called once each job just after ending the event loop  ------------
00411 void 
00412 CentralityProducer::endJob() {
00413 }
00414 }
00415 
00416 //define this as a plug-in
00417 DEFINE_FWK_MODULE(reco::CentralityProducer);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoHI/HiCentralityAlgos/src/CentralityProducer.cc
