Changes to jet pileup subtraction

RecoHI/HiJetAlgos/interface/MultipleAlgoIterator.h

@@ -5,22 +5,19 @@
 
 class MultipleAlgoIterator : public PileUpSubtractor {
 public:
-  MultipleAlgoIterator(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
-      : PileUpSubtractor(iConfig, std::move(iC)),
-        sumRecHits_(iConfig.getParameter<bool>("sumRecHits")),
-        dropZeroTowers_(iConfig.getUntrackedParameter<bool>("dropZeroTowers", true)) {
-    ;
-  }
+  MultipleAlgoIterator(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC);
   void offsetCorrectJets() override;
   void rescaleRMS(double s);
   double getEt(const reco::CandidatePtr& in) const;
   double getEta(const reco::CandidatePtr& in) const;
   void calculatePedestal(std::vector<fastjet::PseudoJet> const& coll) override;
   void subtractPedestal(std::vector<fastjet::PseudoJet>& coll) override;
+  void calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput) override;
 
+private:
+  double minimumTowersFraction_;
   bool sumRecHits_;
   bool dropZeroTowers_;
-  ~MultipleAlgoIterator() override { ; }
 };
 
 #endif

RecoHI/HiJetAlgos/plugins/ParticleTowerProducer.cc

@@ -0,0 +1,197 @@
+// -*- C++ -*-
+//
+// Package:    RecoHI/HiJetAlgos
+// Class:      ParticleTowerProducer
+//
+/**\class ParticleTowerProducer RecoHI/HiJetAlgos/plugins/ParticleTowerProducer.cc
+ Description: [one line class summary]
+ Implementation:
+     [Notes on implementation]
+*/
+//
+// Original Author:  Yetkin Yilmaz,32 4-A08,+41227673039,
+//         Created:  Thu Jan 20 19:53:58 CET 2011
+//
+//
+
+#include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"
+#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
+#include "DataFormats/Candidate/interface/Particle.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "DataFormats/Common/interface/SortedCollection.h"
+#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/src/WorkerMaker.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+#include "DataFormats/DetId/interface/DetId.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include <cmath>
+#include <cstdlib>
+#include <memory>
+#include <string>
+#include <vector>
+#include <map>
+#include <utility>
+
+class ParticleTowerProducer : public edm::stream::EDProducer<> {
+public:
+  explicit ParticleTowerProducer(const edm::ParameterSet&);
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+  int eta2ieta(double eta) const;
+  int phi2iphi(double phi, int ieta) const;
+  double ieta2eta(int ieta) const;
+  double iphi2phi(int iphi, int ieta) const;
+  // ----------member data ---------------------------
+
+  edm::EDGetTokenT<reco::PFCandidateCollection> src_;
+  const bool useHF_;
+
+  // tower edges from fast sim, used starting at index 30 for the HF
+  static constexpr int ietaMax = 42;
+  static constexpr double etaedge[ietaMax] = {
+      0.000, 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131,
+      1.218, 1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650,
+      2.853, 3.000, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889, 5.191};
+};
+//
+// constructors and destructor
+//
+ParticleTowerProducer::ParticleTowerProducer(const edm::ParameterSet& iConfig)
+    : src_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("src"))),
+      useHF_(iConfig.getParameter<bool>("useHF")) {
+  produces<CaloTowerCollection>();
+}
+
+//
+// member functions
+//
+
+// ------------ method called to produce the data  ------------
+void ParticleTowerProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  using namespace edm;
+
+  typedef std::pair<int, int> EtaPhi;
+  typedef std::map<EtaPhi, double> EtaPhiMap;
+  EtaPhiMap towers_;
+  towers_.clear();
+
+  auto const& inputs = iEvent.get(src_);
+  for (auto const& particle : inputs) {
+    double eta = particle.eta();
+
+    int ieta = eta2ieta(eta);
+    int iphi = phi2iphi(particle.phi(), ieta);
+
+    if (!useHF_ && abs(ieta) > 29)
+      continue;
+
+    EtaPhi ep(ieta, iphi);
+    towers_[ep] += particle.et();
+  }
+
+  auto prod = std::make_unique<CaloTowerCollection>();
+  prod->reserve(towers_.size());
+  for (auto const& tower : towers_) {
+    EtaPhi ep = tower.first;
+    double et = tower.second;
+
+    int ieta = ep.first;
+    int iphi = ep.second;
+
+    CaloTowerDetId newTowerId(ieta, iphi);  // totally dummy id
+
+    // currently sets et = pt, mass to zero
+    // pt, eta, phi, mass
+    reco::Particle::PolarLorentzVector p4(et, ieta2eta(ieta), iphi2phi(iphi, ieta), 0.);
+    GlobalPoint point(p4.x(), p4.y(), p4.z());
+    prod->emplace_back(newTowerId, p4.e(), 0, 0, 0, 0, p4, point, point);
+  }
+
+  iEvent.put(std::move(prod));
+}
+
+// Taken from FastSimulation/CalorimeterProperties/src/HCALProperties.cc
+// Note this returns an abs(ieta)
+int ParticleTowerProducer::eta2ieta(double eta) const {
+  // binary search in the array of towers eta edges
+
+  int ieta = 0;
+  while (fabs(eta) > etaedge[ieta] && ieta < ietaMax - 1) {
+    ++ieta;
+  }
+
+  if (eta < 0)
+    ieta = -ieta;
+  return ieta;
+}
+
+int ParticleTowerProducer::phi2iphi(double phi, int ieta) const {
+  phi = angle0to2pi::make0To2pi(phi);
+  int nphi = 72;
+  int n = 1;
+  if (abs(ieta) > 20)
+    n = 2;
+  if (abs(ieta) >= 40)
+    n = 4;
+
+  int iphi = (int)std::ceil(phi / 2.0 / M_PI * nphi / n);
+
+  iphi = n * (iphi - 1) + 1;
+
+  return iphi;
+}
+
+double ParticleTowerProducer::iphi2phi(int iphi, int ieta) const {
+  double phi = 0;
+  int nphi = 72;
+
+  int n = 1;
+  if (abs(ieta) > 20)
+    n = 2;
+  if (abs(ieta) >= 40)
+    n = 4;
+
+  int myphi = (iphi - 1) / n + 1;
+
+  phi = 2. * M_PI * (myphi - 0.5) / nphi * n;
+  while (phi > M_PI)
+    phi -= 2. * M_PI;
+
+  return phi;
+}
+
+double ParticleTowerProducer::ieta2eta(int ieta) const {
+  int sign = 1;
+  if (ieta < 0) {
+    sign = -1;
+    ieta = -ieta;
+  }
+
+  double eta = sign * (etaedge[ieta] + etaedge[ieta - 1]) / 2.;
+  return eta;
+}
+
+void ParticleTowerProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  // particleTowerProducer
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("src", edm::InputTag("particleFlow"));
+  desc.add<bool>("useHF", true);
+  descriptions.add("particleTowerProducer", desc);
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(ParticleTowerProducer);

