FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: JetDefinition.cc 2712 2011-11-16 22:38:44Z 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/JetDefinition.hh" 00030 #include "fastjet/Error.hh" 00031 #include "fastjet/CompositeJetStructure.hh" 00032 #include<sstream> 00033 00034 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00035 00036 using namespace std; 00037 00038 const double JetDefinition::max_allowable_R = 1000.0; 00039 00040 //---------------------------------------------------------------------- 00041 // [NB: implementation was getting complex, so in 2.4-devel moved it 00042 // from .hh to .cc] 00043 JetDefinition::JetDefinition(JetAlgorithm jet_algorithm_in, 00044 double R_in, 00045 Strategy strategy_in, 00046 RecombinationScheme recomb_scheme_in, 00047 int nparameters) : 00048 _jet_algorithm(jet_algorithm_in), _Rparam(R_in), _strategy(strategy_in) { 00049 00050 // set R parameter or ensure its sensibleness, as appropriate 00051 if (_jet_algorithm == ee_kt_algorithm) { 00052 _Rparam = 4.0; // introduce a fictional R that ensures that 00053 // our clustering sequence will not produce 00054 // "beam" jets except when only a single particle remains. 00055 // Any value > 2 would have done here 00056 } else { 00057 // We maintain some limit on R because particles with pt=0, m=0 00058 // can have rapidities O(100000) and one doesn't want the 00059 // clustering to start including them as if their rapidities were 00060 // physical. 00061 if (R_in > max_allowable_R) { 00062 ostringstream oss; 00063 oss << "Requested R = " << R_in << " for jet definition is larger than max_allowable_R = " << max_allowable_R; 00064 throw Error(oss.str()); 00065 } 00066 } 00067 00068 // cross-check the number of parameters that were declared in setting up the 00069 // algorithm (passed internally from the public constructors) 00070 switch (jet_algorithm_in) { 00071 case ee_kt_algorithm: 00072 if (nparameters != 0) { 00073 ostringstream oss; 00074 oss << "ee_kt_algorithm should be constructed with 0 parameters but was called with " 00075 << nparameters << " parameter(s)\n"; 00076 throw Error(oss.str()); 00077 } 00078 break; 00079 case genkt_algorithm: 00080 case ee_genkt_algorithm: 00081 if (nparameters != 2) { 00082 ostringstream oss; 00083 oss << "(ee_)genkt_algorithm should be constructed with 2 parameters but was called with " 00084 << nparameters << " parameter(s)\n"; 00085 throw Error(oss.str()); 00086 } 00087 break; 00088 default: 00089 if (nparameters != 1) { 00090 ostringstream oss; 00091 oss << "The jet algorithm you requested (" 00092 << jet_algorithm_in << ") should be constructed with 1 parameter but was called with " 00093 << nparameters << " parameter(s)\n"; 00094 throw Error(oss.str()); 00095 } 00096 } 00097 00098 // make sure the strategy requested is sensible 00099 assert (_strategy != plugin_strategy); 00100 00101 _plugin = NULL; 00102 set_recombination_scheme(recomb_scheme_in); 00103 set_extra_param(0.0); // make sure it's defined 00104 } 00105 00106 00107 //---------------------------------------------------------------------- 00108 string JetDefinition::description() const { 00109 ostringstream name; 00110 if (jet_algorithm() == plugin_algorithm) { 00111 return plugin()->description(); 00112 } else if (jet_algorithm() == kt_algorithm) { 00113 name << "Longitudinally invariant kt algorithm with R = " << R(); 00114 name << " and " << recombiner()->description(); 00115 } else if (jet_algorithm() == cambridge_algorithm) { 00116 name << "Longitudinally invariant Cambridge/Aachen algorithm with R = " 00117 << R() ; 00118 name << " and " << recombiner()->description(); 00119 } else if (jet_algorithm() == antikt_algorithm) { 00120 name << "Longitudinally invariant anti-kt algorithm with R = " 00121 << R() ; 00122 name << " and " << recombiner()->description(); 00123 } else if (jet_algorithm() == genkt_algorithm) { 00124 name << "Longitudinally invariant generalised kt algorithm with R = " 00125 << R() << ", p = " << extra_param(); 00126 name << " and " << recombiner()->description(); 00127 } else if (jet_algorithm() == cambridge_for_passive_algorithm) { 00128 name << "Longitudinally invariant Cambridge/Aachen algorithm with R = " 00129 << R() << "and a special hack whereby particles with kt < " 00130 << extra_param() << "are treated as passive ghosts"; 00131 } else if (jet_algorithm() == ee_kt_algorithm) { 00132 name << "e+e- kt (Durham) algorithm (NB: no R)"; 00133 name << " with " << recombiner()->description(); 00134 } else if (jet_algorithm() == ee_genkt_algorithm) { 00135 name << "e+e- generalised kt algorithm with R = " 00136 << R() << ", p = " << extra_param(); 00137 name << " and " << recombiner()->description(); 00138 } else if (jet_algorithm() == undefined_jet_algorithm) { 00139 name << "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ; 00140 } else { 00141 throw Error("JetDefinition::description(): unrecognized jet_algorithm"); 00142 } 00143 return name.str(); 00144 } 00145 00146 00147 void JetDefinition::set_recombination_scheme( 00148 RecombinationScheme recomb_scheme) { 00149 _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme); 00150 00151 // do not forget to delete the existing recombiner if needed 00152 if (_recombiner_shared()) _recombiner_shared.reset(); 00153 00154 _recombiner = 0; 00155 } 00156 00157 00158 // returns true if the current jet definitions shares the same 00159 // recombiner as teh one passed as an argument 00160 bool JetDefinition::has_same_recombiner(const JetDefinition &other_jd) const{ 00161 // first make sure that they have the same recombination scheme 00162 const RecombinationScheme & scheme = recombination_scheme(); 00163 if (other_jd.recombination_scheme() != scheme) return false; 00164 00165 // if the scheme is "external", also check that they ahve the same 00166 // recombiner 00167 return (scheme != external_scheme) 00168 || (recombiner() == other_jd.recombiner()); 00169 } 00170 00171 /// allows to let the JetDefinition handle the deletion of the 00172 /// recombiner when it is no longer used 00173 void JetDefinition::delete_recombiner_when_unused(){ 00174 if (_recombiner == 0){ 00175 throw Error("tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme"); 00176 } 00177 00178 _recombiner_shared.reset(_recombiner); 00179 } 00180 00181 /// allows to let the JetDefinition handle the deletion of the 00182 /// plugin when it is no longer used 00183 void JetDefinition::delete_plugin_when_unused(){ 00184 if (_plugin == 0){ 00185 throw Error("tried to call JetDefinition::delete_plugin_when_unused() for a JetDefinition without a plugin"); 00186 } 00187 00188 _plugin_shared.reset(_plugin); 00189 } 00190 00191 00192 00193 string JetDefinition::DefaultRecombiner::description() const { 00194 switch(_recomb_scheme) { 00195 case E_scheme: 00196 return "E scheme recombination"; 00197 case pt_scheme: 00198 return "pt scheme recombination"; 00199 case pt2_scheme: 00200 return "pt2 scheme recombination"; 00201 case Et_scheme: 00202 return "Et scheme recombination"; 00203 case Et2_scheme: 00204 return "Et2 scheme recombination"; 00205 case BIpt_scheme: 00206 return "boost-invariant pt scheme recombination"; 00207 case BIpt2_scheme: 00208 return "boost-invariant pt2 scheme recombination"; 00209 default: 00210 ostringstream err; 00211 err << "DefaultRecombiner: unrecognized recombination scheme " 00212 << _recomb_scheme; 00213 throw Error(err.str()); 00214 } 00215 } 00216 00217 00218 void JetDefinition::DefaultRecombiner::recombine( 00219 const PseudoJet & pa, const PseudoJet & pb, 00220 PseudoJet & pab) const { 00221 00222 double weighta, weightb; 00223 00224 switch(_recomb_scheme) { 00225 case E_scheme: 00226 // a call to reset turns out to be somewhat more efficient 00227 // than a sum and assignment 00228 //pab = pa + pb; 00229 pab.reset(pa.px()+pb.px(), 00230 pa.py()+pb.py(), 00231 pa.pz()+pb.pz(), 00232 pa.E ()+pb.E ()); 00233 return; 00234 // all remaining schemes are massless recombinations and locally 00235 // we just set weights, while the hard work is done below... 00236 case pt_scheme: 00237 case Et_scheme: 00238 case BIpt_scheme: 00239 weighta = pa.perp(); 00240 weightb = pb.perp(); 00241 break; 00242 case pt2_scheme: 00243 case Et2_scheme: 00244 case BIpt2_scheme: 00245 weighta = pa.perp2(); 00246 weightb = pb.perp2(); 00247 break; 00248 default: 00249 ostringstream err; 00250 err << "DefaultRecombiner: unrecognized recombination scheme " 00251 << _recomb_scheme; 00252 throw Error(err.str()); 00253 } 00254 00255 double perp_ab = pa.perp() + pb.perp(); 00256 if (perp_ab != 0.0) { // weights also non-zero... 00257 double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb); 00258 00259 // take care with periodicity in phi... 00260 double phi_a = pa.phi(), phi_b = pb.phi(); 00261 if (phi_a - phi_b > pi) phi_b += twopi; 00262 if (phi_a - phi_b < -pi) phi_b -= twopi; 00263 double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb); 00264 00265 // this is much more efficient... 00266 pab.reset_PtYPhiM(perp_ab,y_ab,phi_ab); 00267 // pab = PseudoJet(perp_ab*cos(phi_ab), 00268 // perp_ab*sin(phi_ab), 00269 // perp_ab*sinh(y_ab), 00270 // perp_ab*cosh(y_ab)); 00271 } else { // weights are zero 00272 //pab = PseudoJet(0.0,0.0,0.0,0.0); 00273 pab.reset(0.0, 0.0, 0.0, 0.0); 00274 } 00275 } 00276 00277 00278 void JetDefinition::DefaultRecombiner::preprocess(PseudoJet & p) const { 00279 switch(_recomb_scheme) { 00280 case E_scheme: 00281 case BIpt_scheme: 00282 case BIpt2_scheme: 00283 break; 00284 case pt_scheme: 00285 case pt2_scheme: 00286 { 00287 // these schemes (as in the ktjet implementation) need massless 00288 // initial 4-vectors with essentially E=|p|. 00289 double newE = sqrt(p.perp2()+p.pz()*p.pz()); 00290 p.reset_momentum(p.px(), p.py(), p.pz(), newE); 00291 // FJ2.x version 00292 // int user_index = p.user_index(); 00293 // p = PseudoJet(p.px(), p.py(), p.pz(), newE); 00294 // p.set_user_index(user_index); 00295 } 00296 break; 00297 case Et_scheme: 00298 case Et2_scheme: 00299 { 00300 // these schemes (as in the ktjet implementation) need massless 00301 // initial 4-vectors with essentially E=|p|. 00302 double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz()); 00303 p.reset_momentum(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E()); 00304 // FJ2.x version 00305 // int user_index = p.user_index(); 00306 // p = PseudoJet(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E()); 00307 // p.set_user_index(user_index); 00308 } 00309 break; 00310 default: 00311 ostringstream err; 00312 err << "DefaultRecombiner: unrecognized recombination scheme " 00313 << _recomb_scheme; 00314 throw Error(err.str()); 00315 } 00316 } 00317 00318 void JetDefinition::Plugin::set_ghost_separation_scale(double /*scale*/) const { 00319 throw Error("set_ghost_separation_scale not supported"); 00320 } 00321 00322 00323 00324 //------------------------------------------------------------------------------- 00325 // helper functions to build a jet made of pieces 00326 // 00327 // This is the extended version with support for a user-defined 00328 // recombination-scheme 00329 // ------------------------------------------------------------------------------- 00330 00331 // build a "CompositeJet" from the vector of its pieces 00332 // 00333 // the user passes the reciombination scheme used to "sum" the pieces. 00334 PseudoJet join(const vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner){ 00335 // compute the total momentum 00336 //-------------------------------------------------- 00337 PseudoJet result; // automatically initialised to 0 00338 if (pieces.size()>0){ 00339 result = pieces[0]; 00340 for (unsigned int i=1; i<pieces.size(); i++) 00341 recombiner.plus_equal(result, pieces[i]); 00342 } 00343 00344 // attach a CompositeJetStructure to the result 00345 //-------------------------------------------------- 00346 CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces, &recombiner); 00347 00348 result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct)); 00349 00350 return result; 00351 } 00352 00353 // build a "CompositeJet" from a single PseudoJet 00354 PseudoJet join(const PseudoJet & j1, 00355 const JetDefinition::Recombiner & recombiner){ 00356 return join(vector<PseudoJet>(1,j1), recombiner); 00357 } 00358 00359 // build a "CompositeJet" from two PseudoJet 00360 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, 00361 const JetDefinition::Recombiner & recombiner){ 00362 vector<PseudoJet> pieces; 00363 pieces.push_back(j1); 00364 pieces.push_back(j2); 00365 return join(pieces, recombiner); 00366 } 00367 00368 // build a "CompositeJet" from 3 PseudoJet 00369 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, 00370 const JetDefinition::Recombiner & recombiner){ 00371 vector<PseudoJet> pieces; 00372 pieces.push_back(j1); 00373 pieces.push_back(j2); 00374 pieces.push_back(j3); 00375 return join(pieces, recombiner); 00376 } 00377 00378 // build a "CompositeJet" from 4 PseudoJet 00379 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4, 00380 const JetDefinition::Recombiner & recombiner){ 00381 vector<PseudoJet> pieces; 00382 pieces.push_back(j1); 00383 pieces.push_back(j2); 00384 pieces.push_back(j3); 00385 pieces.push_back(j4); 00386 return join(pieces, recombiner); 00387 } 00388 00389 00390 00391 00392 FASTJET_END_NAMESPACE