JetAlgorithms/src/CATopJetAlgorithm.cc

0001 // Original author: Brock Tweedie (JHU)
0002 // Ported to CMSSW by: Sal Rappoccio (JHU)
0003
0004 #include "RecoJets/JetAlgorithms/interface/CATopJetAlgorithm.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "DataFormats/Math/interface/deltaPhi.h"
0010 #include "DataFormats/Math/interface/angle.h"
0011
0012 using namespace std;
0013 using namespace reco;
0014 using namespace edm;
0015
0016 //  Run the algorithm
0017 //  ------------------
0018 void CATopJetAlgorithm::run(const vector<fastjet::PseudoJet>& cell_particles,
0019                             vector<fastjet::PseudoJet>& hardjetsOutput,
0020                             std::shared_ptr<fastjet::ClusterSequence>& fjClusterSeq) {
0021   if (verbose_)
0022     cout << "Welcome to CATopSubJetAlgorithm::run" << endl;
0023
0024   // Sum Et of the event
0025   double sumEt = 0.;
0026
0027   //make a list of input objects ordered by ET and calculate sum et
0028   // list of fastjet pseudojet constituents
0029   for (unsigned i = 0; i < cell_particles.size(); ++i) {
0030     sumEt += cell_particles[i].perp();
0031   }
0032
0033   // Determine which bin we are in for et clustering
0034   int sumEtBinId = -1;
0035   for (unsigned int i = 0; i < sumEtBins_.size(); ++i) {
0036     if (sumEt > sumEtBins_[i])
0037       sumEtBinId = i;
0038   }
0039   if (verbose_)
0040     cout << "Using sumEt = " << sumEt << ", bin = " << sumEtBinId << endl;
0041
0042   // If the sum et is too low, exit
0043   if (sumEtBinId < 0) {
0044     return;
0045   }
0046
0047   // empty 4-vector
0048   fastjet::PseudoJet blankJetA(0, 0, 0, 0);
0049   blankJetA.set_user_index(-1);
0050   const fastjet::PseudoJet blankJet = blankJetA;
0051
0052   // Define adjacency variables which depend on which sumEtBin we are in
0053   double deltarcut = deltarBins_[sumEtBinId];
0054   double nCellMin = nCellBins_[sumEtBinId];
0055
0056   if (verbose_)
0057     cout << "useAdjacency_ = " << useAdjacency_ << endl;
0058   if (verbose_ && useAdjacency_ == 0)
0059     cout << "No Adjacency" << endl;
0060   if (verbose_ && useAdjacency_ == 1)
0061     cout << "using deltar adjacency" << endl;
0062   if (verbose_ && useAdjacency_ == 2)
0063     cout << "using modified adjacency" << endl;
0064   if (verbose_ && useAdjacency_ == 3)
0065     cout << "using calorimeter nearest neighbor based adjacency" << endl;
0066   if (verbose_ && useAdjacency_ == 1)
0067     cout << "Using deltarcut = " << deltarcut << endl;
0068   if (verbose_ && useAdjacency_ == 3)
0069     cout << "Using nCellMin = " << nCellMin << endl;
0070
0071   if (verbose_)
0072     cout << "About to do jet clustering in CA" << endl;
0073   // run the jet clustering
0074
0075   //cluster the jets with the jet definition jetDef:
0076   // run algorithm
0077   // std::shared_ptr<fastjet::ClusterSequence> fjClusterSeq;
0078   // if ( !doAreaFastjet_ ) {
0079   //   fjClusterSeq = std::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequence( cell_particles, jetDef ) );
0080   // } else if (voronoiRfact_ <= 0) {
0081   //   fjClusterSeq = std::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequenceArea( cell_particles, jetDef , *fjActiveArea_ ) );
0082   // } else {
0083   //   fjClusterSeq = std::shared_ptr<fastjet::ClusterSequence>( new fastjet::ClusterSequenceVoronoiArea( cell_particles, jetDef , fastjet::VoronoiAreaSpec(voronoiRfact_) ) );
0084   // }
0085
0086   if (verbose_)
0087     cout << "Getting inclusive jets" << endl;
0088   // Get the transient inclusive jets
0089   vector<fastjet::PseudoJet> inclusiveJets = fjClusterSeq->inclusive_jets(ptMin_);
0090
0091   if (verbose_)
0092     cout << "Getting central jets" << endl;
0093   // Find the transient central jets
0094   vector<fastjet::PseudoJet> centralJets;
0095   for (unsigned int i = 0; i < inclusiveJets.size(); i++) {
0096     if (inclusiveJets[i].perp() > ptMin_ && fabs(inclusiveJets[i].rapidity()) < centralEtaCut_) {
0097       centralJets.push_back(inclusiveJets[i]);
0098     }
0099   }
0100   // Sort the transient central jets in Et
0101   sort(centralJets.begin(), centralJets.end(), greaterByEtPseudoJet);
0102
0103   // These will store the 4-vectors of each hard jet
0104   vector<math::XYZTLorentzVector> p4_hardJets;
0105
0106   // These will store the indices of each subjet that
0107   // are present in each jet
0108   vector<vector<int> > indices(centralJets.size());
0109
0110   // Loop over central jets, attempt to find substructure
0111   vector<fastjet::PseudoJet>::iterator jetIt = centralJets.begin(), centralJetsEnd = centralJets.end();
0112   if (verbose_)
0113     cout << "Loop over jets" << endl;
0114   int i = 0;
0115   for (; jetIt != centralJetsEnd; ++jetIt) {
0116     if (verbose_)
0117       cout << "\nJet " << i << endl;
0118     i++;
0119     fastjet::PseudoJet localJet = *jetIt;
0120
0121     // Get the 4-vector for this jet
0122     p4_hardJets.push_back(math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e()));
0123
0124     // jet decomposition.  try to find 3 or 4 hard, well-localized subjets, characteristic of a boosted top.
0125     if (verbose_)
0126       cout << "local jet pt = " << localJet.perp() << endl;
0127     if (verbose_)
0128       cout << "deltap = " << ptFracBins_[sumEtBinId] << endl;
0129
0130     double ptHard = ptFracBins_[sumEtBinId] * localJet.perp();
0131     vector<fastjet::PseudoJet> leftoversAll;
0132
0133     // stage 1:  primary decomposition.  look for when the jet declusters into two hard subjets
0134     if (verbose_)
0135       cout << "Doing decomposition 1" << endl;
0136     fastjet::PseudoJet ja, jb;
0137     vector<fastjet::PseudoJet> leftovers1;
0138     bool hardBreak1 =
0139         decomposeJet(localJet, *fjClusterSeq, cell_particles, ptHard, nCellMin, deltarcut, ja, jb, leftovers1);
0140     leftoversAll.insert(leftoversAll.end(), leftovers1.begin(), leftovers1.end());
0141
0142     // stage 2:  secondary decomposition.  look for when the hard subjets found above further decluster into two hard sub-subjets
0143     //
0144     // ja -> jaa+jab ?
0145     if (verbose_)
0146       cout << "Doing decomposition 2. ja->jaa+jab?" << endl;
0147     fastjet::PseudoJet jaa, jab;
0148     vector<fastjet::PseudoJet> leftovers2a;
0149     bool hardBreak2a = false;
0150     if (hardBreak1)
0151       hardBreak2a = decomposeJet(ja, *fjClusterSeq, cell_particles, ptHard, nCellMin, deltarcut, jaa, jab, leftovers2a);
0152     leftoversAll.insert(leftoversAll.end(), leftovers2a.begin(), leftovers2a.end());
0153     // jb -> jba+jbb ?
0154     if (verbose_)
0155       cout << "Doing decomposition 2. ja->jba+jbb?" << endl;
0156     fastjet::PseudoJet jba, jbb;
0157     vector<fastjet::PseudoJet> leftovers2b;
0158     bool hardBreak2b = false;
0159     if (hardBreak1)
0160       hardBreak2b = decomposeJet(jb, *fjClusterSeq, cell_particles, ptHard, nCellMin, deltarcut, jba, jbb, leftovers2b);
0161     leftoversAll.insert(leftoversAll.end(), leftovers2b.begin(), leftovers2b.end());
0162
0163     // NOTE:  it might be good to consider some checks for whether these subjets can be further decomposed.  e.g., the above procedure leaves
0164     //        open the possibility of "subjets" that actually consist of two or more distinct hard clusters.  however, this kind of thing
0165     //        is a rarity for the simulations so far considered.
0166
0167     // proceed if one or both of the above hard subjets successfully decomposed
0168     if (verbose_)
0169       cout << "Done with decomposition" << endl;
0170
0171     if (verbose_)
0172       cout << "hardBreak1 = " << hardBreak1 << endl;
0173     if (verbose_)
0174       cout << "hardBreak2a = " << hardBreak2a << endl;
0175     if (verbose_)
0176       cout << "hardBreak2b = " << hardBreak2b << endl;
0177
0178     fastjet::PseudoJet hardA = blankJet, hardB = blankJet, hardC = blankJet, hardD = blankJet;
0179     if (!hardBreak1) {
0180       hardA = localJet;
0181       hardB = blankJet;
0182       hardC = blankJet;
0183       hardD = blankJet;
0184       if (verbose_)
0185         cout << "Hardbreak failed. Save subjet1=localJet" << endl;
0186     }
0187     if (hardBreak1 && !hardBreak2a && !hardBreak2b) {
0188       hardA = ja;
0189       hardB = jb;
0190       hardC = blankJet;
0191       hardD = blankJet;
0192       if (verbose_)
0193         cout << "First decomposition succeeded, both second decompositions failed. Save subjet1=ja subjet2=jb" << endl;
0194     }
0195     if (hardBreak1 && hardBreak2a && !hardBreak2b) {
0196       hardA = jaa;
0197       hardB = jab;
0198       hardC = jb;
0199       hardD = blankJet;
0200       if (verbose_)
0201         cout << "First decomposition succeeded, ja split succesfully, jb did not split. Save subjet1=jaa subjet2=jab "
0202                 "subjet3=jb"
0203              << endl;
0204     }
0205     if (hardBreak1 && !hardBreak2a && hardBreak2b) {
0206       hardA = jba;
0207       hardB = jbb;
0208       hardC = ja;
0209       hardD = blankJet;
0210       if (verbose_)
0211         cout << "First decomposition succeeded, jb split succesfully, ja did not split. Save subjet1=jba subjet2=jbb "
0212                 "subjet3=ja"
0213              << endl;
0214     }
0215     if (hardBreak1 && hardBreak2a && hardBreak2b) {
0216       hardA = jaa;
0217       hardB = jab;
0218       hardC = jba;
0219       hardD = jbb;
0220       if (verbose_)
0221         cout << "First decomposition and both secondary decompositions succeeded. Save subjet1=jaa subjet2=jab "
0222                 "subjet3=jba subjet4=jbb"
0223              << endl;
0224     }
0225
0226     // check if we are left with >= 3 hard subjets
0227     fastjet::PseudoJet subjet1 = blankJet;
0228     fastjet::PseudoJet subjet2 = blankJet;
0229     fastjet::PseudoJet subjet3 = blankJet;
0230     fastjet::PseudoJet subjet4 = blankJet;
0231     subjet1 = hardA;
0232     subjet2 = hardB;
0233     subjet3 = hardC;
0234     subjet4 = hardD;
0235
0236     // record the hard subjets
0237     vector<fastjet::PseudoJet> hardSubjets;
0238
0239     if (verbose_) {
0240       std::cout << "HardA : user_index = " << hardA.user_index() << ", (Pt,Y,Phi,M) = (" << hardA.pt() << ", "
0241                 << hardA.rapidity() << ", " << hardA.phi() << ", " << hardA.m() << ")" << std::endl;
0242
0243       std::cout << "HardB : user_index = " << hardB.user_index() << ", (Pt,Y,Phi,M) = (" << hardB.pt() << ", "
0244                 << hardB.rapidity() << ", " << hardB.phi() << ", " << hardB.m() << ")" << std::endl;
0245
0246       std::cout << "HardC : user_index = " << hardC.user_index() << ", (Pt,Y,Phi,M) = (" << hardC.pt() << ", "
0247                 << hardC.rapidity() << ", " << hardC.phi() << ", " << hardC.m() << ")" << std::endl;
0248
0249       std::cout << "HardD : user_index = " << hardD.user_index() << ", (Pt,Y,Phi,M) = (" << hardD.pt() << ", "
0250                 << hardD.rapidity() << ", " << hardD.phi() << ", " << hardD.m() << ")" << std::endl;
0251     }
0252
0253     // Check to see if any subjects are counted amongst the "hard" subjets from previous
0254     // lines. NOTE: In Fastjet 3.0, the default "user_index" changed from 0 to -1, so
0255     // this can no longer be used as a designator for the veto of "blankJet" subjets,
0256     // and now switch to pt > some small value.
0257     if (subjet1.pt() > 0.0001)
0258       hardSubjets.push_back(subjet1);
0259     if (subjet2.pt() > 0.0001)
0260       hardSubjets.push_back(subjet2);
0261     if (subjet3.pt() > 0.0001)
0262       hardSubjets.push_back(subjet3);
0263     if (subjet4.pt() > 0.0001)
0264       hardSubjets.push_back(subjet4);
0265     sort(hardSubjets.begin(), hardSubjets.end(), greaterByEtPseudoJet);
0266
0267     // Use new fastjet functionality to create a Pseudojet from constituents
0268     fastjet::PseudoJet candidate = join(hardSubjets);
0269     // Reset the jet's 4-vector to the "ungroomed" value
0270     candidate.reset_momentum(jetIt->px(), jetIt->py(), jetIt->pz(), jetIt->e());
0271
0272     if (verbose_) {
0273       std::cout << "Final top-jet candidate: (Pt,Y,Phi,M) = (" << candidate.pt() << ", " << candidate.rapidity() << ", "
0274                 << candidate.phi() << ", " << candidate.m() << ")" << std::endl;
0275       std::vector<fastjet::PseudoJet> pieces = candidate.pieces();
0276       std::cout << "Number of pieces = " << pieces.size() << std::endl;
0277       for (std::vector<fastjet::PseudoJet>::const_iterator ibegin = pieces.begin(), iend = pieces.end(), i = ibegin;
0278            i != iend;
0279            ++i) {
0280         std::cout << "   Piece : " << i - ibegin << ", (Pt,Y,Phi,M) = (" << i->pt() << ", " << i->rapidity() << ", "
0281                   << i->phi() << ", " << i->m() << ")" << std::endl;
0282       }
0283     }
0284     // Add to the list
0285     hardjetsOutput.push_back(candidate);
0286   }
0287 }
0288
0289 //-----------------------------------------------------------------------
0290 // determine whether two clusters (made of calorimeter towers) are living on "adjacent" cells.  if they are, then
0291 // we probably shouldn't consider them to be independent objects!
0292 //
0293 // From Sal: Ignoring genjet case
0294 //
0295 bool CATopJetAlgorithm::adjacentCells(const fastjet::PseudoJet& jet1,
0296                                       const fastjet::PseudoJet& jet2,
0297                                       const vector<fastjet::PseudoJet>& cell_particles,
0298                                       const fastjet::ClusterSequence& theClusterSequence,
0299                                       double nCellMin) const {
0300   double eta1 = jet1.rapidity();
0301   double phi1 = jet1.phi();
0302   double eta2 = jet2.rapidity();
0303   double phi2 = jet2.phi();
0304
0305   double deta = abs(eta2 - eta1) / 0.087;
0306   double dphi = fabs(reco::deltaPhi(phi2, phi1)) / 0.087;
0307
0308   return ((deta + dphi) <= nCellMin);
0309 }
0310
0311 //-------------------------------------------------------------------------
0312 // attempt to decompose a jet into "hard" subjets, where hardness is set by ptHard
0313 //
0314 bool CATopJetAlgorithm::decomposeJet(const fastjet::PseudoJet& theJet,
0315                                      const fastjet::ClusterSequence& theClusterSequence,
0316                                      const vector<fastjet::PseudoJet>& cell_particles,
0317                                      double ptHard,
0318                                      double nCellMin,
0319                                      double deltarcut,
0320                                      fastjet::PseudoJet& ja,
0321                                      fastjet::PseudoJet& jb,
0322                                      vector<fastjet::PseudoJet>& leftovers) const {
0323   bool goodBreak;
0324   fastjet::PseudoJet j = theJet;
0325   double InputObjectPt = j.perp();
0326   if (verbose_)
0327     cout << "Input Object Pt = " << InputObjectPt << endl;
0328   if (verbose_)
0329     cout << "ptHard = " << ptHard << endl;
0330   leftovers.clear();
0331   if (verbose_)
0332     cout << "start while loop" << endl;
0333
0334   while (true) {  // watch out for infinite loop!
0335     goodBreak = theClusterSequence.has_parents(j, ja, jb);
0336     if (!goodBreak) {
0337       if (verbose_)
0338         cout << "bad break. this is one cell. can't decluster anymore." << endl;
0339       break;  // this is one cell, can't decluster anymore
0340     }
0341
0342     if (verbose_)
0343       cout << "good break. ja Pt = " << ja.perp() << " jb Pt = " << jb.perp() << endl;
0344
0345     /// Adjacency Requirement ///
0346
0347     // check if clusters are adjacent using a constant deltar adjacency.
0348     double clusters_deltar = fabs(ja.eta() - jb.eta()) + fabs(deltaPhi(ja.phi(), jb.phi()));
0349
0350     if (verbose_ && useAdjacency_ == 1)
0351       cout << "clusters_deltar = " << clusters_deltar << endl;
0352     if (verbose_ && useAdjacency_ == 1)
0353       cout << "deltar cut = " << deltarcut << endl;
0354
0355     if (useAdjacency_ == 1 && clusters_deltar < deltarcut) {
0356       if (verbose_)
0357         cout << "clusters too close. consant adj. break." << endl;
0358       break;
0359     }
0360
0361     // Check if clusters are adjacent using a DeltaR adjacency which is a function of pT.
0362     double clusters_deltaR = deltaR(ja.rapidity(), ja.phi(), jb.rapidity(), jb.phi());
0363
0364     if (verbose_ && useAdjacency_ == 2)
0365       cout << "clusters_deltaR = " << clusters_deltaR << endl;
0366     if (verbose_ && useAdjacency_ == 2)
0367       cout << "0.4-0.0004*InputObjectPt = " << 0.4 - 0.0004 * InputObjectPt << endl;
0368
0369     if (useAdjacency_ == 2 && clusters_deltaR < 0.4 - 0.0004 * InputObjectPt) {
0370       if (verbose_)
0371         cout << "clusters too close. modified adj. break." << endl;
0372       break;
0373     }
0374
0375     // Check if clusters are adjacent in the calorimeter.
0376     if (useAdjacency_ == 3 && adjacentCells(ja, jb, cell_particles, theClusterSequence, nCellMin)) {
0377       if (verbose_)
0378         cout << "clusters too close in the calorimeter. calorimeter adj. break." << endl;
0379       break;  // the clusters are "adjacent" in the calorimeter => shouldn't have decomposed
0380     }
0381
0382     if (verbose_)
0383       cout << "clusters pass distance cut" << endl;
0384
0385     /// Pt Fraction Requirement ///
0386
0387     if (verbose_)
0388       cout << "ptHard = " << ptHard << endl;
0389
0390     if (ja.perp() < ptHard && jb.perp() < ptHard) {
0391       if (verbose_)
0392         cout << "two soft clusters. dead end" << endl;
0393       break;  // broke into two soft clusters, dead end
0394     }
0395
0396     if (ja.perp() > ptHard && jb.perp() > ptHard) {
0397       if (verbose_)
0398         cout << "two hard clusters. done" << endl;
0399       return true;  // broke into two hard clusters, we're done!
0400     }
0401
0402     else if (ja.perp() > jb.perp()) {  // broke into one hard and one soft, ditch the soft one and try again
0403       if (verbose_)
0404         cout << "ja hard jb soft. try to split hard. j = ja" << endl;
0405       j = ja;
0406       vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(jb);
0407       leftovers.insert(leftovers.end(), particles.begin(), particles.end());
0408     } else {
0409       if (verbose_)
0410         cout << "ja hard jb soft. try to split hard. j = jb" << endl;
0411       j = jb;
0412       vector<fastjet::PseudoJet> particles = theClusterSequence.constituents(ja);
0413       leftovers.insert(leftovers.end(), particles.begin(), particles.end());
0414     }
0415   }
0416
0417   if (verbose_)
0418     cout << "did not decluster." << endl;  // did not decluster into hard subjets
0419
0420   ja.reset(0, 0, 0, 0);
0421   jb.reset(0, 0, 0, 0);
0422   leftovers.clear();
0423   return false;
0424 }