JetPlusTrackCorrector.h

Go to the documentation of this file.

#ifndef RecoJets_JetPlusTrack_JetPlusTrackCorrector_h
#define RecoJets_JetPlusTrack_JetPlusTrackCorrector_h

#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/Event.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/ValueMap.h"
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
#include "DataFormats/JetReco/interface/JetTracksAssociation.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

#include "boost/range/iterator_range.hpp"
#include <sstream>
#include <string>

namespace edm {
  class Event;
  class EventSetup;
  class ParameterSet;
}

// --------------------------------------------------------
// -------------------- Helper classes --------------------
// --------------------------------------------------------


namespace jpt {
  
  class Map {

  public:

    Map( std::string, bool verbose = false );
    Map();
    ~Map();
    
    uint32_t nEtaBins() const;
    uint32_t nPtBins() const;
    
    double eta( uint32_t ) const;
    double pt( uint32_t ) const;

    uint32_t etaBin( double eta ) const;
    uint32_t ptBin( double pt ) const;
    
    double value( uint32_t eta_bin, uint32_t pt_bin ) const;

    double binCenterEta( uint32_t ) const;
    double binCenterPt( uint32_t ) const;
    
    void clear();
    void print( std::stringstream& ss ) const;

  private:

    class Element {
    public:
      Element() : ieta_(0), ipt_(0), eta_(0.), pt_(0.), val_(0.) {;} 
      uint32_t ieta_;
      uint32_t ipt_;
      double eta_;
      double pt_;
      double val_;
    };
    
    typedef std::vector<double> VDouble;
    typedef std::vector<VDouble> VVDouble;
    
    std::vector<double> eta_;
    std::vector<double> pt_;
    VVDouble data_;
    
  };

  inline uint32_t Map::nEtaBins() const { return eta_.size(); }
  inline uint32_t Map::nPtBins() const { return pt_.size(); }
  
  class Efficiency {

  public:

    Efficiency( const jpt::Map& response,
        const jpt::Map& efficiency,
        const jpt::Map& leakage );
      
    typedef std::pair<uint16_t,double> Pair;
    
    uint16_t nTrks( uint32_t eta_bin, uint32_t pt_bin ) const;
    
    double inConeCorr( uint32_t eta_bin, uint32_t pt_bin ) const;
    double outOfConeCorr( uint32_t eta_bin, uint32_t pt_bin ) const;
    
    uint32_t nEtaBins() const;
    uint32_t nPtBins() const;

    uint32_t size() const;
    bool empty() const;
    
    void addE( uint32_t eta_bin, uint32_t pt_bin, double energy );
    void reset();
    
    void print() const;
    
  private:
    
    Efficiency();

    double sumE( uint32_t eta_bin, uint32_t pt_bin ) const;
    double meanE( uint32_t eta_bin, uint32_t pt_bin ) const;

    bool check( uint32_t eta_bin, uint32_t pt_bin, std::string name = "check" ) const;

    typedef std::vector<Pair> VPair;
    typedef std::vector<VPair> VVPair;
    VVPair data_;

    const jpt::Map& response_;
    const jpt::Map& efficiency_;
    const jpt::Map& leakage_;
    
  };
  
  inline uint32_t Efficiency::nEtaBins() const { return response_.nEtaBins(); }
  inline uint32_t Efficiency::nPtBins() const { return response_.nPtBins(); }
  inline uint32_t Efficiency::size() const { return data_.size(); }
  inline bool Efficiency::empty() const { return data_.empty(); }

  class JetTracks {
  public:
    JetTracks();
    ~JetTracks();
    void clear();
    reco::TrackRefVector vertex_;
    reco::TrackRefVector caloFace_;
  };

  class MatchedTracks {
  public:
    MatchedTracks();
    ~MatchedTracks();
    void clear();
    reco::TrackRefVector inVertexInCalo_;
    reco::TrackRefVector outOfVertexInCalo_;
    reco::TrackRefVector inVertexOutOfCalo_; 
  };
  
//class Sum {
//public:
//  Sum() { theResponseOfChargedWithEff = 0.; theResponseOfChargedFull = 0.;}
//  void set( double a, double b) const {theResponseOfChargedWithEff=a;theResponseOfChargedFull=b;}
//  double theResponseOfChargedWithEff;
//  double theResponseOfChargedFull;
//};

}


