fastjet 2.4.5
ClusterSequenceActiveAreaExplicitGhosts.cc
Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: ClusterSequenceActiveAreaExplicitGhosts.cc 1571 2009-05-25 15:32:40Z salam $
00003 //
00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
00030 
00031 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
00032 #include<limits>
00033 
00034 using namespace std;
00035 
00036 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00037 
00038 // save some typing
00039 typedef ClusterSequenceActiveAreaExplicitGhosts ClustSeqActAreaEG;
00040 
00041 
00042 LimitedWarning ClustSeqActAreaEG::_warnings;
00043 
00044 //----------------------------------------------------------------------
00046 void ClustSeqActAreaEG::_add_ghosts (
00047                          const GhostedAreaSpec & ghost_spec) {
00048 
00049   // add the ghosts to the jets
00050   ghost_spec.add_ghosts(_jets);
00051 
00052   // now add labelling...
00053   for (unsigned i = _initial_hard_n; i < _jets.size(); i++) {
00054     //_jets[i].set_user_index(1);
00055     _is_pure_ghost.push_back(true);
00056   }
00057 
00058   // and record some info from the ghost_spec
00059   _ghost_area = ghost_spec.actual_ghost_area();
00060   _n_ghosts   = ghost_spec.n_ghosts();
00061 }
00062 
00063 
00064 //----------------------------------------------------------------------
00065 // return the area of a jet
00066 double ClustSeqActAreaEG::area (const PseudoJet & jet) const {
00067   return _areas[jet.cluster_hist_index()];
00068 }
00069 
00070 
00071 //----------------------------------------------------------------------
00072 // return the total area
00073 double ClustSeqActAreaEG::total_area () const {
00074   return _n_ghosts * _ghost_area;
00075 }
00076 
00077 
00078 //----------------------------------------------------------------------
00079 // return the extended area of a jet
00080 PseudoJet ClustSeqActAreaEG::area_4vector (const PseudoJet & jet) const {
00081   return _area_4vectors[jet.cluster_hist_index()];
00082 }
00083 
00084 //----------------------------------------------------------------------
00085 bool ClustSeqActAreaEG::is_pure_ghost(const PseudoJet & jet) const 
00086 {
00087   return _is_pure_ghost[jet.cluster_hist_index()];
00088 }
00089 
00090 //----------------------------------------------------------------------
00091 bool ClustSeqActAreaEG::is_pure_ghost(int hist_ix) const 
00092 {
00093   return hist_ix >= 0 ? _is_pure_ghost[hist_ix] : false;
00094 }
00095 
00096 //----------------------------------------------------------------------
00097 double ClustSeqActAreaEG::empty_area(const RangeDefinition & range) const {
00098   vector<PseudoJet> unclust = unclustered_particles();
00099   double area = 0.0;
00100   for (unsigned iu = 0; iu < unclust.size();  iu++) {
00101     if (is_pure_ghost(unclust[iu]) && range.is_in_range(unclust[iu])) {
00102       area += _ghost_area;
00103     }
00104   }
00105   return area;
00106 }
00107 
00108 //======================================================================
00109 // sort out the areas
00110 void ClustSeqActAreaEG::_post_process() {
00111 
00112   // first check for danger signals.
00113   // Establish largest ghost transverse momentum
00114   _max_ghost_perp2 = 0.0;
00115   for (int i = 0; i < _initial_n; i++) {
00116     if (_is_pure_ghost[i] && _jets[i].perp2() > _max_ghost_perp2) 
00117       _max_ghost_perp2 = _jets[i].perp2();
00118   }
00119 
00120   // now find out if any of the particles are close to danger
00121   double danger_ratio = numeric_limits<double>::epsilon();
00122   danger_ratio = danger_ratio * danger_ratio;
00123   _has_dangerous_particles = false;
00124   for (int i = 0; i < _initial_n; i++) {
00125     if (!_is_pure_ghost[i] && 
00126         danger_ratio * _jets[i].perp2() <=  _max_ghost_perp2) {
00127       _has_dangerous_particles = true;
00128       break;
00129     }
00130   }
00131 
00132   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");
00133 
00134   // sort out sizes
00135   _areas.resize(_history.size());
00136   _area_4vectors.resize(_history.size());
00137   _is_pure_ghost.resize(_history.size());
00138   
00139   // First set up areas for the initial particles (ghost=_ghost_area,
00140   // real particles = 0); recall that _initial_n here is the number of
00141   // particles including ghosts
00142   for (int i = 0; i < _initial_n; i++) {
00143     if (_is_pure_ghost[i]) {
00144       _areas[i] = _ghost_area;
00145       // normalise pt to be _ghost_area (NB we make use of fact that
00146       // for initial particles, jet and clust_hist index are the same).
00147       _area_4vectors[i] = (_ghost_area/_jets[i].perp()) * _jets[i];
00148     } else {
00149       _areas[i] = 0;
00150       _area_4vectors[i] = PseudoJet(0.0,0.0,0.0,0.0);
00151     }
00152   }
00153   
00154   // next follow the branching through and set up the areas 
00155   // and ghost-nature at each step of the clustering (rather than
00156   // each jet).
00157   for (unsigned i = _initial_n; i < _history.size(); i++) {
00158     if (_history[i].parent2 == BeamJet) {
00159       _is_pure_ghost[i]  = _is_pure_ghost[_history[i].parent1];
00160       _areas[i]          = _areas[_history[i].parent1];
00161       _area_4vectors[i] = _area_4vectors[_history[i].parent1];
00162     } else {
00163       _is_pure_ghost[i]  = _is_pure_ghost[_history[i].parent1] && 
00164                            _is_pure_ghost[_history[i].parent2]   ;
00165       _areas[i]          = _areas[_history[i].parent1] + 
00166                            _areas[_history[i].parent2]  ;                          
00167       _jet_def.recombiner()->recombine(_area_4vectors[_history[i].parent1], 
00168                            _area_4vectors[_history[i].parent2],
00169                            _area_4vectors[i]);  
00170 //      _area_4vectors[i] = _area_4vectors[_history[i].parent1] + 
00171 //                          _area_4vectors[_history[i].parent2]  ;
00172     }
00173 
00174   }
00175   
00176 }
00177 
00178 FASTJET_END_NAMESPACE
00179 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines