29 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
34 FASTJET_BEGIN_NAMESPACE
37 typedef ClusterSequenceActiveAreaExplicitGhosts ClustSeqActAreaEG;
40 LimitedWarning ClustSeqActAreaEG::_warnings;
44 void ClustSeqActAreaEG::_add_ghosts (
45 const GhostedAreaSpec & ghost_spec) {
48 ghost_spec.add_ghosts(_jets);
51 for (
unsigned i = _initial_hard_n; i < _jets.size(); i++) {
53 _is_pure_ghost.push_back(
true);
57 _ghost_area = ghost_spec.actual_ghost_area();
58 _n_ghosts = ghost_spec.n_ghosts();
64 double ClustSeqActAreaEG::area (
const PseudoJet & jet)
const {
71 double ClustSeqActAreaEG::total_area ()
const {
72 return _n_ghosts * _ghost_area;
83 bool ClustSeqActAreaEG::is_pure_ghost(
const PseudoJet & jet)
const
89 bool ClustSeqActAreaEG::is_pure_ghost(
int hist_ix)
const
91 return hist_ix >= 0 ? _is_pure_ghost[hist_ix] :
false;
95 double ClustSeqActAreaEG::empty_area(
const Selector & selector)
const {
98 throw Error(
"ClusterSequenceActiveAreaExplicitGhosts: empty area can only be computed from selectors applying jet by jet");
101 vector<PseudoJet> unclust = unclustered_particles();
102 double area_local = 0.0;
103 for (
unsigned iu = 0; iu < unclust.size(); iu++) {
104 if (is_pure_ghost(unclust[iu]) && selector.
pass(unclust[iu])) {
105 area_local += _ghost_area;
113 void ClustSeqActAreaEG::_post_process() {
117 _max_ghost_perp2 = 0.0;
118 for (
int i = 0; i < _initial_n; i++) {
119 if (_is_pure_ghost[i] && _jets[i].perp2() > _max_ghost_perp2)
120 _max_ghost_perp2 = _jets[i].perp2();
124 double danger_ratio = numeric_limits<double>::epsilon();
125 danger_ratio = danger_ratio * danger_ratio;
126 _has_dangerous_particles =
false;
127 for (
int i = 0; i < _initial_n; i++) {
128 if (!_is_pure_ghost[i] &&
129 danger_ratio * _jets[i].perp2() <= _max_ghost_perp2) {
130 _has_dangerous_particles =
true;
135 if (_has_dangerous_particles) _warnings.warn(
"ClusterSequenceActiveAreaExplicitGhosts: \n ghosts not sufficiently soft wrt some of the input particles\n a common cause is (unphysical?) input particles with pt=0 but finite rapidity");
138 _areas.resize(_history.size());
139 _area_4vectors.resize(_history.size());
140 _is_pure_ghost.resize(_history.size());
158 for (
int i = 0; i < _initial_n; i++) {
159 if (_is_pure_ghost[i]) {
160 _areas[i] = _ghost_area;
169 _area_4vectors[i].reset_momentum(_jets[i]);
170 _area_4vectors[i] *= (_ghost_area/_jets[i].perp());
173 _area_4vectors[i] = PseudoJet(0.0,0.0,0.0,0.0);
180 for (
unsigned i = _initial_n; i < _history.size(); i++) {
181 if (_history[i].parent2 == BeamJet) {
182 _is_pure_ghost[i] = _is_pure_ghost[_history[i].parent1];
183 _areas[i] = _areas[_history[i].parent1];
184 _area_4vectors[i] = _area_4vectors[_history[i].parent1];
186 _is_pure_ghost[i] = _is_pure_ghost[_history[i].parent1] &&
187 _is_pure_ghost[_history[i].parent2] ;
188 _areas[i] = _areas[_history[i].parent1] +
189 _areas[_history[i].parent2] ;
190 _jet_def.recombiner()->recombine(_area_4vectors[_history[i].parent1],
191 _area_4vectors[_history[i].parent2],
201 FASTJET_END_NAMESPACE