FastJet 3.0.2
CASubJetTagger.cc
00001 //STARTHEADER
00002 // $Id: CASubJetTagger.cc 2689 2011-11-14 14:51:06Z soyez $
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/CASubJetTagger.hh>
00030 #include <fastjet/ClusterSequence.hh>
00031 
00032 #include <algorithm>
00033 #include <cmath>
00034 #include <sstream>
00035 
00036 using namespace std;
00037 
00038 FASTJET_BEGIN_NAMESPACE
00039 
00040 
00041 LimitedWarning CASubJetTagger::_non_ca_warnings;
00042 
00043 // the tagger's description
00044 //----------------------------------------------------------------------
00045 string CASubJetTagger::description() const{
00046   ostringstream oss;
00047   oss << "CASubJetTagger with z_threshold=" << _z_threshold ;
00048   if (_absolute_z_cut) oss << " (defined wrt original jet)";
00049   oss << " and scale choice ";
00050   switch (_scale_choice) {
00051   case kt2_distance:         oss << "kt2_distance";         break;
00052   case jade_distance:        oss << "jade_distance";        break;
00053   case jade2_distance:       oss << "jade2_distance";       break;
00054   case plain_distance:       oss << "plain_distance";       break;
00055   case mass_drop_distance:   oss << "mass_drop_distance";   break;
00056   case dot_product_distance: oss << "dot_product_distance"; break;
00057   default:
00058     throw Error("unrecognized scale choice");
00059   }
00060 
00061   return oss.str();
00062 }
00063 
00064 // run the tagger on the given cs/jet
00065 // returns the tagged PseudoJet if successful, 0 otherwise
00066 //----------------------------------------------------------------------
00067 PseudoJet CASubJetTagger::result(const PseudoJet & jet) const{
00068   // make sure that the jet results from a Cambridge/Aachen clustering
00069   if (jet.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm)
00070     _non_ca_warnings.warn("CASubJetTagger should only be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk");
00071 
00072   // recurse in the jet to find the max distance
00073   JetAux aux;
00074   aux.jet          = PseudoJet();
00075   aux.aux_distance = -numeric_limits<double>::max();
00076   aux.delta_r      = 0.0;
00077   aux.z            = 1.0;
00078   _recurse_through_jet(jet, aux, jet); // last arg remains original jet
00079 
00080   // create the result and its associated structure
00081   PseudoJet result_local = aux.jet;
00082 
00083   // the tagger is considered to have failed if aux has never been set
00084   // (in which case it will not have parents).
00085   if (result_local == PseudoJet()) return result_local;
00086 
00087   // otherwise sort out the structure
00088   CASubJetTaggerStructure * s = new CASubJetTaggerStructure(result_local);
00089 //  s->_original_jet = jet;
00090   s->_scale_choice = _scale_choice;
00091   s->_distance     = aux.aux_distance;
00092   s->_absolute_z   = _absolute_z_cut;
00093   s->_z            = aux.z;
00094 
00095   result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
00096 
00097   return result_local;
00098 }
00099 
00100 
00101 ///----------------------------------------------------------------------
00102 /// work through the jet, establishing a distance at each branching
00103 inline void CASubJetTagger::_recurse_through_jet(const PseudoJet & jet, JetAux &aux, const PseudoJet & original_jet) const {
00104 
00105   PseudoJet parent1, parent2;
00106   if (! jet.has_parents(parent1, parent2)) return;
00107 
00108   /// make sure the objects are not _too_ close together
00109   if (parent1.squared_distance(parent2) < _dr2_min) return;
00110 
00111   // distance
00112   double dist=0.0;
00113   switch (_scale_choice) {
00114   case kt2_distance:
00115     // a standard (LI) kt distance
00116     dist = parent1.kt_distance(parent2);
00117     break;
00118   case jade_distance:
00119     // something a bit like a mass: pti ptj Delta R_ij^2
00120     dist = parent1.perp()*parent2.perp()*parent1.squared_distance(parent2);
00121     break;
00122   case jade2_distance:
00123     // something a bit like a mass*deltaR^2: pti ptj Delta R_ij^4
00124     dist = parent1.perp()*parent2.perp()*pow(parent1.squared_distance(parent2),2);
00125     break;
00126   case plain_distance:
00127     // Delta R_ij^2
00128     dist = parent1.squared_distance(parent2);
00129     break;
00130   case mass_drop_distance:
00131     // Delta R_ij^2
00132     dist = jet.m() - std::max(parent1.m(),parent2.m());
00133     break;
00134   case dot_product_distance:
00135     // parent1 . parent2 
00136     // ( = jet.m2() - parent1.m2() - parent2.m() in a 
00137     // 4-vector recombination scheme) 
00138     dist = dot_product(parent1, parent2);
00139     break;
00140   default:
00141     throw Error("unrecognized scale choice");
00142   }
00143 
00144   // check the z cut
00145   bool zcut1 = true;
00146   bool zcut2 = true;
00147   double z2 = 0.0;
00148 
00149   // not very efficient -- sort out later
00150   if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2);
00151 
00152   if (_absolute_z_cut) {
00153     z2    = parent2.perp() / original_jet.perp();
00154     zcut1 = parent1.perp() / original_jet.perp() >= _z_threshold;
00155   } else {
00156     z2    = parent2.perp()/(parent1.perp()+parent2.perp());
00157   }
00158   zcut2 = z2 >= _z_threshold;
00159 
00160   if (zcut1 && zcut2){
00161     if (dist > aux.aux_distance){
00162       aux.jet          = jet;
00163       aux.aux_distance = dist;
00164       aux.delta_r      = sqrt(parent1.squared_distance(parent2));
00165       aux.z            = z2; // the softest
00166     }
00167   }    
00168 
00169   if (zcut1) _recurse_through_jet(parent1, aux, original_jet);
00170   if (zcut2) _recurse_through_jet(parent2, aux, original_jet);
00171 }
00172 
00173 FASTJET_END_NAMESPACE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends