Fork me on GitHub

source: git/external/fastjet/GhostedAreaSpec.cc@ 0e6b14e

ImprovedOutputFile Timing
Last change on this file since 0e6b14e was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 7.3 KB
Line 
1//FJSTARTHEADER
2// $Id: GhostedAreaSpec.cc 4354 2018-04-22 07:12:37Z salam $
3//
4// Copyright (c) 2005-2018, 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
40BasicRandom<double> GhostedAreaSpec::_random_generator;
41LimitedWarning GhostedAreaSpec::_warn_fj2_placement_deprecated;
42
43/// explicit constructor
44GhostedAreaSpec::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
78void GhostedAreaSpec::set_fj2_placement(bool val) {
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...
87void GhostedAreaSpec::_initialize() {
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
114void 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
163string GhostedAreaSpec::description() const {
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
179FASTJET_END_NAMESPACE
180
Note: See TracBrowser for help on using the repository browser.