Loading...
Searching...
No Matches
GhostedAreaSpec.cc
//FJSTARTHEADER
// $Id$
//
// Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development. They are described in the original FastJet paper,
// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
// FastJet as part of work towards a scientific publication, please
// quote the version you use and include a citation to the manual and
// optionally also to hep-ph/0512210.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//FJENDHEADER
#include "fastjet/GhostedAreaSpec.hh"
#include "fastjet/Error.hh"
#include<iostream>
#include<sstream>
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
// in order to keep thread-safety, have an independent random
// generator for each thread
//#ifdef FASTJET_HAVE_LIMITED_THREAD_SAFETY
// thread_local
//#endif
BasicRandom<double> GhostedAreaSpec::_random_generator;
LimitedWarning GhostedAreaSpec::_warn_fj2_placement_deprecated;
LimitedWarning GhostedAreaSpec::_warn_fixed_last_seeds_nrepeat_gt_1;
/// explicit constructor
GhostedAreaSpec::GhostedAreaSpec(
const Selector & selector,
int repeat_in ,
double ghost_area_in ,
double grid_scatter_in ,
double pt_scatter_in ,
double mean_ghost_pt_in ,
BasicRandom<double> *user_random_generator):
_repeat(repeat_in),
_ghost_area(ghost_area_in),
_grid_scatter(grid_scatter_in),
_pt_scatter(pt_scatter_in),
_mean_ghost_pt(mean_ghost_pt_in),
_fj2_placement(false),
_selector(selector),
_actual_ghost_area(-1.0),
_user_random_generator(user_random_generator)
{
// check the selector has the properties needed -- an area and
// applicability jet-by-jet (the latter follows automatically from
// the former?)
if (!_selector.has_finite_area()) throw Error("To construct a GhostedAreaSpec with a Selector, the selector must have a finite area");
if (!_selector.applies_jet_by_jet()) throw Error("To construct a GhostedAreaSpec with a Selector, the selector must apply jet-by-jet");
// get the internal rapidity extent from the selector
double ghost_maxrap_local, ghost_minrap_local;
_selector.get_rapidity_extent(ghost_minrap_local, ghost_maxrap_local);
_ghost_maxrap = 0.5*(ghost_maxrap_local - ghost_minrap_local);
_ghost_rap_offset = 0.5*(ghost_maxrap_local + ghost_minrap_local);
_initialize();
}
//======================================================================
// sets fj2 ghost placement
void GhostedAreaSpec::set_fj2_placement(bool val) {
_fj2_placement = val; _initialize();
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.");
}
//======================================================================
/// sets the detailed parameters for the ghosts (which may not be quite
/// the same as those requested -- this is in order for things to fit
/// in nicely into 2pi etc...
void GhostedAreaSpec::_initialize() {
// add on area-measuring dummy particles
_drap = sqrt(_ghost_area);
_dphi = _drap;
if (_fj2_placement) {
_nphi = int(ceil(twopi/_dphi)); _dphi = twopi/_nphi;
_nrap = int(ceil(_ghost_maxrap/_drap)); _drap = _ghost_maxrap / _nrap;
_actual_ghost_area = _dphi * _drap;
_n_ghosts = (2*_nrap+1)*_nphi;
} else {
// for FJ3, update the ghost placement as follows
// - use nearest int rather than ceiling in determining number of
// phi and rapidity locations, because this is more stable when
// the user is trying to get an exact number based on the area
// - rather than placing ghosts up to maximum rapidity
_nphi = int(twopi/_dphi + 0.5); _dphi = twopi/_nphi;
_nrap = int(_ghost_maxrap/_drap + 0.5); _drap = _ghost_maxrap / _nrap;
_actual_ghost_area = _dphi * _drap;
_n_ghosts = (2*_nrap)*_nphi;
}
// checkpoint the status of the random number generator.
checkpoint_random();
//_random_generator.info(cerr);
_fixed_seed=vector<int>();
}
//----------------------------------------------------------------------
/// adds the ghost 4-momenta to the vector of PseudoJet's
void GhostedAreaSpec::add_ghosts(vector<PseudoJet> & event) const {
// determine the number of ghosts in rapidity
double rap_offset;
int nrap_upper;
if (_fj2_placement) {
rap_offset = 0.0;
nrap_upper = _nrap;
} else {
rap_offset = 0.5;
nrap_upper = _nrap-1;
}
// generate all the random numbers
// Note that we have 3 numbers per ghost
unsigned int n_random = (nrap_upper+_nrap+1)*_nphi*3;
double * all_random = new double[n_random];
if (_fixed_seed.size()){
if (_repeat > 1) _warn_fixed_last_seeds_nrepeat_gt_1
.warn("Using fixed seeds (or accessing last used seeds) not sensible with repeat>1");
// take a copy of the random generator and use that to generate
// things with fixed seeds
BasicRandom<double> local_rand = generator_at_own_risk();
local_rand.set_status(_fixed_seed);
_last_used_seed = _fixed_seed;
local_rand(n_random, all_random);
// // update the seeds so next time we get the next set
// local_rand.get_status(_fixed_seeds);
} else {
// make sure we use a lock and get the seeds used for this
// generation
_our_rand(n_random, all_random, _last_used_seed);
}
// a counter for where we are in the random array
unsigned int random_counter=0;
// add momenta for ghosts
for (int irap = -_nrap; irap <= nrap_upper; irap++) {
for (int iphi = 0; iphi < _nphi; iphi++) {
// include random offsets for all quantities
//----------------------------------------------
// NB: in FJ2 we'd exchanged the px and py components relative to a
// standard definition of phi; to preserve the same areas as fj2
// we now generate a "phi_fj2", and then convert to a standard phi
//double phi_fj2 = (iphi+0.5) * _dphi + _dphi*(_our_rand()-0.5)*_grid_scatter;
double phi_fj2 = (iphi+0.5) * _dphi + _dphi*(all_random[random_counter++]-0.5)*_grid_scatter;
double phi;
if (_fj2_placement) phi = 0.5*pi - phi_fj2;
else phi = phi_fj2;
//double rap = (irap+rap_offset) * _drap + _drap*(_our_rand()-0.5)*_grid_scatter
// + _ghost_rap_offset ;
double rap = (irap+rap_offset) * _drap + _drap*(all_random[random_counter++]-0.5)*_grid_scatter
+ _ghost_rap_offset ;
//double pt = _mean_ghost_pt*(1+(_our_rand()-0.5)*_pt_scatter);
double pt = _mean_ghost_pt*(1+(all_random[random_counter++]-0.5)*_pt_scatter);
double exprap = exp(+rap);
double pminus = pt/exprap;
double pplus = pt*exprap;
double px = pt*cos(phi);
double py = pt*sin(phi);
PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
// this call fills in the PseudoJet's cached rap,phi information,
// based on pre-existing knowledge. Watch out: if you get the hint
// wrong nobody will tell you, but you will certainly mess up
// your results.
mom.set_cached_rap_phi(rap,phi);
// if we have an active selector and the particle does not pass the
// selection condition, move on to the next momentum
if (_selector.worker().get() && !_selector.pass(mom)) continue;
event.push_back(mom);
}
}
// release memory
delete[] all_random;
// safety check
assert(random_counter==n_random);
}
string GhostedAreaSpec::description() const {
ostringstream ostr;
ostr << "ghosts of area " << actual_ghost_area()
<< " (had requested " << ghost_area() << ")";
if (_selector.worker().get())
ostr << ", placed according to selector (" << _selector.description() << ")";
else
ostr << ", placed up to y = " << ghost_maxrap() ;
ostr << ", scattered wrt to perfect grid by (rel) " << grid_scatter()
<< ", mean_ghost_pt = " << mean_ghost_pt()
<< ", rel pt_scatter = " << pt_scatter()
<< ", n repetitions of ghost distributions = " << repeat();
return ostr.str();
}
FASTJET_END_NAMESPACE
fastjet::BasicRandom
Definition BasicRandom.hh:49
fastjet::Error
base class corresponding to errors that can be thrown by FastJet
Definition Error.hh:52
fastjet::GhostedAreaSpec::description
std::string description() const
for a summary
Definition GhostedAreaSpec.cc:207
fastjet::GhostedAreaSpec::add_ghosts
void add_ghosts(std::vector< PseudoJet > &) const
push a set of ghost 4-momenta onto the back of the vector of PseudoJets
Definition GhostedAreaSpec.cc:123
fastjet::GhostedAreaSpec::set_fj2_placement
void set_fj2_placement(bool val)
if val is true, set ghost placement as it was in FastJet 2.X.
Definition GhostedAreaSpec.cc:86
fastjet::GhostedAreaSpec::_initialize
void _initialize()
does the initialization of actual ghost parameters
Definition GhostedAreaSpec.cc:95
fastjet::GhostedAreaSpec::GhostedAreaSpec
GhostedAreaSpec()
default constructor
Definition GhostedAreaSpec.hh:106
fastjet::GhostedAreaSpec::generator_at_own_risk
BasicRandom< double > & generator_at_own_risk() const
very deprecated public access to the generator itself
Definition GhostedAreaSpec.hh:310
fastjet::PseudoJet
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition PseudoJet.hh:68
fastjet::PseudoJet::set_cached_rap_phi
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:445
fastjet::Selector
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition Selector.hh:149