PileUpSubtractor.h

Go to the documentation of this file.

#ifndef __PUSubtractor__
#define __PUSubtractor__

#include <vector>

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/GhostedAreaSpec.hh"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"

class PileUpSubtractor {
public:
  typedef std::shared_ptr<fastjet::ClusterSequence> ClusterSequencePtr;
  typedef std::shared_ptr<fastjet::GhostedAreaSpec> ActiveAreaSpecPtr;
  typedef std::shared_ptr<fastjet::RangeDefinition> RangeDefPtr;
  typedef std::shared_ptr<fastjet::JetDefinition> JetDefPtr;

  PileUpSubtractor(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC);
  virtual ~PileUpSubtractor() { ; }

  virtual void setDefinition(JetDefPtr const& jetDef);
  virtual void reset(std::vector<edm::Ptr<reco::Candidate> >& input,
                     std::vector<fastjet::PseudoJet>& towers,
                     std::vector<fastjet::PseudoJet>& output);
  virtual void setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup);
  virtual void calculatePedestal(std::vector<fastjet::PseudoJet> const& coll);
  virtual void subtractPedestal(std::vector<fastjet::PseudoJet>& coll);
  virtual void calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput);
  virtual void offsetCorrectJets();
  virtual double getMeanAtTower(const reco::CandidatePtr& in) const;
  virtual double getSigmaAtTower(const reco::CandidatePtr& in) const;
  virtual double getPileUpAtTower(const reco::CandidatePtr& in) const;
  virtual double getPileUpEnergy(int ijet) const { return jetOffset_[ijet]; }
  virtual double getCone(double cone, double eta, double phi, double& et, double& pu);
  int getN(const reco::CandidatePtr& in) const;
  int getNwithJets(const reco::CandidatePtr& in) const;

  int ieta(const reco::CandidatePtr& in) const;
  int iphi(const reco::CandidatePtr& in) const;

protected:
  // From jet producer
  JetDefPtr fjJetDefinition_;                         // fastjet jet definition
  ClusterSequencePtr fjClusterSeq_;                   // fastjet cluster sequence
  std::vector<edm::Ptr<reco::Candidate> >* inputs_;   // input candidates
  std::vector<fastjet::PseudoJet>* fjInputs_;         // fastjet inputs
  std::vector<fastjet::PseudoJet>* fjJets_;           // fastjet jets
  std::vector<fastjet::PseudoJet> fjOriginalInputs_;  // to back-up unsubtracted fastjet inputs

  // PU subtraction parameters
  bool reRunAlgo_;
  bool doAreaFastjet_;
  bool doRhoFastjet_;
  double jetPtMin_;
  double puPtMin_;

  double ghostEtaMax;
  int activeAreaRepeats;
  double ghostArea;

  double nSigmaPU_;                     // number of sigma for pileup
  double radiusPU_;                     // pileup radius
  ActiveAreaSpecPtr fjActiveArea_;      // fastjet active area definition
  CaloGeometry const* geo_;             // geometry
  int ietamax_;                         // maximum eta in geometry
  int ietamin_;                         // minimum eta in geometry
  std::vector<HcalDetId> allgeomid_;    // all det ids in the geometry
  std::map<int, int> geomtowers_;       // map of geometry towers to det id
  std::map<int, int> ntowersWithJets_;  // number of towers with jets
  std::map<int, double> esigma_;        // energy sigma
  std::map<int, double> emean_;         // energy mean

  std::vector<double> jetOffset_;
};

#include "FWCore/PluginManager/interface/PluginFactory.h"
namespace edm {
  class ParameterSet;
  class EventSetup;
  class ConsumesCollector;
}  // namespace edm
typedef edmplugin::PluginFactory<PileUpSubtractor*(const edm::ParameterSet&, edm::ConsumesCollector&&)>
    PileUpSubtractorFactory;

#endif