FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: JHTopTagger.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/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_local = join<JHTopTaggerStructure>(W, non_W, *rec); 00128 JHTopTaggerStructure *s = (JHTopTaggerStructure*) result_local.structure_non_const_ptr(); 00129 s->_cos_theta_w = _cos_theta_W(result_local); 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_local) || ! _W_selector.pass(W) 00139 ) { 00140 result_local *= 0.0; 00141 } 00142 00143 return result_local; 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_local; 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_local.push_back(p1); 00184 result_local.push_back(p2); 00185 break; 00186 } 00187 return result_local; 00188 } 00189 00190 00191 00192 00193 FASTJET_END_NAMESPACE 00194