CATopJetAlgorithm Class Reference

#include <CATopJetAlgorithm.h>

Public Member Functions
	CATopJetAlgorithm (const 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)
void	run (const std::vector< fastjet::PseudoJet > &cell_particles, std::vector< fastjet::PseudoJet > &hardjetsOutput, boost::shared_ptr< fastjet::ClusterSequence > &fjClusterSeq)
	Find the ProtoJets from the collection of input Candidates. More...

Private Member Functions
bool	adjacentCells (const fastjet::PseudoJet &jet1, const fastjet::PseudoJet &jet2, const std::vector< fastjet::PseudoJet > &cell_particles, const fastjet::ClusterSequence &theClusterSequence, double nCellMin) const
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

Private Attributes
int	algorithm_
double	centralEtaCut_
std::vector< double >	deltarBins_
double	etFrac_
std::string	jetType_
edm::InputTag	mSrc_
std::vector< double >	nCellBins_
std::vector< double >	ptFracBins_
double	ptMin_
std::vector< double >	rBins_
double	seedThreshold_
std::vector< double >	sumEtBins_
double	sumEtEtaCut_
int	useAdjacency_
bool	useMaxTower_
bool	verbose_

Detailed Description

Definition at line 46 of file CATopJetAlgorithm.h.

Constructor & Destructor Documentation

CATopJetAlgorithm::CATopJetAlgorithm ( const 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  )

inline

Constructor

Definition at line 50 of file CATopJetAlgorithm.h.

References run().

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

      { }

Member Function Documentation

bool CATopJetAlgorithm::adjacentCells ( const fastjet::PseudoJet &  jet1, const fastjet::PseudoJet &  jet2, const std::vector< fastjet::PseudoJet > &  cell_particles, const fastjet::ClusterSequence &  theClusterSequence, double  nCellMin  ) const

private

Definition at line 268 of file CATopJetAlgorithm.cc.

References funct::abs(), and reco::deltaPhi().

                                                              {
    
    
    double eta1 = jet1.rapidity();
    double phi1 = jet1.phi();
    double eta2 = jet2.rapidity();
    double phi2 = jet2.phi();
    
    double deta = abs(eta2 - eta1) / 0.087;
    double dphi = fabs( reco::deltaPhi(phi2,phi1) ) / 0.087;
    
    return ( ( deta + dphi ) <= nCellMin );
}

bool CATopJetAlgorithm::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

private

Adjacency Requirement ///

Pt Fraction Requirement ///

Definition at line 290 of file CATopJetAlgorithm.cc.

References gather_cfg::cout, hiPixelPairStep_cff::deltaPhi, deltaR(), and HadronAndPartonSelector_cfi::particles.

                                                                                 {
    
    bool goodBreak;
    fastjet::PseudoJet j = theJet;
    double InputObjectPt = j.perp();
    if ( verbose_ )cout<<"Input Object Pt = "<<InputObjectPt<<endl;
    if ( verbose_ )cout<<"ptHard = "<<ptHard<<endl;
    leftovers.clear();
    if ( verbose_ )cout<<"start while loop"<<endl;
    
    while (true) {                                                      // watch out for infinite loop!
        goodBreak = theClusterSequence.has_parents(j,ja,jb);
        if (!goodBreak){
            if ( verbose_ )cout<<"bad break. this is one cell. can't decluster anymore."<<endl;
            break;         // this is one cell, can't decluster anymore
        }
        
        if ( verbose_ )cout<<"good break. ja Pt = "<<ja.perp()<<" jb Pt = "<<jb.perp()<<endl;
        
        
        // check if clusters are adjacent using a constant deltar adjacency.
        double clusters_deltar=fabs(ja.eta()-jb.eta())+fabs(deltaPhi(ja.phi(),jb.phi()));
        
        if ( verbose_  && useAdjacency_ ==1)cout<<"clusters_deltar = "<<clusters_deltar<<endl;
        if ( verbose_  && useAdjacency_ ==1)cout<<"deltar cut = "<<deltarcut<<endl;
        
        if ( useAdjacency_==1 && clusters_deltar < deltarcut){
            if ( verbose_ )cout<<"clusters too close. consant adj. break."<<endl;
            break;
        } 
        
        // Check if clusters are adjacent using a DeltaR adjacency which is a function of pT.
        double clusters_deltaR=deltaR( ja.rapidity(), ja.phi(), jb.rapidity(), jb.phi() );
        
        if ( verbose_  && useAdjacency_ ==2)cout<<"clusters_deltaR = "<<clusters_deltaR<<endl;
        if ( verbose_  && useAdjacency_ ==2)cout<<"0.4-0.0004*InputObjectPt = "<<0.4-0.0004*InputObjectPt<<endl;
        
        if ( useAdjacency_==2 && clusters_deltaR < 0.4-0.0004*InputObjectPt)
        {
            if ( verbose_ )cout<<"clusters too close. modified adj. break."<<endl;
            break;
        } 

        // Check if clusters are adjacent in the calorimeter. 
        if ( useAdjacency_==3 &&  adjacentCells(ja,jb,cell_particles,theClusterSequence,nCellMin) ){                  
            if ( verbose_ )cout<<"clusters too close in the calorimeter. calorimeter adj. break."<<endl;
            break;         // the clusters are "adjacent" in the calorimeter => shouldn't have decomposed
        }
                
        if ( verbose_ )cout<<"clusters pass distance cut"<<endl;
        
        
        if ( verbose_ )cout<<"ptHard = "<<ptHard<<endl;
        
        if (ja.perp() < ptHard && jb.perp() < ptHard){
            if ( verbose_ )cout<<"two soft clusters. dead end"<<endl;
            break;         // broke into two soft clusters, dead end
        }
        
        if (ja.perp() > ptHard && jb.perp() > ptHard){
            if ( verbose_ )cout<<"two hard clusters. done"<<endl;
            return true;   // broke into two hard clusters, we're done!
        }
        
        else if (ja.perp() > jb.perp()) {                              // broke into one hard and one soft, ditch the soft one and try again
            if ( verbose_ )cout<<"ja hard jb soft. try to split hard. j = ja"<<endl; 
            j = ja;
            vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(jb);
            leftovers.insert(leftovers.end(),particles.begin(),particles.end());
        }
        else {
            if ( verbose_ )cout<<"ja hard jb soft. try to split hard. j = jb"<<endl; 
            j = jb;
            vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(ja);
            leftovers.insert(leftovers.end(),particles.begin(),particles.end());
        }
    }
    
    if ( verbose_ )cout<<"did not decluster."<<endl;  // did not decluster into hard subjets
    
    ja.reset(0,0,0,0);
    jb.reset(0,0,0,0);
    leftovers.clear();
    return false;
}

void CATopJetAlgorithm::run	(	const std::vector< fastjet::PseudoJet > &	cell_particles,
		std::vector< fastjet::PseudoJet > &	hardjetsOutput,
		boost::shared_ptr< fastjet::ClusterSequence > &	fjClusterSeq
	)

Find the ProtoJets from the collection of input Candidates.

Definition at line 19 of file CATopJetAlgorithm.cc.

References gather_cfg::cout, mps_fire::i, join(), perp(), csvLumiCalc::pieces, and objects.METAnalyzer::sumEt.

Referenced by CATopJetAlgorithm().

