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