FastJet 3.0.5
GhostedAreaSpec.cc
00001 //STARTHEADER
00002 // $Id: GhostedAreaSpec.cc 2728 2011-11-20 14:18:59Z 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/GhostedAreaSpec.hh"
00030 #include "fastjet/Error.hh"
00031 #include<iostream>
00032 #include<sstream>
00033 
00034 using namespace std;
00035 
00036 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00037 
00038 BasicRandom<double> GhostedAreaSpec::_random_generator;
00039 LimitedWarning GhostedAreaSpec::_warn_fj2_placement_deprecated;
00040 
00041 /// explicit constructor
00042 GhostedAreaSpec::GhostedAreaSpec(
00043                            const Selector & selector,
00044                            int    repeat_in        ,
00045                            double ghost_area_in    ,   
00046                            double grid_scatter_in  , 
00047                            double pt_scatter_in    ,   
00048                            double mean_ghost_pt_in 
00049                           ): 
00050     _repeat(repeat_in), 
00051     _ghost_area(ghost_area_in), 
00052     _grid_scatter(grid_scatter_in),  
00053     _pt_scatter(pt_scatter_in), 
00054     _mean_ghost_pt(mean_ghost_pt_in),
00055     _fj2_placement(false),
00056     _selector(selector),
00057     _actual_ghost_area(-1.0)
00058   {
00059     // check the selector has the properties needed -- an area and
00060     // applicability jet-by-jet (the latter follows automatically from
00061     // the former?)
00062     if (!_selector.has_finite_area()) throw Error("To construct a GhostedAreaSpec with a Selector, the selector must have a finite area");
00063     if (!_selector.applies_jet_by_jet()) throw Error("To construct a GhostedAreaSpec with a Selector, the selector must apply jet-by-jet");
00064     // get the internal rapidity extent from the selector
00065     double ghost_maxrap_local, ghost_minrap_local;
00066     _selector.get_rapidity_extent(ghost_minrap_local, ghost_maxrap_local);
00067     _ghost_maxrap     = 0.5*(ghost_maxrap_local - ghost_minrap_local); 
00068     _ghost_rap_offset = 0.5*(ghost_maxrap_local + ghost_minrap_local);
00069     
00070     _initialize();
00071   
00072 }
00073 
00074 //======================================================================
00075 // sets fj2 ghost placement
00076 void GhostedAreaSpec::set_fj2_placement(bool val) {
00077   _fj2_placement  = val; _initialize();
00078   if (val) _warn_fj2_placement_deprecated.warn("FJ2 placement of ghosts can lead to systematic edge effects in area evaluation and is deprecated. Prefer new (default) FJ3 placement.");
00079 }
00080 
00081 //======================================================================
00082 /// sets the detailed parameters for the ghosts (which may not be quite
00083 /// the same as those requested -- this is in order for things to fit
00084 /// in nicely into 2pi etc...
00085 void GhostedAreaSpec::_initialize() {
00086   // add on area-measuring dummy particles
00087   _drap = sqrt(_ghost_area);
00088   _dphi = _drap;
00089   if (_fj2_placement) {
00090     _nphi = int(ceil(twopi/_dphi)); _dphi = twopi/_nphi;
00091     _nrap = int(ceil(_ghost_maxrap/_drap)); _drap = _ghost_maxrap / _nrap;
00092     _actual_ghost_area = _dphi * _drap;
00093     _n_ghosts   = (2*_nrap+1)*_nphi;
00094   } else {
00095     // for FJ3, update the ghost placement as follows
00096     // - use nearest int rather than ceiling in determining number of
00097     //   phi and rapidity locations, because this is more stable when
00098     //   the user is trying to get an exact number based on the area
00099     // - rather than placing ghosts up to maximum rapidity
00100     _nphi = int(twopi/_dphi + 0.5); _dphi = twopi/_nphi;
00101     _nrap = int(_ghost_maxrap/_drap + 0.5); _drap = _ghost_maxrap / _nrap;
00102     _actual_ghost_area = _dphi * _drap;
00103     _n_ghosts   = (2*_nrap)*_nphi;
00104   }
00105   // checkpoint the status of the random number generator.
00106   checkpoint_random();
00107   //_random_generator.info(cerr);
00108 }
00109 
00110 //----------------------------------------------------------------------
00111 /// adds the ghost 4-momenta to the vector of PseudoJet's
00112 void GhostedAreaSpec::add_ghosts(vector<PseudoJet> & event) const {
00113 
00114   double rap_offset;
00115   int nrap_upper;
00116   if (_fj2_placement) {
00117     rap_offset  = 0.0;
00118     nrap_upper  = _nrap;
00119   } else {
00120     rap_offset  = 0.5;
00121     nrap_upper  = _nrap-1;
00122   }
00123 
00124   // add momenta for ghosts
00125   for (int irap = -_nrap; irap <= nrap_upper; irap++) {
00126     for (int iphi = 0; iphi < _nphi; iphi++) {
00127      
00128       // include random offsets for all quantities
00129       //----------------------------------------------
00130       // NB: in FJ2 we'd exchanged the px and py components relative to a
00131       // standard definition of phi; to preserve the same areas as fj2
00132       // we now generate a "phi_fj2", and then convert to a standard phi
00133       double phi_fj2 = (iphi+0.5) * _dphi + _dphi*(_our_rand()-0.5)*_grid_scatter;
00134       double phi;
00135       if (_fj2_placement) phi = 0.5*pi - phi_fj2;
00136       else                phi = phi_fj2;
00137       double rap = (irap+rap_offset) * _drap + _drap*(_our_rand()-0.5)*_grid_scatter
00138                                                          + _ghost_rap_offset ;
00139       double pt = _mean_ghost_pt*(1+(_our_rand()-0.5)*_pt_scatter);
00140 
00141       double exprap = exp(+rap);
00142       double pminus = pt/exprap;
00143       double pplus  = pt*exprap;
00144       double px = pt*cos(phi);
00145       double py = pt*sin(phi);
00146       PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00147       // this call fills in the PseudoJet's cached rap,phi information,
00148       // based on pre-existing knowledge. Watch out: if you get the hint
00149       // wrong nobody will tell you, but you will certainly mess up
00150       // your results.
00151       mom.set_cached_rap_phi(rap,phi);
00152 
00153       // if we have an active selector and the particle does not pass the 
00154       // selection condition, move on to the next momentum
00155       if (_selector.worker().get() && !_selector.pass(mom)) continue;
00156       event.push_back(mom);
00157     }
00158   }
00159 }
00160 
00161 string GhostedAreaSpec::description() const {
00162 
00163   ostringstream ostr;
00164   ostr << "ghosts of area " << actual_ghost_area() 
00165        << " (had requested " << ghost_area() << ")";
00166   if (_selector.worker().get()) 
00167     ostr << ", placed according to selector (" << _selector.description() << ")";
00168   else
00169     ostr << ", placed up to y = " << ghost_maxrap() ;
00170   ostr << ", scattered wrt to perfect grid by (rel) " << grid_scatter() 
00171        << ", mean_ghost_pt = " << mean_ghost_pt()
00172        << ", rel pt_scatter =  " << pt_scatter()
00173        << ", n repetitions of ghost distributions =  " << repeat();
00174   return ostr.str();
00175 }
00176 
00177 FASTJET_END_NAMESPACE
00178 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends