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