JetDefinition.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: JetDefinition.cc 958 2007-11-28 17:43:48Z cacciari $
00003 //
00004 // Copyright (c) 2005-2007, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
00030 
00031 #include "fastjet/JetDefinition.hh"
00032 #include "fastjet/Error.hh"
00033 #include<sstream>
00034 
00035 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00036 
00037 using namespace std;
00038 
00039 string JetDefinition::description() const {
00040   ostringstream name;
00041   if (jet_algorithm() == plugin_algorithm) {
00042     return plugin()->description();
00043   } else if (jet_algorithm() == kt_algorithm) {
00044     name << "Longitudinally invariant kt algorithm with R = " << R();
00045     name << " and " << recombiner()->description();
00046   } else if (jet_algorithm() == cambridge_algorithm) {
00047     name << "Longitudinally invariant Cambridge/Aachen algorithm with R = " 
00048          << R() ;
00049     name << " and " << recombiner()->description();
00050   } else if (jet_algorithm() == antikt_algorithm) {
00051     name << "Longitudinally invariant anti-kt algorithm with R = " 
00052          << R() ;
00053     name << " and " << recombiner()->description();
00054   } else if (jet_algorithm() == cambridge_for_passive_algorithm) {
00055     name << "Longitudinally invariant Cambridge/Aachen algorithm with R = " 
00056          << R() << "and a special hack whereby particles with kt < " 
00057          << extra_param() << "are treated as passive ghosts";
00058   } else {
00059     throw Error("Unrecognized jet_finder");
00060   }
00061   return name.str();
00062 }
00063 
00064 
00065 void JetDefinition::set_recombination_scheme(
00066                                RecombinationScheme recomb_scheme) {
00067   _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme);
00068   _recombiner = 0;
00069 }
00070 
00071 
00072 string JetDefinition::DefaultRecombiner::description() const {
00073   switch(_recomb_scheme) {
00074   case E_scheme:
00075     return "E scheme recombination";
00076   case pt_scheme:
00077     return "pt scheme recombination";
00078   case pt2_scheme:
00079     return "pt2 scheme recombination";
00080   case Et_scheme:
00081     return "Et scheme recombination";
00082   case Et2_scheme:
00083     return "Et2 scheme recombination";
00084   case BIpt_scheme:
00085     return "boost-invariant pt scheme recombination";
00086   case BIpt2_scheme:
00087     return "boost-invariant pt2 scheme recombination";
00088   default:
00089     ostringstream err;
00090     err << "DefaultRecombiner: unrecognized recombination scheme " 
00091         << _recomb_scheme;
00092     throw Error(err.str());
00093   }
00094 }
00095 
00096 
00097 void JetDefinition::DefaultRecombiner::recombine(
00098            const PseudoJet & pa, const PseudoJet & pb,
00099            PseudoJet & pab) const {
00100   
00101   double weighta, weightb;
00102 
00103   switch(_recomb_scheme) {
00104   case E_scheme:
00105     pab = pa + pb; 
00106     pab.set_user_index(0);
00107     return;
00108   // all remaining schemes are massless recombinations and locally
00109   // we just set weights, while the hard work is done below...
00110   case pt_scheme:
00111   case Et_scheme:
00112   case BIpt_scheme:
00113     weighta = pa.perp(); 
00114     weightb = pb.perp();
00115     break;
00116   case pt2_scheme:
00117   case Et2_scheme:
00118   case BIpt2_scheme:
00119     weighta = pa.perp2(); 
00120     weightb = pb.perp2();
00121     break;
00122   default:
00123     ostringstream err;
00124     err << "DefaultRecombiner: unrecognized recombination scheme " 
00125         << _recomb_scheme;
00126     throw Error(err.str());
00127   }
00128 
00129   double perp_ab = pa.perp() + pb.perp();
00130   if (perp_ab != 0.0) { // weights also non-zero...
00131     double y_ab    = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
00132     
00133     // take care with periodicity in phi...
00134     double phi_a = pa.phi(), phi_b = pb.phi();
00135     if (phi_a - phi_b > pi)  phi_b += twopi;
00136     if (phi_a - phi_b < -pi) phi_b -= twopi;
00137     double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
00138     
00139     pab = PseudoJet(perp_ab*cos(phi_ab),
00140                     perp_ab*sin(phi_ab),
00141                     perp_ab*sinh(y_ab),
00142                     perp_ab*cosh(y_ab));
00143   } else { // weights are zero
00144     pab = PseudoJet(0.0,0.0,0.0,0.0);
00145   }
00146   pab.set_user_index(0);
00147 }
00148 
00149 
00150 void JetDefinition::DefaultRecombiner::preprocess(PseudoJet & p) const {
00151   switch(_recomb_scheme) {
00152   case E_scheme:
00153   case BIpt_scheme:
00154   case BIpt2_scheme:
00155     break;
00156   case pt_scheme:
00157   case pt2_scheme:
00158     {
00159       // these schemes (as in the ktjet implementation) need massless
00160       // initial 4-vectors with essentially E=|p|.
00161       double newE = sqrt(p.perp2()+p.pz()*p.pz());
00162       int    user_index = p.user_index();
00163       p = PseudoJet(p.px(), p.py(), p.pz(), newE);
00164       p.set_user_index(user_index);
00165     }
00166     break;
00167   case Et_scheme:
00168   case Et2_scheme:
00169     {
00170       // these schemes (as in the ktjet implementation) need massless
00171       // initial 4-vectors with essentially E=|p|.
00172       double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
00173       int    user_index = p.user_index();
00174       p = PseudoJet(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
00175       p.set_user_index(user_index);
00176     }
00177     break;
00178   default:
00179     ostringstream err;
00180     err << "DefaultRecombiner: unrecognized recombination scheme " 
00181         << _recomb_scheme;
00182     throw Error(err.str());
00183   }
00184 }
00185 
00186 void JetDefinition::Plugin::set_ghost_separation_scale(double scale) const {
00187       throw Error("set_ghost_separation_scale not supported");
00188 }
00189  
00190 
00191 FASTJET_END_NAMESPACE

Generated on Tue Dec 18 17:05:03 2007 for fastjet by  doxygen 1.5.2