{
    if ( verbose_ ) cout << "Welcome to CATopSubJetAlgorithm::run" << endl;
    
    // Sum Et of the event
    double sumEt = 0.;
    
    //make a list of input objects ordered by ET and calculate sum et
    // list of fastjet pseudojet constituents
    for (unsigned i = 0; i < cell_particles.size(); ++i) {
        sumEt += cell_particles[i].perp();
    }
    
    // Determine which bin we are in for et clustering
    int sumEtBinId = -1;
    for ( unsigned int i = 0; i < sumEtBins_.size(); ++i ) {
        if ( sumEt > sumEtBins_[i] ) sumEtBinId = i;
    }
    if ( verbose_ ) cout << "Using sumEt = " << sumEt << ", bin = " << sumEtBinId << endl;
    
    // If the sum et is too low, exit
    if ( sumEtBinId < 0 ) {    
        return;
    }
    
    // empty 4-vector
    fastjet::PseudoJet blankJetA(0,0,0,0);
    blankJetA.set_user_index(-1);
    const fastjet::PseudoJet blankJet = blankJetA;
    
    // Define adjacency variables which depend on which sumEtBin we are in
    double deltarcut = deltarBins_[sumEtBinId];
    double nCellMin = nCellBins_[sumEtBinId];
    
    if ( verbose_ )cout<<"useAdjacency_ = "<<useAdjacency_<<endl;
    if ( verbose_ && useAdjacency_==0)cout<<"No Adjacency"<<endl;
    if ( verbose_ && useAdjacency_==1)cout<<"using deltar adjacency"<<endl;
    if ( verbose_ && useAdjacency_==2)cout<<"using modified adjacency"<<endl;
    if ( verbose_ && useAdjacency_==3)cout<<"using calorimeter nearest neighbor based adjacency"<<endl;
    if ( verbose_ && useAdjacency_==1)cout << "Using deltarcut = " << deltarcut << endl;
    if ( verbose_ && useAdjacency_==3)cout << "Using nCellMin = " << nCellMin << endl;
    
    if ( verbose_ ) cout << "About to do jet clustering in CA" << endl;
    // run the jet clustering

    //cluster the jets with the jet definition jetDef:
    // run algorithm
    // boost::shared_ptr<fastjet::ClusterSequence> fjClusterSeq;
    // if ( !doAreaFastjet_ ) {
    //   fjClusterSeq = boost::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequence( cell_particles, jetDef ) );
    // } else if (voronoiRfact_ <= 0) {
    //   fjClusterSeq = boost::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequenceArea( cell_particles, jetDef , *fjActiveArea_ ) );
    // } else {
    //   fjClusterSeq = boost::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequenceVoronoiArea( cell_particles, jetDef , fastjet::VoronoiAreaSpec(voronoiRfact_) ) );
    // }
    
    if ( verbose_ ) cout << "Getting inclusive jets" << endl;
    // Get the transient inclusive jets
    vector<fastjet::PseudoJet> inclusiveJets = fjClusterSeq->inclusive_jets(ptMin_);
    
    if ( verbose_ ) cout << "Getting central jets" << endl;
    // Find the transient central jets
    vector<fastjet::PseudoJet> centralJets;
    for (unsigned int i = 0; i < inclusiveJets.size(); i++) {
        
        if (inclusiveJets[i].perp() > ptMin_ && fabs(inclusiveJets[i].rapidity()) < centralEtaCut_) {
            centralJets.push_back(inclusiveJets[i]);
        }
    }
    // Sort the transient central jets in Et
    GreaterByEtPseudoJet compEt;
    sort( centralJets.begin(), centralJets.end(), compEt );
    
    // These will store the 4-vectors of each hard jet
    vector<math::XYZTLorentzVector> p4_hardJets;
    
    // These will store the indices of each subjet that 
    // are present in each jet
    vector<vector<int> > indices( centralJets.size() );
    
    // Loop over central jets, attempt to find substructure
    vector<fastjet::PseudoJet>::iterator jetIt = centralJets.begin(),
      centralJetsEnd = centralJets.end();
    if ( verbose_ )cout<<"Loop over jets"<<endl;
    int i=0;
    for ( ; jetIt != centralJetsEnd; ++jetIt ) {
        if ( verbose_ )cout<<"\nJet "<<i<<endl;
        i++;
        fastjet::PseudoJet localJet = *jetIt;
        
        // Get the 4-vector for this jet
        p4_hardJets.push_back( math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e() ));
        
        // jet decomposition.  try to find 3 or 4 hard, well-localized subjets, characteristic of a boosted top.
        if ( verbose_ )cout<<"local jet pt = "<<localJet.perp()<<endl;
        if ( verbose_ )cout<<"deltap = "<<ptFracBins_[sumEtBinId]<<endl;
        
        double ptHard = ptFracBins_[sumEtBinId]*localJet.perp();
        vector<fastjet::PseudoJet> leftoversAll;
        
        // stage 1:  primary decomposition.  look for when the jet declusters into two hard subjets
        if ( verbose_ ) cout << "Doing decomposition 1" << endl;
        fastjet::PseudoJet ja, jb;
        vector<fastjet::PseudoJet> leftovers1;
        bool hardBreak1 = decomposeJet(localJet,*fjClusterSeq,cell_particles,ptHard,nCellMin,deltarcut,ja,jb,leftovers1);
        leftoversAll.insert(leftoversAll.end(),leftovers1.begin(),leftovers1.end());
        
        // stage 2:  secondary decomposition.  look for when the hard subjets found above further decluster into two hard sub-subjets
        //
        // ja -> jaa+jab ?
        if ( verbose_ ) cout << "Doing decomposition 2. ja->jaa+jab?" << endl;
        fastjet::PseudoJet jaa, jab;
        vector<fastjet::PseudoJet> leftovers2a;
        bool hardBreak2a = false;
        if (hardBreak1)  hardBreak2a = decomposeJet(ja,*fjClusterSeq,cell_particles,ptHard,nCellMin,deltarcut,jaa,jab,leftovers2a);
        leftoversAll.insert(leftoversAll.end(),leftovers2a.begin(),leftovers2a.end());
        // jb -> jba+jbb ?
        if ( verbose_ ) cout << "Doing decomposition 2. ja->jba+jbb?" << endl;
        fastjet::PseudoJet jba, jbb;
        vector<fastjet::PseudoJet> leftovers2b;
        bool hardBreak2b = false;
        if (hardBreak1)  hardBreak2b = decomposeJet(jb,*fjClusterSeq,cell_particles,ptHard,nCellMin,deltarcut,jba,jbb,leftovers2b);
        leftoversAll.insert(leftoversAll.end(),leftovers2b.begin(),leftovers2b.end());
        
        // NOTE:  it might be good to consider some checks for whether these subjets can be further decomposed.  e.g., the above procedure leaves
        //        open the possibility of "subjets" that actually consist of two or more distinct hard clusters.  however, this kind of thing
        //        is a rarity for the simulations so far considered.
        
        // proceed if one or both of the above hard subjets successfully decomposed
        if ( verbose_ ) cout << "Done with decomposition" << endl;

        if ( verbose_ ) cout<<"hardBreak1 = "<<hardBreak1<<endl;
        if ( verbose_ ) cout<<"hardBreak2a = "<<hardBreak2a<<endl;
        if ( verbose_ ) cout<<"hardBreak2b = "<<hardBreak2b<<endl;
            
        fastjet::PseudoJet hardA = blankJet, hardB = blankJet, hardC = blankJet, hardD = blankJet;
        if (!hardBreak1) { 
          hardA = localJet;  
          hardB = blankJet;  
          hardC = blankJet; 
          hardD = blankJet; 
          if(verbose_)cout<<"Hardbreak failed. Save subjet1=localJet"<<endl;
        } 
        if (hardBreak1 && !hardBreak2a && !hardBreak2b) { 
          hardA = ja;  
          hardB = jb;  
          hardC = blankJet; 
          hardD = blankJet; 
          if(verbose_)cout<<"First decomposition succeeded, both second decompositions failed. Save subjet1=ja subjet2=jb"<<endl;
        }
        if (hardBreak1 && hardBreak2a && !hardBreak2b) { 
          hardA = jaa; 
          hardB = jab; 
          hardC = jb; 
          hardD = blankJet; 
          if(verbose_)cout<<"First decomposition succeeded, ja split succesfully, jb did not split. Save subjet1=jaa subjet2=jab subjet3=jb"<<endl;
        }
        if (hardBreak1 && !hardBreak2a &&  hardBreak2b) { 
          hardA = jba; 
          hardB = jbb; 
          hardC = ja; 
          hardD = blankJet; 
          if(verbose_)cout<<"First decomposition succeeded, jb split succesfully, ja did not split. Save subjet1=jba subjet2=jbb subjet3=ja"<<endl;
        }
        if (hardBreak1 && hardBreak2a &&  hardBreak2b) { 
          hardA = jaa; 
          hardB = jab; 
          hardC = jba; 
          hardD = jbb; 
          if(verbose_)cout<<"First decomposition and both secondary decompositions succeeded. Save subjet1=jaa subjet2=jab subjet3=jba subjet4=jbb"<<endl;
        }
        
        // check if we are left with >= 3 hard subjets
        fastjet::PseudoJet subjet1 = blankJet;
        fastjet::PseudoJet subjet2 = blankJet;
        fastjet::PseudoJet subjet3 = blankJet;
        fastjet::PseudoJet subjet4 = blankJet;
        subjet1 = hardA; subjet2 = hardB; subjet3 = hardC; subjet4 = hardD;
        
        // record the hard subjets
        vector<fastjet::PseudoJet> hardSubjets;

        if ( verbose_ ) {
          std::cout << "HardA : user_index = " << hardA.user_index() << ", (Pt,Y,Phi,M) = (" 
                << hardA.pt() << ", " << hardA.rapidity() << ", " 
                << hardA.phi() << ", " << hardA.m() << ")" << std::endl;

          std::cout << "HardB : user_index = " << hardB.user_index() << ", (Pt,Y,Phi,M) = (" 
                << hardB.pt() << ", " << hardB.rapidity() << ", " 
                << hardB.phi() << ", " << hardB.m() << ")" << std::endl;

          std::cout << "HardC : user_index = " << hardC.user_index() << ", (Pt,Y,Phi,M) = (" 
                << hardC.pt() << ", " << hardC.rapidity() << ", " 
                << hardC.phi() << ", " << hardC.m() << ")" << std::endl;

          std::cout << "HardD : user_index = " << hardD.user_index() << ", (Pt,Y,Phi,M) = (" 
                << hardD.pt() << ", " << hardD.rapidity() << ", " 
                << hardD.phi() << ", " << hardD.m() << ")" << std::endl;
        }

        // Check to see if any subjects are counted amongst the "hard" subjets from previous
        // lines. NOTE: In Fastjet 3.0, the default "user_index" changed from 0 to -1, so
        // this can no longer be used as a designator for the veto of "blankJet" subjets,
        // and now switch to pt > some small value. 
        if ( subjet1.pt() > 0.0001 )
            hardSubjets.push_back(subjet1);
        if ( subjet2.pt() > 0.0001 )
            hardSubjets.push_back(subjet2);
        if ( subjet3.pt() > 0.0001 )
            hardSubjets.push_back(subjet3);
        if ( subjet4.pt() > 0.0001 )
            hardSubjets.push_back(subjet4);
        sort(hardSubjets.begin(), hardSubjets.end(), compEt );

        // Use new fastjet functionality to create a Pseudojet from constituents
        fastjet::PseudoJet candidate = join(hardSubjets);
        // Reset the jet's 4-vector to the "ungroomed" value
        candidate.reset_momentum( jetIt->px(), jetIt->py(), jetIt->pz(), jetIt->e() );

        if ( verbose_ ) {
          std::cout << "Final top-jet candidate: (Pt,Y,Phi,M) = (" 
                << candidate.pt() << ", " << candidate.rapidity() << ", " 
                << candidate.phi() << ", " << candidate.m() << ")" << std::endl;
          std::vector<fastjet::PseudoJet> pieces = candidate.pieces();
          std::cout << "Number of pieces = " << pieces.size() << std::endl;
          for ( std::vector<fastjet::PseudoJet>::const_iterator ibegin = pieces.begin(), iend = pieces.end(), i = ibegin;
            i != iend; ++i ) {
            std::cout << "   Piece : " << i - ibegin << ", (Pt,Y,Phi,M) = (" 
                  << i->pt() << ", " << i->rapidity() << ", " 
                  << i->phi() << ", " << i->m() << ")" << std::endl;
          }
        }
        // Add to the list
        hardjetsOutput.push_back( candidate );      
    }
}

Member Data Documentation

int CATopJetAlgorithm::algorithm_

private

Definition at line 93 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::centralEtaCut_

private

Definition at line 99 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::deltarBins_

private

Definition at line 104 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::etFrac_

private

Definition at line 110 of file CATopJetAlgorithm.h.

std::string CATopJetAlgorithm::jetType_

private

Definition at line 111 of file CATopJetAlgorithm.h.

edm::InputTag CATopJetAlgorithm::mSrc_

private

Definition at line 91 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::nCellBins_

private

Definition at line 105 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::ptFracBins_

private

Definition at line 103 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::ptMin_

private

Definition at line 100 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::rBins_

private

Definition at line 102 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::seedThreshold_

private

Definition at line 107 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::sumEtBins_

private

Definition at line 101 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::sumEtEtaCut_

private

Definition at line 109 of file CATopJetAlgorithm.h.

int CATopJetAlgorithm::useAdjacency_

private

Definition at line 94 of file CATopJetAlgorithm.h.

bool CATopJetAlgorithm::useMaxTower_

private

Definition at line 108 of file CATopJetAlgorithm.h.

bool CATopJetAlgorithm::verbose_

private

Definition at line 92 of file CATopJetAlgorithm.h.