[5] | 1 | // -*- C++ -*- |
---|
| 2 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 3 | // File: area.h // |
---|
| 4 | // Description: header file for the computation of jet area // |
---|
| 5 | // This file is part of the SISCone project. // |
---|
| 6 | // For more details, see http://projects.hepforge.org/siscone // |
---|
| 7 | // // |
---|
| 8 | // Copyright (c) 2006 Gavin Salam and Gregory Soyez // |
---|
| 9 | // // |
---|
| 10 | // This program is free software; you can redistribute it and/or modify // |
---|
| 11 | // it under the terms of the GNU General Public License as published by // |
---|
| 12 | // the Free Software Foundation; either version 2 of the License, or // |
---|
| 13 | // (at your option) any later version. // |
---|
| 14 | // // |
---|
| 15 | // This program is distributed in the hope that it will be useful, // |
---|
| 16 | // but WITHOUT ANY WARRANTY; without even the implied warranty of // |
---|
| 17 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // |
---|
| 18 | // GNU General Public License for more details. // |
---|
| 19 | // // |
---|
| 20 | // You should have received a copy of the GNU General Public License // |
---|
| 21 | // along with this program; if not, write to the Free Software // |
---|
| 22 | // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // |
---|
| 23 | // // |
---|
| 24 | // $Revision:: 859 $// |
---|
| 25 | // $Date:: 2012-11-28 02:49:23 +0100 (Wed, 28 Nov 2012) $// |
---|
| 26 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 27 | |
---|
| 28 | #include "area.h" |
---|
| 29 | #include "momentum.h" |
---|
| 30 | #include <stdlib.h> |
---|
| 31 | #include <iostream> |
---|
| 32 | |
---|
| 33 | namespace siscone{ |
---|
| 34 | using namespace std; |
---|
| 35 | |
---|
| 36 | /******************************************************* |
---|
| 37 | * Cjet_area implementation * |
---|
| 38 | * real Jet information, including its area(s) * |
---|
| 39 | * * |
---|
| 40 | * This class contains information for one single jet. * |
---|
| 41 | * That is, first, its momentum carrying information * |
---|
| 42 | * about its centre and pT, and second, its particle * |
---|
| 43 | * contents (from CJeT). * |
---|
| 44 | * Compared to the Cjet class, it also includes the * |
---|
| 45 | * passive and active areas of the jet computed using * |
---|
| 46 | * the Carea class. * |
---|
| 47 | *******************************************************/ |
---|
| 48 | |
---|
| 49 | // default ctor |
---|
| 50 | //-------------- |
---|
| 51 | Cjet_area::Cjet_area(){ |
---|
| 52 | active_area = passive_area = 0.0; |
---|
| 53 | } |
---|
| 54 | |
---|
| 55 | // jet-initiated ctor |
---|
| 56 | //------------------- |
---|
| 57 | Cjet_area::Cjet_area(Cjet &j){ |
---|
| 58 | v = j.v; |
---|
| 59 | n = j.n; |
---|
| 60 | contents = j.contents; |
---|
| 61 | |
---|
| 62 | pass = j.pass; |
---|
| 63 | |
---|
| 64 | pt_tilde = j.pt_tilde; |
---|
| 65 | sm_var2 = j.sm_var2; |
---|
| 66 | |
---|
| 67 | active_area = passive_area = 0.0; |
---|
| 68 | } |
---|
| 69 | |
---|
| 70 | // default dtor |
---|
| 71 | //-------------- |
---|
| 72 | Cjet_area::~Cjet_area(){ |
---|
| 73 | |
---|
| 74 | } |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | /****************************************************************** |
---|
| 78 | * Csiscone_area implementation * |
---|
| 79 | * class for the computation of jet areas. * |
---|
| 80 | * * |
---|
| 81 | * This is the class user should use whenever you want to compute * |
---|
| 82 | * the jet area (passive and active). * |
---|
| 83 | * It uses the SISCone algorithm to perform the jet analysis. * |
---|
| 84 | ******************************************************************/ |
---|
| 85 | |
---|
| 86 | // default ctor |
---|
| 87 | //------------- |
---|
| 88 | Carea::Carea(){ |
---|
| 89 | grid_size = 60; // 3600 particles added |
---|
| 90 | grid_eta_max = 6.0; // maybe having twice more points in eta than in phi should be nice |
---|
| 91 | grid_shift = 0.5; // 50% "shacking" |
---|
| 92 | |
---|
| 93 | pt_soft = 1e-100; |
---|
| 94 | pt_shift = 0.05; |
---|
| 95 | pt_soft_min = 1e-90; |
---|
| 96 | } |
---|
| 97 | |
---|
| 98 | // default dtor |
---|
| 99 | //------------- |
---|
| 100 | Carea::~Carea(){ |
---|
| 101 | |
---|
| 102 | } |
---|
| 103 | |
---|
| 104 | /* |
---|
| 105 | * compute the jet areas from a given particle set. |
---|
| 106 | * The parameters of this method are the ones which control the jet clustering alghorithm. |
---|
| 107 | * Note that the pt_min is not allowed here soince the jet-area determination involves soft |
---|
| 108 | * particles/jets and thus is used internally. |
---|
| 109 | * - _particles list of particles |
---|
| 110 | * - _radius cone radius |
---|
| 111 | * - _f shared energy threshold for splitting&merging |
---|
| 112 | * - _n_pass_max maximum number of passes (0=full search, the default) |
---|
| 113 | * - _split_merge_scale the scale choice for the split-merge procedure |
---|
| 114 | * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation. |
---|
| 115 | * SM_Et is IR safe but not boost invariant and not implemented(!) |
---|
| 116 | * SM_mt is IR safe for hadronic events, but not for decays of two |
---|
| 117 | * back-to-back particles of identical mass |
---|
| 118 | * SM_pttilde |
---|
| 119 | * is always IR safe, and also boost invariant (default) |
---|
| 120 | * - _hard_only when this is set on, only hard jets are computed |
---|
| 121 | * and not the purely ghosted jets (default: false) |
---|
| 122 | * return the jets together with their areas |
---|
| 123 | * The return value is the number of jets (including pure-ghost ones if they are included) |
---|
| 124 | ********************************************************************************************/ |
---|
| 125 | int Carea::compute_areas(std::vector<Cmomentum> &_particles, double _radius, double _f, |
---|
| 126 | int _n_pass_max, Esplit_merge_scale _split_merge_scale, |
---|
| 127 | bool _hard_only){ |
---|
| 128 | |
---|
| 129 | vector<Cmomentum> all_particles; |
---|
| 130 | |
---|
| 131 | // put "hardest cut-off" if needed |
---|
| 132 | // this avoids computation of ghosted jets when not required and |
---|
| 133 | // significantly shortens the SM. |
---|
| 134 | if (_hard_only){ |
---|
| 135 | SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min; |
---|
| 136 | } |
---|
| 137 | |
---|
| 138 | // clear potential previous runs |
---|
| 139 | jet_areas.clear(); |
---|
| 140 | |
---|
| 141 | // put initial set of particles in the list |
---|
| 142 | int n_hard = _particles.size(); |
---|
| 143 | all_particles = _particles; |
---|
| 144 | |
---|
| 145 | // build the set of ghost particles and add them to the particle list |
---|
| 146 | int i,j; |
---|
| 147 | double eta_g,phi_g,pt_g; |
---|
| 148 | |
---|
| 149 | for (i=0;i<grid_size;i++){ |
---|
| 150 | for (j=0;j<grid_size;j++){ |
---|
| 151 | 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); |
---|
| 152 | phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); |
---|
| 153 | pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0)))); |
---|
| 154 | 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))); |
---|
| 155 | } |
---|
| 156 | } |
---|
| 157 | |
---|
| 158 | // run clustering with all particles. |
---|
| 159 | // the split-merge here dynamically accounts for the purely soft jets. |
---|
| 160 | // we therefore end up with the active area for the jets |
---|
| 161 | int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, 0.0, _split_merge_scale); |
---|
| 162 | |
---|
| 163 | // save jets in the Cjet_area structure |
---|
| 164 | // and determine their size |
---|
| 165 | // jet contents is ordered by increasing index of the initial |
---|
| 166 | // particles. Hence, we look for the first particle with index >= n_hard |
---|
| 167 | // and deduce the number of ghosts in the jet, hence the area. |
---|
| 168 | double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); |
---|
| 169 | |
---|
| 170 | for (i=0;i<(int) jets.size();i++){ |
---|
| 171 | jet_areas.push_back(jets[i]); |
---|
| 172 | j=0; |
---|
| 173 | while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++; |
---|
| 174 | jet_areas[i].active_area = (jets[i].n-j)*area_factor; |
---|
| 175 | } |
---|
| 176 | |
---|
| 177 | // determine passive jet area |
---|
| 178 | // for that onem we use the pt_min cut in order to remove purely |
---|
| 179 | // soft jets from the SM procedure |
---|
| 180 | recompute_jets(_f, pt_soft_min); |
---|
| 181 | |
---|
| 182 | // for the area computation, we assume the jete order is the same! |
---|
| 183 | for (i=0;i<(int) jets.size();i++){ |
---|
| 184 | j=0; |
---|
| 185 | while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++; |
---|
| 186 | jet_areas[i].passive_area = (jets[i].n-j)*area_factor; |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | // Note: |
---|
| 190 | // there surely remain purely soft je at the end |
---|
| 191 | // their active area is 0 by default (see ctor) |
---|
| 192 | |
---|
| 193 | jets.clear(); |
---|
| 194 | |
---|
| 195 | return n_jets; |
---|
| 196 | } |
---|
| 197 | |
---|
| 198 | /* |
---|
| 199 | * compute the passive jet areas from a given particle set. |
---|
| 200 | * The parameters of this method are the ones which control the jet clustering alghorithm. |
---|
| 201 | * Note that the pt_min is not allowed here soince the jet-area determination involves soft |
---|
| 202 | * particles/jets and thus is used internally. |
---|
| 203 | * - _particles list of particles |
---|
| 204 | * - _radius cone radius |
---|
| 205 | * - _f shared energy threshold for splitting&merging |
---|
| 206 | * - _n_pass_max maximum number of passes (0=full search, the default) |
---|
| 207 | * - _split_merge_scale the scale choice for the split-merge procedure |
---|
| 208 | * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation. |
---|
| 209 | * SM_Et is IR safe but not boost invariant and not implemented(!) |
---|
| 210 | * SM_mt is IR safe for hadronic events, but not for decays of two |
---|
| 211 | * back-to-back particles of identical mass |
---|
| 212 | * SM_pttilde |
---|
| 213 | * is always IR safe, and also boost invariant (default) |
---|
| 214 | * return the jets together with their passive areas |
---|
| 215 | * The return value is the number of jets |
---|
| 216 | ********************************************************************************************/ |
---|
| 217 | int Carea::compute_passive_areas(std::vector<Cmomentum> &_particles, double _radius, double _f, |
---|
| 218 | int _n_pass_max, Esplit_merge_scale _split_merge_scale){ |
---|
| 219 | |
---|
| 220 | vector<Cmomentum> all_particles; |
---|
| 221 | |
---|
| 222 | // in the case of passive area, we do not need |
---|
| 223 | // to put the ghosts in the stable-cone search |
---|
| 224 | // (they do no influence the list of stable cones) |
---|
| 225 | // Here's how it goes... |
---|
| 226 | stable_cone_soft_pt2_cutoff = pt_soft_min*pt_soft_min; |
---|
| 227 | |
---|
| 228 | // clear potential previous runs |
---|
| 229 | jet_areas.clear(); |
---|
| 230 | |
---|
| 231 | // put initial set of particles in the list |
---|
| 232 | int n_hard = _particles.size(); |
---|
| 233 | all_particles = _particles; |
---|
| 234 | |
---|
| 235 | // build the set of ghost particles and add them to the particle list |
---|
| 236 | int i,j; |
---|
| 237 | double eta_g,phi_g,pt_g; |
---|
| 238 | |
---|
| 239 | for (i=0;i<grid_size;i++){ |
---|
| 240 | for (j=0;j<grid_size;j++){ |
---|
| 241 | 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); |
---|
| 242 | phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); |
---|
| 243 | pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0)))); |
---|
| 244 | 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))); |
---|
| 245 | } |
---|
| 246 | } |
---|
| 247 | |
---|
| 248 | // determine passive jet area |
---|
| 249 | // for that onem we use the pt_min cut in order to remove purely |
---|
| 250 | // soft jets from the SM procedure |
---|
| 251 | int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, pt_soft_min, _split_merge_scale); |
---|
| 252 | |
---|
| 253 | // save jets in the Cjet_area structure |
---|
| 254 | // and determine their size |
---|
| 255 | // jet contents is ordered by increasing index of the initial |
---|
| 256 | // particles. Hence, we look for the first particle with index >= n_hard |
---|
| 257 | // and deduce the number of ghosts in the jet, hence the area. |
---|
| 258 | double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); |
---|
| 259 | for (i=0;i<(int) jets.size();i++){ |
---|
| 260 | j=0; |
---|
| 261 | while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++; |
---|
| 262 | jet_areas[i].passive_area = (jets[i].n-j)*area_factor; |
---|
| 263 | } |
---|
| 264 | |
---|
| 265 | jets.clear(); |
---|
| 266 | |
---|
| 267 | return n_jets; |
---|
| 268 | } |
---|
| 269 | |
---|
| 270 | /* |
---|
| 271 | * compute the active jet areas from a given particle set. |
---|
| 272 | * The parameters of this method are the ones which control the jet clustering alghorithm. |
---|
| 273 | * Note that the pt_min is not allowed here soince the jet-area determination involves soft |
---|
| 274 | * particles/jets and thus is used internally. |
---|
| 275 | * - _particles list of particles |
---|
| 276 | * - _radius cone radius |
---|
| 277 | * - _f shared energy threshold for splitting&merging |
---|
| 278 | * - _n_pass_max maximum number of passes (0=full search, the default) |
---|
| 279 | * - _split_merge_scale the scale choice for the split-merge procedure |
---|
| 280 | * NOTE: SM_pt leads to IR unsafety for some events with momentum conservation. |
---|
| 281 | * SM_Et is IR safe but not boost invariant and not implemented(!) |
---|
| 282 | * SM_mt is IR safe for hadronic events, but not for decays of two |
---|
| 283 | * back-to-back particles of identical mass |
---|
| 284 | * SM_pttilde |
---|
| 285 | * is always IR safe, and also boost invariant (default) |
---|
| 286 | * - _hard_only when this is set on, only hard jets are computed |
---|
| 287 | * and not the purely ghosted jets (default: false) |
---|
| 288 | * return the jets together with their active areas |
---|
| 289 | * The return value is the number of jets (including pure-ghost ones if they are included) |
---|
| 290 | ********************************************************************************************/ |
---|
| 291 | int Carea::compute_active_areas(std::vector<Cmomentum> &_particles, double _radius, double _f, |
---|
| 292 | int _n_pass_max, Esplit_merge_scale _split_merge_scale, |
---|
| 293 | bool _hard_only){ |
---|
| 294 | |
---|
| 295 | vector<Cmomentum> all_particles; |
---|
| 296 | |
---|
| 297 | // put "hardest cut-off" if needed |
---|
| 298 | // this avoids computation of ghosted jets when not required and |
---|
| 299 | // significantly shortens the SM. |
---|
| 300 | if (_hard_only){ |
---|
| 301 | SM_var2_hardest_cut_off = pt_soft_min*pt_soft_min; |
---|
| 302 | } |
---|
| 303 | |
---|
| 304 | // clear potential previous runs |
---|
| 305 | jet_areas.clear(); |
---|
| 306 | |
---|
| 307 | // put initial set of particles in the list |
---|
| 308 | int n_hard = _particles.size(); |
---|
| 309 | all_particles = _particles; |
---|
| 310 | |
---|
| 311 | // build the set of ghost particles and add them to the particle list |
---|
| 312 | int i,j; |
---|
| 313 | double eta_g,phi_g,pt_g; |
---|
| 314 | |
---|
| 315 | for (i=0;i<grid_size;i++){ |
---|
| 316 | for (j=0;j<grid_size;j++){ |
---|
| 317 | 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); |
---|
| 318 | phi_g = M_PI *(-1.0+2.0*(j+0.5+grid_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0))))/grid_size); |
---|
| 319 | pt_g = pt_soft*(1.0+pt_shift*(-1.0+2.0*(rand()/(RAND_MAX+1.0)))); |
---|
| 320 | 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))); |
---|
| 321 | } |
---|
| 322 | } |
---|
| 323 | |
---|
| 324 | // run clustering with all particles. |
---|
| 325 | // the split-merge here dynamically accounts for the purely soft jets. |
---|
| 326 | // we therefore end up with the active area for the jets |
---|
| 327 | int n_jets = compute_jets(all_particles, _radius, _f, _n_pass_max, 0.0, _split_merge_scale); |
---|
| 328 | |
---|
| 329 | // save jets in the Cjet_area structure |
---|
| 330 | // and determine their size |
---|
| 331 | // jet contents is ordered by increasing index of the initial |
---|
| 332 | // particles. Hence, we look for the first particle with index >= n_hard |
---|
| 333 | // and deduce the number of ghosts in the jet, hence the area. |
---|
| 334 | double area_factor = (2.0*grid_eta_max/grid_size)*(twopi/grid_size); |
---|
| 335 | |
---|
| 336 | for (i=0;i<(int) jets.size();i++){ |
---|
| 337 | jet_areas.push_back(jets[i]); |
---|
| 338 | j=0; |
---|
| 339 | while ((j<jets[i].n) && (jets[i].contents[j]<n_hard)) j++; |
---|
| 340 | jet_areas[i].active_area = (jets[i].n-j)*area_factor; |
---|
| 341 | } |
---|
| 342 | |
---|
| 343 | jets.clear(); |
---|
| 344 | |
---|
| 345 | return n_jets; |
---|
| 346 | } |
---|
| 347 | |
---|
| 348 | } |
---|