00001 //STARTHEADER 00002 // $Id: ClusterSequenceActiveAreaExplicitGhosts.cc 957 2007-11-22 18:58:45Z 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 & area_spec) { 00048 00049 // add the ghosts to the jets 00050 area_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 area_spec 00059 _ghost_area = area_spec.actual_ghost_area(); 00060 _n_ghosts = area_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