// -------------------------------------------------------
// -------------------- JPT algorithm --------------------
// -------------------------------------------------------

class JetPlusTrackCorrector {

  // ---------- Public interface ----------
  
 public: 

  JetPlusTrackCorrector( const edm::ParameterSet& iPS, edm::ConsumesCollector&& iC );

  virtual ~JetPlusTrackCorrector();

  // Typedefs for 4-momentum
  typedef math::XYZTLorentzVector       P4;
  typedef math::PtEtaPhiELorentzVectorD PtEtaPhiE;
  typedef math::PtEtaPhiMLorentzVectorD PtEtaPhiM;
  
  double correction( const reco::Jet&, const reco::Jet&, const edm::Event&, const edm::EventSetup&, P4&,
             jpt::MatchedTracks &pions,
             jpt::MatchedTracks &muons,
             jpt::MatchedTracks &elecs,
             bool &validMatches) ;
  
  double correction( const reco::Jet&, const reco::Jet&, const edm::Event&, const edm::EventSetup&, 
             jpt::MatchedTracks &pions,
             jpt::MatchedTracks &muons,
             jpt::MatchedTracks &elecs,
             bool &validMatches) ;
  
  double correction( const reco::Jet& ) const;

  double correction( const P4& ) const;

  
//  double correctAA(  const reco::Jet&, const reco::TrackRefVector&, double&, const reco::TrackRefVector&,const reco::TrackRefVector&, double&) const;
  double correctAA(  const reco::Jet&, const reco::TrackRefVector&, double&, const reco::TrackRefVector&,const reco::TrackRefVector&,
                                                           double,
                                                           const reco::TrackRefVector&) const;

  
  bool eventRequired() const;
  
  bool vectorialCorrection() const;
  
  // ---------- Extended interface ----------


  double getResponseOfChargedWithEff() {return theResponseOfChargedWithEff;}
  double getResponseOfChargedWithoutEff() {return theResponseOfChargedWithoutEff;}
  double getSumPtWithEff() {return theSumPtWithEff;}
  double getSumPtWithoutEff() {return theSumPtWithoutEff;}
  double getSumEnergyWithEff() {return theSumEnergyWithEff;}
  double getSumEnergyWithoutEff() {return theSumEnergyWithoutEff;}
  double getSumPtForBeta() {return theSumPtForBeta;}


  bool canCorrect( const reco::Jet& ) const;
  
  bool matchTracks( const reco::Jet&, 
            const edm::Event&, 
            const edm::EventSetup&,
            jpt::MatchedTracks& pions, 
            jpt::MatchedTracks& muons, 
            jpt::MatchedTracks& elecs ) ;
  
  P4 pionCorrection( const P4& jet, const jpt::MatchedTracks& pions ) ;
  
  P4 muonCorrection( const P4& jet, const jpt::MatchedTracks& muons ) ;
  
  P4 elecCorrection( const P4& jet, const jpt::MatchedTracks& elecs ) const;

   
//  double theResponseOfChargedWithEff;
//  double theResponseOfChargedFull;
//  double setResp(double a) const {return a;}
//  double getChargedResponsesWithEff() { return theResponseOfChargedWithEff;}
//  double getChargedResponsesFull() { return theResponseOfChargedFull;}
  
  // ---------- Protected interface ----------

 protected: 

  // Some useful typedefs
  typedef reco::MuonCollection RecoMuons;
  typedef reco::GsfElectronCollection RecoElectrons;
  typedef edm::ValueMap<float> RecoElectronIds;
  typedef reco::JetTracksAssociation::Container JetTracksAssociations;
  typedef reco::TrackRefVector TrackRefs;

  virtual bool jetTrackAssociation( const reco::Jet&, 
                    const edm::Event&, 
                    const edm::EventSetup&,
                    jpt::JetTracks& ) const;
  
