| 1 | //STARTHEADER
|
|---|
| 2 | // $Id$
|
|---|
| 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/ClusterSequenceActiveAreaExplicitGhosts.hh"
|
|---|
| 30 | #include<limits>
|
|---|
| 31 |
|
|---|
| 32 | using namespace std;
|
|---|
| 33 |
|
|---|
| 34 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 35 |
|
|---|
| 36 | // save some typing
|
|---|
| 37 | typedef ClusterSequenceActiveAreaExplicitGhosts ClustSeqActAreaEG;
|
|---|
| 38 |
|
|---|
| 39 |
|
|---|
| 40 | LimitedWarning ClustSeqActAreaEG::_warnings;
|
|---|
| 41 |
|
|---|
| 42 | //----------------------------------------------------------------------
|
|---|
| 43 | ///
|
|---|
| 44 | void ClustSeqActAreaEG::_add_ghosts (
|
|---|
| 45 | const GhostedAreaSpec & ghost_spec) {
|
|---|
| 46 |
|
|---|
| 47 | // add the ghosts to the jets
|
|---|
| 48 | ghost_spec.add_ghosts(_jets);
|
|---|
| 49 |
|
|---|
| 50 | // now add labelling...
|
|---|
| 51 | for (unsigned i = _initial_hard_n; i < _jets.size(); i++) {
|
|---|
| 52 | //_jets[i].set_user_index(1);
|
|---|
| 53 | _is_pure_ghost.push_back(true);
|
|---|
| 54 | }
|
|---|
| 55 |
|
|---|
| 56 | // and record some info from the ghost_spec
|
|---|
| 57 | _ghost_area = ghost_spec.actual_ghost_area();
|
|---|
| 58 | _n_ghosts = ghost_spec.n_ghosts();
|
|---|
| 59 | }
|
|---|
| 60 |
|
|---|
| 61 |
|
|---|
| 62 | //----------------------------------------------------------------------
|
|---|
| 63 | // return the area of a jet
|
|---|
| 64 | double ClustSeqActAreaEG::area (const PseudoJet & jet) const {
|
|---|
| 65 | return _areas[jet.cluster_hist_index()];
|
|---|
| 66 | }
|
|---|
| 67 |
|
|---|
| 68 |
|
|---|
| 69 | //----------------------------------------------------------------------
|
|---|
| 70 | // return the total area
|
|---|
| 71 | double ClustSeqActAreaEG::total_area () const {
|
|---|
| 72 | return _n_ghosts * _ghost_area;
|
|---|
| 73 | }
|
|---|
| 74 |
|
|---|
| 75 |
|
|---|
| 76 | //----------------------------------------------------------------------
|
|---|
| 77 | // return the extended area of a jet
|
|---|
| 78 | PseudoJet ClustSeqActAreaEG::area_4vector (const PseudoJet & jet) const {
|
|---|
| 79 | return _area_4vectors[jet.cluster_hist_index()];
|
|---|
| 80 | }
|
|---|
| 81 |
|
|---|
| 82 | //----------------------------------------------------------------------
|
|---|
| 83 | bool ClustSeqActAreaEG::is_pure_ghost(const PseudoJet & jet) const
|
|---|
| 84 | {
|
|---|
| 85 | return _is_pure_ghost[jet.cluster_hist_index()];
|
|---|
| 86 | }
|
|---|
| 87 |
|
|---|
| 88 | //----------------------------------------------------------------------
|
|---|
| 89 | bool ClustSeqActAreaEG::is_pure_ghost(int hist_ix) const
|
|---|
| 90 | {
|
|---|
| 91 | return hist_ix >= 0 ? _is_pure_ghost[hist_ix] : false;
|
|---|
| 92 | }
|
|---|
| 93 |
|
|---|
| 94 | //----------------------------------------------------------------------
|
|---|
| 95 | double ClustSeqActAreaEG::empty_area(const Selector & selector) const {
|
|---|
| 96 | // make sure that the selector applies jet by jet
|
|---|
| 97 | if (! selector.applies_jet_by_jet()){
|
|---|
| 98 | throw Error("ClusterSequenceActiveAreaExplicitGhosts: empty area can only be computed from selectors applying jet by jet");
|
|---|
| 99 | }
|
|---|
| 100 |
|
|---|
| 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;
|
|---|
| 106 | }
|
|---|
| 107 | }
|
|---|
| 108 | return area_local;
|
|---|
| 109 | }
|
|---|
| 110 |
|
|---|
| 111 | //======================================================================
|
|---|
| 112 | // sort out the areas
|
|---|
| 113 | void ClustSeqActAreaEG::_post_process() {
|
|---|
| 114 |
|
|---|
| 115 | // first check for danger signals.
|
|---|
| 116 | // Establish largest ghost transverse momentum
|
|---|
| 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();
|
|---|
| 121 | }
|
|---|
| 122 |
|
|---|
| 123 | // now find out if any of the particles are close to danger
|
|---|
| 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;
|
|---|
| 131 | break;
|
|---|
| 132 | }
|
|---|
| 133 | }
|
|---|
| 134 |
|
|---|
| 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");
|
|---|
| 136 |
|
|---|
| 137 | // sort out sizes
|
|---|
| 138 | _areas.resize(_history.size());
|
|---|
| 139 | _area_4vectors.resize(_history.size());
|
|---|
| 140 | _is_pure_ghost.resize(_history.size());
|
|---|
| 141 |
|
|---|
| 142 | // copy(_jets.begin(), _jets.begin()+_initial_n, _area_4vectors.begin());
|
|---|
| 143 | // for (int i = 0; i < _initial_n; i++) {
|
|---|
| 144 | // if (_is_pure_ghost[i]) {
|
|---|
| 145 | // _areas[i] = _ghost_area;
|
|---|
| 146 | // // normalise pt to be _ghost_area (NB we make use of fact that
|
|---|
| 147 | // // for initial particles, jet and clust_hist index are the same).
|
|---|
| 148 | // _area_4vectors[i] *= (_ghost_area/_jets[i].perp());
|
|---|
| 149 | // } else {
|
|---|
| 150 | // _areas[i] = 0;
|
|---|
| 151 | // _area_4vectors[i].reset(0,0,0,0);
|
|---|
| 152 | // }
|
|---|
| 153 | // }
|
|---|
| 154 |
|
|---|
| 155 | // First set up areas for the initial particles (ghost=_ghost_area,
|
|---|
| 156 | // real particles = 0); recall that _initial_n here is the number of
|
|---|
| 157 | // particles including ghosts
|
|---|
| 158 | for (int i = 0; i < _initial_n; i++) {
|
|---|
| 159 | if (_is_pure_ghost[i]) {
|
|---|
| 160 | _areas[i] = _ghost_area;
|
|---|
| 161 | // normalise pt to be _ghost_area (NB we make use of fact that
|
|---|
| 162 | // for initial particles, jet and clust_hist index are the same).
|
|---|
| 163 | //_area_4vectors[i] = (_ghost_area/_jets[i].perp()) * _jets[i];
|
|---|
| 164 |
|
|---|
| 165 | // NB: we use reset_momentum here, to ensure that the area 4
|
|---|
| 166 | // vectors do not acquire any structure (the structure would not
|
|---|
| 167 | // be meaningful for an area, and it messes up the use count (->
|
|---|
| 168 | // memory leaks if the user call delete_self_when_unused).
|
|---|
| 169 | _area_4vectors[i].reset_momentum(_jets[i]);
|
|---|
| 170 | _area_4vectors[i] *= (_ghost_area/_jets[i].perp());
|
|---|
| 171 | } else {
|
|---|
| 172 | _areas[i] = 0;
|
|---|
| 173 | _area_4vectors[i] = PseudoJet(0.0,0.0,0.0,0.0);
|
|---|
| 174 | }
|
|---|
| 175 | }
|
|---|
| 176 |
|
|---|
| 177 | // next follow the branching through and set up the areas
|
|---|
| 178 | // and ghost-nature at each step of the clustering (rather than
|
|---|
| 179 | // each jet).
|
|---|
| 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];
|
|---|
| 185 | } else {
|
|---|
| 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],
|
|---|
| 192 | _area_4vectors[i]);
|
|---|
| 193 | // _area_4vectors[i] = _area_4vectors[_history[i].parent1] +
|
|---|
| 194 | // _area_4vectors[_history[i].parent2] ;
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 | }
|
|---|
| 198 |
|
|---|
| 199 | }
|
|---|
| 200 |
|
|---|
| 201 | FASTJET_END_NAMESPACE
|
|---|
| 202 |
|
|---|
