HiCentralityAlgos/plugins/CentralityProducer.cc

0001 //
0002 // Original Author:  Yetkin Yilmaz, Young Soo Park
0003 //         Created:  Wed Jun 11 15:31:41 CEST 2008
0004 //
0005 //
0006
0007 // system include files
0008 #include <memory>
0009 #include <iostream>
0010
0011 // user include files
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/global/EDProducer.h"
0014 #include "FWCore/Utilities/interface/StreamID.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019 #include "FWCore/Framework/interface/ESHandle.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022 #include "FWCore/Utilities/interface/ESGetToken.h"
0023
0024 #include "DataFormats/Candidate/interface/Candidate.h"
0025
0026 #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
0027 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0028 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0029
0030 #include "DataFormats/Common/interface/Ref.h"
0031 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0032 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0033 #include "DataFormats/TrackReco/interface/Track.h"
0034 #include "DataFormats/VertexReco/interface/Vertex.h"
0035 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0036 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0037
0038 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0039 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0041 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0042 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0043
0044 using namespace std;
0045 using namespace edm;
0046 using namespace reco;
0047
0048 //
0049 // class declaration
0050 //
0051
0052 namespace reco {
0053   class CentralityProducer : public edm::global::EDProducer<> {
0054   public:
0055     explicit CentralityProducer(const edm::ParameterSet&);
0056     ~CentralityProducer() override = default;
0057     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0058
0059   private:
0060     void beginJob() override;
0061     void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0062     void endJob() override;
0063
0064     // ----------member data ---------------------------
0065
0066     const bool produceHFhits_;
0067     const bool produceHFtowers_;
0068     const bool produceEcalhits_;
0069     const bool produceZDChits_;
0070     const bool lowGainZDC_;
0071     const bool produceETmidRap_;
0072     const bool producePixelhits_;
0073     const bool doPixelCut_;
0074     const bool produceTracks_;
0075     const bool producePixelTracks_;
0076     const bool producePF_;
0077
0078     const double midRapidityRange_;
0079     const double trackPtCut_;
0080     const double trackEtaCut_;
0081     const double hfEtaCut_;
0082
0083     const bool reuseAny_;
0084     const bool useQuality_;
0085     const reco::TrackBase::TrackQuality trackQuality_;
0086
0087     const edm::EDGetTokenT<HFRecHitCollection> srcHFhits_;
0088     const edm::EDGetTokenT<CaloTowerCollection> srcTowers_;
0089     const edm::EDGetTokenT<EcalRecHitCollection> srcEBhits_;
0090     const edm::EDGetTokenT<EcalRecHitCollection> srcEEhits_;
0091     const edm::EDGetTokenT<ZDCRecHitCollection> srcZDChits_;
0092     const edm::EDGetTokenT<SiPixelRecHitCollection> srcPixelhits_;
0093     const edm::EDGetTokenT<TrackCollection> srcTracks_;
0094     const edm::EDGetTokenT<TrackCollection> srcPixelTracks_;
0095     const edm::EDGetTokenT<reco::CandidateView> srcPF_;
0096     const edm::EDGetTokenT<VertexCollection> srcVertex_;
0097     const edm::EDGetTokenT<Centrality> reuseTag_;
0098
0099     const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeom_;
0100     const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeom_;
0101     const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopo_;
0102   };
0103
0104   //
0105   // constants, enums and typedefs
0106   //
0107
0108   //
0109   // static data member definitions
0110   //
0111
0112   //
0113   // constructors and destructor
0114   //
0115   CentralityProducer::CentralityProducer(const edm::ParameterSet& iConfig)
0116       : produceHFhits_(iConfig.getParameter<bool>("produceHFhits")),
0117         produceHFtowers_(iConfig.getParameter<bool>("produceHFtowers")),
0118         produceEcalhits_(iConfig.getParameter<bool>("produceEcalhits")),
0119         produceZDChits_(iConfig.getParameter<bool>("produceZDChits")),
0120         lowGainZDC_(iConfig.getParameter<bool>("lowGainZDC")),
0121         produceETmidRap_(iConfig.getParameter<bool>("produceETmidRapidity")),
0122         producePixelhits_(iConfig.getParameter<bool>("producePixelhits")),
0123         doPixelCut_(iConfig.getParameter<bool>("doPixelCut")),
0124         produceTracks_(iConfig.getParameter<bool>("produceTracks")),
0125         producePixelTracks_(iConfig.getParameter<bool>("producePixelTracks")),
0126         producePF_(iConfig.getParameter<bool>("producePF")),
0127         midRapidityRange_(iConfig.getParameter<double>("midRapidityRange")),
0128         trackPtCut_(iConfig.getParameter<double>("trackPtCut")),
0129         trackEtaCut_(iConfig.getParameter<double>("trackEtaCut")),
0130         hfEtaCut_(iConfig.getParameter<double>("hfEtaCut")),
0131         reuseAny_(iConfig.getParameter<bool>("reUseCentrality")),
0132         useQuality_(iConfig.getParameter<bool>("useQuality")),
0133         trackQuality_(TrackBase::qualityByName(iConfig.getParameter<std::string>("trackQuality"))),
0134         srcHFhits_(produceHFhits_ ? consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("srcHFhits"))
0135                                   : edm::EDGetTokenT<HFRecHitCollection>()),
0136         srcTowers_((produceHFtowers_ || produceETmidRap_)
0137                        ? consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("srcTowers"))
0138                        : edm::EDGetTokenT<CaloTowerCollection>()),
0139         srcEBhits_(produceEcalhits_ ? consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("srcEBhits"))
0140                                     : edm::EDGetTokenT<EcalRecHitCollection>()),
0141         srcEEhits_(produceEcalhits_ ? consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("srcEEhits"))
0142                                     : edm::EDGetTokenT<EcalRecHitCollection>()),
0143         srcZDChits_(produceZDChits_ ? consumes<ZDCRecHitCollection>(iConfig.getParameter<edm::InputTag>("srcZDChits"))
0144                                     : edm::EDGetTokenT<ZDCRecHitCollection>()),
0145         srcPixelhits_(producePixelhits_
0146                           ? consumes<SiPixelRecHitCollection>(iConfig.getParameter<edm::InputTag>("srcPixelhits"))
0147                           : edm::EDGetTokenT<SiPixelRecHitCollection>()),
0148         srcTracks_(produceTracks_ ? consumes<TrackCollection>(iConfig.getParameter<edm::InputTag>("srcTracks"))
0149                                   : edm::EDGetTokenT<TrackCollection>()),
0150         srcPixelTracks_(producePixelTracks_
0151                             ? consumes<TrackCollection>(iConfig.getParameter<edm::InputTag>("srcPixelTracks"))
0152                             : edm::EDGetTokenT<TrackCollection>()),
0153         srcPF_(producePF_ ? consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("srcPF"))
0154                           : edm::EDGetTokenT<reco::CandidateView>()),
0155         srcVertex_((produceTracks_ || producePixelTracks_)
0156                        ? consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("srcVertex"))
0157                        : edm::EDGetTokenT<VertexCollection>()),
0158         reuseTag_(reuseAny_ ? consumes<Centrality>(iConfig.getParameter<edm::InputTag>("srcReUse"))
0159                             : edm::EDGetTokenT<Centrality>()),
0160         caloGeom_(produceEcalhits_ ? esConsumes<CaloGeometry, CaloGeometryRecord>()
0161                                    : edm::ESGetToken<CaloGeometry, CaloGeometryRecord>()),
0162         trackerGeom_(producePixelhits_ ? esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()
0163                                        : edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord>()),
0164         trackerTopo_(producePixelhits_ ? esConsumes<TrackerTopology, TrackerTopologyRcd>()
0165                                        : edm::ESGetToken<TrackerTopology, TrackerTopologyRcd>()) {
0166     produces<Centrality>();
0167   }
0168
0169   //
0170   // member functions
0171   //
0172
0173   // ------------ method called to produce the data  ------------
0174   void CentralityProducer::produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0175     auto creco = std::make_unique<Centrality>();
0176     Handle<Centrality> inputCentrality;
0177     if (reuseAny_)
0178       iEvent.getByToken(reuseTag_, inputCentrality);
0179
0180     if (produceHFhits_) {
0181       creco->etHFhitSumPlus_ = 0;
0182       creco->etHFhitSumMinus_ = 0;
0183       Handle<HFRecHitCollection> hits;
0184       iEvent.getByToken(srcHFhits_, hits);
0185       for (size_t ihit = 0; ihit < hits->size(); ++ihit) {
0186         const HFRecHit& rechit = (*hits)[ihit];
0187         if (rechit.id().ieta() > 0)
0188           creco->etHFhitSumPlus_ += rechit.energy();
0189         if (rechit.id().ieta() < 0)
0190           creco->etHFhitSumMinus_ += rechit.energy();
0191       }
0192     } else {
0193       if (reuseAny_) {
0194         creco->etHFhitSumMinus_ = inputCentrality->EtHFhitSumMinus();
0195         creco->etHFhitSumPlus_ = inputCentrality->EtHFhitSumPlus();
0196       }
0197     }
0198
0199     if (produceHFtowers_ || produceETmidRap_) {
0200       creco->etHFtowerSumPlus_ = 0;
0201       creco->etHFtowerSumMinus_ = 0;
0202       creco->etHFtowerSumECutPlus_ = 0;
0203       creco->etHFtowerSumECutMinus_ = 0;
0204       creco->etMidRapiditySum_ = 0;
0205
0206       Handle<CaloTowerCollection> towers;
0207       iEvent.getByToken(srcTowers_, towers);
0208
0209       for (size_t i = 0; i < towers->size(); ++i) {
0210         const CaloTower& tower = (*towers)[i];
0211         double eta = tower.eta();
0212         if (produceHFtowers_) {
0213           bool isHF = tower.ietaAbs() > 29;
0214           if (isHF && eta > 0) {
0215             creco->etHFtowerSumPlus_ += tower.pt();
0216             if (tower.energy() > 1.5)
0217               creco->etHFtowerSumECutPlus_ += tower.pt();
0218             if (eta > hfEtaCut_)
0219               creco->etHFtruncatedPlus_ += tower.pt();
0220           } else if (isHF && eta < 0) {
0221             creco->etHFtowerSumMinus_ += tower.pt();
0222             if (tower.energy() > 1.5)
0223               creco->etHFtowerSumECutMinus_ += tower.pt();
0224             if (eta < -hfEtaCut_)
0225               creco->etHFtruncatedMinus_ += tower.pt();
0226           }
0227         } else {
0228           if (reuseAny_) {
0229             creco->etHFtowerSumMinus_ = inputCentrality->EtHFtowerSumMinus();
0230             creco->etHFtowerSumPlus_ = inputCentrality->EtHFtowerSumPlus();
0231             creco->etHFtowerSumECutMinus_ = inputCentrality->EtHFtowerSumECutMinus();
0232             creco->etHFtowerSumECutPlus_ = inputCentrality->EtHFtowerSumECutPlus();
0233             creco->etHFtruncatedMinus_ = inputCentrality->EtHFtruncatedMinus();
0234             creco->etHFtruncatedPlus_ = inputCentrality->EtHFtruncatedPlus();
0235           }
0236         }
0237         if (produceETmidRap_) {
0238           if (std::abs(eta) < midRapidityRange_)
0239             creco->etMidRapiditySum_ += tower.pt() / (midRapidityRange_ * 2.);
0240         } else if (reuseAny_)
0241           creco->etMidRapiditySum_ = inputCentrality->EtMidRapiditySum();
0242       }
0243     } else {
0244       if (reuseAny_) {
0245         creco->etHFtowerSumMinus_ = inputCentrality->EtHFtowerSumMinus();
0246         creco->etHFtowerSumPlus_ = inputCentrality->EtHFtowerSumPlus();
0247         creco->etHFtowerSumECutMinus_ = inputCentrality->EtHFtowerSumECutMinus();
0248         creco->etHFtowerSumECutPlus_ = inputCentrality->EtHFtowerSumECutPlus();
0249         creco->etMidRapiditySum_ = inputCentrality->EtMidRapiditySum();
0250       }
0251     }
0252
0253     if (producePF_) {
0254       creco->etPFhfSumPlus_ = 0;
0255       creco->etPFhfSumMinus_ = 0;
0256       for (const auto& pf : iEvent.get(srcPF_)) {
0257         if (pf.pdgId() != 1 && pf.pdgId() != 2)
0258           continue;
0259         if (pf.eta() > 0)
0260           creco->etPFhfSumPlus_ += pf.pt();
0261         else
0262           creco->etPFhfSumMinus_ += pf.pt();
0263       }
0264     } else if (reuseAny_) {
0265       creco->etPFhfSumMinus_ = inputCentrality->EtPFhfSumMinus();
0266       creco->etPFhfSumPlus_ = inputCentrality->EtPFhfSumPlus();
0267     }
0268
0269     if (produceEcalhits_) {
0270       creco->etEESumPlus_ = 0;
0271       creco->etEESumMinus_ = 0;
0272       creco->etEBSum_ = 0;
0273
0274       Handle<EcalRecHitCollection> ebHits;
0275       Handle<EcalRecHitCollection> eeHits;
0276       iEvent.getByToken(srcEBhits_, ebHits);
0277       iEvent.getByToken(srcEEhits_, eeHits);
0278
0279       edm::ESHandle<CaloGeometry> cGeo = iSetup.getHandle(caloGeom_);
0280       for (unsigned int i = 0; i < ebHits->size(); ++i) {
0281         const EcalRecHit& hit = (*ebHits)[i];
0282         const GlobalPoint& pos = cGeo->getPosition(hit.id());
0283         double et = hit.energy() * (pos.perp() / pos.mag());
0284         creco->etEBSum_ += et;
0285       }
0286
0287       for (unsigned int i = 0; i < eeHits->size(); ++i) {
0288         const EcalRecHit& hit = (*eeHits)[i];
0289         const GlobalPoint& pos = cGeo->getPosition(hit.id());
0290         double et = hit.energy() * (pos.perp() / pos.mag());
0291         if (pos.z() > 0) {
0292           creco->etEESumPlus_ += et;
0293         } else {
0294           creco->etEESumMinus_ += et;
0295         }
0296       }
0297     } else {
0298       if (reuseAny_) {
0299         creco->etEESumMinus_ = inputCentrality->EtEESumMinus();
0300         creco->etEESumPlus_ = inputCentrality->EtEESumPlus();
0301         creco->etEBSum_ = inputCentrality->EtEBSum();
0302       }
0303     }
0304
0305     if (producePixelhits_) {
0306       edm::ESHandle<TrackerGeometry> tGeo = iSetup.getHandle(trackerGeom_);
0307       edm::ESHandle<TrackerTopology> topo = iSetup.getHandle(trackerTopo_);
0308       creco->pixelMultiplicity_ = 0;
0309       const SiPixelRecHitCollection* rechits;
0310       Handle<SiPixelRecHitCollection> rchts;
0311       iEvent.getByToken(srcPixelhits_, rchts);
0312       rechits = rchts.product();
0313       int nPixel = 0;
0314       int nPixel_plus = 0;
0315       int nPixel_minus = 0;
0316       for (SiPixelRecHitCollection::const_iterator it = rechits->begin(); it != rechits->end(); it++) {
0317         SiPixelRecHitCollection::DetSet hits = *it;
0318         DetId detId = DetId(hits.detId());
0319         SiPixelRecHitCollection::const_iterator recHitMatch = rechits->find(detId);
0320         const SiPixelRecHitCollection::DetSet recHitRange = *recHitMatch;
0321         for (SiPixelRecHitCollection::DetSet::const_iterator recHitIterator = recHitRange.begin();
0322              recHitIterator != recHitRange.end();
0323              ++recHitIterator) {
0324           // add selection if needed, now all hits.
0325           const SiPixelRecHit* recHit = &(*recHitIterator);
0326           const PixelGeomDetUnit* pixelLayer =
0327               dynamic_cast<const PixelGeomDetUnit*>(tGeo->idToDet(recHit->geographicalId()));
0328           GlobalPoint gpos = pixelLayer->toGlobal(recHit->localPosition());
0329           math::XYZVector rechitPos(gpos.x(), gpos.y(), gpos.z());
0330           double eta = rechitPos.eta();
0331           int clusterSize = recHit->cluster()->size();
0332           unsigned layer = topo->layer(detId);
0333           if (doPixelCut_) {
0334             if (detId.det() == DetId::Tracker && detId.subdetId() == PixelSubdetector::PixelBarrel) {
0335               double abeta = std::abs(eta);
0336               if (layer == 1) {
0337                 if (18 * abeta - 40 > clusterSize)
0338                   continue;
0339               } else if (layer == 2) {
0340                 if (6 * abeta - 7.2 > clusterSize)
0341                   continue;
0342               } else if (layer == 3 || layer == 4) {
0343                 if (4 * abeta - 2.4 > clusterSize)
0344                   continue;
0345               }
0346             }
0347           }
0348           nPixel++;
0349           if (eta >= 0)
0350             nPixel_plus++;
0351           else if (eta < 0)
0352             nPixel_minus++;
0353         }
0354       }
0355       creco->pixelMultiplicity_ = nPixel;
0356       creco->pixelMultiplicityPlus_ = nPixel_plus;
0357       creco->pixelMultiplicityMinus_ = nPixel_minus;
0358     } else {
0359       if (reuseAny_) {
0360         creco->pixelMultiplicity_ = inputCentrality->multiplicityPixel();
0361         creco->pixelMultiplicityPlus_ = inputCentrality->multiplicityPixelPlus();
0362         creco->pixelMultiplicityMinus_ = inputCentrality->multiplicityPixelMinus();
0363       }
0364     }
0365
0366     if (produceTracks_) {
0367       double vx = -999.;
0368       double vy = -999.;
0369       double vz = -999.;
0370       double vxError = -999.;
0371       double vyError = -999.;
0372       double vzError = -999.;
0373       math::XYZVector vtxPos(0, 0, 0);
0374
0375       Handle<VertexCollection> recoVertices;
0376       iEvent.getByToken(srcVertex_, recoVertices);
0377       unsigned int daughter = 0;
0378       int greatestvtx = 0;
0379
0380       for (unsigned int i = 0; i < recoVertices->size(); ++i) {
0381         daughter = (*recoVertices)[i].tracksSize();
0382         if (daughter > (*recoVertices)[greatestvtx].tracksSize())
0383           greatestvtx = i;
0384       }
0385
0386       if (!recoVertices->empty()) {
0387         vx = (*recoVertices)[greatestvtx].position().x();
0388         vy = (*recoVertices)[greatestvtx].position().y();
0389         vz = (*recoVertices)[greatestvtx].position().z();
0390         vxError = (*recoVertices)[greatestvtx].xError();
0391         vyError = (*recoVertices)[greatestvtx].yError();
0392         vzError = (*recoVertices)[greatestvtx].zError();
0393       }
0394
0395       vtxPos = math::XYZVector(vx, vy, vz);
0396
0397       Handle<TrackCollection> tracks;
0398       iEvent.getByToken(srcTracks_, tracks);
0399       int nTracks = 0;
0400
0401       double trackCounter = 0;
0402       double trackCounterEta = 0;
0403       double trackCounterEtaPt = 0;
0404
0405       for (unsigned int i = 0; i < tracks->size(); ++i) {
0406         const Track& track = (*tracks)[i];
0407         if (useQuality_ && !track.quality(trackQuality_))
0408           continue;
0409
0410         if (track.pt() > trackPtCut_)
0411           trackCounter++;
0412         if (std::abs(track.eta()) < trackEtaCut_) {
0413           trackCounterEta++;
0414           if (track.pt() > trackPtCut_)
0415             trackCounterEtaPt++;
0416         }
0417
0418         math::XYZPoint v1(vx, vy, vz);
0419         double dz = track.dz(v1);
0420         double dzsigma2 = track.dzError() * track.dzError() + vzError * vzError;
0421         double dxy = track.dxy(v1);
0422         double dxysigma2 = track.dxyError() * track.dxyError() + vxError * vyError;
0423
0424         const double pterrcut = 0.1;
0425         const double dzrelcut = 3.0;
0426         const double dxyrelcut = 3.0;
0427
0428         if (track.quality(trackQuality_) && track.pt() > 0.4 && std::abs(track.eta()) < 2.4 &&
0429             track.ptError() / track.pt() < pterrcut && dz * dz < dzrelcut * dzrelcut * dzsigma2 &&
0430             dxy * dxy < dxyrelcut * dxyrelcut * dxysigma2) {
0431           nTracks++;
0432         }
0433       }
0434
0435       creco->trackMultiplicity_ = nTracks;
0436       creco->ntracksPtCut_ = trackCounter;
0437       creco->ntracksEtaCut_ = trackCounterEta;
0438       creco->ntracksEtaPtCut_ = trackCounterEtaPt;
0439
0440     } else {
0441       if (reuseAny_) {
0442         creco->trackMultiplicity_ = inputCentrality->Ntracks();
0443         creco->ntracksPtCut_ = inputCentrality->NtracksPtCut();
0444         creco->ntracksEtaCut_ = inputCentrality->NtracksEtaCut();
0445         creco->ntracksEtaPtCut_ = inputCentrality->NtracksEtaPtCut();
0446       }
0447     }
0448
0449     if (producePixelTracks_) {
0450       Handle<TrackCollection> pixeltracks;
0451       iEvent.getByToken(srcPixelTracks_, pixeltracks);
0452       int nPixelTracks = pixeltracks->size();
0453       int nPixelTracksPlus = 0;
0454       int nPixelTracksMinus = 0;
0455
0456       for (auto const& track : *pixeltracks) {
0457         if (track.eta() < 0)
0458           nPixelTracksMinus++;
0459         else
0460           nPixelTracksPlus++;
0461       }
0462       creco->nPixelTracks_ = nPixelTracks;
0463       creco->nPixelTracksPlus_ = nPixelTracksPlus;
0464       creco->nPixelTracksMinus_ = nPixelTracksMinus;
0465     } else {
0466       if (reuseAny_) {
0467         creco->nPixelTracks_ = inputCentrality->NpixelTracks();
0468         creco->nPixelTracksPlus_ = inputCentrality->NpixelTracksPlus();
0469         creco->nPixelTracksMinus_ = inputCentrality->NpixelTracksMinus();
0470       }
0471     }
0472
0473     if (produceZDChits_) {
0474       creco->zdcSumPlus_ = 0;
0475       creco->zdcSumMinus_ = 0;
0476
0477       Handle<ZDCRecHitCollection> hits;
0478       bool zdcAvailable = iEvent.getByToken(srcZDChits_, hits);
0479       if (zdcAvailable) {
0480         for (size_t ihit = 0; ihit < hits->size(); ++ihit) {
0481           const ZDCRecHit& rechit = (*hits)[ihit];
0482           if (rechit.id().zside() > 0) {
0483             if (lowGainZDC_) {
0484               creco->zdcSumPlus_ += rechit.lowGainEnergy();
0485             } else {
0486               creco->zdcSumPlus_ += rechit.energy();
0487             }
0488           }
0489           if (rechit.id().zside() < 0) {
0490             if (lowGainZDC_) {
0491               creco->zdcSumMinus_ += rechit.lowGainEnergy();
0492             } else {
0493               creco->zdcSumMinus_ += rechit.energy();
0494             }
0495           }
0496         }
0497       } else {
0498         creco->zdcSumPlus_ = -9;
0499         creco->zdcSumMinus_ = -9;
0500       }
0501     } else {
0502       if (reuseAny_) {
0503         creco->zdcSumMinus_ = inputCentrality->zdcSumMinus();
0504         creco->zdcSumPlus_ = inputCentrality->zdcSumPlus();
0505       }
0506     }
0507
0508     iEvent.put(std::move(creco));
0509   }
0510
0511   void CentralityProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0512     edm::ParameterSetDescription desc;
0513     desc.add<bool>("produceHFhits", true);
0514     desc.add<bool>("produceHFtowers", true);
0515     desc.add<bool>("produceEcalhits", true);
0516     desc.add<bool>("produceZDChits", true);
0517     desc.add<bool>("produceETmidRapidity", true);
0518     desc.add<bool>("producePixelhits", true);
0519     desc.add<bool>("produceTracks", true);
0520     desc.add<bool>("producePixelTracks", true);
0521     desc.add<bool>("producePF", true);
0522     desc.add<bool>("reUseCentrality", false);
0523     desc.add<edm::InputTag>("srcHFhits", edm::InputTag("hfreco"));
0524     desc.add<edm::InputTag>("srcTowers", edm::InputTag("towerMaker"));
0525     desc.add<edm::InputTag>("srcEBhits", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0526     desc.add<edm::InputTag>("srcEEhits", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0527     desc.add<edm::InputTag>("srcZDChits", edm::InputTag("zdcreco"));
0528     desc.add<edm::InputTag>("srcPixelhits", edm::InputTag("siPixelRecHits"));
0529     desc.add<edm::InputTag>("srcTracks", edm::InputTag("hiGeneralTracks"));
0530     desc.add<edm::InputTag>("srcVertex", edm::InputTag("hiSelectedVertex"));
0531     desc.add<edm::InputTag>("srcReUse", edm::InputTag("hiCentrality"));
0532     desc.add<edm::InputTag>("srcPixelTracks", edm::InputTag("hiPixel3PrimTracks"));
0533     desc.add<edm::InputTag>("srcPF", edm::InputTag("particleFlow"));
0534     desc.add<bool>("doPixelCut", true);
0535     desc.add<bool>("useQuality", true);
0536     desc.add<string>("trackQuality", "highPurity");
0537     desc.add<double>("trackEtaCut", 2);
0538     desc.add<double>("trackPtCut", 1);
0539     desc.add<double>("hfEtaCut", 4)->setComment("hf above the absolute value of this cut is used");
0540     desc.add<double>("midRapidityRange", 1);
0541     desc.add<bool>("lowGainZDC", true);
0542
0543     descriptions.addDefault(desc);
0544   }
0545
0546   // ------------ method called once each job just before starting event loop  ------------
0547   void CentralityProducer::beginJob() {}
0548
0549   // ------------ method called once each job just after ending the event loop  ------------
0550   void CentralityProducer::endJob() {}
0551
0552   //define this as a plug-in
0553   DEFINE_FWK_MODULE(CentralityProducer);
0554
0555 }  // namespace reco