  bool jtaUsingEventData( const reco::Jet&, 
              const edm::Event&, 
              jpt::JetTracks& ) const;
  
  virtual void matchTracks( const jpt::JetTracks&,
                const edm::Event&, 
                jpt::MatchedTracks& pions, 
                jpt::MatchedTracks& muons, 
                jpt::MatchedTracks& elecs ) ;

  P4 pionCorrection( const P4& jet, 
             const TrackRefs& pions, 
             jpt::Efficiency&,
             bool in_cone_at_vertex,
             bool in_cone_at_calo_face ) ; 

  P4 muonCorrection( const P4& jet, 
             const TrackRefs& muons, 
             bool in_cone_at_vertex,
             bool in_cone_at_calo_face ) ;
  
  P4 elecCorrection( const P4& jet, 
             const TrackRefs& elecs, 
             bool in_cone_at_vertex,
             bool in_cone_at_calo_face ) ;

  P4 jetDirFromTracks( const P4& jet, 
               const jpt::MatchedTracks& pions,
               const jpt::MatchedTracks& muons,
               const jpt::MatchedTracks& elecs ) const;
  
  P4 calculateCorr( const P4& jet, 
            const TrackRefs&, 
            jpt::Efficiency&,
            bool in_cone_at_vertex,
            bool in_cone_at_calo_face,
            double mass,
            bool is_pion,
            double mip ) ;
  
  P4 pionEfficiency( const P4& jet, 
             const jpt::Efficiency&,
             bool in_cone_at_calo_face ) ;
  
  double checkScale( const P4& jet, P4& corrected ) const;
  
  bool getMuons( const edm::Event&, edm::Handle<RecoMuons>& ) const;

  bool getElectrons( const edm::Event&, 
             edm::Handle<RecoElectrons>&, 
             edm::Handle<RecoElectronIds>& ) const;
  
  bool matchMuons( TrackRefs::const_iterator&,
           const edm::Handle<RecoMuons>& ) const;
  
  bool matchElectrons( TrackRefs::const_iterator&,
               const edm::Handle<RecoElectrons>&, 
               const edm::Handle<RecoElectronIds>& ) const;
  
  bool failTrackQuality( TrackRefs::const_iterator& ) const;

  bool findTrack( const jpt::JetTracks&, 
          TrackRefs::const_iterator&,
          TrackRefs::iterator& ) const;

  bool findTrack( const jpt::MatchedTracks& pions, 
          const jpt::MatchedTracks& muons,
          const jpt::MatchedTracks& electrons,
          TrackRefs::const_iterator& ) const;

  bool tracksInCalo( const jpt::MatchedTracks& pions, 
             const jpt::MatchedTracks& muons,
             const jpt::MatchedTracks& elecs ) const;
  
  void rebuildJta( const reco::Jet&, 
           const JetTracksAssociations&, 
           TrackRefs& included,
           TrackRefs& excluded ) const;
  
  void excludeJta( const reco::Jet&, 
           const JetTracksAssociations&, 
           TrackRefs& included,
           const TrackRefs& excluded ) const;
  
  const jpt::Map& responseMap() const;
  const jpt::Map& efficiencyMap() const;
  const jpt::Map& leakageMap() const;
  
  JetPlusTrackCorrector() {;}

  // ---------- Protected member data ----------

 protected:
  
  // Some general configuration
  bool verbose_;
  bool vectorial_;
  bool vecResponse_;
  bool useInConeTracks_;
  bool useOutOfConeTracks_;
  bool useOutOfVertexTracks_;
  bool usePions_;
  bool useEff_;
  bool useMuons_;
  bool useElecs_;
  bool useTrackQuality_;
  
  // Jet-track association
  edm::InputTag jetTracksAtVertex_;
  edm::InputTag jetTracksAtCalo_;
  int jetSplitMerge_;
  edm::InputTag srcPVs_;
  double ptErrorQuality_;
  double dzVertexCut_;
  reco::Particle::Point vertex_;

  // Muons and electrons
  edm::InputTag muons_;
  edm::InputTag electrons_; 
  edm::InputTag electronIds_;
  
  // Filter tracks by quality
  reco::TrackBase::TrackQuality trackQuality_;

  // Response and efficiency maps  
  const jpt::Map response_;
  const jpt::Map efficiency_;
  const jpt::Map leakage_;

  // Mass    
  double pionMass_;
  double muonMass_;
  double elecMass_;

  // Jet-related
  double maxEta_;
  float theResponseOfChargedWithEff;
  float theResponseOfChargedWithoutEff;
  float theSumPtWithEff;
  float theSumPtWithoutEff;
  float theSumEnergyWithEff;
  float theSumEnergyWithoutEff;
  float theSumPtForBeta;
  jpt::Efficiency not_used{response_,efficiency_,leakage_};

  // ---------- Private data members ----------

 private:
  edm::EDGetTokenT<reco::JetTracksAssociation::Container> input_jetTracksAtVertex_token_;
  edm::EDGetTokenT<reco::JetTracksAssociation::Container> input_jetTracksAtCalo_token_;
  edm::EDGetTokenT<RecoMuons> inut_reco_muons_token_;
  edm::EDGetTokenT<reco::VertexCollection> input_pvCollection_token_;
  edm::EDGetTokenT<RecoElectrons> input_reco_elecs_token_;     
  edm::EDGetTokenT<RecoElectronIds> input_reco_elec_ids_token_;

};

// ---------- Inline methods ----------

inline double JetPlusTrackCorrector::correction( const reco::Jet& fJet, const reco::Jet& fJetcalo,
                         const edm::Event& event,
                         const edm::EventSetup& setup,
                         jpt::MatchedTracks &pions,
                         jpt::MatchedTracks &muons,
                         jpt::MatchedTracks &elecs,
                         bool &validMatches) {
  P4 not_used_for_scalar_correction;
  return correction( fJet, fJetcalo, event, setup, not_used_for_scalar_correction,pions,muons,elecs,validMatches );
}

inline bool JetPlusTrackCorrector::eventRequired() const { return true; }
inline bool JetPlusTrackCorrector::vectorialCorrection() const { return vectorial_; }
inline bool JetPlusTrackCorrector::canCorrect( const reco::Jet& jet ) const { return ( fabs( jet.eta() ) <= maxEta_ ); }

inline JetPlusTrackCorrector::P4 JetPlusTrackCorrector::pionCorrection( const P4& jet, 
                                    const TrackRefs& pions, 
                                    jpt::Efficiency& eff,
                                    bool in_cone_at_vertex,
                                    bool in_cone_at_calo_face ) {
  return calculateCorr( jet, pions, eff, in_cone_at_vertex, in_cone_at_calo_face, pionMass_, true, -1. );
}

inline JetPlusTrackCorrector::P4 JetPlusTrackCorrector::muonCorrection( const P4& jet, 
                                    const TrackRefs& muons, 
                                    bool in_cone_at_vertex,
                                    bool in_cone_at_calo_face ) {
  return calculateCorr( jet, muons, not_used, in_cone_at_vertex, in_cone_at_calo_face, muonMass_, false, 2. );
} 

inline JetPlusTrackCorrector::P4 JetPlusTrackCorrector::elecCorrection( const P4& jet, 
                                    const TrackRefs& elecs, 
                                    bool in_cone_at_vertex,
                                    bool in_cone_at_calo_face ) {
  return calculateCorr( jet, elecs, not_used, in_cone_at_vertex, in_cone_at_calo_face, elecMass_, false, 0. ); 
} 

inline double JetPlusTrackCorrector::checkScale( const P4& jet, P4& corrected ) const {
  if ( jet.energy() > 0. && ( corrected.energy() / jet.energy() ) < 0. ) { 
    corrected = jet; 
  }
  return corrected.energy() / jet.energy();
}

inline const jpt::Map& JetPlusTrackCorrector::responseMap() const { return response_; }
inline const jpt::Map& JetPlusTrackCorrector::efficiencyMap() const { return efficiency_; }
inline const jpt::Map& JetPlusTrackCorrector::leakageMap() const { return leakage_; }

#endif // RecoJets_JetPlusTracks_JetPlusTrackCorrector_h