FastJet 3.0.2
|
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