CATopJetAlgorithm.h

Go to the documentation of this file.

#ifndef RecoJets_JetAlgorithms_CATopJetAlgorithm_h
#define RecoJets_JetAlgorithms_CATopJetAlgorithm_h



/* *********************************************************
 * \class CATopJetAlgorithm
 * Jet producer to produce top jets using the C-A algorithm to break
 * jets into subjets as described here:
 * "Top-tagging: A Method for Identifying Boosted Hadronic Tops"
 * David E. Kaplan, Keith Rehermann, Matthew D. Schwartz, Brock Tweedie
 * arXiv:0806.0848v1 [hep-ph] 
 *
 ************************************************************/

#include <vector>
#include <list>
#include <functional>
#include <TMath.h>
#include <iostream>

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/JetReco/interface/BasicJet.h"
#include "DataFormats/JetReco/interface/BasicJetCollection.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "RecoJets/JetAlgorithms/interface/JetAlgoHelper.h"
#include "RecoJets/JetAlgorithms/interface/CompoundPseudoJet.h"
#include "DataFormats/Candidate/interface/LeafCandidate.h"
#include "FWCore/Framework/interface/Event.h"

#include "RecoJets/JetAlgorithms/interface/CompoundPseudoJet.h"


#include <fastjet/JetDefinition.hh>
#include <fastjet/PseudoJet.hh>
#include <fastjet/ClusterSequence.hh>
#include <fastjet/ActiveAreaSpec.hh>
#include <fastjet/ClusterSequenceArea.hh>

class CATopJetAlgorithm{
 public:
  CATopJetAlgorithm(edm::InputTag mSrc,
            bool   verbose,
            int algorithm,
            int useAdjacency,
            double centralEtaCut,            
            double ptMin,                    
            const std::vector<double> & sumEtBins,      
            const std::vector<double> & rBins,       
            const std::vector<double> & ptFracBins,  
            const std::vector<double> & deltarBins,
            const std::vector<double> & nCellBins,
            double seedThreshold,            
            bool   useMaxTower,
            double sumEtEtaCut,
            double etFrac,
            boost::shared_ptr<fastjet::JetDefinition> fjJetDefinition,
            bool doAreaFastjet,
            boost::shared_ptr<fastjet::ActiveAreaSpec> fjActiveArea,
            double voronoiRfact) :
    mSrc_          (mSrc          ),
    verbose_       (verbose       ),
    algorithm_     (algorithm     ),
    useAdjacency_  (useAdjacency  ),
    centralEtaCut_ (centralEtaCut ), 
    ptMin_         (ptMin         ),         
    sumEtBins_     (sumEtBins     ),        
    rBins_         (rBins         ),         
    ptFracBins_    (ptFracBins    ),  
    deltarBins_    (deltarBins    ),
    nCellBins_     (nCellBins     ),
    seedThreshold_ (seedThreshold ), 
    useMaxTower_   (useMaxTower   ),
    sumEtEtaCut_   (sumEtEtaCut   ),   
    etFrac_        (etFrac        ),
    fjJetDefinition_(fjJetDefinition),
    doAreaFastjet_ (doAreaFastjet),
    fjActiveArea_  (fjActiveArea),
    voronoiRfact_  (voronoiRfact)
      { }

    void run( const std::vector<fastjet::PseudoJet> & cell_particles, 
          std::vector<CompoundPseudoJet> & hardjetsOutput 
          );

 private:

  edm::InputTag       mSrc_;                    //<! calo tower input source
  bool                verbose_;                 //<!
  int                 algorithm_;               //<! 0 = KT, 1 = CA, 2 = anti-KT
  int                 useAdjacency_;            //<! choose adjacency requirement:
                                                //<!    0 = no adjacency
                                                //<!    1 = deltar adjacency 
                                                //<!    2 = modified adjacency
                                                //<!    3 = calotower neirest neigbor based adjacency (untested)
  double              centralEtaCut_;           //<! eta for defining "central" jets                                     
  double              ptMin_;                   //<! lower pt cut on which jets to reco                                  
  std::vector<double> sumEtBins_;               //<! sumEt bins over which cuts vary. vector={bin 0 lower bound, bin 1 lower bound, ...}   
  std::vector<double> rBins_;                   //<! Jet distance paramter R. R values depend on sumEt bins.         
  std::vector<double> ptFracBins_;              //<! deltap = fraction of full jet pt for a subjet to be consider "hard". 
  std::vector<double> deltarBins_;              //<! Applicable only if useAdjacency=1. deltar adjacency values for each sumEtBin
  std::vector<double> nCellBins_;               //<! Applicable only if useAdjacency=3. number of cells apart for two subjets to be considered "independent"
  // NOT USED:
  double              seedThreshold_;           //<! calo tower seed threshold - NOT USED 
  bool                useMaxTower_;             //<! use max tower for jet adjacency criterion, false is to use the centroid - NOT USED
  double              sumEtEtaCut_;             //<! eta for event SumEt - NOT USED                                 
  double              etFrac_;                  //<! fraction of event sumEt / 2 for a jet to be considered "hard" - NOT USED 
  std::string         jetType_;                 //<! CaloJets or GenJets - NOT USED
  boost::shared_ptr<fastjet::JetDefinition> fjJetDefinition_; //<! jet definition to use
  bool                doAreaFastjet_; //<! whether or not to use the fastjet area
  boost::shared_ptr<fastjet::ActiveAreaSpec> fjActiveArea_; //<! fastjet area spec
  double              voronoiRfact_;  //<! fastjet voronoi area R factor

  // Decide if the two jets are in adjacent cells    
  bool adjacentCells(const fastjet::PseudoJet & jet1, const fastjet::PseudoJet & jet2, 
             const std::vector<fastjet::PseudoJet> & cell_particles,
             const fastjet::ClusterSequence & theClusterSequence,
             double nCellMin ) const;


  // Attempt to break up one "hard" jet into two "soft" jets

  bool decomposeJet(const fastjet::PseudoJet & theJet, 
            const fastjet::ClusterSequence & theClusterSequence, 
            const std::vector<fastjet::PseudoJet> & cell_particles,
            double ptHard, double nCellMin, double deltarcut,
            fastjet::PseudoJet & ja, fastjet::PseudoJet & jb, 
            std::vector<fastjet::PseudoJet> & leftovers) const;

};



#endif