FastJet 3.0.2
|
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