source: HiSusy/trunk/Delphes-3.0.0/external/fastjet/GhostedAreaSpec.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 7.1 KB
Line 
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
34using namespace std;
35
36FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
37
38BasicRandom<double> GhostedAreaSpec::_random_generator;
39LimitedWarning GhostedAreaSpec::_warn_fj2_placement_deprecated;
40
41/// explicit constructor
42GhostedAreaSpec::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
76void 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...
85void 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
112void 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
161string 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
177FASTJET_END_NAMESPACE
178
Note: See TracBrowser for help on using the repository browser.