1 | //STARTHEADER |
---|
2 | // $Id: GhostedAreaSpec.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
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 |
---|
76 | void GhostedAreaSpec::set_fj2_placement(bool val) { |
---|
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... |
---|
85 | void GhostedAreaSpec::_initialize() { |
---|
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 | |
---|
161 | string GhostedAreaSpec::description() const { |
---|
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 | |
---|