FastJet  3.0.6
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
GhostedAreaSpec.cc
1 //STARTHEADER
2 // $Id: GhostedAreaSpec.cc 2728 2011-11-20 14:18:59Z salam $
3 //
4 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #include "fastjet/GhostedAreaSpec.hh"
30 #include "fastjet/Error.hh"
31 #include<iostream>
32 #include<sstream>
33 
34 using namespace std;
35 
36 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37 
38 BasicRandom<double> GhostedAreaSpec::_random_generator;
39 LimitedWarning GhostedAreaSpec::_warn_fj2_placement_deprecated;
40 
41 /// explicit constructor
42 GhostedAreaSpec::GhostedAreaSpec(
43  const Selector & selector,
44  int repeat_in ,
45  double ghost_area_in ,
46  double grid_scatter_in ,
47  double pt_scatter_in ,
48  double mean_ghost_pt_in
49  ):
50  _repeat(repeat_in),
51  _ghost_area(ghost_area_in),
52  _grid_scatter(grid_scatter_in),
53  _pt_scatter(pt_scatter_in),
54  _mean_ghost_pt(mean_ghost_pt_in),
55  _fj2_placement(false),
56  _selector(selector),
57  _actual_ghost_area(-1.0)
58  {
59  // check the selector has the properties needed -- an area and
60  // applicability jet-by-jet (the latter follows automatically from
61  // the former?)
62  if (!_selector.has_finite_area()) throw Error("To construct a GhostedAreaSpec with a Selector, the selector must have a finite area");
63  if (!_selector.applies_jet_by_jet()) throw Error("To construct a GhostedAreaSpec with a Selector, the selector must apply jet-by-jet");
64  // get the internal rapidity extent from the selector
65  double ghost_maxrap_local, ghost_minrap_local;
66  _selector.get_rapidity_extent(ghost_minrap_local, ghost_maxrap_local);
67  _ghost_maxrap = 0.5*(ghost_maxrap_local - ghost_minrap_local);
68  _ghost_rap_offset = 0.5*(ghost_maxrap_local + ghost_minrap_local);
69 
70  _initialize();
71 
72 }
73 
74 //======================================================================
75 // sets fj2 ghost placement
77  _fj2_placement = val; _initialize();
78  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.");
79 }
80 
81 //======================================================================
82 /// sets the detailed parameters for the ghosts (which may not be quite
83 /// the same as those requested -- this is in order for things to fit
84 /// in nicely into 2pi etc...
86  // add on area-measuring dummy particles
87  _drap = sqrt(_ghost_area);
88  _dphi = _drap;
89  if (_fj2_placement) {
90  _nphi = int(ceil(twopi/_dphi)); _dphi = twopi/_nphi;
91  _nrap = int(ceil(_ghost_maxrap/_drap)); _drap = _ghost_maxrap / _nrap;
92  _actual_ghost_area = _dphi * _drap;
93  _n_ghosts = (2*_nrap+1)*_nphi;
94  } else {
95  // for FJ3, update the ghost placement as follows
96  // - use nearest int rather than ceiling in determining number of
97  // phi and rapidity locations, because this is more stable when
98  // the user is trying to get an exact number based on the area
99  // - rather than placing ghosts up to maximum rapidity
100  _nphi = int(twopi/_dphi + 0.5); _dphi = twopi/_nphi;
101  _nrap = int(_ghost_maxrap/_drap + 0.5); _drap = _ghost_maxrap / _nrap;
102  _actual_ghost_area = _dphi * _drap;
103  _n_ghosts = (2*_nrap)*_nphi;
104  }
105  // checkpoint the status of the random number generator.
106  checkpoint_random();
107  //_random_generator.info(cerr);
108 }
109 
110 //----------------------------------------------------------------------
111 /// adds the ghost 4-momenta to the vector of PseudoJet's
112 void GhostedAreaSpec::add_ghosts(vector<PseudoJet> & event) const {
113 
114  double rap_offset;
115  int nrap_upper;
116  if (_fj2_placement) {
117  rap_offset = 0.0;
118  nrap_upper = _nrap;
119  } else {
120  rap_offset = 0.5;
121  nrap_upper = _nrap-1;
122  }
123 
124  // add momenta for ghosts
125  for (int irap = -_nrap; irap <= nrap_upper; irap++) {
126  for (int iphi = 0; iphi < _nphi; iphi++) {
127 
128  // include random offsets for all quantities
129  //----------------------------------------------
130  // NB: in FJ2 we'd exchanged the px and py components relative to a
131  // standard definition of phi; to preserve the same areas as fj2
132  // we now generate a "phi_fj2", and then convert to a standard phi
133  double phi_fj2 = (iphi+0.5) * _dphi + _dphi*(_our_rand()-0.5)*_grid_scatter;
134  double phi;
135  if (_fj2_placement) phi = 0.5*pi - phi_fj2;
136  else phi = phi_fj2;
137  double rap = (irap+rap_offset) * _drap + _drap*(_our_rand()-0.5)*_grid_scatter
138  + _ghost_rap_offset ;
139  double pt = _mean_ghost_pt*(1+(_our_rand()-0.5)*_pt_scatter);
140 
141  double exprap = exp(+rap);
142  double pminus = pt/exprap;
143  double pplus = pt*exprap;
144  double px = pt*cos(phi);
145  double py = pt*sin(phi);
146  PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
147  // this call fills in the PseudoJet's cached rap,phi information,
148  // based on pre-existing knowledge. Watch out: if you get the hint
149  // wrong nobody will tell you, but you will certainly mess up
150  // your results.
151  mom.set_cached_rap_phi(rap,phi);
152 
153  // if we have an active selector and the particle does not pass the
154  // selection condition, move on to the next momentum
155  if (_selector.worker().get() && !_selector.pass(mom)) continue;
156  event.push_back(mom);
157  }
158  }
159 }
160 
162 
163  ostringstream ostr;
164  ostr << "ghosts of area " << actual_ghost_area()
165  << " (had requested " << ghost_area() << ")";
166  if (_selector.worker().get())
167  ostr << ", placed according to selector (" << _selector.description() << ")";
168  else
169  ostr << ", placed up to y = " << ghost_maxrap() ;
170  ostr << ", scattered wrt to perfect grid by (rel) " << grid_scatter()
171  << ", mean_ghost_pt = " << mean_ghost_pt()
172  << ", rel pt_scatter = " << pt_scatter()
173  << ", n repetitions of ghost distributions = " << repeat();
174  return ostr.str();
175 }
176 
177 FASTJET_END_NAMESPACE
178 
void set_cached_rap_phi(double rap, double phi)
in some cases when setting a 4-momentum, the user/program knows what rapidity and azimuth are associa...
Definition: PseudoJet.cc:326
void get_rapidity_extent(double &rapmin, double &rapmax) const
returns the rapidity range for which it may return &quot;true&quot;
Definition: Selector.hh:213
bool applies_jet_by_jet() const
returns true if this can be applied jet by jet
Definition: Selector.hh:196
void warn(const std::string &warning)
outputs a warning to standard error (or the user&#39;s default warning stream if set) ...
void _initialize()
does the initialization of actual ghost parameters
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:41
bool has_finite_area() const
returns true if it has a meaningful and finite area (i.e.
Definition: Selector.hh:231
std::string description() const
returns a textual description of the selector
Definition: Selector.hh:218
const SharedPtr< SelectorWorker > & worker() const
returns a (reference to) the underlying worker&#39;s shared pointer
Definition: Selector.hh:258
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
void set_fj2_placement(bool val)
if val is true, set ghost placement as it was in FastJet 2.X.
std::string description() const
for a summary
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
Definition: Selector.hh:171
void add_ghosts(std::vector< PseudoJet > &) const
push a set of ghost 4-momenta onto the back of the vector of PseudoJets
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:65