FastJet 3.0.4
JetDefinition.cc
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends