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, std::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 41 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 45 of file CATopJetAlgorithm.h.

      : 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 295 of file CATopJetAlgorithm.cc.

References funct::abs(), reco::deltaPhi(), HLT_FULL_cff::eta1, and HLT_FULL_cff::eta2.

                                                             {
  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 314 of file CATopJetAlgorithm.cc.

References gather_cfg::cout, srCondWrite_cfg::deltaPhi, HLT_FULL_cff::deltaR, and dqmiolumiharvest::j.

                                                                                {
  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,
		std::shared_ptr< fastjet::ClusterSequence > &	fjClusterSeq
	)

Find the ProtoJets from the collection of input Candidates.

Definition at line 18 of file CATopJetAlgorithm.cc.

References gather_cfg::cout, greaterByEtPseudoJet(), mps_fire::i, dqmdumpme::indices, join(), perp(), csv2json::pieces, and ptMin_().

                                                                                 {
  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
  // std::shared_ptr<fastjet::ClusterSequence> fjClusterSeq;
  // if ( !doAreaFastjet_ ) {
  //   fjClusterSeq = std::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequence( cell_particles, jetDef ) );
  // } else if (voronoiRfact_ <= 0) {
  //   fjClusterSeq = std::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequenceArea( cell_particles, jetDef , *fjActiveArea_ ) );
  // } else {
  //   fjClusterSeq = std::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
  sort(centralJets.begin(), centralJets.end(), greaterByEtPseudoJet);

  // 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(), greaterByEtPseudoJet);

    // 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 86 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::centralEtaCut_

private

Definition at line 92 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::deltarBins_

private

Definition at line 98 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::etFrac_

private

Definition at line 105 of file CATopJetAlgorithm.h.

std::string CATopJetAlgorithm::jetType_

private

Definition at line 106 of file CATopJetAlgorithm.h.

edm::InputTag CATopJetAlgorithm::mSrc_

private

Definition at line 84 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::nCellBins_

private

Definition at line 100 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::ptFracBins_

private

Definition at line 97 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::ptMin_

private

Definition at line 93 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::rBins_

private

Definition at line 96 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::seedThreshold_

private

Definition at line 102 of file CATopJetAlgorithm.h.

std::vector<double> CATopJetAlgorithm::sumEtBins_

private

Definition at line 95 of file CATopJetAlgorithm.h.

double CATopJetAlgorithm::sumEtEtaCut_

private

Definition at line 104 of file CATopJetAlgorithm.h.

int CATopJetAlgorithm::useAdjacency_

private

Definition at line 87 of file CATopJetAlgorithm.h.

bool CATopJetAlgorithm::useMaxTower_

private

Definition at line 103 of file CATopJetAlgorithm.h.

bool CATopJetAlgorithm::verbose_

private

Definition at line 85 of file CATopJetAlgorithm.h.