FastJet 3.0.4
ClusterSequenceActiveAreaExplicitGhosts.cc
00001 //STARTHEADER
00002 // $Id: ClusterSequenceActiveAreaExplicitGhosts.cc 2687 2011-11-14 11:17:51Z soyez $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
00030 #include<limits>
00031 
00032 using namespace std;
00033 
00034 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00035 
00036 // save some typing
00037 typedef ClusterSequenceActiveAreaExplicitGhosts ClustSeqActAreaEG;
00038 
00039 
00040 LimitedWarning ClustSeqActAreaEG::_warnings;
00041 
00042 //----------------------------------------------------------------------
00043 ///
00044 void ClustSeqActAreaEG::_add_ghosts (
00045                          const GhostedAreaSpec & ghost_spec) {
00046 
00047   // add the ghosts to the jets
00048   ghost_spec.add_ghosts(_jets);
00049 
00050   // now add labelling...
00051   for (unsigned i = _initial_hard_n; i < _jets.size(); i++) {
00052     //_jets[i].set_user_index(1);
00053     _is_pure_ghost.push_back(true);
00054   }
00055 
00056   // and record some info from the ghost_spec
00057   _ghost_area = ghost_spec.actual_ghost_area();
00058   _n_ghosts   = ghost_spec.n_ghosts();
00059 }
00060 
00061 
00062 //----------------------------------------------------------------------
00063 // return the area of a jet
00064 double ClustSeqActAreaEG::area (const PseudoJet & jet) const {
00065   return _areas[jet.cluster_hist_index()];
00066 }
00067 
00068 
00069 //----------------------------------------------------------------------
00070 // return the total area
00071 double ClustSeqActAreaEG::total_area () const {
00072   return _n_ghosts * _ghost_area;
00073 }
00074 
00075 
00076 //----------------------------------------------------------------------
00077 // return the extended area of a jet
00078 PseudoJet ClustSeqActAreaEG::area_4vector (const PseudoJet & jet) const {
00079   return _area_4vectors[jet.cluster_hist_index()];
00080 }
00081 
00082 //----------------------------------------------------------------------
00083 bool ClustSeqActAreaEG::is_pure_ghost(const PseudoJet & jet) const 
00084 {
00085   return _is_pure_ghost[jet.cluster_hist_index()];
00086 }
00087 
00088 //----------------------------------------------------------------------
00089 bool ClustSeqActAreaEG::is_pure_ghost(int hist_ix) const 
00090 {
00091   return hist_ix >= 0 ? _is_pure_ghost[hist_ix] : false;
00092 }
00093 
00094 //----------------------------------------------------------------------
00095 double ClustSeqActAreaEG::empty_area(const Selector & selector) const {
00096   // make sure that the selector applies jet by jet
00097   if (! selector.applies_jet_by_jet()){
00098     throw Error("ClusterSequenceActiveAreaExplicitGhosts: empty area can only be computed from selectors applying jet by jet");
00099   }
00100 
00101   vector<PseudoJet> unclust = unclustered_particles();
00102   double area_local = 0.0;
00103   for (unsigned iu = 0; iu < unclust.size();  iu++) {
00104     if (is_pure_ghost(unclust[iu]) && selector.pass(unclust[iu])) {
00105       area_local += _ghost_area;
00106     }
00107   }
00108   return area_local;
00109 }
00110 
00111 //======================================================================
00112 // sort out the areas
00113 void ClustSeqActAreaEG::_post_process() {
00114 
00115   // first check for danger signals.
00116   // Establish largest ghost transverse momentum
00117   _max_ghost_perp2 = 0.0;
00118   for (int i = 0; i < _initial_n; i++) {
00119     if (_is_pure_ghost[i] && _jets[i].perp2() > _max_ghost_perp2) 
00120       _max_ghost_perp2 = _jets[i].perp2();
00121   }
00122 
00123   // now find out if any of the particles are close to danger
00124   double danger_ratio = numeric_limits<double>::epsilon();
00125   danger_ratio = danger_ratio * danger_ratio;
00126   _has_dangerous_particles = false;
00127   for (int i = 0; i < _initial_n; i++) {
00128     if (!_is_pure_ghost[i] && 
00129         danger_ratio * _jets[i].perp2() <=  _max_ghost_perp2) {
00130       _has_dangerous_particles = true;
00131       break;
00132     }
00133   }
00134 
00135   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");
00136 
00137   // sort out sizes
00138   _areas.resize(_history.size());
00139   _area_4vectors.resize(_history.size());
00140   _is_pure_ghost.resize(_history.size());
00141   
00142 //   copy(_jets.begin(), _jets.begin()+_initial_n, _area_4vectors.begin());
00143 //   for (int i = 0; i < _initial_n; i++) {
00144 //     if (_is_pure_ghost[i]) {
00145 //       _areas[i] = _ghost_area;
00146 //       // normalise pt to be _ghost_area (NB we make use of fact that
00147 //       // for initial particles, jet and clust_hist index are the same).
00148 //       _area_4vectors[i] *= (_ghost_area/_jets[i].perp());
00149 //     } else {
00150 //       _areas[i] = 0;
00151 //       _area_4vectors[i].reset(0,0,0,0);
00152 //     }
00153 //   }
00154   
00155   // First set up areas for the initial particles (ghost=_ghost_area,
00156   // real particles = 0); recall that _initial_n here is the number of
00157   // particles including ghosts
00158   for (int i = 0; i < _initial_n; i++) {
00159     if (_is_pure_ghost[i]) {
00160       _areas[i] = _ghost_area;
00161       // normalise pt to be _ghost_area (NB we make use of fact that
00162       // for initial particles, jet and clust_hist index are the same).
00163       //_area_4vectors[i] = (_ghost_area/_jets[i].perp()) * _jets[i];
00164       
00165       // NB: we use reset_momentum here, to ensure that the area 4
00166       // vectors do not acquire any structure (the structure would not
00167       // be meaningful for an area, and it messes up the use count (->
00168       // memory leaks if the user call delete_self_when_unused).
00169       _area_4vectors[i].reset_momentum(_jets[i]);
00170       _area_4vectors[i] *= (_ghost_area/_jets[i].perp());
00171     } else {
00172       _areas[i] = 0;
00173       _area_4vectors[i] = PseudoJet(0.0,0.0,0.0,0.0);
00174     }
00175   }
00176   
00177   // next follow the branching through and set up the areas 
00178   // and ghost-nature at each step of the clustering (rather than
00179   // each jet).
00180   for (unsigned i = _initial_n; i < _history.size(); i++) {
00181     if (_history[i].parent2 == BeamJet) {
00182       _is_pure_ghost[i]  = _is_pure_ghost[_history[i].parent1];
00183       _areas[i]          = _areas[_history[i].parent1];
00184       _area_4vectors[i] = _area_4vectors[_history[i].parent1];
00185     } else {
00186       _is_pure_ghost[i]  = _is_pure_ghost[_history[i].parent1] && 
00187                            _is_pure_ghost[_history[i].parent2]   ;
00188       _areas[i]          = _areas[_history[i].parent1] + 
00189                            _areas[_history[i].parent2]  ;                          
00190       _jet_def.recombiner()->recombine(_area_4vectors[_history[i].parent1], 
00191                            _area_4vectors[_history[i].parent2],
00192                            _area_4vectors[i]);  
00193 //      _area_4vectors[i] = _area_4vectors[_history[i].parent1] + 
00194 //                          _area_4vectors[_history[i].parent2]  ;
00195     }
00196 
00197   }
00198   
00199 }
00200 
00201 FASTJET_END_NAMESPACE
00202 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends