Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:26 

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // SubjetFilterAlgorithm
0004 // ---------------------
0005 //
0006 // This algorithm implements the Subjet/Filter jet reconstruction
0007 // which was first proposed here: http://arXiv.org/abs/0802.2470
0008 //
0009 // The implementation is largely based on fastjet_boosted_higgs.cc provided
0010 // with the fastjet package.
0011 //
0012 // see: https://twiki.cern.ch/twiki/bin/view/CMS/SWGuideSubjetFilterJetProducer
0013 //
0014 //                       David Lopes-Pegna           <david.lopes-pegna@cern.ch>
0015 //            25/11/2009 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
0016 ////////////////////////////////////////////////////////////////////////////////
0017 
0018 #include "RecoJets/JetAlgorithms/interface/SubjetFilterAlgorithm.h"
0019 
0020 #include <fastjet/ClusterSequenceArea.hh>
0021 
0022 #include <iostream>
0023 #include <iomanip>
0024 #include <sstream>
0025 #include <cmath>
0026 
0027 using namespace std;
0028 
0029 ostream& operator<<(ostream& ostr, fastjet::PseudoJet& jet);
0030 
0031 ////////////////////////////////////////////////////////////////////////////////
0032 // construction / destruction
0033 ////////////////////////////////////////////////////////////////////////////////
0034 
0035 //______________________________________________________________________________
0036 SubjetFilterAlgorithm::SubjetFilterAlgorithm(const std::string& moduleLabel,
0037                                              const std::string& jetAlgorithm,
0038                                              unsigned nFatMax,
0039                                              double rParam,
0040                                              double rFilt,
0041                                              double jetPtMin,
0042                                              double massDropCut,
0043                                              double asymmCut,
0044                                              bool asymmCutLater,
0045                                              bool doAreaFastjet,
0046                                              double ghostEtaMax,
0047                                              int activeAreaRepeats,
0048                                              double ghostArea,
0049                                              bool verbose)
0050     : moduleLabel_(moduleLabel),
0051       jetAlgorithm_(jetAlgorithm),
0052       nFatMax_(nFatMax),
0053       rParam_(rParam),
0054       rFilt_(rFilt),
0055       jetPtMin_(jetPtMin),
0056       massDropCut_(massDropCut),
0057       asymmCut2_(asymmCut * asymmCut),
0058       asymmCutLater_(asymmCutLater),
0059       doAreaFastjet_(doAreaFastjet),
0060       ghostEtaMax_(ghostEtaMax),
0061       activeAreaRepeats_(activeAreaRepeats),
0062       ghostArea_(ghostArea),
0063       verbose_(verbose),
0064       nevents_(0),
0065       ntotal_(0),
0066       nfound_(0),
0067       fjJetDef_(nullptr),
0068       fjAreaDef_(nullptr) {
0069   // FASTJET JET DEFINITION
0070   if (jetAlgorithm == "CambridgeAachen" || jetAlgorithm_ == "ca")
0071     fjJetDef_ = new fastjet::JetDefinition(fastjet::cambridge_algorithm, rParam_);
0072   else if (jetAlgorithm == "AntiKt" || jetAlgorithm_ == "ak")
0073     fjJetDef_ = new fastjet::JetDefinition(fastjet::antikt_algorithm, rParam_);
0074   else if (jetAlgorithm == "Kt" || jetAlgorithm_ == "kt")
0075     fjJetDef_ = new fastjet::JetDefinition(fastjet::kt_algorithm, rParam_);
0076   else
0077     throw cms::Exception("InvalidJetAlgo") << "Jet Algorithm for SubjetFilterAlgorithm is invalid: " << jetAlgorithm_
0078                                            << ", use (ca|CambridgeAachen)|(Kt|kt)|(AntiKt|ak)" << endl;
0079 
0080   // FASTJET AREA DEFINITION
0081   if (doAreaFastjet_) {
0082     fastjet::GhostedAreaSpec ghostSpec(ghostEtaMax_, activeAreaRepeats_, ghostArea_);
0083     fjAreaDef_ = new fastjet::AreaDefinition(fastjet::active_area_explicit_ghosts, ghostSpec);
0084   }
0085 }
0086 
0087 //______________________________________________________________________________
0088 SubjetFilterAlgorithm::~SubjetFilterAlgorithm() {
0089   if (nullptr != fjJetDef_)
0090     delete fjJetDef_;
0091   if (nullptr != fjAreaDef_)
0092     delete fjAreaDef_;
0093 }
0094 
0095 ////////////////////////////////////////////////////////////////////////////////
0096 // implementation of member functions
0097 ////////////////////////////////////////////////////////////////////////////////
0098 
0099 //______________________________________________________________________________
0100 void SubjetFilterAlgorithm::run(const std::vector<fastjet::PseudoJet>& fjInputs,
0101                                 std::vector<CompoundPseudoJet>& fjJets,
0102                                 const edm::EventSetup& iSetup) {
0103   nevents_++;
0104 
0105   if (verbose_)
0106     cout << endl << nevents_ << ". EVENT" << endl;
0107 
0108   fastjet::ClusterSequence* cs = (doAreaFastjet_) ? new fastjet::ClusterSequenceArea(fjInputs, *fjJetDef_, *fjAreaDef_)
0109                                                   : new fastjet::ClusterSequence(fjInputs, *fjJetDef_);
0110 
0111   vector<fastjet::PseudoJet> fjFatJets = fastjet::sorted_by_pt(cs->inclusive_jets(jetPtMin_));
0112 
0113   size_t nFat = (nFatMax_ == 0) ? fjFatJets.size() : std::min(fjFatJets.size(), (size_t)nFatMax_);
0114 
0115   for (size_t iFat = 0; iFat < nFat; iFat++) {
0116     if (verbose_)
0117       cout << endl << iFat << ". FATJET: " << fjFatJets[iFat] << endl;
0118 
0119     fastjet::PseudoJet fjFatJet = fjFatJets[iFat];
0120     fastjet::PseudoJet fjCurrentJet(fjFatJet);
0121     fastjet::PseudoJet fjSubJet1, fjSubJet2;
0122     bool hadSubJets;
0123 
0124     vector<CompoundPseudoSubJet> subJets;
0125 
0126     // FIND SUBJETS PASSING MASSDROP [AND ASYMMETRY] CUT(S)
0127     while ((hadSubJets = cs->has_parents(fjCurrentJet, fjSubJet1, fjSubJet2))) {
0128       if (fjSubJet1.m() < fjSubJet2.m())
0129         swap(fjSubJet1, fjSubJet2);
0130 
0131       if (verbose_)
0132         cout << "SUBJET CANDIDATES: " << fjSubJet1 << endl << "                   " << fjSubJet2 << endl;
0133       if (verbose_)
0134         cout << "md=" << fjSubJet1.m() / fjCurrentJet.m()
0135              << " y=" << fjSubJet1.kt_distance(fjSubJet2) / fjCurrentJet.m2() << endl;
0136 
0137       if (fjSubJet1.m() < massDropCut_ * fjCurrentJet.m() &&
0138           (asymmCutLater_ || fjSubJet1.kt_distance(fjSubJet2) > asymmCut2_ * fjCurrentJet.m2())) {
0139         break;
0140       } else {
0141         fjCurrentJet = fjSubJet1;
0142       }
0143     }
0144 
0145     if (!hadSubJets) {
0146       if (verbose_)
0147         cout << "FAILED TO SELECT SUBJETS" << endl;
0148     }
0149     // FOUND TWO GOOD SUBJETS PASSING MASSDROP CUT
0150     else {
0151       if (verbose_)
0152         cout << "SUBJETS selected" << endl;
0153 
0154       if (asymmCutLater_ && fjSubJet1.kt_distance(fjSubJet2) <= asymmCut2_ * fjCurrentJet.m2()) {
0155         if (verbose_)
0156           cout << "FAILED y-cut" << endl;
0157       }
0158       // PASSED ASYMMETRY (Y) CUT
0159       else {
0160         if (verbose_)
0161           cout << "PASSED y cut" << endl;
0162 
0163         vector<fastjet::PseudoJet> fjFilterJets;
0164         double Rbb = std::sqrt(fjSubJet1.squared_distance(fjSubJet2));
0165         double Rfilt = std::min(0.5 * Rbb, rFilt_);
0166         double dcut = Rfilt * Rfilt / rParam_ / rParam_;
0167         fjFilterJets = fastjet::sorted_by_pt(cs->exclusive_subjets(fjCurrentJet, dcut));
0168 
0169         if (verbose_) {
0170           cout << "Rbb=" << Rbb << ", Rfilt=" << Rfilt << endl;
0171           cout << "FILTER JETS: " << flush;
0172           for (size_t i = 0; i < fjFilterJets.size(); i++) {
0173             if (i > 0)
0174               cout << "             " << flush;
0175             cout << fjFilterJets[i] << endl;
0176           }
0177         }
0178 
0179         vector<fastjet::PseudoJet> fjSubJets;
0180         fjSubJets.push_back(fjSubJet1);
0181         fjSubJets.push_back(fjSubJet2);
0182         size_t nFilter = fjFilterJets.size();
0183         for (size_t iFilter = 0; iFilter < nFilter; iFilter++)
0184           fjSubJets.push_back(fjFilterJets[iFilter]);
0185 
0186         for (size_t iSub = 0; iSub < fjSubJets.size(); iSub++) {
0187           vector<fastjet::PseudoJet> fjConstituents = cs->constituents(fjSubJets[iSub]);
0188           vector<int> constituents;
0189           for (size_t iConst = 0; iConst < fjConstituents.size(); iConst++) {
0190             int userIndex = fjConstituents[iConst].user_index();
0191             if (userIndex >= 0)
0192               constituents.push_back(userIndex);
0193           }
0194 
0195           double subJetArea = (doAreaFastjet_) ? ((fastjet::ClusterSequenceArea*)cs)->area(fjSubJets[iSub]) : 0.0;
0196 
0197           if (iSub < 2 || !constituents.empty())
0198             subJets.push_back(CompoundPseudoSubJet(fjSubJets[iSub], subJetArea, constituents));
0199         }
0200 
0201       }  // PASSED Y CUT
0202 
0203     }  // PASSED MASSDROP CUT
0204 
0205     if (verbose_)
0206       cout << "write fatjet with " << subJets.size() << " sub+filter jets" << endl;
0207 
0208     double fatJetArea = (doAreaFastjet_) ? ((fastjet::ClusterSequenceArea*)cs)->area(fjFatJet) : 0.0;
0209 
0210     fjJets.push_back(CompoundPseudoJet(fjFatJet, fatJetArea, subJets));
0211 
0212     ntotal_++;
0213     if (subJets.size() > 3)
0214       nfound_++;
0215 
0216   }  // LOOP OVER FATJETS
0217 
0218   if (verbose_)
0219     cout << endl << fjJets.size() << " FATJETS written\n" << endl;
0220 
0221   delete cs;
0222 
0223   return;
0224 }
0225 
0226 //______________________________________________________________________________
0227 string SubjetFilterAlgorithm::summary() const {
0228   double eff = (ntotal_ > 0) ? nfound_ / (double)ntotal_ : 0;
0229   std::stringstream ss;
0230   ss << "************************************************************\n"
0231      << "* " << moduleLabel_ << " (SubjetFilterAlgorithm) SUMMARY:\n"
0232      << "************************************************************\n"
0233      << "nevents = " << nevents_ << endl
0234      << "ntotal  = " << ntotal_ << endl
0235      << "nfound  = " << nfound_ << endl
0236      << "eff     = " << eff << endl
0237      << "************************************************************\n";
0238   return ss.str();
0239 }
0240 
0241 /// does the actual work for printing out a jet
0242 ostream& operator<<(ostream& ostr, fastjet::PseudoJet& jet) {
0243   ostr << "pt=" << setw(10) << jet.perp() << " eta=" << setw(6) << jet.eta() << " m=" << setw(10) << jet.m();
0244   return ostr;
0245 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!