RecoHI/HiJetAlgos/plugins/plugins.cc

@@ -13,8 +13,6 @@ DEFINE_EDM_PLUGIN(PileUpSubtractorFactory, MultipleAlgoIterator, "MultipleAlgoIt
 DEFINE_EDM_PLUGIN(PileUpSubtractorFactory, ParametrizedSubtractor, "ParametrizedSubtractor");
 DEFINE_EDM_PLUGIN(PileUpSubtractorFactory, ReflectedIterator, "ReflectedIterator");
 
-#include "RecoHI/HiJetAlgos/interface/ParticleTowerProducer.h"
-DEFINE_FWK_MODULE(ParticleTowerProducer);
 #include "RecoHI/HiJetAlgos/interface/HiGenCleaner.h"
 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
 #include "DataFormats/JetReco/interface/GenJetCollection.h"

RecoHI/HiJetAlgos/python/HiRecoPFJets_cff.py

@@ -5,10 +5,8 @@
 from RecoHI.HiJetAlgos.HiPFJetParameters_cff import *
 
 #pseudo towers for noise suppression background subtraction
-PFTowers = cms.EDProducer("ParticleTowerProducer",
-                          src = cms.InputTag("particleFlow"),
-                          useHF = cms.bool(False)
-                          )
+import RecoHI.HiJetAlgos.particleTowerProducer_cfi as _mod
+PFTowers = _mod.particleTowerProducer.clone(useHF = False)
 
 #dummy sequence to speed-up reconstruction in pp_on_AA era
 pfNoPileUpJMEHI = cms.EDFilter('GenericPFCandidateSelector',

RecoHI/HiJetAlgos/src/MultipleAlgoIterator.cc

@@ -3,10 +3,19 @@
 #include "RecoHI/HiJetAlgos/interface/MultipleAlgoIterator.h"
 #include "DataFormats/Candidate/interface/CandidateFwd.h"
 #include "DataFormats/Candidate/interface/Candidate.h"
+#include "DataFormats/Math/interface/deltaR.h"
 
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 using namespace std;
 
+MultipleAlgoIterator::MultipleAlgoIterator(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
+    : PileUpSubtractor(iConfig, std::move(iC)),
+      minimumTowersFraction_(iConfig.getParameter<double>("minimumTowersFraction")),
+      sumRecHits_(iConfig.getParameter<bool>("sumRecHits")),
+      dropZeroTowers_(iConfig.getUntrackedParameter<bool>("dropZeroTowers", true)) {
+  LogDebug("PileUpSubtractor") << "LIMITING THE MINIMUM TOWERS FRACTION TO : " << minimumTowersFraction_ << endl;
+}
+
 void MultipleAlgoIterator::rescaleRMS(double s) {
   for (std::map<int, double>::iterator iter = esigma_.begin(); iter != esigma_.end(); ++iter) {
     iter->second = s * (iter->second);
@@ -160,6 +169,59 @@ void MultipleAlgoIterator::calculatePedestal(vector<fastjet::PseudoJet> const& c
   }
 }
 
+void MultipleAlgoIterator::calculateOrphanInput(vector<fastjet::PseudoJet>& orphanInput) {
+  LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
+
+  (*fjInputs_) = fjOriginalInputs_;
+
+  vector<int> jettowers;                   // vector of towers indexed by "user_index"
+  vector<pair<int, int> > excludedTowers;  // vector of excluded ieta, iphi values
+
+  for (auto const& pseudojetTMP : *fjJets_) {
+    if (pseudojetTMP.perp() < puPtMin_)
+      continue;
+
+    // find towers within radiusPU_ of this jet
+    for (auto const im : allgeomid_) {
+      double dr = reco::deltaR(geo_->getPosition(im), pseudojetTMP);
+      auto exclude = find(excludedTowers.begin(), excludedTowers.end(), pair<int, int>(im.ieta(), im.iphi()));
+      if (dr < radiusPU_ && exclude == excludedTowers.end() &&
+          (geomtowers_[im.ieta()] - ntowersWithJets_[im.ieta()]) > minimumTowersFraction_ * (geomtowers_[im.ieta()])) {
+        ntowersWithJets_[im.ieta()]++;
+
+        excludedTowers.push_back(pair<int, int>(im.ieta(), im.iphi()));
+      }
+    }
+
+    for (auto const& it : *fjInputs_) {
+      int index = it.user_index();
+      int ie = ieta((*inputs_)[index]);
+      int ip = iphi((*inputs_)[index]);
+      auto exclude = find(excludedTowers.begin(), excludedTowers.end(), pair<int, int>(ie, ip));
+      if (exclude != excludedTowers.end()) {
+        jettowers.push_back(index);
+      }
+    }  // initial input collection
+
+  }  // pseudojets
+
+  //
+  // Create a new collections from the towers not included in jets
+  //
+
+  for (auto const& it : *fjInputs_) {
+    int index = it.user_index();
+    vector<int>::const_iterator itjet = find(jettowers.begin(), jettowers.end(), index);
+    if (itjet == jettowers.end()) {
+      const reco::CandidatePtr& originalTower = (*inputs_)[index];
+      fastjet::PseudoJet orphan(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
+      orphan.set_user_index(index);
+
+      orphanInput.push_back(orphan);
+    }
+  }
+}
+
 double MultipleAlgoIterator::getEt(const reco::CandidatePtr& in) const {
   const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
   const GlobalPoint& pos = geo_->getPosition(ctc->id());