JetAlgorithms/interface/CMSTopTagger.h

0001 // Original JHTopTagger code 2011 Matteo Cacciari, Gregory Soyez,
0002 //  and Gavin Salam.
0003 // Modifications to make it CMSTopTagger 2011 Christopher
0004 //  Vermilion.
0005 //
0006 //----------------------------------------------------------------------
0007 //  This file is free software; you can redistribute it and/or modify
0008 //  it under the terms of the GNU General Public License as published by
0009 //  the Free Software Foundation; either version 3 of the License, or
0010 //  (at your option) any later version.
0011 //
0012 //  This file is distributed in the hope that it will be useful,
0013 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0015 //  GNU General Public License for more details.
0016 //
0017 //  The GNU General Public License is available at
0018 //  http://www.gnu.org/licenses/gpl.html or you can write to the Free Software
0019 //  Foundation, Inc.:
0020 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
0021 //----------------------------------------------------------------------
0022
0023 #ifndef __CMSTOPTAGGER_HH__
0024 #define __CMSTOPTAGGER_HH__
0025
0026 #include <fastjet/tools/JHTopTagger.hh>
0027 #include <fastjet/CompositeJetStructure.hh>
0028 #include <fastjet/LimitedWarning.hh>
0029 #include <fastjet/Error.hh>
0030 #include <fastjet/JetDefinition.hh>
0031 #include <fastjet/ClusterSequence.hh>
0032
0033 #include <sstream>
0034 #include <limits>
0035
0036 FASTJET_BEGIN_NAMESPACE
0037
0038 /// An implementation of the "CMS Top Tagger", as described in CMS-PAS-JME-10-013,
0039 /// which is based on the "Johns Hopkins" top tagger (arXiv:0806.0848; Kaplan,
0040 /// Rehermann, Schwartz and Tweedie).  An earlier version of this tagger was
0041 /// described in CMS-PAS-JME-09-001; both implementations can be used via
0042 /// different values passed to the constructor (see below).
0043 ///
0044 /// This code is based on JHTopTagger.{hh,cc}, part of the FastJet package,
0045 /// written by Matteo Cacciari, Gavin P. Salam and Gregory Soyez, and released
0046 /// under the GNU General Public License.
0047 ///
0048 /// The CMS top tagger is extremely similar to JH; the difference is in the
0049 /// choice of parameters and in how the W is identified.  Accordingly I reuse
0050 /// the JHTopTaggerStructure.
0051 class CMSTopTaggerStructure;
0052
0053 class CMSTopTagger : public TopTaggerBase {
0054 public:
0055   /// The parameters are the following:
0056   ///  \param delta_p        fractional pt cut imposed on the subjets
0057   ///                         (computed as a fraction of the original jet)
0058   ///  \param delta_r        minimum distance between 2 subjets
0059   ///                         (computed as sqrt((y1-y2)^2+(phi1-phi2)^2))
0060   ///  \param A              the actual DeltaR cut is (delta_r - A * pT_child)
0061   ///
0062   /// The default values of these parameters are taken from
0063   /// CMS-PAS-JME-10-013.  For the older tagger described in CMS-PAS-JME-09-001,
0064   /// use delta_p=0.05, delta_r=0.0, A=0.0
0065   CMSTopTagger(double delta_p = 0.05, double delta_r = 0.4, double A = 0.0004);
0066
0067   /// returns a textual description of the tagger
0068   std::string description() const override;
0069
0070   /// runs the tagger on the given jet and
0071   /// returns the tagged PseudoJet if successful, or a PseudoJet==0 otherwise
0072   /// (standard access is through operator()).
0073   ///  \param jet   the PseudoJet to tag
0074   PseudoJet result(const PseudoJet& jet) const override;
0075
0076   // the type of the associated structure
0077   typedef CMSTopTaggerStructure StructureType;
0078
0079 protected:
0080   /// runs the Johns Hopkins decomposition procedure
0081   std::vector<PseudoJet> _split_once(const PseudoJet& jet_to_split, const PseudoJet& reference_jet) const;
0082
0083   /// find the indices corresponding to the minimum mass pairing in subjets
0084   /// only considers the hardest 3
0085   void _find_min_mass(const std::vector<PseudoJet>& subjets, int& i, int& j) const;
0086
0087   double _delta_p, _delta_r, _A;
0088   mutable LimitedWarning _warnings_nonca;
0089 };
0090
0091 /// Basically just a copy of JHTopTaggerStructure, but this way CMSTopTagger can
0092 /// be a friend.
0093 class CMSTopTaggerStructure : public JHTopTaggerStructure {
0094 public:
0095   CMSTopTaggerStructure(const std::vector<PseudoJet>& pieces, const JetDefinition::Recombiner* recombiner = nullptr)
0096       : JHTopTaggerStructure(pieces, recombiner) {}
0097
0098 protected:
0099   friend class CMSTopTagger;
0100 };
0101
0102 //----------------------------------------------------------------------
0103 inline CMSTopTagger::CMSTopTagger(double delta_p, double delta_r, double A)
0104     : _delta_p(delta_p), _delta_r(delta_r), _A(A) {}
0105
0106 //------------------------------------------------------------------------
0107 // description of the tagger
0108 inline std::string CMSTopTagger::description() const {
0109   std::ostringstream oss;
0110   oss << "CMSTopTagger with delta_p=" << _delta_p << ", delta_r=" << _delta_r << ", and A=" << _A;
0111   oss << description_of_selectors();
0112   return oss.str();
0113 }
0114
0115 //------------------------------------------------------------------------
0116 // returns the tagged PseudoJet if successful, 0 otherwise
0117 //  - jet   the PseudoJet to tag
0118 inline PseudoJet CMSTopTagger::result(const PseudoJet& jet) const {
0119   // make sure that there is a "regular" cluster sequence associated
0120   // with the jet. Note that we also check it is valid (to avoid a
0121   // more criptic error later on)
0122   if (!jet.has_valid_cluster_sequence()) {
0123     throw Error("CMSTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
0124   }
0125
0126   // warn if the jet has not been clustered with a Cambridge/Aachen
0127   // algorithm
0128   if (jet.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm)
0129     _warnings_nonca.warn(
0130         "CMSTopTagger should only be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms "
0131         "at your own risk.");
0132
0133   // do the first splitting
0134   std::vector<PseudoJet> split0 = _split_once(jet, jet);
0135   if (split0.empty())
0136     return PseudoJet();
0137
0138   // now try a second splitting on each of the resulting objects
0139   std::vector<PseudoJet> subjets;
0140   for (unsigned i = 0; i < 2; i++) {
0141     std::vector<PseudoJet> split1 = _split_once(split0[i], jet);
0142     if (!split1.empty()) {
0143       subjets.push_back(split1[0]);
0144       subjets.push_back(split1[1]);
0145     } else {
0146       subjets.push_back(split0[i]);
0147     }
0148   }
0149
0150   // make sure things make sense
0151   if (subjets.size() < 3)
0152     return PseudoJet();
0153
0154   // now find the pair of subjets with minimum mass (only taking hardest three)
0155   int ii = -1, jj = -1;
0156   _find_min_mass(subjets, ii, jj);
0157
0158   // order the subjets in the following order:
0159   //  - hardest of the W subjets
0160   //  - softest of the W subjets
0161   //  - hardest of the remaining subjets
0162   //  - softest of the remaining subjets (if any)
0163   if (ii > 0)
0164     std::swap(subjets[ii], subjets[0]);
0165   if (jj > 1)
0166     std::swap(subjets[jj], subjets[1]);
0167   if (subjets[0].perp2() < subjets[1].perp2())
0168     std::swap(subjets[0], subjets[1]);
0169   if ((subjets.size() > 3) && (subjets[2].perp2() < subjets[3].perp2()))
0170     std::swap(subjets[2], subjets[3]);
0171
0172   // create the result and its structure
0173   const JetDefinition::Recombiner* rec = jet.associated_cluster_sequence()->jet_def().recombiner();
0174
0175   PseudoJet W = join(subjets[0], subjets[1], *rec);
0176   PseudoJet non_W;
0177   if (subjets.size() > 3) {
0178     non_W = join(subjets[2], subjets[3], *rec);
0179   } else {
0180     non_W = join(subjets[2], *rec);
0181   }
0182   PseudoJet result = join<CMSTopTaggerStructure>(W, non_W, *rec);
0183   CMSTopTaggerStructure* s = (CMSTopTaggerStructure*)result.structure_non_const_ptr();
0184   s->_cos_theta_w = _cos_theta_W(result);
0185
0186   // Note that we could perhaps ensure this cut before constructing
0187   // the result structure but this has the advantage that the top
0188   // 4-vector is already available and does not have to de re-computed
0189   if (!_top_selector.pass(result) || !_W_selector.pass(W)) {
0190     result *= 0.0;
0191   }
0192
0193   result = join(
0194       subjets);  //Added by J. Pilot to combine the (up to 4) subjets identified in the decomposition instead of just the W and non_W components
0195
0196   return result;
0197 }
0198
0199 // runs the Johns Hopkins decomposition procedure
0200 inline std::vector<PseudoJet> CMSTopTagger::_split_once(const PseudoJet& jet_to_split,
0201                                                         const PseudoJet& reference_jet) const {
0202   PseudoJet this_jet = jet_to_split;
0203   PseudoJet p1, p2;
0204   std::vector<PseudoJet> result;
0205   while (this_jet.has_parents(p1, p2)) {
0206     if (p2.perp2() > p1.perp2())
0207       std::swap(p1, p2);  // order with hardness
0208     if (p1.perp() < _delta_p * reference_jet.perp())
0209       break;  // harder is too soft wrt original jet
0210     double DR = p1.delta_R(p2);
0211     if (DR < _delta_r - _A * this_jet.perp())
0212       break;  // distance is too small
0213     if (p2.perp() < _delta_p * reference_jet.perp()) {
0214       this_jet = p1;  // softer is too soft wrt original, so ignore it
0215       continue;
0216     }
0217     //result.push_back(this_jet);
0218     result.push_back(p1);
0219     result.push_back(p2);
0220     break;
0221   }
0222   return result;
0223 }
0224
0225 // find the indices corresponding to the minimum mass pairing in subjets
0226 inline void CMSTopTagger::_find_min_mass(const std::vector<PseudoJet>& subjets, int& i, int& j) const {
0227   assert(subjets.size() > 1 && subjets.size() < 5);
0228
0229   // if four subjets found, only consider three hardest
0230   unsigned softest = 5;  // invalid value
0231   if (subjets.size() == 4) {
0232     double min_pt = std::numeric_limits<double>::max();
0233     ;
0234     for (unsigned ii = 0; ii < subjets.size(); ++ii) {
0235       if (subjets[ii].perp() < min_pt) {
0236         min_pt = subjets[ii].perp();
0237         softest = ii;
0238       }
0239     }
0240   }
0241
0242   double min_mass = std::numeric_limits<double>::max();
0243   for (unsigned ii = 0; ii + 1 < subjets.size(); ++ii) {  // don't do size()-1: (unsigned(0)-1 != -1) !!
0244     if (ii == softest)
0245       continue;
0246     for (unsigned jj = ii + 1; jj < subjets.size(); ++jj) {
0247       if (jj == softest)
0248         continue;
0249       if ((subjets[ii] + subjets[jj]).m() < min_mass) {
0250         min_mass = (subjets[ii] + subjets[jj]).m();
0251         i = ii;
0252         j = jj;
0253       }
0254     }
0255   }
0256 }
0257
0258 FASTJET_END_NAMESPACE
0259
0260 #endif  // __CMSTOPTAGGER_HH__