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