|
SISCone
2.0.5
|
00001 // -*- C++ -*- 00003 // File: area.h // 00004 // Description: header file for the computation of jet area // 00005 // This file is part of the SISCone project. // 00006 // For more details, see http://projects.hepforge.org/siscone // 00007 // // 00008 // Copyright (c) 2006 Gavin Salam and Gregory Soyez // 00009 // // 00010 // This program is free software; you can redistribute it and/or modify // 00011 // it under the terms of the GNU General Public License as published by // 00012 // the Free Software Foundation; either version 2 of the License, or // 00013 // (at your option) any later version. // 00014 // // 00015 // This program is distributed in the hope that it will be useful, // 00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of // 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // 00018 // GNU General Public License for more details. // 00019 // // 00020 // You should have received a copy of the GNU General Public License // 00021 // along with this program; if not, write to the Free Software // 00022 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // 00023 // // 00024 // $Revision:: 149 $// 00025 // $Date:: 2007-03-15 00:13:58 +0100 (Thu, 15 Mar 2007) $// 00027 00028 #include "area.h" 00029 #include "momentum.h" 00030 #include <stdlib.h> 00031 #include <iostream> 00032 00033 namespace siscone{ 00034 using namespace std; 00035 00036 /******************************************************* 00037 * Cjet_area implementation * 00038 * real Jet information, including its area(s) * 00039 * * 00040 * This class contains information for one single jet. * 00041 * That is, first, its momentum carrying information * 00042 * about its centre and pT, and second, its particle * 00043 * contents (from CJeT). * 00044 * Compared to the Cjet class, it also includes the * 00045 * passive and active areas of the jet computed using * 00046 * the Carea class. * 00047 *******************************************************/ 00048 00049 // default ctor 00050 //-------------- 00051 Cjet_area::Cjet_area(){ 00052 active_area = passive_area = 0.0; 00053 } 00054 00055 // jet-initiated ctor 00056 //------------------- 00057 Cjet_area::Cjet_area(Cjet &j){ 00058 v = j.v; 00059 n = j.n; 00060 contents = j.contents; 00061 00062 pass = j.pass; 00063 00064 pt_tilde = j.pt_tilde; 00065 sm_var2 = j.sm_var2; 00066 00067 active_area = passive_area = 0.0; 00068 } 00069 00070 // default dtor 00071 //-------------- 00072 Cjet_area::~Cjet_area(){ 00073 00074 } 00075 00076 00077 /****************************************************************** 00078 * Csiscone_area implementation * 00079 * class for the computation of jet areas. * 00080 * * 00081 * This is the class user should use whenever you want to compute * 00082 * the jet area (passive and active). * 00083 * It uses the SISCone algorithm to perform the jet analysis. * 00084 ******************************************************************/ 00085 00086 // default ctor 00087 //------------- 00088 Carea::Carea(){ 00089 grid_size = 60; // 3600 particles added 00090 grid_eta_max = 6.0; // maybe having twice more points in eta than in phi should be nice 00091 grid_shift = 0.5; // 50% "shacking" 00092 00093 pt_soft = 1e-100; 00094 pt_shift = 0.05; 00095 pt_soft_min = 1e-90; 00096 } 00097 00098 // default dtor 00099 //------------- 00100 Carea::~Carea(){ 00101 00102 } 00103 00104 /* 00105 * compute the jet areas from a given particle set. 00106 * The parameters of this method are the ones which control the jet clustering alghorithm. 00107 * Note that the pt_min is not allowed here soince the jet-area determination involves soft 00108 * particles/jets and thus is used internally. 00109 * - _particles list of particles 00110 * - _radius cone radius 00111 * - _f shared energy threshold for splitting&merging 00112 * - _n_pass_max maximum number of passes (0=full search, the default) 00113 * - _split_merge_scale the scale choice for the split-merge procedure 00114 * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation. 00115 * SM_Et is IR safe but not boost invariant and not implemented(!) 00116 * SM_mt is IR safe for hadronic events, but not for decays of two 00117 * back-to-back particles of identical mass 00118 * SM_pttilde 00119 * is always IR safe, and also boost invariant (default) 00120 * - _hard_only when this is set on, only hard jets are computed 00121 * and not the purely ghosted jets (default: false) 00122 * return the jets together with their areas 00123 * The return value is the number of jets (including pure-ghost ones if they are included) 00124 ********************************************************************************************/ 00125 int Carea::compute_areas(std::vector<Cmomentum> &_particles, double _radius, double _f, 00126 int _n_pass_max, Esplit_merge_scale _split_merge_scale, 00127 bool _hard_only){ 00128 00129 vector<Cmomentum> all_particles; 00130 00131 // put "hardest cut-off" if needed 00132 // this avoids computation of ghosted jets when not required and 00133 // significantly shortens the SM. 00134 if (_hard_only){ 00135 SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min; 00136 } 00137 00138 // clear potential previous runs 00139 jet_areas.clear(); 00140 00141 // put initial set of particles in the list 00142 int n_hard = _particles.size(); 00143 all_particles = _particles; 00144 00145 // build the set of ghost particles and add them to the particle list 00146 int i,j; 00147 double eta_g,phi_g,pt_g; 00148 00149 for (i=0;i<grid_size;i++){ 00150 for (j=0;j<grid_size;j++){ 00151 eta_g = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); 00152 phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); 00153 pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0)))); 00154 all_particles.push_back(Cmomentum(pt_g*cos(phi_g),pt_g*sin(phi_g),pt_g*sinh(eta_g),pt_g*cosh(eta_g))); 00155 } 00156 } 00157 00158 // run clustering with all particles. 00159 // the split-merge here dynamically accounts for the purely soft jets. 00160 // we therefore end up with the active area for the jets 00161 int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, 0.0, _split_merge_scale); 00162 00163 // save jets in the Cjet_area structure 00164 // and determine their size 00165 // jet contents is ordered by increasing index of the initial 00166 // particles. Hence, we look for the first particle with index >= n_hard 00167 // and deduce the number of ghosts in the jet, hence the area. 00168 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); 00169 00170 for (i=0;i<(int) jets.size();i++){ 00171 jet_areas.push_back(jets[i]); 00172 j=0; 00173 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++; 00174 jet_areas[i].active_area = (jets[i].n-j)*area_factor; 00175 } 00176 00177 // determine passive jet area 00178 // for that onem we use the pt_min cut in order to remove purely 00179 // soft jets from the SM procedure 00180 recompute_jets(_f, pt_soft_min); 00181 00182 // for the area computation, we assume the jete order is the same! 00183 for (i=0;i<(int) jets.size();i++){ 00184 j=0; 00185 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++; 00186 jet_areas[i].passive_area = (jets[i].n-j)*area_factor; 00187 } 00188 00189 // Note: 00190 // there surely remain purely soft je at the end 00191 // their active area is 0 by default (see ctor) 00192 00193 jets.clear(); 00194 00195 return n_jets; 00196 } 00197 00198 /* 00199 * compute the passive jet areas from a given particle set. 00200 * The parameters of this method are the ones which control the jet clustering alghorithm. 00201 * Note that the pt_min is not allowed here soince the jet-area determination involves soft 00202 * particles/jets and thus is used internally. 00203 * - _particles list of particles 00204 * - _radius cone radius 00205 * - _f shared energy threshold for splitting&merging 00206 * - _n_pass_max maximum number of passes (0=full search, the default) 00207 * - _split_merge_scale the scale choice for the split-merge procedure 00208 * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation. 00209 * SM_Et is IR safe but not boost invariant and not implemented(!) 00210 * SM_mt is IR safe for hadronic events, but not for decays of two 00211 * back-to-back particles of identical mass 00212 * SM_pttilde 00213 * is always IR safe, and also boost invariant (default) 00214 * return the jets together with their passive areas 00215 * The return value is the number of jets 00216 ********************************************************************************************/ 00217 int Carea::compute_passive_areas(std::vector<Cmomentum> &_particles, double _radius, double _f, 00218 int _n_pass_max, Esplit_merge_scale _split_merge_scale){ 00219 00220 vector<Cmomentum> all_particles; 00221 00222 // in the case of passive area, we do not need 00223 // to put the ghosts in the stable-cone search 00224 // (they do no influence the list of stable cones) 00225 // Here's how it goes... 00226 stable_cone_soft_pt2_cutoff = pt_soft_min*pt_soft_min; 00227 00228 // clear potential previous runs 00229 jet_areas.clear(); 00230 00231 // put initial set of particles in the list 00232 int n_hard = _particles.size(); 00233 all_particles = _particles; 00234 00235 // build the set of ghost particles and add them to the particle list 00236 int i,j; 00237 double eta_g,phi_g,pt_g; 00238 00239 for (i=0;i<grid_size;i++){ 00240 for (j=0;j<grid_size;j++){ 00241 eta_g = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); 00242 phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); 00243 pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0)))); 00244 all_particles.push_back(Cmomentum(pt_g*cos(phi_g),pt_g*sin(phi_g),pt_g*sinh(eta_g),pt_g*cosh(eta_g))); 00245 } 00246 } 00247 00248 // determine passive jet area 00249 // for that onem we use the pt_min cut in order to remove purely 00250 // soft jets from the SM procedure 00251 int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, pt_soft_min, _split_merge_scale); 00252 00253 // save jets in the Cjet_area structure 00254 // and determine their size 00255 // jet contents is ordered by increasing index of the initial 00256 // particles. Hence, we look for the first particle with index >= n_hard 00257 // and deduce the number of ghosts in the jet, hence the area. 00258 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); 00259 for (i=0;i<(int) jets.size();i++){ 00260 j=0; 00261 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++; 00262 jet_areas[i].passive_area = (jets[i].n-j)*area_factor; 00263 } 00264 00265 jets.clear(); 00266 00267 return n_jets; 00268 } 00269 00270 /* 00271 * compute the active jet areas from a given particle set. 00272 * The parameters of this method are the ones which control the jet clustering alghorithm. 00273 * Note that the pt_min is not allowed here soince the jet-area determination involves soft 00274 * particles/jets and thus is used internally. 00275 * - _particles list of particles 00276 * - _radius cone radius 00277 * - _f shared energy threshold for splitting&merging 00278 * - _n_pass_max maximum number of passes (0=full search, the default) 00279 * - _split_merge_scale the scale choice for the split-merge procedure 00280 * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation. 00281 * SM_Et is IR safe but not boost invariant and not implemented(!) 00282 * SM_mt is IR safe for hadronic events, but not for decays of two 00283 * back-to-back particles of identical mass 00284 * SM_pttilde 00285 * is always IR safe, and also boost invariant (default) 00286 * - _hard_only when this is set on, only hard jets are computed 00287 * and not the purely ghosted jets (default: false) 00288 * return the jets together with their active areas 00289 * The return value is the number of jets (including pure-ghost ones if they are included) 00290 ********************************************************************************************/ 00291 int Carea::compute_active_areas(std::vector<Cmomentum> &_particles, double _radius, double _f, 00292 int _n_pass_max, Esplit_merge_scale _split_merge_scale, 00293 bool _hard_only){ 00294 00295 vector<Cmomentum> all_particles; 00296 00297 // put "hardest cut-off" if needed 00298 // this avoids computation of ghosted jets when not required and 00299 // significantly shortens the SM. 00300 if (_hard_only){ 00301 SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min; 00302 } 00303 00304 // clear potential previous runs 00305 jet_areas.clear(); 00306 00307 // put initial set of particles in the list 00308 int n_hard = _particles.size(); 00309 all_particles = _particles; 00310 00311 // build the set of ghost particles and add them to the particle list 00312 int i,j; 00313 double eta_g,phi_g,pt_g; 00314 00315 for (i=0;i<grid_size;i++){ 00316 for (j=0;j<grid_size;j++){ 00317 eta_g = grid_eta_max*(-1.0+2.0*(i+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); 00318 phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); 00319 pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0)))); 00320 all_particles.push_back(Cmomentum(pt_g*cos(phi_g),pt_g*sin(phi_g),pt_g*sinh(eta_g),pt_g*cosh(eta_g))); 00321 } 00322 } 00323 00324 // run clustering with all particles. 00325 // the split-merge here dynamically accounts for the purely soft jets. 00326 // we therefore end up with the active area for the jets 00327 int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, 0.0, _split_merge_scale); 00328 00329 // save jets in the Cjet_area structure 00330 // and determine their size 00331 // jet contents is ordered by increasing index of the initial 00332 // particles. Hence, we look for the first particle with index >= n_hard 00333 // and deduce the number of ghosts in the jet, hence the area. 00334 double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); 00335 00336 for (i=0;i<(int) jets.size();i++){ 00337 jet_areas.push_back(jets[i]); 00338 j=0; 00339 while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++; 00340 jet_areas[i].active_area = (jets[i].n-j)*area_factor; 00341 } 00342 00343 jets.clear(); 00344 00345 return n_jets; 00346 } 00347 00348 }