[5] | 1 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2 | // File: split_merge.cpp // |
---|
| 3 | // Description: source file for splitting/merging (contains the CJet class) // |
---|
| 4 | // This file is part of the SISCone project. // |
---|
| 5 | // For more details, see http://projects.hepforge.org/siscone // |
---|
| 6 | // // |
---|
| 7 | // Copyright (c) 2006 Gavin Salam and Gregory Soyez // |
---|
| 8 | // // |
---|
| 9 | // This program is free software; you can redistribute it and/or modify // |
---|
| 10 | // it under the terms of the GNU General Public License as published by // |
---|
| 11 | // the Free Software Foundation; either version 2 of the License, or // |
---|
| 12 | // (at your option) any later version. // |
---|
| 13 | // // |
---|
| 14 | // This program is distributed in the hope that it will be useful, // |
---|
| 15 | // but WITHOUT ANY WARRANTY; without even the implied warranty of // |
---|
| 16 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // |
---|
| 17 | // GNU General Public License for more details. // |
---|
| 18 | // // |
---|
| 19 | // You should have received a copy of the GNU General Public License // |
---|
| 20 | // along with this program; if not, write to the Free Software // |
---|
| 21 | // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // |
---|
| 22 | // // |
---|
| 23 | // $Revision:: 859 $// |
---|
| 24 | // $Date:: 2012-11-28 02:49:23 +0100 (Wed, 28 Nov 2012) $// |
---|
| 25 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 26 | |
---|
| 27 | #include "split_merge.h" |
---|
| 28 | #include "siscone_error.h" |
---|
| 29 | #include "momentum.h" |
---|
| 30 | #include <math.h> |
---|
| 31 | #include <limits> // for max |
---|
| 32 | #include <iostream> |
---|
| 33 | #include <algorithm> |
---|
| 34 | #include <sstream> |
---|
| 35 | #include <cassert> |
---|
| 36 | |
---|
| 37 | namespace siscone{ |
---|
| 38 | |
---|
| 39 | using namespace std; |
---|
| 40 | |
---|
| 41 | /******************************************************** |
---|
| 42 | * class Cjet implementation * |
---|
| 43 | * real Jet information. * |
---|
| 44 | * This class contains information for one single jet. * |
---|
| 45 | * That is, first, its momentum carrying information * |
---|
| 46 | * about its centre and pT, and second, its particle * |
---|
| 47 | * contents * |
---|
| 48 | ********************************************************/ |
---|
| 49 | // default ctor |
---|
| 50 | //-------------- |
---|
| 51 | Cjet::Cjet(){ |
---|
| 52 | n = 0; |
---|
| 53 | v = Cmomentum(); |
---|
| 54 | pt_tilde = 0.0; |
---|
| 55 | sm_var2 = 0.0; |
---|
| 56 | } |
---|
| 57 | |
---|
| 58 | // default dtor |
---|
| 59 | //-------------- |
---|
| 60 | Cjet::~Cjet(){ |
---|
| 61 | |
---|
| 62 | } |
---|
| 63 | |
---|
| 64 | // ordering of jets in pt (e.g. used in final jets ordering) |
---|
| 65 | //----------------------------------------------------------- |
---|
| 66 | bool jets_pt_less(const Cjet &j1, const Cjet &j2){ |
---|
| 67 | return j1.v.perp2() > j2.v.perp2(); |
---|
| 68 | } |
---|
| 69 | |
---|
| 70 | |
---|
| 71 | /******************************************************** |
---|
| 72 | * Csplit_merge_ptcomparison implementation * |
---|
| 73 | * This deals with the ordering of the jets candidates * |
---|
| 74 | ********************************************************/ |
---|
| 75 | |
---|
| 76 | // odering of two jets |
---|
| 77 | // The variable of the ordering is pt or mt |
---|
| 78 | // depending on 'split_merge_scale' choice |
---|
| 79 | // |
---|
| 80 | // with EPSILON_SPLITMERGE defined, this attempts to identify |
---|
| 81 | // delicate cases where two jets have identical momenta except for |
---|
| 82 | // softish particles -- the difference of pt's may not be correctly |
---|
| 83 | // identified normally and so lead to problems for the fate of the |
---|
| 84 | // softish particle. |
---|
| 85 | // |
---|
| 86 | // NB: there is a potential issue in momentum-conserving events, |
---|
| 87 | // whereby the harder of two jets becomes ill-defined when a soft |
---|
| 88 | // particle is emitted --- this may have a knock-on effect on |
---|
| 89 | // subsequent split-merge steps in cases with sufficiently large R |
---|
| 90 | // (but we don't know what the limit is...) |
---|
| 91 | //------------------------------------------------------------------ |
---|
| 92 | bool Csplit_merge_ptcomparison::operator ()(const Cjet &jet1, const Cjet &jet2) const{ |
---|
| 93 | double q1, q2; |
---|
| 94 | |
---|
| 95 | // compute the value for comparison for both jets |
---|
| 96 | // This depends on the choice of variable (mt is the default) |
---|
| 97 | q1 = jet1.sm_var2; |
---|
| 98 | q2 = jet2.sm_var2; |
---|
| 99 | |
---|
| 100 | bool res = q1 > q2; |
---|
| 101 | |
---|
| 102 | // if we enable the refined version of the comparison (see defines.h), |
---|
| 103 | // we compute the difference more precisely when the two jets are very |
---|
| 104 | // close in the ordering variable. |
---|
| 105 | #ifdef EPSILON_SPLITMERGE |
---|
| 106 | if ( (fabs(q1-q2) < EPSILON_SPLITMERGE*max(q1,q2)) && |
---|
| 107 | (jet1.v.ref != jet2.v.ref) ) { |
---|
| 108 | // get the momentum of the difference |
---|
| 109 | Cmomentum difference; |
---|
| 110 | double pt_tilde_difference; |
---|
| 111 | get_difference(jet1,jet2,&difference,&pt_tilde_difference); |
---|
| 112 | |
---|
| 113 | // use the following relation: pt1^2 - pt2^2 = (pt1+pt2)*(pt1-pt2) |
---|
| 114 | double qdiff; |
---|
| 115 | Cmomentum sum = jet1.v ; |
---|
| 116 | sum += jet2.v; |
---|
| 117 | double pt_tilde_sum = jet1.pt_tilde + jet2.pt_tilde; |
---|
| 118 | |
---|
| 119 | // depending on the choice of ordering variable, set the result |
---|
| 120 | switch (split_merge_scale){ |
---|
| 121 | case SM_mt: |
---|
| 122 | qdiff = sum.E*difference.E - sum.pz*difference.pz; |
---|
| 123 | break; |
---|
| 124 | case SM_pt: |
---|
| 125 | qdiff = sum.px*difference.px + sum.py*difference.py; |
---|
| 126 | break; |
---|
| 127 | case SM_pttilde: |
---|
| 128 | qdiff = pt_tilde_sum*pt_tilde_difference; |
---|
| 129 | break; |
---|
| 130 | case SM_Et: |
---|
| 131 | // diff = E^2 (dpt^2 pz^2- pt^2 dpz^2) |
---|
| 132 | // + dE^2 (pt^2+pz^2) pt2^2 |
---|
| 133 | // where, unless explicitely specified the variables |
---|
| 134 | // refer to the first jet or differences jet1-jet2. |
---|
| 135 | qdiff = jet1.v.E*jet1.v.E* |
---|
| 136 | ((sum.px*difference.px + sum.py*difference.py)*jet1.v.pz*jet1.v.pz |
---|
| 137 | -jet1.v.perp2()*sum.pz*difference.pz) |
---|
| 138 | +sum.E*difference.E*(jet1.v.perp2()+jet1.v.pz*jet1.v.pz)*jet2.v.perp2(); |
---|
| 139 | break; |
---|
| 140 | default: |
---|
| 141 | throw Csiscone_error("Unsupported split-merge scale choice: " |
---|
| 142 | + SM_scale_name()); |
---|
| 143 | } |
---|
| 144 | res = qdiff > 0; |
---|
| 145 | } |
---|
| 146 | #endif // EPSILON_SPLITMERGE |
---|
| 147 | |
---|
| 148 | return res; |
---|
| 149 | } |
---|
| 150 | |
---|
| 151 | |
---|
| 152 | /// return a name for the sm scale choice |
---|
| 153 | /// NB: used internally and also by fastjet |
---|
| 154 | std::string split_merge_scale_name(Esplit_merge_scale sms) { |
---|
| 155 | switch(sms) { |
---|
| 156 | case SM_pt: |
---|
| 157 | return "pt (IR unsafe)"; |
---|
| 158 | case SM_Et: |
---|
| 159 | return "Et (boost dep.)"; |
---|
| 160 | case SM_mt: |
---|
| 161 | return "mt (IR safe except for pairs of identical decayed heavy particles)"; |
---|
| 162 | case SM_pttilde: |
---|
| 163 | return "pttilde (scalar sum of pt's)"; |
---|
| 164 | default: |
---|
| 165 | return "[SM scale without a name]"; |
---|
| 166 | } |
---|
| 167 | } |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | // get the difference between 2 jets |
---|
| 171 | // - j1 first jet |
---|
| 172 | // - j2 second jet |
---|
| 173 | // - v jet1-jet2 |
---|
| 174 | // - pt_tilde jet1-jet2 pt_tilde |
---|
| 175 | // return true if overlapping, false if disjoint |
---|
| 176 | //----------------------------------------------- |
---|
| 177 | void Csplit_merge_ptcomparison::get_difference(const Cjet &j1, const Cjet &j2, Cmomentum *v, double *pt_tilde) const { |
---|
| 178 | int i1,i2; |
---|
| 179 | |
---|
| 180 | // initialise |
---|
| 181 | i1=i2=0; |
---|
| 182 | *v = Cmomentum(); |
---|
| 183 | *pt_tilde = 0.0; |
---|
| 184 | |
---|
| 185 | // compute overlap |
---|
| 186 | // at the same time, we store union in indices |
---|
| 187 | do{ |
---|
| 188 | if (j1.contents[i1]==j2.contents[i2]) { |
---|
| 189 | i1++; |
---|
| 190 | i2++; |
---|
| 191 | } else if (j1.contents[i1]<j2.contents[i2]){ |
---|
| 192 | (*v) += (*particles)[j1.contents[i1]]; |
---|
| 193 | (*pt_tilde) += (*pt)[j1.contents[i1]]; |
---|
| 194 | i1++; |
---|
| 195 | } else if (j1.contents[i1]>j2.contents[i2]){ |
---|
| 196 | (*v) -= (*particles)[j2.contents[i2]]; |
---|
| 197 | (*pt_tilde) -= (*pt)[j2.contents[i2]]; |
---|
| 198 | i2++; |
---|
| 199 | } else { |
---|
| 200 | throw Csiscone_error("get_non_overlap reached part it should never have seen..."); |
---|
| 201 | } |
---|
| 202 | } while ((i1<j1.n) && (i2<j2.n)); |
---|
| 203 | |
---|
| 204 | // deal with particles at the end of the list... |
---|
| 205 | while (i1 < j1.n) { |
---|
| 206 | (*v) += (*particles)[j1.contents[i1]]; |
---|
| 207 | (*pt_tilde) += (*pt)[j1.contents[i1]]; |
---|
| 208 | i1++; |
---|
| 209 | } |
---|
| 210 | while (i2 < j2.n) { |
---|
| 211 | (*v) -= (*particles)[j2.contents[i2]]; |
---|
| 212 | (*pt_tilde) -= (*pt)[j2.contents[i2]]; |
---|
| 213 | i2++; |
---|
| 214 | } |
---|
| 215 | } |
---|
| 216 | |
---|
| 217 | |
---|
| 218 | /******************************************************** |
---|
| 219 | * class Csplit_merge implementation * |
---|
| 220 | * Class used to split and merge jets. * |
---|
| 221 | ********************************************************/ |
---|
| 222 | // default ctor |
---|
| 223 | //-------------- |
---|
| 224 | Csplit_merge::Csplit_merge(){ |
---|
| 225 | merge_identical_protocones = false; |
---|
| 226 | #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES |
---|
| 227 | #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE |
---|
| 228 | merge_identical_protocones = true; |
---|
| 229 | #endif |
---|
| 230 | #endif |
---|
| 231 | indices = NULL; |
---|
| 232 | |
---|
| 233 | // ensure that ptcomparison points to our set of particles (though params not correct) |
---|
| 234 | ptcomparison.particles = &particles; |
---|
| 235 | ptcomparison.pt = &pt; |
---|
| 236 | candidates.reset(new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison)); |
---|
| 237 | |
---|
| 238 | // no hardest cut (col-unsafe) |
---|
| 239 | SM_var2_hardest_cut_off = -1.0; |
---|
| 240 | |
---|
| 241 | // no pt cutoff for the particles to put in p_uncol_hard |
---|
| 242 | stable_cone_soft_pt2_cutoff = -1.0; |
---|
| 243 | |
---|
| 244 | // no pt-weighted splitting |
---|
| 245 | use_pt_weighted_splitting = false; |
---|
| 246 | } |
---|
| 247 | |
---|
| 248 | |
---|
| 249 | // default dtor |
---|
| 250 | //-------------- |
---|
| 251 | Csplit_merge::~Csplit_merge(){ |
---|
| 252 | full_clear(); |
---|
| 253 | } |
---|
| 254 | |
---|
| 255 | |
---|
| 256 | // initialisation function |
---|
| 257 | // - _particles list of particles |
---|
| 258 | // - protocones list of protocones (initial jet candidates) |
---|
| 259 | // - R2 cone radius (squared) |
---|
| 260 | // - ptmin minimal pT allowed for jets |
---|
| 261 | //------------------------------------------------------------- |
---|
| 262 | int Csplit_merge::init(vector<Cmomentum> & /*_particles*/, vector<Cmomentum> *protocones, double R2, double ptmin){ |
---|
| 263 | // browse protocones |
---|
| 264 | return add_protocones(protocones, R2, ptmin); |
---|
| 265 | } |
---|
| 266 | |
---|
| 267 | |
---|
| 268 | // initialisation function for particle list |
---|
| 269 | // - _particles list of particles |
---|
| 270 | //------------------------------------------------------------- |
---|
| 271 | int Csplit_merge::init_particles(vector<Cmomentum> &_particles){ |
---|
| 272 | full_clear(); |
---|
| 273 | |
---|
| 274 | // compute the list of particles |
---|
| 275 | // here, we do not need to throw away particles |
---|
| 276 | // with infinite rapidity (colinear with the beam) |
---|
| 277 | particles = _particles; |
---|
| 278 | n = particles.size(); |
---|
| 279 | |
---|
| 280 | // build the vector of particles' pt |
---|
| 281 | pt.resize(n); |
---|
| 282 | for (int i=0;i<n;i++) |
---|
| 283 | pt[i] = particles[i].perp(); |
---|
| 284 | |
---|
| 285 | // ensure that ptcomparison points to our set of particles (though params not correct) |
---|
| 286 | ptcomparison.particles = &particles; |
---|
| 287 | ptcomparison.pt = &pt; |
---|
| 288 | |
---|
| 289 | // set up the list of particles left. |
---|
| 290 | init_pleft(); |
---|
| 291 | |
---|
| 292 | indices = new int[n]; |
---|
| 293 | |
---|
| 294 | return 0; |
---|
| 295 | } |
---|
| 296 | |
---|
| 297 | |
---|
| 298 | // build initial list of remaining particles |
---|
| 299 | //------------------------------------------ |
---|
| 300 | int Csplit_merge::init_pleft(){ |
---|
| 301 | // at this level, we only rule out particles with |
---|
| 302 | // infinite rapidity |
---|
| 303 | // in the parent particle list, index contain the run |
---|
| 304 | // at which particles are puts in jets: |
---|
| 305 | // - -1 means infinity rapidity |
---|
| 306 | // - 0 means not included |
---|
| 307 | // - i mean included at run i |
---|
| 308 | int i,j; |
---|
| 309 | |
---|
| 310 | // copy particles removing the ones with infinite rapidity |
---|
| 311 | // determine min,max eta |
---|
| 312 | j=0; |
---|
| 313 | double eta_min=0.0; /// for the Ceta_phi_range static member! |
---|
| 314 | double eta_max=0.0; /// for the Ceta_phi_range static member! |
---|
| 315 | p_remain.clear(); |
---|
| 316 | for (i=0;i<n;i++){ |
---|
| 317 | // set ref for checkxor |
---|
| 318 | particles[i].ref.randomize(); |
---|
| 319 | |
---|
| 320 | // check if rapidity is not infinite or ill-defined |
---|
| 321 | if (fabs(particles[i].pz) < (particles[i].E)){ |
---|
| 322 | p_remain.push_back(particles[i]); |
---|
| 323 | // set up parent index for tracability |
---|
| 324 | p_remain[j].parent_index = i; |
---|
| 325 | // left particles are marked with a 1 |
---|
| 326 | // IMPORTANT NOTE: the meaning of index in p_remain is |
---|
| 327 | // somehow different as in the initial particle list. |
---|
| 328 | // here, within one pass, we use the index to track whether |
---|
| 329 | // a particle is included in the current pass (index set to 0 |
---|
| 330 | // in add_protocones) or still remain (index still 1) |
---|
| 331 | p_remain[j].index = 1; |
---|
| 332 | |
---|
| 333 | j++; |
---|
| 334 | // set up parent-particle index |
---|
| 335 | particles[i].index = 0; |
---|
| 336 | |
---|
| 337 | eta_min = min(eta_min, particles[i].eta); |
---|
| 338 | eta_max = max(eta_max, particles[i].eta); |
---|
| 339 | } else { |
---|
| 340 | particles[i].index = -1; |
---|
| 341 | } |
---|
| 342 | } |
---|
| 343 | n_left = p_remain.size(); |
---|
| 344 | n_pass = 0; |
---|
| 345 | |
---|
| 346 | Ceta_phi_range epr; |
---|
| 347 | epr.eta_min = eta_min-0.01; |
---|
| 348 | epr.eta_max = eta_max+0.01; |
---|
| 349 | |
---|
| 350 | merge_collinear_and_remove_soft(); |
---|
| 351 | |
---|
| 352 | return 0; |
---|
| 353 | } |
---|
| 354 | |
---|
| 355 | |
---|
| 356 | // partial clearance |
---|
| 357 | // we want to keep particle list and indices |
---|
| 358 | // for future usage, so do not clear it ! |
---|
| 359 | // this is done in full_clear |
---|
| 360 | //---------------------------------------- |
---|
| 361 | int Csplit_merge::partial_clear(){ |
---|
| 362 | // release jets |
---|
| 363 | |
---|
| 364 | // set up the auto_ptr for the multiset with the _current_ state of |
---|
| 365 | // ptcomparison (which may have changed since we constructed the |
---|
| 366 | // class) |
---|
| 367 | candidates.reset(new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison)); |
---|
| 368 | |
---|
| 369 | // start off with huge number |
---|
| 370 | most_ambiguous_split = numeric_limits<double>::max(); |
---|
| 371 | |
---|
| 372 | jets.clear(); |
---|
| 373 | #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES |
---|
| 374 | if (merge_identical_protocones) |
---|
| 375 | cand_refs.clear(); |
---|
| 376 | #endif |
---|
| 377 | |
---|
| 378 | p_remain.clear(); |
---|
| 379 | |
---|
| 380 | return 0; |
---|
| 381 | } |
---|
| 382 | |
---|
| 383 | |
---|
| 384 | // full clearance |
---|
| 385 | //---------------- |
---|
| 386 | int Csplit_merge::full_clear(){ |
---|
| 387 | partial_clear(); |
---|
| 388 | |
---|
| 389 | // clear previously allocated memory |
---|
| 390 | if (indices != NULL){ |
---|
| 391 | delete[] indices; |
---|
| 392 | } |
---|
| 393 | particles.clear(); |
---|
| 394 | |
---|
| 395 | return 0; |
---|
| 396 | } |
---|
| 397 | |
---|
| 398 | |
---|
| 399 | // build the list 'p_uncol_hard' from p_remain by clustering collinear particles |
---|
| 400 | // note that thins in only used for stable-cone detection |
---|
| 401 | // so the parent_index field is unnecessary |
---|
| 402 | //------------------------------------------------------------------------- |
---|
| 403 | int Csplit_merge::merge_collinear_and_remove_soft(){ |
---|
| 404 | int i,j; |
---|
| 405 | vector<Cmomentum> p_sorted; |
---|
| 406 | bool collinear; |
---|
| 407 | double dphi; |
---|
| 408 | |
---|
| 409 | p_uncol_hard.clear(); |
---|
| 410 | |
---|
| 411 | // we first sort the particles according to their rapidity |
---|
| 412 | for (i=0;i<n_left;i++) |
---|
| 413 | p_sorted.push_back(p_remain[i]); |
---|
| 414 | sort(p_sorted.begin(), p_sorted.end(), momentum_eta_less); |
---|
| 415 | |
---|
| 416 | // then we cluster particles looping over the particles in the following way |
---|
| 417 | // if (a particle i has same eta-phi a one after (j)) |
---|
| 418 | // then add momentum i to j |
---|
| 419 | // else add i to the p_uncol_hard list |
---|
| 420 | i = 0; |
---|
| 421 | while (i<n_left){ |
---|
| 422 | // check if the particle passes the stable_cone_soft_pt2_cutoff |
---|
| 423 | if (p_sorted[i].perp2()<stable_cone_soft_pt2_cutoff) { |
---|
| 424 | i++; |
---|
| 425 | continue; |
---|
| 426 | } |
---|
| 427 | |
---|
| 428 | // check if there is(are) particle(s) with the 'same' eta |
---|
| 429 | collinear = false; |
---|
| 430 | j=i+1; |
---|
| 431 | while ((j<n_left) && (fabs(p_sorted[j].eta-p_sorted[i].eta)<EPSILON_COLLINEAR) && (!collinear)){ |
---|
| 432 | dphi = fabs(p_sorted[j].phi-p_sorted[i].phi); |
---|
| 433 | if (dphi>M_PI) dphi = twopi-dphi; |
---|
| 434 | if (dphi<EPSILON_COLLINEAR){ |
---|
| 435 | // i is collinear with j; add the momentum (i) to the momentum (j) |
---|
| 436 | p_sorted[j] += p_sorted[i]; |
---|
| 437 | // set collinearity test to true |
---|
| 438 | collinear = true; |
---|
| 439 | } |
---|
| 440 | j++; |
---|
| 441 | } |
---|
| 442 | // if no collinearity detected, add the particle to our list |
---|
| 443 | if (!collinear) |
---|
| 444 | p_uncol_hard.push_back(p_sorted[i]); |
---|
| 445 | i++; |
---|
| 446 | } |
---|
| 447 | |
---|
| 448 | return 0; |
---|
| 449 | } |
---|
| 450 | |
---|
| 451 | |
---|
| 452 | // add a list of protocones |
---|
| 453 | // - protocones list of protocones (initial jet candidates) |
---|
| 454 | // - R2 cone radius (squared) |
---|
| 455 | // - ptmin minimal pT allowed for jets |
---|
| 456 | //------------------------------------------------------------- |
---|
| 457 | int Csplit_merge::add_protocones(vector<Cmomentum> *protocones, double R2, double ptmin){ |
---|
| 458 | int i; |
---|
| 459 | Cmomentum *c; |
---|
| 460 | Cmomentum *v; |
---|
| 461 | double eta, phi; |
---|
| 462 | double dx, dy; |
---|
| 463 | double R; |
---|
| 464 | Cjet jet; |
---|
| 465 | |
---|
| 466 | if (protocones->size()==0) |
---|
| 467 | return 1; |
---|
| 468 | |
---|
| 469 | pt_min2 = ptmin*ptmin; |
---|
| 470 | R = sqrt(R2); |
---|
| 471 | |
---|
| 472 | // browse protocones |
---|
| 473 | // for each of them, build the list of particles in them |
---|
| 474 | for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){ |
---|
| 475 | // initialise variables |
---|
| 476 | c = &(*p_it); |
---|
| 477 | |
---|
| 478 | // note: cones have been tested => their (eta,phi) coordinates are computed |
---|
| 479 | eta = c->eta; |
---|
| 480 | phi = c->phi; |
---|
| 481 | |
---|
| 482 | // browse particles to create cone contents |
---|
| 483 | // note that jet is always initialised with default values at this level |
---|
| 484 | jet.v = Cmomentum(); |
---|
| 485 | jet.pt_tilde=0; |
---|
| 486 | jet.contents.clear(); |
---|
| 487 | for (i=0;i<n_left;i++){ |
---|
| 488 | v = &(p_remain[i]); |
---|
| 489 | // for security, avoid including particles with infinite rapidity) |
---|
| 490 | // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING |
---|
| 491 | //if (fabs(v->pz)!=v->E){ |
---|
| 492 | dx = eta - v->eta; |
---|
| 493 | dy = fabs(phi - v->phi); |
---|
| 494 | if (dy>M_PI) |
---|
| 495 | dy -= twopi; |
---|
| 496 | if (dx*dx+dy*dy<R2){ |
---|
| 497 | jet.contents.push_back(v->parent_index); |
---|
| 498 | jet.v+= *v; |
---|
| 499 | jet.pt_tilde+= pt[v->parent_index]; |
---|
| 500 | v->index=0; |
---|
| 501 | } |
---|
| 502 | } |
---|
| 503 | jet.n=jet.contents.size(); |
---|
| 504 | |
---|
| 505 | // set the momentum in protocones |
---|
| 506 | // (it was only known through eta and phi up to now) |
---|
| 507 | *c = jet.v; |
---|
| 508 | c->eta = eta; // restore exact original coords |
---|
| 509 | c->phi = phi; // to avoid rounding error inconsistencies |
---|
| 510 | |
---|
| 511 | // set the jet range |
---|
| 512 | jet.range=Ceta_phi_range(eta,phi,R); |
---|
| 513 | |
---|
| 514 | #ifdef DEBUG_SPLIT_MERGE |
---|
| 515 | cout << "adding jet: "; |
---|
| 516 | for (int i2=0;i2<jet.n;i2++) |
---|
| 517 | cout << jet.contents[i2] << " "; |
---|
| 518 | cout << endl; |
---|
| 519 | #endif |
---|
| 520 | |
---|
| 521 | // add it to the list of jets |
---|
| 522 | insert(jet); |
---|
| 523 | } |
---|
| 524 | |
---|
| 525 | // update list of included particles |
---|
| 526 | n_pass++; |
---|
| 527 | |
---|
| 528 | #ifdef DEBUG_SPLIT_MERGE |
---|
| 529 | cout << "remaining particles: "; |
---|
| 530 | #endif |
---|
| 531 | int j=0; |
---|
| 532 | for (i=0;i<n_left;i++){ |
---|
| 533 | if (p_remain[i].index){ |
---|
| 534 | // copy particle |
---|
| 535 | p_remain[j]=p_remain[i]; |
---|
| 536 | p_remain[j].parent_index = p_remain[i].parent_index; |
---|
| 537 | p_remain[j].index=1; |
---|
| 538 | // set run in initial list |
---|
| 539 | particles[p_remain[j].parent_index].index = n_pass; |
---|
| 540 | #ifdef DEBUG_SPLIT_MERGE |
---|
| 541 | cout << p_remain[j].parent_index << " "; |
---|
| 542 | #endif |
---|
| 543 | j++; |
---|
| 544 | } |
---|
| 545 | } |
---|
| 546 | #ifdef DEBUG_SPLIT_MERGE |
---|
| 547 | cout << endl; |
---|
| 548 | #endif |
---|
| 549 | n_left = j; |
---|
| 550 | p_remain.resize(j); |
---|
| 551 | |
---|
| 552 | merge_collinear_and_remove_soft(); |
---|
| 553 | |
---|
| 554 | return 0; |
---|
| 555 | } |
---|
| 556 | |
---|
| 557 | |
---|
| 558 | /* |
---|
| 559 | * really do the splitting and merging |
---|
| 560 | * At the end, the vector jets is filled with the jets found. |
---|
| 561 | * the 'contents' field of each jets contains the indices |
---|
| 562 | * of the particles included in that jet. |
---|
| 563 | * - overlap_tshold threshold for splitting/merging transition |
---|
| 564 | * - ptmin minimal pT allowed for jets |
---|
| 565 | * return the number of jets is returned |
---|
| 566 | ******************************************************************/ |
---|
| 567 | int Csplit_merge::perform(double overlap_tshold, double ptmin){ |
---|
| 568 | // iterators for the 2 jets |
---|
| 569 | cjet_iterator j1; |
---|
| 570 | cjet_iterator j2; |
---|
| 571 | |
---|
| 572 | pt_min2 = ptmin*ptmin; |
---|
| 573 | |
---|
| 574 | if (candidates->size()==0) |
---|
| 575 | return 0; |
---|
| 576 | |
---|
| 577 | if (overlap_tshold>=1.0 || overlap_tshold <= 0) { |
---|
| 578 | ostringstream message; |
---|
| 579 | message << "Illegal value for overlap_tshold, f = " << overlap_tshold; |
---|
| 580 | message << " (legal values are 0<f<1)"; |
---|
| 581 | throw Csiscone_error(message.str()); |
---|
| 582 | } |
---|
| 583 | |
---|
| 584 | // overlap (the contents of this variable depends on the choice for |
---|
| 585 | // the split--merge variable.) |
---|
| 586 | // Note that the square of the ovelap is used |
---|
| 587 | double overlap2; |
---|
| 588 | |
---|
| 589 | // avoid to compute tshold*tshold at each overlap |
---|
| 590 | double overlap_tshold2 = overlap_tshold*overlap_tshold; |
---|
| 591 | |
---|
| 592 | do{ |
---|
| 593 | if (candidates->size()>0){ |
---|
| 594 | // browse for the first jet |
---|
| 595 | j1 = candidates->begin(); |
---|
| 596 | |
---|
| 597 | // if hardest jet does not pass threshold then nothing else will |
---|
| 598 | // either so one stops the split merge. |
---|
| 599 | if (j1->sm_var2<SM_var2_hardest_cut_off) {break;} |
---|
| 600 | |
---|
| 601 | // browse for the second jet |
---|
| 602 | j2 = j1; |
---|
| 603 | j2++; |
---|
| 604 | int j2_relindex = 1; // used only in ifdef, but costs little so keep it outside |
---|
| 605 | |
---|
| 606 | while (j2 != candidates->end()){ |
---|
| 607 | #ifdef DEBUG_SPLIT_MERGE |
---|
| 608 | show(); |
---|
| 609 | #endif |
---|
| 610 | // check overlapping |
---|
| 611 | if (get_overlap(*j1, *j2, &overlap2)){ |
---|
| 612 | // check if overlapping energy passes threshold |
---|
| 613 | // Note that this depends on the ordering variable |
---|
| 614 | #ifdef DEBUG_SPLIT_MERGE |
---|
| 615 | cout << "overlap between cdt 1 and cdt " << j2_relindex+1 << " with overlap " |
---|
| 616 | << sqrt(overlap2/j2->sm_var2) << endl<<endl; |
---|
| 617 | #endif |
---|
| 618 | if (overlap2<overlap_tshold2*j2->sm_var2){ |
---|
| 619 | // split jets |
---|
| 620 | split(j1, j2); |
---|
| 621 | |
---|
| 622 | // update iterators |
---|
| 623 | j2 = j1 = candidates->begin(); |
---|
| 624 | j2_relindex = 0; |
---|
| 625 | } else { |
---|
| 626 | // merge jets |
---|
| 627 | merge(j1, j2); |
---|
| 628 | |
---|
| 629 | // update iterators |
---|
| 630 | j2 = j1 = candidates->begin(); |
---|
| 631 | j2_relindex = 0; |
---|
| 632 | } |
---|
| 633 | } |
---|
| 634 | // watch out: split/merge might have caused new jets with pt < |
---|
| 635 | // ptmin to disappear, so the total number of jets may |
---|
| 636 | // have changed by more than expected and j2 might already by |
---|
| 637 | // the end of the candidates list... |
---|
| 638 | j2_relindex++; |
---|
| 639 | if (j2 != candidates->end()) j2++; |
---|
| 640 | } // end of loop on the second jet |
---|
| 641 | |
---|
| 642 | if (j1 != candidates->end()) { |
---|
| 643 | // all "second jet" passed without overlapping |
---|
| 644 | // (otherwise we won't leave the j2 loop) |
---|
| 645 | // => store jet 1 as real jet |
---|
| 646 | jets.push_back(*j1); |
---|
| 647 | jets[jets.size()-1].v.build_etaphi(); |
---|
| 648 | // a bug where the contents has non-zero size has been cropping |
---|
| 649 | // up in many contexts -- so catch it! |
---|
| 650 | assert(j1->contents.size() > 0); |
---|
| 651 | jets[jets.size()-1].pass = particles[j1->contents[0]].index; |
---|
| 652 | #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES |
---|
| 653 | cand_refs.erase(j1->v.ref); |
---|
| 654 | #endif |
---|
| 655 | candidates->erase(j1); |
---|
| 656 | |
---|
| 657 | //// test that the hardest jet pass the potential cut-off |
---|
| 658 | //if ((candidates->size()!=0) && |
---|
| 659 | // (candidates->begin()->sm_var2<SM_var2_hardest_cut_off)){ |
---|
| 660 | // candidates->clear(); |
---|
| 661 | //} |
---|
| 662 | } |
---|
| 663 | } |
---|
| 664 | } while (candidates->size()>0); |
---|
| 665 | |
---|
| 666 | // sort jets by pT |
---|
| 667 | sort(jets.begin(), jets.end(), jets_pt_less); |
---|
| 668 | #ifdef DEBUG_SPLIT_MERGE |
---|
| 669 | show(); |
---|
| 670 | #endif |
---|
| 671 | |
---|
| 672 | return jets.size(); |
---|
| 673 | } |
---|
| 674 | |
---|
| 675 | |
---|
| 676 | |
---|
| 677 | // save the event on disk |
---|
| 678 | // - flux stream used to save jet contents |
---|
| 679 | //-------------------------------------------- |
---|
| 680 | int Csplit_merge::save_contents(FILE *flux){ |
---|
| 681 | jet_iterator it_j; |
---|
| 682 | Cjet *j1; |
---|
| 683 | int i1, i2; |
---|
| 684 | |
---|
| 685 | fprintf(flux, "# %d jets found\n", (int) jets.size()); |
---|
| 686 | fprintf(flux, "# columns are: eta, phi, pt and number of particles for each jet\n"); |
---|
| 687 | for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){ |
---|
| 688 | j1 = &(*it_j); |
---|
| 689 | j1->v.build_etaphi(); |
---|
| 690 | fprintf(flux, "%f\t%f\t%e\t%d\n", |
---|
| 691 | j1->v.eta, j1->v.phi, j1->v.perp(), j1->n); |
---|
| 692 | } |
---|
| 693 | |
---|
| 694 | fprintf(flux, "# jet contents\n"); |
---|
| 695 | fprintf(flux, "# columns are: eta, phi, pt, particle index and jet number\n"); |
---|
| 696 | for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){ |
---|
| 697 | j1 = &(*it_j); |
---|
| 698 | for (i2=0;i2<j1->n;i2++) |
---|
| 699 | fprintf(flux, "%f\t%f\t%e\t%d\t%d\n", |
---|
| 700 | particles[j1->contents[i2]].eta, particles[j1->contents[i2]].phi, |
---|
| 701 | particles[j1->contents[i2]].perp(), j1->contents[i2], i1); |
---|
| 702 | } |
---|
| 703 | |
---|
| 704 | return 0; |
---|
| 705 | } |
---|
| 706 | |
---|
| 707 | |
---|
| 708 | // show current jets/candidate status |
---|
| 709 | //------------------------------------ |
---|
| 710 | int Csplit_merge::show(){ |
---|
| 711 | jet_iterator it_j; |
---|
| 712 | cjet_iterator it_c; |
---|
| 713 | Cjet *j; |
---|
| 714 | const Cjet *c; |
---|
| 715 | int i1, i2; |
---|
| 716 | |
---|
| 717 | for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){ |
---|
| 718 | j = &(*it_j); |
---|
| 719 | fprintf(stdout, "jet %2d: %e\t%e\t%e\t%e\t", i1+1, |
---|
| 720 | j->v.px, j->v.py, j->v.pz, j->v.E); |
---|
| 721 | for (i2=0;i2<j->n;i2++) |
---|
| 722 | fprintf(stdout, "%d ", j->contents[i2]); |
---|
| 723 | fprintf(stdout, "\n"); |
---|
| 724 | } |
---|
| 725 | |
---|
| 726 | for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){ |
---|
| 727 | c = &(*it_c); |
---|
| 728 | fprintf(stdout, "cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1, |
---|
| 729 | c->v.px, c->v.py, c->v.pz, c->v.E, sqrt(c->sm_var2)); |
---|
| 730 | for (i2=0;i2<c->n;i2++) |
---|
| 731 | fprintf(stdout, "%d ", c->contents[i2]); |
---|
| 732 | fprintf(stdout, "\n"); |
---|
| 733 | } |
---|
| 734 | |
---|
| 735 | fprintf(stdout, "\n"); |
---|
| 736 | return 0; |
---|
| 737 | } |
---|
| 738 | |
---|
| 739 | |
---|
| 740 | // get the overlap between 2 jets |
---|
| 741 | // - j1 first jet |
---|
| 742 | // - j2 second jet |
---|
| 743 | // - overlap2 returned overlap^2 (determined by the choice of SM variable) |
---|
| 744 | // return true if overlapping, false if disjoint |
---|
| 745 | //--------------------------------------------------------------------- |
---|
| 746 | bool Csplit_merge::get_overlap(const Cjet &j1, const Cjet &j2, double *overlap2){ |
---|
| 747 | // check if ranges overlap |
---|
| 748 | if (!is_range_overlap(j1.range,j2.range)) |
---|
| 749 | return false; |
---|
| 750 | |
---|
| 751 | int i1,i2; |
---|
| 752 | bool is_overlap; |
---|
| 753 | |
---|
| 754 | // initialise |
---|
| 755 | i1=i2=idx_size=0; |
---|
| 756 | is_overlap = false; |
---|
| 757 | Cmomentum v; |
---|
| 758 | double pt_tilde=0.0; |
---|
| 759 | |
---|
| 760 | // compute overlap |
---|
| 761 | // at the same time, we store union in indices |
---|
| 762 | do{ |
---|
| 763 | if (j1.contents[i1]<j2.contents[i2]){ |
---|
| 764 | indices[idx_size] = j1.contents[i1]; |
---|
| 765 | i1++; |
---|
| 766 | } else if (j1.contents[i1]>j2.contents[i2]){ |
---|
| 767 | indices[idx_size] = j2.contents[i2]; |
---|
| 768 | i2++; |
---|
| 769 | } else { // (j1.contents[i1]==j2.contents[i2]) |
---|
| 770 | v += particles[j1.contents[i1]]; |
---|
| 771 | pt_tilde += pt[j1.contents[i1]]; |
---|
| 772 | indices[idx_size] = j1.contents[i1]; |
---|
| 773 | i1++; |
---|
| 774 | i2++; |
---|
| 775 | is_overlap = true; |
---|
| 776 | } |
---|
| 777 | idx_size++; |
---|
| 778 | } while ((i1<j1.n) && (i2<j2.n)); |
---|
| 779 | |
---|
| 780 | // finish computing union |
---|
| 781 | // (only needed if overlap !) |
---|
| 782 | if (is_overlap){ |
---|
| 783 | while (i1<j1.n){ |
---|
| 784 | indices[idx_size] = j1.contents[i1]; |
---|
| 785 | i1++; |
---|
| 786 | idx_size++; |
---|
| 787 | } |
---|
| 788 | while (i2<j2.n){ |
---|
| 789 | indices[idx_size] = j2.contents[i2]; |
---|
| 790 | i2++; |
---|
| 791 | idx_size++; |
---|
| 792 | } |
---|
| 793 | } |
---|
| 794 | |
---|
| 795 | // assign the overlapping var as return variable |
---|
| 796 | (*overlap2) = get_sm_var2(v, pt_tilde); |
---|
| 797 | |
---|
| 798 | return is_overlap; |
---|
| 799 | } |
---|
| 800 | |
---|
| 801 | |
---|
| 802 | |
---|
| 803 | // split the two given jet. |
---|
| 804 | // during this procedure, the jets j1 & j2 are replaced |
---|
| 805 | // by 2 new jets. Common particles are associted to the |
---|
| 806 | // closest initial jet. |
---|
| 807 | // - it_j1 iterator of the first jet in 'candidates' |
---|
| 808 | // - it_j2 iterator of the second jet in 'candidates' |
---|
| 809 | // - j1 first jet (Cjet instance) |
---|
| 810 | // - j2 second jet (Cjet instance) |
---|
| 811 | // return true on success, false on error |
---|
| 812 | //////////////////////////////////////////////////////////////// |
---|
| 813 | bool Csplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){ |
---|
| 814 | int i1, i2; |
---|
| 815 | Cjet jet1, jet2; |
---|
| 816 | double eta1, phi1, pt1_weight, eta2, phi2, pt2_weight; |
---|
| 817 | double dx1, dy1, dx2, dy2; |
---|
| 818 | Cmomentum tmp; |
---|
| 819 | Cmomentum *v; |
---|
| 820 | |
---|
| 821 | // shorthand to avoid having "->" all over the place |
---|
| 822 | const Cjet & j1 = * it_j1; |
---|
| 823 | const Cjet & j2 = * it_j2; |
---|
| 824 | |
---|
| 825 | i1=i2=0; |
---|
| 826 | jet2.v = jet1.v = Cmomentum(); |
---|
| 827 | jet2.pt_tilde = jet1.pt_tilde = 0.0; |
---|
| 828 | |
---|
| 829 | // compute centroids |
---|
| 830 | // When use_pt_weighted_splitting is activated, the |
---|
| 831 | // "geometrical" distance is weighted by the inverse |
---|
| 832 | // of the pt of the protojet |
---|
| 833 | // This is stored in pt{1,2}_weight |
---|
| 834 | tmp = j1.v; |
---|
| 835 | tmp.build_etaphi(); |
---|
| 836 | eta1 = tmp.eta; |
---|
| 837 | phi1 = tmp.phi; |
---|
| 838 | pt1_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0; |
---|
| 839 | |
---|
| 840 | tmp = j2.v; |
---|
| 841 | tmp.build_etaphi(); |
---|
| 842 | eta2 = tmp.eta; |
---|
| 843 | phi2 = tmp.phi; |
---|
| 844 | pt2_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0; |
---|
| 845 | |
---|
| 846 | jet1.v = jet2.v = Cmomentum(); |
---|
| 847 | |
---|
| 848 | // compute jet splitting |
---|
| 849 | do{ |
---|
| 850 | if (j1.contents[i1]<j2.contents[i2]){ |
---|
| 851 | // particle i1 belong only to jet 1 |
---|
| 852 | v = &(particles[j1.contents[i1]]); |
---|
| 853 | jet1.contents.push_back(j1.contents[i1]); |
---|
| 854 | jet1.v += *v; |
---|
| 855 | jet1.pt_tilde += pt[j1.contents[i1]]; |
---|
| 856 | i1++; |
---|
| 857 | jet1.range.add_particle(v->eta,v->phi); |
---|
| 858 | } else if (j1.contents[i1]>j2.contents[i2]){ |
---|
| 859 | // particle i2 belong only to jet 2 |
---|
| 860 | v = &(particles[j2.contents[i2]]); |
---|
| 861 | jet2.contents.push_back(j2.contents[i2]); |
---|
| 862 | jet2.v += *v; |
---|
| 863 | jet2.pt_tilde += pt[j2.contents[i2]]; |
---|
| 864 | i2++; |
---|
| 865 | jet2.range.add_particle(v->eta,v->phi); |
---|
| 866 | } else { // (j1.contents[i1]==j2.contents[i2]) |
---|
| 867 | // common particle, decide which is the closest centre |
---|
| 868 | v = &(particles[j1.contents[i1]]); |
---|
| 869 | |
---|
| 870 | // distance w.r.t. centroid 1 |
---|
| 871 | dx1 = eta1 - v->eta; |
---|
| 872 | dy1 = fabs(phi1 - v->phi); |
---|
| 873 | if (dy1>M_PI) |
---|
| 874 | dy1 -= twopi; |
---|
| 875 | |
---|
| 876 | // distance w.r.t. centroid 2 |
---|
| 877 | dx2 = eta2 - v->eta; |
---|
| 878 | dy2 = fabs(phi2 - v->phi); |
---|
| 879 | if (dy2>M_PI) |
---|
| 880 | dy2 -= twopi; |
---|
| 881 | |
---|
| 882 | //? what when == ? |
---|
| 883 | // When use_pt_weighted_splitting is activated, the |
---|
| 884 | // "geometrical" distance is weighted by the inverse |
---|
| 885 | // of the pt of the protojet |
---|
| 886 | double d1sq = (dx1*dx1+dy1*dy1)*pt1_weight; |
---|
| 887 | double d2sq = (dx2*dx2+dy2*dy2)*pt2_weight; |
---|
| 888 | // do bookkeeping on most ambiguous split |
---|
| 889 | if (fabs(d1sq-d2sq) < most_ambiguous_split) |
---|
| 890 | most_ambiguous_split = fabs(d1sq-d2sq); |
---|
| 891 | |
---|
| 892 | if (d1sq<d2sq){ |
---|
| 893 | // particle i1 belong only to jet 1 |
---|
| 894 | jet1.contents.push_back(j1.contents[i1]); |
---|
| 895 | jet1.v += *v; |
---|
| 896 | jet1.pt_tilde += pt[j1.contents[i1]]; |
---|
| 897 | jet1.range.add_particle(v->eta,v->phi); |
---|
| 898 | } else { |
---|
| 899 | // particle i2 belong only to jet 2 |
---|
| 900 | jet2.contents.push_back(j2.contents[i2]); |
---|
| 901 | jet2.v += *v; |
---|
| 902 | jet2.pt_tilde += pt[j2.contents[i2]]; |
---|
| 903 | jet2.range.add_particle(v->eta,v->phi); |
---|
| 904 | } |
---|
| 905 | |
---|
| 906 | i1++; |
---|
| 907 | i2++; |
---|
| 908 | } |
---|
| 909 | } while ((i1<j1.n) && (i2<j2.n)); |
---|
| 910 | |
---|
| 911 | while (i1<j1.n){ |
---|
| 912 | v = &(particles[j1.contents[i1]]); |
---|
| 913 | jet1.contents.push_back(j1.contents[i1]); |
---|
| 914 | jet1.v += *v; |
---|
| 915 | jet1.pt_tilde += pt[j1.contents[i1]]; |
---|
| 916 | i1++; |
---|
| 917 | jet1.range.add_particle(v->eta,v->phi); |
---|
| 918 | } |
---|
| 919 | while (i2<j2.n){ |
---|
| 920 | v = &(particles[j2.contents[i2]]); |
---|
| 921 | jet2.contents.push_back(j2.contents[i2]); |
---|
| 922 | jet2.v += *v; |
---|
| 923 | jet2.pt_tilde += pt[j2.contents[i2]]; |
---|
| 924 | i2++; |
---|
| 925 | jet2.range.add_particle(v->eta,v->phi); |
---|
| 926 | } |
---|
| 927 | |
---|
| 928 | // finalise jets |
---|
| 929 | jet1.n = jet1.contents.size(); |
---|
| 930 | jet2.n = jet2.contents.size(); |
---|
| 931 | |
---|
| 932 | //jet1.range = j1.range; |
---|
| 933 | //jet2.range = j2.range; |
---|
| 934 | |
---|
| 935 | // remove previous jets |
---|
| 936 | #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES |
---|
| 937 | cand_refs.erase(j1.v.ref); |
---|
| 938 | cand_refs.erase(j2.v.ref); |
---|
| 939 | #endif |
---|
| 940 | candidates->erase(it_j1); |
---|
| 941 | candidates->erase(it_j2); |
---|
| 942 | |
---|
| 943 | // reinsert new ones |
---|
| 944 | insert(jet1); |
---|
| 945 | insert(jet2); |
---|
| 946 | |
---|
| 947 | return true; |
---|
| 948 | } |
---|
| 949 | |
---|
| 950 | // merge the two given jet. |
---|
| 951 | // during this procedure, the jets j1 & j2 are replaced |
---|
| 952 | // by 1 single jets containing both of them. |
---|
| 953 | // - it_j1 iterator of the first jet in 'candidates' |
---|
| 954 | // - it_j2 iterator of the second jet in 'candidates' |
---|
| 955 | // return true on success, false on error |
---|
| 956 | //////////////////////////////////////////////////////////////// |
---|
| 957 | bool Csplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){ |
---|
| 958 | Cjet jet; |
---|
| 959 | int i; |
---|
| 960 | |
---|
| 961 | // build new jet |
---|
| 962 | // note: particles within j1 & j2 have already been stored in indices |
---|
| 963 | for (i=0;i<idx_size;i++){ |
---|
| 964 | jet.contents.push_back(indices[i]); |
---|
| 965 | jet.v += particles[indices[i]]; |
---|
| 966 | jet.pt_tilde += pt[indices[i]]; |
---|
| 967 | } |
---|
| 968 | jet.n = jet.contents.size(); |
---|
| 969 | |
---|
| 970 | // deal with ranges |
---|
| 971 | jet.range = range_union(it_j1->range, it_j2->range); |
---|
| 972 | |
---|
| 973 | // remove old candidates |
---|
| 974 | #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES |
---|
| 975 | if (merge_identical_protocones){ |
---|
| 976 | cand_refs.erase(it_j1->v.ref); |
---|
| 977 | cand_refs.erase(it_j2->v.ref); |
---|
| 978 | } |
---|
| 979 | #endif |
---|
| 980 | candidates->erase(it_j1); |
---|
| 981 | candidates->erase(it_j2); |
---|
| 982 | |
---|
| 983 | // reinsert new candidate |
---|
| 984 | insert(jet); |
---|
| 985 | |
---|
| 986 | return true; |
---|
| 987 | } |
---|
| 988 | |
---|
| 989 | /** |
---|
| 990 | * Check whether or not a jet has to be inserted in the |
---|
| 991 | * list of protojets. If it has, set its sm_variable and |
---|
| 992 | * insert it to the list of protojets. |
---|
| 993 | */ |
---|
| 994 | bool Csplit_merge::insert(Cjet &jet){ |
---|
| 995 | |
---|
| 996 | // eventually check that no other candidate are present with the |
---|
| 997 | // same cone contents. We recall that this automatic merging of |
---|
| 998 | // identical protocones can lead to infrared-unsafe situations. |
---|
| 999 | #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES |
---|
| 1000 | if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second)) |
---|
| 1001 | return false; |
---|
| 1002 | #endif |
---|
| 1003 | |
---|
| 1004 | // check that the protojet has large enough pt |
---|
| 1005 | if (jet.v.perp2()<pt_min2) |
---|
| 1006 | return false; |
---|
| 1007 | |
---|
| 1008 | // assign SM variable |
---|
| 1009 | jet.sm_var2 = get_sm_var2(jet.v, jet.pt_tilde); |
---|
| 1010 | |
---|
| 1011 | // insert the jet. |
---|
| 1012 | candidates->insert(jet); |
---|
| 1013 | |
---|
| 1014 | return true; |
---|
| 1015 | } |
---|
| 1016 | |
---|
| 1017 | /** |
---|
| 1018 | * given a 4-momentum and its associated pT, return the |
---|
| 1019 | * variable that has to be used for SM |
---|
| 1020 | * \param v 4 momentum of the protojet |
---|
| 1021 | * \param pt_tilde pt_tilde of the protojet |
---|
| 1022 | */ |
---|
| 1023 | double Csplit_merge::get_sm_var2(Cmomentum &v, double &pt_tilde){ |
---|
| 1024 | switch(ptcomparison.split_merge_scale) { |
---|
| 1025 | case SM_pt: return v.perp2(); |
---|
| 1026 | case SM_mt: return v.perpmass2(); |
---|
| 1027 | case SM_pttilde: return pt_tilde*pt_tilde; |
---|
| 1028 | case SM_Et: return v.Et2(); |
---|
| 1029 | default: |
---|
| 1030 | throw Csiscone_error("Unsupported split-merge scale choice: " |
---|
| 1031 | + ptcomparison.SM_scale_name()); |
---|
| 1032 | } |
---|
| 1033 | |
---|
| 1034 | return 0.0; |
---|
| 1035 | } |
---|
| 1036 | |
---|
| 1037 | } |
---|