Add PF HF sum energy to reco::Centrality

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<std::vector<reco::Vertex> > vertex;
     iEvent.getByToken(vertexToken, vertex);
     h_vertex_x->Fill(vertex->begin()->x());

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;

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<std::vector<reco::Vertex> > vertex;
     iEvent.getByToken(vertexToken, vertex);
     h_vertex_x->Fill(vertex->begin()->x());

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;

DQMOffline/JetMET/src/JetAnalyzer_HeavyIons.cc

@@ -642,16 +642,16 @@ void JetAnalyzer_HeavyIons::analyze(const edm::Event &mEvent, const edm::EventSe
   edm::Handle<reco::Centrality> cent;
   mEvent.getByToken(centralityToken, cent);  //_centralitytag comes from the cfg
 
+  if (!cent.isValid())
+    return;
+
   mHF->Fill(cent->EtHFtowerSum());
   Float_t HF_energy = cent->EtHFtowerSum();
 
   //for later when centrality gets added to RelVal
   //edm::Handle<int> cbin;
   //mEvent.getByToken(centralityBinToken, cbin);
 
-  if (!cent.isValid())
-    return;
-
   /*int hibin = -999;
   if(cbin.isValid()){
     hibin = *cbin;

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_;

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),

DataFormats/HeavyIonEvent/src/classes_def.xml

@@ -3,7 +3,8 @@
  <version ClassVersion="11" checksum="393681411"/>
 </class>
 <class name="edm::Wrapper<reco::EvtPlane>"/>
-<class name="reco::Centrality" ClassVersion="11">
+<class name="reco::Centrality" ClassVersion="12">
+ <version ClassVersion="12" checksum="926006939"/>
  <version ClassVersion="11" checksum="688182903"/>
  <version ClassVersion="10" checksum="3770933457"/>
 </class>

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();
 

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<SiPixelRecHitCollection> srcPixelhits_;
     const edm::EDGetTokenT<TrackCollection> srcTracks_;
     const edm::EDGetTokenT<TrackCollection> srcPixelTracks_;
+    const edm::EDGetTokenT<reco::CandidateView> srcPF_;
     const edm::EDGetTokenT<VertexCollection> srcVertex_;
     const edm::EDGetTokenT<Centrality> reuseTag_;
 
@@ -121,6 +123,7 @@ namespace reco {
         doPixelCut_(iConfig.getParameter<bool>("doPixelCut")),
         produceTracks_(iConfig.getParameter<bool>("produceTracks")),
         producePixelTracks_(iConfig.getParameter<bool>("producePixelTracks")),
+        producePF_(iConfig.getParameter<bool>("producePF")),
         midRapidityRange_(iConfig.getParameter<double>("midRapidityRange")),
         trackPtCut_(iConfig.getParameter<double>("trackPtCut")),
         trackEtaCut_(iConfig.getParameter<double>("trackEtaCut")),
@@ -147,6 +150,8 @@ namespace reco {
         srcPixelTracks_(producePixelTracks_
                             ? consumes<TrackCollection>(iConfig.getParameter<edm::InputTag>("srcPixelTracks"))
                             : edm::EDGetTokenT<TrackCollection>()),
+        srcPF_(producePF_ ? consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("srcPF"))
+                          : edm::EDGetTokenT<reco::CandidateView>()),
         srcVertex_((produceTracks_ || producePixelTracks_)
                        ? consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("srcVertex"))
                        : edm::EDGetTokenT<VertexCollection>()),
@@ -245,6 +250,22 @@ namespace reco {
       }
     }
 
+    if (producePF_) {
+      creco->etPFhfSumPlus_ = 0;
+      creco->etPFhfSumMinus_ = 0;
+      for (const auto& pf : iEvent.get(srcPF_)) {
+        const auto apdgId = std::abs(pf.pdgId());
+        bool isHF = apdgId == 1 || apdgId == 2;
+        if (isHF && pf.eta() > 0)
+          creco->etPFhfSumPlus_ += pf.pt();
+        else if (isHF && pf.eta() < 0)
+          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<bool>("producePixelhits", true);
     desc.add<bool>("produceTracks", true);
     desc.add<bool>("producePixelTracks", true);
+    desc.add<bool>("producePF", true);
     desc.add<bool>("reUseCentrality", false);
     desc.add<edm::InputTag>("srcHFhits", edm::InputTag("hfreco"));
     desc.add<edm::InputTag>("srcTowers", edm::InputTag("towerMaker"));
@@ -508,6 +530,7 @@ namespace reco {
     desc.add<edm::InputTag>("srcVertex", edm::InputTag("hiSelectedVertex"));
     desc.add<edm::InputTag>("srcReUse", edm::InputTag("hiCentrality"));
     desc.add<edm::InputTag>("srcPixelTracks", edm::InputTag("hiPixel3PrimTracks"));
+    desc.add<edm::InputTag>("srcPF", edm::InputTag("particleFlow"));
     desc.add<bool>("doPixelCut", true);
     desc.add<bool>("useQuality", true);
     desc.add<string>("trackQuality", "highPurity");

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),

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),

RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.cc

@@ -39,7 +39,7 @@ HiFJGridEmptyAreaCalculator::HiFJGridEmptyAreaCalculator(const edm::ParameterSet
   ntotalJet_ = 0;
   nyJet_ = 0;
 
-  pfCandsToken_ = consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"));
+  pfCandsToken_ = consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("pfCandSource"));
   mapEtaToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapEtaEdges"));
   mapRhoToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapToRho"));
   mapRhoMToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapToRhoM"));
@@ -171,11 +171,11 @@ void HiFJGridEmptyAreaCalculator::produce(edm::Event& iEvent, const edm::EventSe
 void HiFJGridEmptyAreaCalculator::calculateGridRho(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
   vector<vector<double>> scalarPt(ny_, vector<double>(nphi_, 0.0));
 
-  edm::Handle<reco::PFCandidateCollection> pfCands;
+  edm::Handle<reco::CandidateView> pfCands;
   iEvent.getByToken(pfCandsToken_, pfCands);
-  const reco::PFCandidateCollection* pfCandidateColl = pfCands.product();
+  const auto& pfCandidateColl = pfCands.product();
   for (unsigned icand = 0; icand < pfCandidateColl->size(); icand++) {
-    const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
+    const auto& pfCandidate = pfCandidateColl->at(icand);
     //use ony the particles within the eta range
     if (pfCandidate.eta() < ymin_ || pfCandidate.eta() > ymax_)
       continue;
@@ -323,7 +323,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 +345,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 +393,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;

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<std::pair<int, int>>& indices);
@@ -106,7 +106,7 @@ class HiFJGridEmptyAreaCalculator : public edm::stream::EDProducer<> {
 
   /// input tokens
   edm::EDGetTokenT<edm::View<reco::Jet>> jetsToken_;
-  edm::EDGetTokenT<reco::PFCandidateCollection> pfCandsToken_;
+  edm::EDGetTokenT<reco::CandidateView> pfCandsToken_;
   edm::EDGetTokenT<std::vector<double>> mapEtaToken_;
   edm::EDGetTokenT<std::vector<double>> mapRhoToken_;
   edm::EDGetTokenT<std::vector<double>> mapRhoMToken_;