FastJet 3.4.3
Loading...
Searching...
No Matches
GhostedAreaSpec.cc
1//FJSTARTHEADER
2// $Id$
3//
4// Copyright (c) 2005-2024, 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
36using namespace std;
37
38FASTJET_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
45BasicRandom<double> GhostedAreaSpec::_random_generator;
46
47LimitedWarning GhostedAreaSpec::_warn_fj2_placement_deprecated;
48LimitedWarning GhostedAreaSpec::_warn_fixed_last_seeds_nrepeat_gt_1;
49
50/// explicit constructor
51GhostedAreaSpec::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
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
123void 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
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
223FASTJET_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:434
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition Selector.hh:149
const SharedPtr< SelectorWorker > & worker() const
returns a (reference to) the underlying worker's shared pointer
Definition Selector.hh:277
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
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