FastJet  3.4.0
GhostedAreaSpec.cc
1 //FJSTARTHEADER
2 // $Id$
3 //
4 // Copyright (c) 2005-2021, 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 // in order to keep thread-safety, have an independent random
41 // generator for each thread
42 //#ifdef FASTJET_HAVE_LIMITED_THREAD_SAFETY
43 // thread_local
44 //#endif
45 BasicRandom<double> GhostedAreaSpec::_random_generator;
46 
47 LimitedWarning GhostedAreaSpec::_warn_fj2_placement_deprecated;
48 LimitedWarning GhostedAreaSpec::_warn_fixed_last_seeds_nrepeat_gt_1;
49 
50 /// explicit constructor
51 GhostedAreaSpec::GhostedAreaSpec(
52  const Selector & selector,
53  int repeat_in ,
54  double ghost_area_in ,
55  double grid_scatter_in ,
56  double pt_scatter_in ,
57  double mean_ghost_pt_in ,
58  BasicRandom<double> *user_random_generator):
59  _repeat(repeat_in),
60  _ghost_area(ghost_area_in),
61  _grid_scatter(grid_scatter_in),
62  _pt_scatter(pt_scatter_in),
63  _mean_ghost_pt(mean_ghost_pt_in),
64  _fj2_placement(false),
65  _selector(selector),
66  _actual_ghost_area(-1.0),
67  _user_random_generator(user_random_generator)
68  {
69  // check the selector has the properties needed -- an area and
70  // applicability jet-by-jet (the latter follows automatically from
71  // the former?)
72  if (!_selector.has_finite_area()) throw Error("To construct a GhostedAreaSpec with a Selector, the selector must have a finite area");
73  if (!_selector.applies_jet_by_jet()) throw Error("To construct a GhostedAreaSpec with a Selector, the selector must apply jet-by-jet");
74  // get the internal rapidity extent from the selector
75  double ghost_maxrap_local, ghost_minrap_local;
76  _selector.get_rapidity_extent(ghost_minrap_local, ghost_maxrap_local);
77  _ghost_maxrap = 0.5*(ghost_maxrap_local - ghost_minrap_local);
78  _ghost_rap_offset = 0.5*(ghost_maxrap_local + ghost_minrap_local);
79 
80  _initialize();
81 
82 }
83 
84 //======================================================================
85 // sets fj2 ghost placement
87  _fj2_placement = val; _initialize();
88  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.");
89 }
90 
91 //======================================================================
92 /// sets the detailed parameters for the ghosts (which may not be quite
93 /// the same as those requested -- this is in order for things to fit
94 /// in nicely into 2pi etc...
96  // add on area-measuring dummy particles
97  _drap = sqrt(_ghost_area);
98  _dphi = _drap;
99  if (_fj2_placement) {
100  _nphi = int(ceil(twopi/_dphi)); _dphi = twopi/_nphi;
101  _nrap = int(ceil(_ghost_maxrap/_drap)); _drap = _ghost_maxrap / _nrap;
102  _actual_ghost_area = _dphi * _drap;
103  _n_ghosts = (2*_nrap+1)*_nphi;
104  } else {
105  // for FJ3, update the ghost placement as follows
106  // - use nearest int rather than ceiling in determining number of
107  // phi and rapidity locations, because this is more stable when
108  // the user is trying to get an exact number based on the area
109  // - rather than placing ghosts up to maximum rapidity
110  _nphi = int(twopi/_dphi + 0.5); _dphi = twopi/_nphi;
111  _nrap = int(_ghost_maxrap/_drap + 0.5); _drap = _ghost_maxrap / _nrap;
112  _actual_ghost_area = _dphi * _drap;
113  _n_ghosts = (2*_nrap)*_nphi;
114  }
115  // checkpoint the status of the random number generator.
116  checkpoint_random();
117  //_random_generator.info(cerr);
118  _fixed_seed=vector<int>();
119 }
120 
121 //----------------------------------------------------------------------
122 /// adds the ghost 4-momenta to the vector of PseudoJet's
123 void GhostedAreaSpec::add_ghosts(vector<PseudoJet> & event) const {
124 
125  // determine the number of ghosts in rapidity
126  double rap_offset;
127  int nrap_upper;
128  if (_fj2_placement) {
129  rap_offset = 0.0;
130  nrap_upper = _nrap;
131  } else {
132  rap_offset = 0.5;
133  nrap_upper = _nrap-1;
134  }
135 
136  // generate all the random numbers
137  // Note that we have 3 numbers per ghost
138  unsigned int n_random = (nrap_upper+_nrap+1)*_nphi*3;
139  double * all_random = new double[n_random];
140  if (_fixed_seed.size()){
141  if (_repeat > 1) _warn_fixed_last_seeds_nrepeat_gt_1
142  .warn("Using fixed seeds (or accessing last used seeds) not sensible with repeat>1");
143  // take a copy of the random generator and use that to generate
144  // things with fixed seeds
145  BasicRandom<double> local_rand = generator_at_own_risk();
146  local_rand.set_status(_fixed_seed);
147  _last_used_seed = _fixed_seed;
148  local_rand(n_random, all_random);
149  // // update the seeds so next time we get the next set
150  // local_rand.get_status(_fixed_seeds);
151  } else {
152  // make sure we use a lock and get the seeds used for this
153  // generation
154  _our_rand(n_random, all_random, _last_used_seed);
155  }
156 
157  // a counter for where we are in the random array
158  unsigned int random_counter=0;
159 
160  // add momenta for ghosts
161  for (int irap = -_nrap; irap <= nrap_upper; irap++) {
162  for (int iphi = 0; iphi < _nphi; iphi++) {
163 
164  // include random offsets for all quantities
165  //----------------------------------------------
166  // NB: in FJ2 we'd exchanged the px and py components relative to a
167  // standard definition of phi; to preserve the same areas as fj2
168  // we now generate a "phi_fj2", and then convert to a standard phi
169  //double phi_fj2 = (iphi+0.5) * _dphi + _dphi*(_our_rand()-0.5)*_grid_scatter;
170  double phi_fj2 = (iphi+0.5) * _dphi + _dphi*(all_random[random_counter++]-0.5)*_grid_scatter;
171  double phi;
172  if (_fj2_placement) phi = 0.5*pi - phi_fj2;
173  else phi = phi_fj2;
174  //double rap = (irap+rap_offset) * _drap + _drap*(_our_rand()-0.5)*_grid_scatter
175  // + _ghost_rap_offset ;
176  double rap = (irap+rap_offset) * _drap + _drap*(all_random[random_counter++]-0.5)*_grid_scatter
177  + _ghost_rap_offset ;
178  //double pt = _mean_ghost_pt*(1+(_our_rand()-0.5)*_pt_scatter);
179  double pt = _mean_ghost_pt*(1+(all_random[random_counter++]-0.5)*_pt_scatter);
180 
181  double exprap = exp(+rap);
182  double pminus = pt/exprap;
183  double pplus = pt*exprap;
184  double px = pt*cos(phi);
185  double py = pt*sin(phi);
186  PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
187  // this call fills in the PseudoJet's cached rap,phi information,
188  // based on pre-existing knowledge. Watch out: if you get the hint
189  // wrong nobody will tell you, but you will certainly mess up
190  // your results.
191  mom.set_cached_rap_phi(rap,phi);
192 
193  // if we have an active selector and the particle does not pass the
194  // selection condition, move on to the next momentum
195  if (_selector.worker().get() && !_selector.pass(mom)) continue;
196  event.push_back(mom);
197  }
198  }
199 
200  // release memory
201  delete[] all_random;
202 
203  // safety check
204  assert(random_counter==n_random);
205 }
206 
208 
209  ostringstream ostr;
210  ostr << "ghosts of area " << actual_ghost_area()
211  << " (had requested " << ghost_area() << ")";
212  if (_selector.worker().get())
213  ostr << ", placed according to selector (" << _selector.description() << ")";
214  else
215  ostr << ", placed up to y = " << ghost_maxrap() ;
216  ostr << ", scattered wrt to perfect grid by (rel) " << grid_scatter()
217  << ", mean_ghost_pt = " << mean_ghost_pt()
218  << ", rel pt_scatter = " << pt_scatter()
219  << ", n repetitions of ghost distributions = " << repeat();
220  return ostr.str();
221 }
222 
223 FASTJET_END_NAMESPACE
224 
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:52
std::string description() const
for a summary
void add_ghosts(std::vector< PseudoJet > &) const
push a set of ghost 4-momenta onto the back of the vector of PseudoJets
void set_fj2_placement(bool val)
if val is true, set ghost placement as it was in FastJet 2.X.
void _initialize()
does the initialization of actual ghost parameters
BasicRandom< double > & generator_at_own_risk() const
very deprecated public access to the generator itself
void warn(const char *warning)
outputs a warning to standard error (or the user's default warning stream if set)
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition: PseudoJet.hh:68
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:445
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:149
bool has_finite_area() const
returns true if it has a meaningful and finite area (i.e.
Definition: Selector.hh:250
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
Definition: Selector.hh:178
bool applies_jet_by_jet() const
returns true if this can be applied jet by jet
Definition: Selector.hh:215
const SharedPtr< SelectorWorker > & worker() const
returns a (reference to) the underlying worker's shared pointer
Definition: Selector.hh:277
std::string description() const
returns a textual description of the selector
Definition: Selector.hh:237
void get_rapidity_extent(double &rapmin, double &rapmax) const
returns the rapidity range for which it may return "true"
Definition: Selector.hh:232