31 #include "fastjet/GhostedAreaSpec.hh"
32 #include "fastjet/Error.hh"
38 FASTJET_BEGIN_NAMESPACE
45 BasicRandom<double> GhostedAreaSpec::_random_generator;
47 LimitedWarning GhostedAreaSpec::_warn_fj2_placement_deprecated;
48 LimitedWarning GhostedAreaSpec::_warn_fixed_last_seeds_nrepeat_gt_1;
51 GhostedAreaSpec::GhostedAreaSpec(
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):
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),
66 _actual_ghost_area(-1.0),
67 _user_random_generator(user_random_generator)
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");
75 double ghost_maxrap_local, ghost_minrap_local;
77 _ghost_maxrap = 0.5*(ghost_maxrap_local - ghost_minrap_local);
78 _ghost_rap_offset = 0.5*(ghost_maxrap_local + ghost_minrap_local);
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.");
97 _drap = sqrt(_ghost_area);
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;
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;
118 _fixed_seed=vector<int>();
128 if (_fj2_placement) {
133 nrap_upper = _nrap-1;
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");
146 local_rand.set_status(_fixed_seed);
147 _last_used_seed = _fixed_seed;
148 local_rand(n_random, all_random);
154 _our_rand(n_random, all_random, _last_used_seed);
158 unsigned int random_counter=0;
161 for (
int irap = -_nrap; irap <= nrap_upper; irap++) {
162 for (
int iphi = 0; iphi < _nphi; iphi++) {
170 double phi_fj2 = (iphi+0.5) * _dphi + _dphi*(all_random[random_counter++]-0.5)*_grid_scatter;
172 if (_fj2_placement) phi = 0.5*pi - phi_fj2;
176 double rap = (irap+rap_offset) * _drap + _drap*(all_random[random_counter++]-0.5)*_grid_scatter
177 + _ghost_rap_offset ;
179 double pt = _mean_ghost_pt*(1+(all_random[random_counter++]-0.5)*_pt_scatter);
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));
195 if (_selector.
worker().get() && !_selector.
pass(mom))
continue;
196 event.push_back(mom);
204 assert(random_counter==n_random);
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() <<
")";
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();
223 FASTJET_END_NAMESPACE