JetAlgorithms/src/CMSInsideOutAlgorithm.cc

0001 #include "RecoJets/JetAlgorithms/interface/CMSInsideOutAlgorithm.h"
0002
0003 #include "RecoJets/JetAlgorithms/interface/CompoundPseudoJet.h"
0004
0005 using namespace std;
0006
0007 void CMSInsideOutAlgorithm::run(const std::vector<fastjet::PseudoJet>& fInput,
0008                                 std::vector<fastjet::PseudoJet>& fOutput) {
0009   //make a list of input objects
0010   list<fastjet::PseudoJet> input;
0011   for (std::vector<fastjet::PseudoJet>::const_iterator candIter = fInput.begin(); candIter != fInput.end();
0012        ++candIter) {
0013     input.push_back(*candIter);
0014   }
0015
0016   while (!input.empty() && input.front().perp() > seedThresholdPt_) {
0017     //get seed eta/phi
0018     double seedEta = input.front().eta();
0019     double seedPhi = input.front().phi();
0020
0021     //find iterators to those objects that are in the max cone size
0022     list<inputListIter> maxCone;
0023     // add seed, then test elements after seed
0024     inputListIter iCand = input.begin();
0025     maxCone.push_back(iCand++);
0026     for (; iCand != input.end(); ++iCand) {
0027       const fastjet::PseudoJet& candidate = *iCand;
0028       if (reco::deltaR2(seedEta, candidate.eta(), seedPhi, candidate.phi()) < maxSizeSquared_)
0029         maxCone.push_back(iCand);
0030     }
0031     //sort objects by increasing DR about the seed  directions
0032     maxCone.sort(ListIteratorLesserByDeltaR(seedEta, seedPhi));
0033     list<inputListIter>::const_iterator position = maxCone.begin();
0034     bool limitReached = false;
0035     double totalET = (**position).perp();
0036     ++position;
0037     while (position != maxCone.end() && !limitReached) {
0038       const fastjet::PseudoJet& theCandidate = **position;
0039       double candidateET = theCandidate.perp() + totalET;
0040       double candDR2 = reco::deltaR2(seedEta, theCandidate.eta(), seedPhi, theCandidate.phi());
0041       if (candDR2 < minSizeSquared_ || candDR2 * candidateET * candidateET < growthParameterSquared_)
0042         totalET = candidateET;
0043       else
0044         limitReached = true;
0045       ++position;
0046     }
0047     //turn this into a final jet
0048     fastjet::PseudoJet final;
0049     for (list<inputListIter>::const_iterator iNewJet = maxCone.begin(); iNewJet != position; ++iNewJet) {
0050       final += **iNewJet;
0051       input.erase(*iNewJet);
0052     }
0053     fOutput.push_back(final);
0054   }  // end loop over seeds
0055   sort(fOutput.begin(), fOutput.end(), greaterByEtPseudoJet);
0056 }