diff --git a/DQM/Physics/src/CentralityDQM.cc b/DQM/Physics/src/CentralityDQM.cc index e626cbcee5048..74f280b71d2de 100644 --- a/DQM/Physics/src/CentralityDQM.cc +++ b/DQM/Physics/src/CentralityDQM.cc @@ -84,6 +84,9 @@ void CentralityDQM::bookHistograms(DQMStore::IBooker& bei, edm::Run const&, edm: h_hiZDC = bei.book1D("h_hiZDC", "h_hiZDC", 600, 0, 6000); h_hiZDCplus = bei.book1D("h_hiZDCplus", "h_hiZDCplus", 600, 0, 6000); h_hiZDCminus = bei.book1D("h_hiZDCminus", "h_hiZDCminus", 600, 0, 6000); + h_hiPF = bei.book1D("h_hiPF", "h_hiPF", 900, 0, 9000); + h_hiPFplus = bei.book1D("h_hiPFplus", "h_hiPFplus", 900, 0, 9000); + h_hiPFminus = bei.book1D("h_hiPFminus", "h_hiPFminus", 900, 0, 9000); h_vertex_x = bei.book1D("h_vertex_x", "h_vertex_x", 400, -4, 4); h_vertex_y = bei.book1D("h_vertex_y", "h_vertex_y", 400, -4, 4); @@ -155,6 +158,10 @@ void CentralityDQM::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe h_hiEB->Fill(cent->EtEBSum()); h_hiET->Fill(cent->EtMidRapiditySum()); + h_hiPF->Fill(cent->EtPFhfSum()); + h_hiPFplus->Fill(cent->EtPFhfSumPlus()); + h_hiPFminus->Fill(cent->EtPFhfSumMinus()); + edm::Handle > vertex; iEvent.getByToken(vertexToken, vertex); h_vertex_x->Fill(vertex->begin()->x()); diff --git a/DQM/Physics/src/CentralityDQM.h b/DQM/Physics/src/CentralityDQM.h index 02a4c90a6dd03..82563cab7429c 100644 --- a/DQM/Physics/src/CentralityDQM.h +++ b/DQM/Physics/src/CentralityDQM.h @@ -76,6 +76,9 @@ class CentralityDQM : public DQMEDAnalyzer { MonitorElement* h_hiZDC; MonitorElement* h_hiZDCplus; MonitorElement* h_hiZDCminus; + MonitorElement* h_hiPF; + MonitorElement* h_hiPFplus; + MonitorElement* h_hiPFminus; MonitorElement* h_vertex_x; MonitorElement* h_vertex_y; diff --git a/DQM/Physics/src/CentralitypADQM.cc b/DQM/Physics/src/CentralitypADQM.cc index 7524faab59868..4d8aa011fae81 100644 --- a/DQM/Physics/src/CentralitypADQM.cc +++ b/DQM/Physics/src/CentralitypADQM.cc @@ -77,6 +77,9 @@ void CentralitypADQM::bookHistograms(DQMStore::IBooker& bei, edm::Run const&, ed h_hiZDC = bei.book1D("h_hiZDC", "h_hiZDC", 600, 0, 6000); h_hiZDCplus = bei.book1D("h_hiZDCplus", "h_hiZDCplus", 600, 0, 6000); h_hiZDCminus = bei.book1D("h_hiZDCminus", "h_hiZDCminus", 600, 0, 6000); + h_hiPF = bei.book1D("h_hiPF", "h_hiPF", 900, 0, 9000); + h_hiPFplus = bei.book1D("h_hiPFplus", "h_hiPFplus", 900, 0, 9000); + h_hiPFminus = bei.book1D("h_hiPFminus", "h_hiPFminus", 900, 0, 9000); h_vertex_x = bei.book1D("h_vertex_x", "h_vertex_x", 400, -4, 4); h_vertex_y = bei.book1D("h_vertex_y", "h_vertex_y", 400, -4, 4); @@ -118,6 +121,10 @@ void CentralitypADQM::analyze(const edm::Event& iEvent, const edm::EventSetup& i h_hiEB->Fill(cent->EtEBSum()); h_hiET->Fill(cent->EtMidRapiditySum()); + h_hiPF->Fill(cent->EtPFhfSum()); + h_hiPFplus->Fill(cent->EtPFhfSumPlus()); + h_hiPFminus->Fill(cent->EtPFhfSumMinus()); + edm::Handle > vertex; iEvent.getByToken(vertexToken, vertex); h_vertex_x->Fill(vertex->begin()->x()); diff --git a/DQM/Physics/src/CentralitypADQM.h b/DQM/Physics/src/CentralitypADQM.h index 0dfd1fd7aba26..ea4842e68fa15 100644 --- a/DQM/Physics/src/CentralitypADQM.h +++ b/DQM/Physics/src/CentralitypADQM.h @@ -73,6 +73,9 @@ class CentralitypADQM : public DQMEDAnalyzer { MonitorElement* h_hiZDC; MonitorElement* h_hiZDCplus; MonitorElement* h_hiZDCminus; + MonitorElement* h_hiPF; + MonitorElement* h_hiPFplus; + MonitorElement* h_hiPFminus; MonitorElement* h_vertex_x; MonitorElement* h_vertex_y; diff --git a/DQMOffline/JetMET/src/JetAnalyzer_HeavyIons.cc b/DQMOffline/JetMET/src/JetAnalyzer_HeavyIons.cc index c5bd7cb970e6a..08d2805e1509c 100644 --- a/DQMOffline/JetMET/src/JetAnalyzer_HeavyIons.cc +++ b/DQMOffline/JetMET/src/JetAnalyzer_HeavyIons.cc @@ -642,6 +642,9 @@ void JetAnalyzer_HeavyIons::analyze(const edm::Event &mEvent, const edm::EventSe edm::Handle cent; mEvent.getByToken(centralityToken, cent); //_centralitytag comes from the cfg + if (!cent.isValid()) + return; + mHF->Fill(cent->EtHFtowerSum()); Float_t HF_energy = cent->EtHFtowerSum(); @@ -649,9 +652,6 @@ void JetAnalyzer_HeavyIons::analyze(const edm::Event &mEvent, const edm::EventSe //edm::Handle cbin; //mEvent.getByToken(centralityBinToken, cbin); - if (!cent.isValid()) - return; - /*int hibin = -999; if(cbin.isValid()){ hibin = *cbin; diff --git a/DataFormats/HeavyIonEvent/interface/Centrality.h b/DataFormats/HeavyIonEvent/interface/Centrality.h index c65f1b9fb377f..903cb72e26f66 100644 --- a/DataFormats/HeavyIonEvent/interface/Centrality.h +++ b/DataFormats/HeavyIonEvent/interface/Centrality.h @@ -55,6 +55,9 @@ namespace reco { double zdcSumPlus() const { return zdcSumPlus_; } double zdcSumMinus() const { return zdcSumMinus_; } double EtMidRapiditySum() const { return etMidRapiditySum_; } + double EtPFhfSum() const { return etPFhfSumPlus_ + etPFhfSumMinus_; } + double EtPFhfSumPlus() const { return etPFhfSumPlus_; } + double EtPFhfSumMinus() const { return etPFhfSumMinus_; } protected: double value_; @@ -84,6 +87,8 @@ namespace reco { double zdcSumPlus_; double zdcSumMinus_; double etMidRapiditySum_; + double etPFhfSumPlus_; + double etPFhfSumMinus_; double ntracksPtCut_; double ntracksEtaCut_; double ntracksEtaPtCut_; diff --git a/DataFormats/HeavyIonEvent/src/Centrality.cc b/DataFormats/HeavyIonEvent/src/Centrality.cc index ce64ad5ade745..15c167020aa4f 100644 --- a/DataFormats/HeavyIonEvent/src/Centrality.cc +++ b/DataFormats/HeavyIonEvent/src/Centrality.cc @@ -31,6 +31,8 @@ Centrality::Centrality(double d, std::string label) zdcSumPlus_(0), zdcSumMinus_(0), etMidRapiditySum_(0), + etPFhfSumPlus_(0), + etPFhfSumMinus_(0), ntracksPtCut_(0), ntracksEtaCut_(0), ntracksEtaPtCut_(0), diff --git a/DataFormats/HeavyIonEvent/src/classes_def.xml b/DataFormats/HeavyIonEvent/src/classes_def.xml index 1ed47600b0800..e7e6481185269 100644 --- a/DataFormats/HeavyIonEvent/src/classes_def.xml +++ b/DataFormats/HeavyIonEvent/src/classes_def.xml @@ -3,7 +3,8 @@ - + + diff --git a/RecoHI/HiCentralityAlgos/plugins/CentralityBinProducer.cc b/RecoHI/HiCentralityAlgos/plugins/CentralityBinProducer.cc index 53a5416e81735..b635dab6bf662 100644 --- a/RecoHI/HiCentralityAlgos/plugins/CentralityBinProducer.cc +++ b/RecoHI/HiCentralityAlgos/plugins/CentralityBinProducer.cc @@ -59,7 +59,10 @@ class CentralityBinProducer : public edm::stream::EDProducer<> { EE = 11, ZDChitsPlus = 12, ZDChitsMinus = 13, - Missing = 14 + PFhf = 14, + PFhfPlus = 15, + PFhfMinus = 16, + Missing = 17 }; public: @@ -131,12 +134,18 @@ CentralityBinProducer::CentralityBinProducer(const edm::ParameterSet& iConfig) : varType_ = ZDChitsPlus; if (centralityVariable_ == "ZDChitsMinus") varType_ = ZDChitsMinus; + if (centralityVariable_ == "PFhf") + varType_ = PFhf; + if (centralityVariable_ == "PFhfPlus") + varType_ = PFhfPlus; + if (centralityVariable_ == "PFhfMinus") + varType_ = PFhfMinus; if (varType_ == Missing) { - std::string errorMessage = "Requested Centrality variable does not exist : " + centralityVariable_ + "\n" + - "Supported variables are: \n" + - "HFtowers HFtowersPlus HFtowersMinus HFtowersTrunc HFtowersPlusTrunc HFtowersMinusTrunc " - "HFhits PixelHits PixelTracks Tracks EB EE" + - "\n"; + std::string errorMessage = + "Requested Centrality variable does not exist : " + centralityVariable_ + "\n" + "Supported variables are: \n" + + "HFtowers HFtowersPlus HFtowersMinus HFtowersTrunc HFtowersPlusTrunc HFtowersMinusTrunc " + "HFhits PixelHits PixelTracks Tracks EB EE ZDChitsPlus ZDChitsMinus PFhf PFhfPlus PFhfMinus" + + "\n"; throw cms::Exception("Configuration", errorMessage); } @@ -207,6 +216,15 @@ void CentralityBinProducer::produce(edm::Event& iEvent, const edm::EventSetup& i case ZDChitsMinus: value = chandle_->zdcSumMinus(); break; + case PFhf: + value = chandle_->EtPFhfSum(); + break; + case PFhfPlus: + value = chandle_->EtPFhfSumPlus(); + break; + case PFhfMinus: + value = chandle_->EtPFhfSumMinus(); + break; default: throw cms::Exception("CentralityBinProducer", "Centrality variable not recognized."); } @@ -236,6 +254,10 @@ void CentralityBinProducer::beginRun(edm::Run const& iRun, const edm::EventSetup varType_ = ZDChitsMinus; if (centralityVariable_ == "ZDChitsMinus") varType_ = ZDChitsPlus; + if (centralityVariable_ == "PFhfPlus") + varType_ = PFhfMinus; + if (centralityVariable_ == "PFhfMinus") + varType_ = PFhfPlus; } prevRun_ = iRun.run(); diff --git a/RecoHI/HiCentralityAlgos/plugins/CentralityProducer.cc b/RecoHI/HiCentralityAlgos/plugins/CentralityProducer.cc index 65f62bda44acd..ad219224affd4 100644 --- a/RecoHI/HiCentralityAlgos/plugins/CentralityProducer.cc +++ b/RecoHI/HiCentralityAlgos/plugins/CentralityProducer.cc @@ -73,6 +73,7 @@ namespace reco { const bool doPixelCut_; const bool produceTracks_; const bool producePixelTracks_; + const bool producePF_; const double midRapidityRange_; const double trackPtCut_; @@ -91,6 +92,7 @@ namespace reco { const edm::EDGetTokenT srcPixelhits_; const edm::EDGetTokenT srcTracks_; const edm::EDGetTokenT srcPixelTracks_; + const edm::EDGetTokenT srcPF_; const edm::EDGetTokenT srcVertex_; const edm::EDGetTokenT reuseTag_; @@ -121,6 +123,7 @@ namespace reco { doPixelCut_(iConfig.getParameter("doPixelCut")), produceTracks_(iConfig.getParameter("produceTracks")), producePixelTracks_(iConfig.getParameter("producePixelTracks")), + producePF_(iConfig.getParameter("producePF")), midRapidityRange_(iConfig.getParameter("midRapidityRange")), trackPtCut_(iConfig.getParameter("trackPtCut")), trackEtaCut_(iConfig.getParameter("trackEtaCut")), @@ -147,6 +150,8 @@ namespace reco { srcPixelTracks_(producePixelTracks_ ? consumes(iConfig.getParameter("srcPixelTracks")) : edm::EDGetTokenT()), + srcPF_(producePF_ ? consumes(iConfig.getParameter("srcPF")) + : edm::EDGetTokenT()), srcVertex_((produceTracks_ || producePixelTracks_) ? consumes(iConfig.getParameter("srcVertex")) : edm::EDGetTokenT()), @@ -245,6 +250,22 @@ namespace reco { } } + if (producePF_) { + creco->etPFhfSumPlus_ = 0; + creco->etPFhfSumMinus_ = 0; + for (const auto& pf : iEvent.get(srcPF_)) { + if (pf.pdgId() != 1 && pf.pdgId() != 2) + continue; + if (pf.eta() > 0) + creco->etPFhfSumPlus_ += pf.pt(); + else + creco->etPFhfSumMinus_ += pf.pt(); + } + } else if (reuseAny_) { + creco->etPFhfSumMinus_ = inputCentrality->EtPFhfSumMinus(); + creco->etPFhfSumPlus_ = inputCentrality->EtPFhfSumPlus(); + } + if (produceEcalhits_) { creco->etEESumPlus_ = 0; creco->etEESumMinus_ = 0; @@ -497,6 +518,7 @@ namespace reco { desc.add("producePixelhits", true); desc.add("produceTracks", true); desc.add("producePixelTracks", true); + desc.add("producePF", true); desc.add("reUseCentrality", false); desc.add("srcHFhits", edm::InputTag("hfreco")); desc.add("srcTowers", edm::InputTag("towerMaker")); @@ -508,6 +530,7 @@ namespace reco { desc.add("srcVertex", edm::InputTag("hiSelectedVertex")); desc.add("srcReUse", edm::InputTag("hiCentrality")); desc.add("srcPixelTracks", edm::InputTag("hiPixel3PrimTracks")); + desc.add("srcPF", edm::InputTag("particleFlow")); desc.add("doPixelCut", true); desc.add("useQuality", true); desc.add("trackQuality", "highPurity"); diff --git a/RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py b/RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py index c3e9d73fc0d5b..2b01d0e3706ef 100644 --- a/RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py +++ b/RecoHI/HiCentralityAlgos/python/HiCentrality_cfi.py @@ -10,6 +10,7 @@ producePixelhits = cms.bool(True), produceTracks = cms.bool(True), producePixelTracks = cms.bool(True), + producePF = cms.bool(True), reUseCentrality = cms.bool(False), srcHFhits = cms.InputTag("hfreco"), @@ -22,6 +23,7 @@ srcVertex= cms.InputTag("hiSelectedVertex"), srcReUse = cms.InputTag("hiCentrality"), srcPixelTracks = cms.InputTag("hiPixel3PrimTracks"), + srcPF = cms.InputTag("particleFlow"), doPixelCut = cms.bool(True), useQuality = cms.bool(True), diff --git a/RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py b/RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py index 76b3be16775b2..15badf74fe346 100644 --- a/RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py +++ b/RecoHI/HiCentralityAlgos/python/pACentrality_cfi.py @@ -10,6 +10,7 @@ producePixelhits = cms.bool(True), produceTracks = cms.bool(True), producePixelTracks = cms.bool(True), + producePF = cms.bool(True), reUseCentrality = cms.bool(False), srcHFhits = cms.InputTag("hfreco"), @@ -22,6 +23,7 @@ srcVertex= cms.InputTag("offlinePrimaryVertices"), srcReUse = cms.InputTag("pACentrality"), srcPixelTracks = cms.InputTag("pixelTracks"), + srcPF = cms.InputTag("particleFlow"), doPixelCut = cms.bool(True), useQuality = cms.bool(True), diff --git a/RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.cc b/RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.cc index b3a73f128d776..03aefef0e4f2e 100644 --- a/RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.cc +++ b/RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.cc @@ -39,7 +39,7 @@ HiFJGridEmptyAreaCalculator::HiFJGridEmptyAreaCalculator(const edm::ParameterSet ntotalJet_ = 0; nyJet_ = 0; - pfCandsToken_ = consumes(iConfig.getParameter("pfCandSource")); + pfCandsToken_ = consumes(iConfig.getParameter("pfCandSource")); mapEtaToken_ = consumes>(iConfig.getParameter("mapEtaEdges")); mapRhoToken_ = consumes>(iConfig.getParameter("mapToRho")); mapRhoMToken_ = consumes>(iConfig.getParameter("mapToRhoM")); @@ -171,11 +171,10 @@ void HiFJGridEmptyAreaCalculator::produce(edm::Event& iEvent, const edm::EventSe void HiFJGridEmptyAreaCalculator::calculateGridRho(const edm::Event& iEvent, const edm::EventSetup& iSetup) { vector> scalarPt(ny_, vector(nphi_, 0.0)); - edm::Handle pfCands; + edm::Handle pfCands; iEvent.getByToken(pfCandsToken_, pfCands); - const reco::PFCandidateCollection* pfCandidateColl = pfCands.product(); - for (unsigned icand = 0; icand < pfCandidateColl->size(); icand++) { - const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand); + const auto& pfCandidateColl = pfCands.product(); + for (const auto& pfCandidate : *pfCandidateColl) { //use ony the particles within the eta range if (pfCandidate.eta() < ymin_ || pfCandidate.eta() > ymax_) continue; @@ -231,17 +230,17 @@ void HiFJGridEmptyAreaCalculator::calculateAreaFractionOfJets(const edm::Event& //calculate jet kt area fraction inside boundary by grid totalInboundArea_ = 0; - for (auto jet = jets->begin(); jet != jets->end(); ++jet) { - if (jet->eta() < etaminJet_ || jet->eta() > etamaxJet_) + for (const auto& jet : *jets) { + if (jet.eta() < etaminJet_ || jet.eta() > etamaxJet_) continue; - double areaKt = jet->jetArea(); - setupGridJet(&*jet); + double areaKt = jet.jetArea(); + setupGridJet(&jet); std::vector> pfIndicesJet; std::vector> pfIndicesJetInbound; int nConstitJet = 0; int nConstitJetInbound = 0; - for (auto daughter : jet->getJetConstituentsQuick()) { + for (const auto& daughter : jet.getJetConstituentsQuick()) { auto pfCandidate = static_cast(daughter); int jeta = tileIndexEtaJet(&*pfCandidate); @@ -323,7 +322,7 @@ void HiFJGridEmptyAreaCalculator::setupGrid(double etamin, double etamax) { //---------------------------------------------------------------------- // retrieve the grid tile index for a given PseudoJet -int HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::PFCandidate* pfCand) { +int HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::Candidate* pfCand) { // directly taking int does not work for values between -1 and 0 // so use floor instead // double iy_double = (p.rap() - _ymin) / _dy; @@ -345,7 +344,7 @@ int HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::PFCandidate* pfCand) { //---------------------------------------------------------------------- // retrieve the grid tile index for a given PseudoJet -int HiFJGridEmptyAreaCalculator::tileIndexEta(const reco::PFCandidate* pfCand) { +int HiFJGridEmptyAreaCalculator::tileIndexEta(const reco::Candidate* pfCand) { // directly taking int does not work for values between -1 and 0 // so use floor instead // double iy_double = (p.rap() - _ymin) / _dy; @@ -393,7 +392,7 @@ void HiFJGridEmptyAreaCalculator::setupGridJet(const reco::Jet* jet) { //---------------------------------------------------------------------- // retrieve the grid tile index for a given PseudoJet -int HiFJGridEmptyAreaCalculator::tileIndexEtaJet(const reco::PFCandidate* pfCand) { +int HiFJGridEmptyAreaCalculator::tileIndexEtaJet(const reco::Candidate* pfCand) { // directly taking int does not work for values between -1 and 0 // so use floor instead // double iy_double = (p.rap() - _ymin) / _dy; @@ -422,9 +421,9 @@ int HiFJGridEmptyAreaCalculator::numJetGridCells(std::vector int highestJEta = lowestJEta; //for each fixed phi value calculate the number of grids in eta - for (unsigned int iconst = 1; iconst < indices.size(); iconst++) { - int jphi = indices[iconst].first; - int jeta = indices[iconst].second; + for (const auto& index : indices) { + int jphi = index.first; + int jeta = index.second; if (jphi == lowestJPhi) { if (jeta < lowestJEta) lowestJEta = jeta; diff --git a/RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.h b/RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.h index 91ad85009a199..735825d5b4181 100644 --- a/RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.h +++ b/RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.h @@ -46,10 +46,10 @@ class HiFJGridEmptyAreaCalculator : public edm::stream::EDProducer<> { void setupGridJet(const reco::Jet* jet); /// retrieve the grid cell index for a given PseudoJet - int tileIndexJet(const reco::PFCandidate* pfCand); - int tileIndexEta(const reco::PFCandidate* pfCand); - int tileIndexEtaJet(const reco::PFCandidate* pfCand); - int tileIndexPhi(const reco::PFCandidate* pfCand); + int tileIndexJet(const reco::Candidate* pfCand); + int tileIndexEta(const reco::Candidate* pfCand); + int tileIndexEtaJet(const reco::Candidate* pfCand); + int tileIndexPhi(const reco::Candidate* pfCand); ///number of grid cells that overlap with jet constituents filling in the in between area int numJetGridCells(std::vector>& indices); @@ -106,7 +106,7 @@ class HiFJGridEmptyAreaCalculator : public edm::stream::EDProducer<> { /// input tokens edm::EDGetTokenT> jetsToken_; - edm::EDGetTokenT pfCandsToken_; + edm::EDGetTokenT pfCandsToken_; edm::EDGetTokenT> mapEtaToken_; edm::EDGetTokenT> mapRhoToken_; edm::EDGetTokenT> mapRhoMToken_;