FastJet 3.0.0
JHTopTagger.cc
00001 //STARTHEADER
00002 // $Id: JHTopTagger.cc 2616 2011-09-30 18:03:40Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 #include <fastjet/tools/JHTopTagger.hh>
00030 #include <fastjet/Error.hh>
00031 #include <fastjet/JetDefinition.hh>
00032 #include <fastjet/ClusterSequence.hh>
00033 #include <sstream>
00034 #include <limits>
00035 
00036 FASTJET_BEGIN_NAMESPACE
00037 
00038 using namespace std;
00039 
00040 //----------------------------------------------------------------------
00041 // JHTopTagger class implementation
00042 //----------------------------------------------------------------------
00043 
00044 LimitedWarning JHTopTagger::_warnings_nonca;
00045 
00046 //------------------------------------------------------------------------
00047 // description of the tagger
00048 string JHTopTagger::description() const{ 
00049   ostringstream oss;
00050   oss << "JHTopTagger with delta_p=" << _delta_p << ", delta_r=" << _delta_r
00051       << ", cos_theta_W_max=" << _cos_theta_W_max
00052       << " and mW = " << _mW;
00053   oss << description_of_selectors();
00054   return oss.str();
00055 }
00056 
00057 //------------------------------------------------------------------------
00058 // returns the tagged PseudoJet if successful, 0 otherwise
00059 //  - jet   the PseudoJet to tag
00060 PseudoJet JHTopTagger::result(const PseudoJet & jet) const{
00061   // make sure that there is a "regular" cluster sequence associated
00062   // with the jet. Note that we also check it is valid (to avoid a
00063   // more criptic error later on)
00064   if (!jet.has_valid_cluster_sequence()){
00065     throw Error("JHTopTagger can only be applied on jets having an associated (and valid) ClusterSequence");
00066   }
00067 
00068   // warn if the jet has not been clustered with a Cambridge/Aachen
00069   // algorithm
00070   if (! jet.validated_cs()->jet_def().jet_algorithm() == cambridge_algorithm)
00071     _warnings_nonca.warn("JHTopTagger should only be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk.");
00072 
00073 
00074   // do the first splitting
00075   vector<PseudoJet> split0 = _split_once(jet, jet);
00076   if (split0.size() == 0) return PseudoJet();
00077 
00078   // now try a second splitting on each of the resulting objects
00079   vector<PseudoJet> subjets;
00080   for (unsigned i = 0; i < 2; i++) {
00081     vector<PseudoJet> split1 = _split_once(split0[i], jet);
00082     if (split1.size() > 0) {
00083       subjets.push_back(split1[0]);
00084       subjets.push_back(split1[1]);
00085     } else {
00086       subjets.push_back(split0[i]);
00087     }
00088   }
00089 
00090   // make sure things make sense
00091   if (subjets.size() < 3) return PseudoJet();
00092 
00093   // now find the pair of objects closest in mass to the W
00094   double dmW_min = numeric_limits<double>::max();
00095   int ii=-1, jj=-1;
00096   for (unsigned i = 0 ; i < subjets.size()-1; i++) {
00097     for (unsigned j = i+1 ; j < subjets.size(); j++) {
00098       double dmW = abs(_mW - (subjets[i]+subjets[j]).m());
00099       if (dmW < dmW_min) {
00100         dmW_min = dmW; ii = i; jj = j;
00101       }
00102     }
00103   }
00104 
00105   // order the subjets in the following order:
00106   //  - hardest of the W subjets
00107   //  - softest of the W subjets
00108   //  - hardest of the remaining subjets
00109   //  - softest of the remaining subjets (if any)
00110   if (ii>0) std::swap(subjets[ii], subjets[0]);
00111   if (jj>1) std::swap(subjets[jj], subjets[1]);
00112   if (subjets[0].perp2() < subjets[1].perp2()) std::swap(subjets[0], subjets[1]);
00113   if ((subjets.size()>3) && (subjets[2].perp2() < subjets[3].perp2())) 
00114     std::swap(subjets[2], subjets[3]);
00115   
00116   // create the result and its structure
00117   const JetDefinition::Recombiner *rec
00118     = jet.associated_cluster_sequence()->jet_def().recombiner();
00119 
00120   PseudoJet W = join(subjets[0], subjets[1], *rec);
00121   PseudoJet non_W;
00122   if (subjets.size()>3) {
00123     non_W = join(subjets[2], subjets[3], *rec);
00124   } else {
00125     non_W = join(subjets[2], *rec);
00126   }
00127   PseudoJet result = join<JHTopTaggerStructure>(W, non_W, *rec);
00128   JHTopTaggerStructure *s = (JHTopTaggerStructure*) result.structure_non_const_ptr();
00129   s->_cos_theta_w = _cos_theta_W(result);
00130 
00131   // if the polarisation angle does not pass the cut, consider that
00132   // the tagging has failed
00133   //
00134   // Note that we could perhaps ensure this cut before constructing
00135   // the result structure but this has the advantage that the top
00136   // 4-vector is already available and does not have to de re-computed
00137   if (s->_cos_theta_w >= _cos_theta_W_max ||
00138       ! _top_selector.pass(result) || ! _W_selector.pass(W)
00139       ) {
00140     result *= 0.0;
00141   }
00142 
00143   return result;
00144 
00145   // // old version
00146   // PseudoJet result = join<JHTopTaggerStructure>(subjets, *rec);
00147   // JHTopTaggerStructure *s = (JHTopTaggerStructure*) result.structure_non_const_ptr();
00148   // //  s->_original_jet = jet;
00149   // s->_W = join(subjets[0], subjets[1], *rec);
00150   // if (subjets.size()>3)
00151   //   s->_non_W = join(subjets[2], subjets[3], *rec);
00152   // else
00153   //   s->_non_W = join(subjets[2], *rec);
00154   // s->_cos_theta_w = _cos_theta_W(result);
00155   // 
00156   // // if the polarisation angle does not pass the cut, consider that
00157   // // the tagging has failed
00158   // //
00159   // // Note that we could perhaps ensure this cut before constructing
00160   // // the result structure but this has the advantage that the top
00161   // // 4-vector is already available and does not have to de re-computed
00162   // if (s->_cos_theta_w >= _cos_theta_W_max)
00163   //   return PseudoJet();
00164   // 
00165   // return result;
00166 }
00167 
00168 // runs the Johns Hopkins decomposition procedure
00169 vector<PseudoJet> JHTopTagger::_split_once(const PseudoJet & jet_to_split,
00170                                            const PseudoJet & reference_jet) const{
00171   PseudoJet this_jet = jet_to_split;
00172   PseudoJet p1, p2;
00173   vector<PseudoJet> result;
00174   while (this_jet.has_parents(p1, p2)) {
00175     if (p2.perp2() > p1.perp2()) std::swap(p1,p2); // order with hardness
00176     if (p1.perp() < _delta_p * reference_jet.perp()) break; // harder is too soft wrt original jet
00177     if ( (abs(p2.rap()-p1.rap()) + abs(p2.delta_phi_to(p1))) < _delta_r) break; // distance is too small
00178     if (p2.perp() < _delta_p * reference_jet.perp()) {
00179       this_jet = p1; // softer is too soft wrt original, so ignore it
00180       continue; 
00181     }
00182     //result.push_back(this_jet);
00183     result.push_back(p1);
00184     result.push_back(p2);
00185     break;
00186   }
00187   return result;
00188 }
00189 
00190 
00191 
00192 
00193 FASTJET_END_NAMESPACE
00194 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends