JetAlgorithms/src/CMSBoostedTauSeedingAlgorithm.cc

0001 // CMSBoostedTauSeedingAlgorithm Package
0002 //
0003 //----------------------------------------------------------------------
0004 // This file is part of FastJet contrib.
0005 //
0006 // It is free software; you can redistribute it and/or modify it under
0007 // the terms of the GNU General Public License as published by the
0008 // Free Software Foundation; either version 2 of the License, or (at
0009 // your option) any later version.
0010 //
0011 // It is distributed in the hope that it will be useful, but WITHOUT
0012 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
0013 // or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
0014 // License for more details.
0015 //
0016 // You should have received a copy of the GNU General Public License
0017 // along with this code. If not, see <http://www.gnu.org/licenses/>.
0018
0019 //CG: From CMSSW 8_0_X cut on eta subjets <3. Anyways taus are cut at 2.3
0020
0021 //----------------------------------------------------------------------
0022
0023 #include "RecoJets/JetAlgorithms/interface/CMSBoostedTauSeedingAlgorithm.h"
0024
0025 FASTJET_BEGIN_NAMESPACE  // defined in fastjet/internal/base.hh
0026
0027     namespace contrib {
0028   /////////////////////////////
0029   // constructor
0030   CMSBoostedTauSeedingAlgorithm::CMSBoostedTauSeedingAlgorithm(double iminPt,
0031                                                                double iminMassDrop,
0032                                                                double imaxMassDrop,
0033                                                                double iminY,
0034                                                                double imaxY,
0035                                                                double iminDeltaR,
0036                                                                double imaxDeltaR,
0037                                                                int maxDepth,
0038                                                                int verbosity)
0039       : ptMin_(iminPt),
0040         muMin_(iminMassDrop),
0041         muMax_(imaxMassDrop),
0042         yMin_(iminY),
0043         yMax_(imaxY),
0044         dRMin_(iminDeltaR),
0045         dRMax_(imaxDeltaR),
0046         maxDepth_(maxDepth),
0047         verbosity_(verbosity) {}
0048
0049   /////////////////////////////
0050   // description
0051   std::string CMSBoostedTauSeedingAlgorithm::description() const {
0052     std::ostringstream oss;
0053     oss << "CMSBoostedTauSeedingAlgorithm algorithm";
0054     return oss.str();
0055   }
0056
0057   /////////////////////////////
0058   void CMSBoostedTauSeedingAlgorithm::dumpSubJetStructure(
0059       const fastjet::PseudoJet& jet, int depth, int maxDepth, const std::string& depth_and_idx_string) const {
0060     if (maxDepth != -1 && depth > maxDepth)
0061       return;
0062     fastjet::PseudoJet subjet1, subjet2;
0063     bool hasSubjets = jet.has_parents(subjet1, subjet2);
0064     if (!hasSubjets)
0065       return;
0066     std::string depth_and_idx_string_subjet1 = depth_and_idx_string;
0067     if (depth_and_idx_string_subjet1.length() > 0)
0068       depth_and_idx_string_subjet1.append(".");
0069     depth_and_idx_string_subjet1.append("0");
0070     for (int iSpace = 0; iSpace < depth; ++iSpace) {
0071       std::cout << " ";
0072     }
0073     std::cout << " jetConstituent #" << depth_and_idx_string_subjet1 << " (depth = " << depth
0074               << "): Pt = " << subjet1.pt() << ","
0075               << " eta = " << subjet1.eta() << ", phi = " << subjet1.phi() << ", mass = " << subjet1.m()
0076               << " (constituents = " << subjet1.constituents().size() << ")" << std::endl;
0077     dumpSubJetStructure(subjet1, depth + 1, maxDepth, depth_and_idx_string_subjet1);
0078     std::string depth_and_idx_string_subjet2 = depth_and_idx_string;
0079     if (depth_and_idx_string_subjet2.length() > 0)
0080       depth_and_idx_string_subjet2.append(".");
0081     depth_and_idx_string_subjet2.append("1");
0082     for (int iSpace = 0; iSpace < depth; ++iSpace) {
0083       std::cout << " ";
0084     }
0085     std::cout << " jetConstituent #" << depth_and_idx_string_subjet2 << " (depth = " << depth
0086               << "): Pt = " << subjet2.pt() << ","
0087               << " eta = " << subjet2.eta() << ", phi = " << subjet2.phi() << ", mass = " << subjet2.m()
0088               << " (constituents = " << subjet2.constituents().size() << ")" << std::endl;
0089     dumpSubJetStructure(subjet2, depth + 1, maxDepth, depth_and_idx_string_subjet2);
0090     for (int iSpace = 0; iSpace < depth; ++iSpace) {
0091       std::cout << " ";
0092     }
0093     double dR = subjet1.delta_R(subjet2);
0094     std::cout << " (mass-drop @ " << depth_and_idx_string << " = " << std::max(subjet1.m(), subjet2.m()) / jet.m()
0095               << ", dR = " << dR << ")" << std::endl;
0096   }
0097
0098   std::pair<PseudoJet, PseudoJet> CMSBoostedTauSeedingAlgorithm::findSubjets(
0099       const PseudoJet& jet, int depth, bool& subjetsFound) const {
0100     const float etaMax_ = 3.0;  // cut on eta subjets <3. Anyway taus are cut at 2.3
0101     if (verbosity_ >= 2) {
0102       std::cout << "<CMSBoostedTauSeedingAlgorithm::findSubjets>:" << std::endl;
0103       std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi()
0104                 << ", mass = " << jet.m() << std::endl;
0105     }
0106
0107     PseudoJet subjet1, subjet2;
0108     bool hasSubjets = jet.has_parents(subjet1, subjet2);
0109     if (hasSubjets && (maxDepth_ == -1 || depth <= maxDepth_)) {
0110       // make subjet1 the more massive jet
0111       if (subjet1.m2() < subjet2.m2()) {
0112         std::swap(subjet1, subjet2);
0113       }
0114       double dR = subjet1.delta_R(subjet2);
0115       double kT = subjet1.kt_distance(subjet2);
0116       double mu = (jet.m() > 0.) ? sqrt(std::max(subjet1.m2(), subjet2.m2()) / jet.m2()) : -1.;
0117       // check if subjets pass selection required for seeding boosted tau reconstruction
0118       if (subjet1.pt() > ptMin_ && fabs(subjet1.eta()) < etaMax_ && subjet2.pt() > ptMin_ &&
0119           fabs(subjet2.eta()) < etaMax_ && dR > dRMin_ && dR < dRMax_ && mu > muMin_ && mu < muMax_ &&
0120           kT < (yMax_ * jet.m2()) && kT > (yMin_ * jet.m2())) {
0121         subjetsFound = true;
0122         return std::make_pair(subjet1, subjet2);
0123       } else if (subjet1.pt() > ptMin_) {
0124         return findSubjets(subjet1, depth + 1, subjetsFound);
0125       } else if (subjet2.pt() > ptMin_) {
0126         return findSubjets(subjet2, depth + 1, subjetsFound);
0127       }
0128     }
0129     subjetsFound = false;
0130     PseudoJet dummy_subjet1, dummy_subjet2;
0131     return std::make_pair(dummy_subjet1, dummy_subjet2);
0132   }
0133
0134   PseudoJet CMSBoostedTauSeedingAlgorithm::result(const PseudoJet& jet) const {
0135     if (verbosity_ >= 1) {
0136       std::cout << "<CMSBoostedTauSeedingAlgorithm::findSubjets>:" << std::endl;
0137       std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi()
0138                 << ", mass = " << jet.m() << std::endl;
0139     }
0140
0141     if (verbosity_ >= 2) {
0142       dumpSubJetStructure(jet, 0, maxDepth_, "");
0143     }
0144
0145     bool subjetsFound = false;
0146     std::pair<PseudoJet, PseudoJet> subjets = findSubjets(jet, 0, subjetsFound);
0147     if (subjetsFound) {
0148       // fill structure for returning result
0149       PseudoJet subjet1 = subjets.first;
0150       PseudoJet subjet2 = subjets.second;
0151       if (verbosity_ >= 1) {
0152         std::cout << "before recombination:" << std::endl;
0153         std::cout << " subjet #1: Pt = " << subjet1.pt() << ", eta = " << subjet1.eta() << ", phi = " << subjet1.phi()
0154                   << ", mass = " << subjet1.m() << std::endl;
0155         std::cout << " subjet #2: Pt = " << subjet2.pt() << ", eta = " << subjet2.eta() << ", phi = " << subjet2.phi()
0156                   << ", mass = " << subjet2.m() << std::endl;
0157       }
0158
0159       const JetDefinition::Recombiner* rec = jet.associated_cluster_sequence()->jet_def().recombiner();
0160       PseudoJet result_local = join(subjet1, subjet2, *rec);
0161       if (verbosity_ >= 1) {
0162         std::cout << "after recombination:" << std::endl;
0163         std::vector<fastjet::PseudoJet> subjets = result_local.pieces();
0164         int idx_subjet = 0;
0165         for (std::vector<fastjet::PseudoJet>::const_iterator subjet = subjets.begin(); subjet != subjets.end();
0166              ++subjet) {
0167           std::cout << " subjet #" << idx_subjet << ": Pt = " << subjet->pt() << ", eta = " << subjet->eta()
0168                     << ", phi = " << subjet->phi() << ", mass = " << subjet->m()
0169                     << " (#constituents = " << subjet->constituents().size() << ")" << std::endl;
0170           std::vector<fastjet::PseudoJet> constituents = subjet->constituents();
0171           int idx_constituent = 0;
0172           for (std::vector<fastjet::PseudoJet>::const_iterator constituent = constituents.begin();
0173                constituent != constituents.end();
0174                ++constituent) {
0175             if (constituent->pt() < 1.e-3)
0176               continue;  // CV: skip ghosts
0177             std::cout << "  constituent #" << idx_constituent << ": Pt = " << constituent->pt()
0178                       << ", eta = " << constituent->eta() << ", phi = " << constituent->phi() << ","
0179                       << " mass = " << constituent->m() << std::endl;
0180             ++idx_constituent;
0181           }
0182           ++idx_subjet;
0183         }
0184       }
0185
0186       CMSBoostedTauSeedingAlgorithmStructure* s = new CMSBoostedTauSeedingAlgorithmStructure(result_local);
0187       //s->_original_jet = jet;
0188       s->_mu = (jet.m2() > 0.) ? sqrt(std::max(subjet1.m2(), subjet2.m2()) / jet.m2()) : 0.;
0189       s->_y = (jet.m2() > 0.) ? subjet1.kt_distance(subjet2) / jet.m2() : 0.;
0190       s->_dR = subjet1.delta_R(subjet2);
0191       s->_pt = subjet2.pt();
0192
0193       result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
0194
0195       return result_local;
0196     } else {
0197       // no subjets for seeding boosted tau reconstruction found, return an empty PseudoJet
0198       if (verbosity_ >= 1) {
0199         std::cout << "No subjets found." << std::endl;
0200       }
0201       return PseudoJet();
0202     }
0203   }
0204
0205 }  // namespace contrib
0206
0207 FASTJET_END_NAMESPACE