FastJet 3.0.2
|
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