Fork me on GitHub

source: git/external/fastjet/plugins/SISCone/split_merge.cc@ 8abab33

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 8abab33 was 273e668, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

upgrade FastJet to version 3.1.0

  • Property mode set to 100644
File size: 34.8 KB
RevLine 
[d7d2da3]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// //
[273e668]23// $Revision:: 370 $//
24// $Date:: 2014-09-04 17:03:15 +0200 (Thu, 04 Sep 2014) $//
[d7d2da3]25///////////////////////////////////////////////////////////////////////////////
26
27#include "split_merge.h"
28#include "siscone_error.h"
29#include "momentum.h"
30#include <limits> // for max
31#include <iostream>
32#include <algorithm>
33#include <sstream>
34#include <cassert>
[273e668]35#include <cmath>
[d7d2da3]36
37namespace siscone{
38
39using 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//--------------
51Cjet::Cjet(){
52 n = 0;
53 v = Cmomentum();
54 pt_tilde = 0.0;
55 sm_var2 = 0.0;
56}
57
58// default dtor
59//--------------
60Cjet::~Cjet(){
61
62}
63
64// ordering of jets in pt (e.g. used in final jets ordering)
65//-----------------------------------------------------------
66bool 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//------------------------------------------------------------------
92bool 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
154std::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//-----------------------------------------------
177void 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//--------------
224Csplit_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
[273e668]231 _user_scale = NULL;
[d7d2da3]232 indices = NULL;
233
234 // ensure that ptcomparison points to our set of particles (though params not correct)
235 ptcomparison.particles = &particles;
236 ptcomparison.pt = &pt;
237 candidates.reset(new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison));
238
239 // no hardest cut (col-unsafe)
[273e668]240 SM_var2_hardest_cut_off = -numeric_limits<double>::max();
[d7d2da3]241
242 // no pt cutoff for the particles to put in p_uncol_hard
243 stable_cone_soft_pt2_cutoff = -1.0;
244
245 // no pt-weighted splitting
246 use_pt_weighted_splitting = false;
247}
248
249
250// default dtor
251//--------------
252Csplit_merge::~Csplit_merge(){
253 full_clear();
254}
255
256
257// initialisation function
258// - _particles list of particles
259// - protocones list of protocones (initial jet candidates)
260// - R2 cone radius (squared)
261// - ptmin minimal pT allowed for jets
262//-------------------------------------------------------------
263int Csplit_merge::init(vector<Cmomentum> & /*_particles*/, vector<Cmomentum> *protocones, double R2, double ptmin){
264 // browse protocones
265 return add_protocones(protocones, R2, ptmin);
266}
267
268
269// initialisation function for particle list
270// - _particles list of particles
271//-------------------------------------------------------------
272int Csplit_merge::init_particles(vector<Cmomentum> &_particles){
273 full_clear();
274
275 // compute the list of particles
276 // here, we do not need to throw away particles
277 // with infinite rapidity (colinear with the beam)
278 particles = _particles;
279 n = particles.size();
280
281 // build the vector of particles' pt
282 pt.resize(n);
283 for (int i=0;i<n;i++)
284 pt[i] = particles[i].perp();
285
286 // ensure that ptcomparison points to our set of particles (though params not correct)
287 ptcomparison.particles = &particles;
288 ptcomparison.pt = &pt;
289
290 // set up the list of particles left.
291 init_pleft();
292
293 indices = new int[n];
294
295 return 0;
296}
297
298
299// build initial list of remaining particles
300//------------------------------------------
301int Csplit_merge::init_pleft(){
302 // at this level, we only rule out particles with
303 // infinite rapidity
304 // in the parent particle list, index contain the run
305 // at which particles are puts in jets:
306 // - -1 means infinity rapidity
307 // - 0 means not included
308 // - i mean included at run i
309 int i,j;
310
311 // copy particles removing the ones with infinite rapidity
312 // determine min,max eta
313 j=0;
314 double eta_min=0.0; /// for the Ceta_phi_range static member!
315 double eta_max=0.0; /// for the Ceta_phi_range static member!
316 p_remain.clear();
317 for (i=0;i<n;i++){
318 // set ref for checkxor
319 particles[i].ref.randomize();
320
321 // check if rapidity is not infinite or ill-defined
322 if (fabs(particles[i].pz) < (particles[i].E)){
323 p_remain.push_back(particles[i]);
324 // set up parent index for tracability
325 p_remain[j].parent_index = i;
326 // left particles are marked with a 1
327 // IMPORTANT NOTE: the meaning of index in p_remain is
328 // somehow different as in the initial particle list.
329 // here, within one pass, we use the index to track whether
330 // a particle is included in the current pass (index set to 0
331 // in add_protocones) or still remain (index still 1)
332 p_remain[j].index = 1;
333
334 j++;
335 // set up parent-particle index
336 particles[i].index = 0;
337
338 eta_min = min(eta_min, particles[i].eta);
339 eta_max = max(eta_max, particles[i].eta);
340 } else {
341 particles[i].index = -1;
342 }
343 }
344 n_left = p_remain.size();
345 n_pass = 0;
346
347 Ceta_phi_range epr;
348 epr.eta_min = eta_min-0.01;
349 epr.eta_max = eta_max+0.01;
350
351 merge_collinear_and_remove_soft();
352
353 return 0;
354}
355
356
357// partial clearance
358// we want to keep particle list and indices
359// for future usage, so do not clear it !
360// this is done in full_clear
361//----------------------------------------
362int Csplit_merge::partial_clear(){
363 // release jets
364
365 // set up the auto_ptr for the multiset with the _current_ state of
366 // ptcomparison (which may have changed since we constructed the
367 // class)
368 candidates.reset(new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison));
369
370 // start off with huge number
371 most_ambiguous_split = numeric_limits<double>::max();
372
373 jets.clear();
374#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
375 if (merge_identical_protocones)
376 cand_refs.clear();
377#endif
378
379 p_remain.clear();
380
381 return 0;
382}
383
384
385// full clearance
386//----------------
387int Csplit_merge::full_clear(){
388 partial_clear();
389
390 // clear previously allocated memory
391 if (indices != NULL){
392 delete[] indices;
393 }
394 particles.clear();
395
396 return 0;
397}
398
399
400// build the list 'p_uncol_hard' from p_remain by clustering collinear particles
401// note that thins in only used for stable-cone detection
402// so the parent_index field is unnecessary
403//-------------------------------------------------------------------------
404int Csplit_merge::merge_collinear_and_remove_soft(){
405 int i,j;
406 vector<Cmomentum> p_sorted;
407 bool collinear;
408 double dphi;
409
410 p_uncol_hard.clear();
411
412 // we first sort the particles according to their rapidity
413 for (i=0;i<n_left;i++)
414 p_sorted.push_back(p_remain[i]);
415 sort(p_sorted.begin(), p_sorted.end(), momentum_eta_less);
416
417 // then we cluster particles looping over the particles in the following way
418 // if (a particle i has same eta-phi a one after (j))
419 // then add momentum i to j
420 // else add i to the p_uncol_hard list
421 i = 0;
422 while (i<n_left){
423 // check if the particle passes the stable_cone_soft_pt2_cutoff
424 if (p_sorted[i].perp2()<stable_cone_soft_pt2_cutoff) {
425 i++;
426 continue;
427 }
428
429 // check if there is(are) particle(s) with the 'same' eta
430 collinear = false;
431 j=i+1;
432 while ((j<n_left) && (fabs(p_sorted[j].eta-p_sorted[i].eta)<EPSILON_COLLINEAR) && (!collinear)){
433 dphi = fabs(p_sorted[j].phi-p_sorted[i].phi);
434 if (dphi>M_PI) dphi = twopi-dphi;
435 if (dphi<EPSILON_COLLINEAR){
436 // i is collinear with j; add the momentum (i) to the momentum (j)
437 p_sorted[j] += p_sorted[i];
438 // set collinearity test to true
439 collinear = true;
440 }
441 j++;
442 }
443 // if no collinearity detected, add the particle to our list
444 if (!collinear)
445 p_uncol_hard.push_back(p_sorted[i]);
446 i++;
447 }
448
449 return 0;
450}
451
452
453// add a list of protocones
454// - protocones list of protocones (initial jet candidates)
455// - R2 cone radius (squared)
456// - ptmin minimal pT allowed for jets
457//-------------------------------------------------------------
458int Csplit_merge::add_protocones(vector<Cmomentum> *protocones, double R2, double ptmin){
459 int i;
460 Cmomentum *c;
461 Cmomentum *v;
462 double eta, phi;
463 double dx, dy;
464 double R;
465 Cjet jet;
466
467 if (protocones->size()==0)
468 return 1;
469
470 pt_min2 = ptmin*ptmin;
471 R = sqrt(R2);
472
473 // browse protocones
474 // for each of them, build the list of particles in them
475 for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
476 // initialise variables
477 c = &(*p_it);
478
479 // note: cones have been tested => their (eta,phi) coordinates are computed
480 eta = c->eta;
481 phi = c->phi;
482
483 // browse particles to create cone contents
484 // note that jet is always initialised with default values at this level
485 jet.v = Cmomentum();
486 jet.pt_tilde=0;
487 jet.contents.clear();
488 for (i=0;i<n_left;i++){
489 v = &(p_remain[i]);
490 // for security, avoid including particles with infinite rapidity)
491 // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING
492 //if (fabs(v->pz)!=v->E){
493 dx = eta - v->eta;
494 dy = fabs(phi - v->phi);
495 if (dy>M_PI)
496 dy -= twopi;
497 if (dx*dx+dy*dy<R2){
498 jet.contents.push_back(v->parent_index);
499 jet.v+= *v;
500 jet.pt_tilde+= pt[v->parent_index];
501 v->index=0;
502 }
503 }
504 jet.n=jet.contents.size();
505
506 // set the momentum in protocones
507 // (it was only known through eta and phi up to now)
508 *c = jet.v;
509 c->eta = eta; // restore exact original coords
510 c->phi = phi; // to avoid rounding error inconsistencies
511
512 // set the jet range
513 jet.range=Ceta_phi_range(eta,phi,R);
514
515#ifdef DEBUG_SPLIT_MERGE
516 cout << "adding jet: ";
517 for (int i2=0;i2<jet.n;i2++)
518 cout << jet.contents[i2] << " ";
519 cout << endl;
520#endif
521
522 // add it to the list of jets
523 insert(jet);
524 }
525
526 // update list of included particles
527 n_pass++;
528
529#ifdef DEBUG_SPLIT_MERGE
530 cout << "remaining particles: ";
531#endif
532 int j=0;
533 for (i=0;i<n_left;i++){
534 if (p_remain[i].index){
535 // copy particle
536 p_remain[j]=p_remain[i];
537 p_remain[j].parent_index = p_remain[i].parent_index;
538 p_remain[j].index=1;
539 // set run in initial list
540 particles[p_remain[j].parent_index].index = n_pass;
541#ifdef DEBUG_SPLIT_MERGE
542 cout << p_remain[j].parent_index << " ";
543#endif
544 j++;
545 }
546 }
547#ifdef DEBUG_SPLIT_MERGE
548 cout << endl;
549#endif
550 n_left = j;
551 p_remain.resize(j);
552
553 merge_collinear_and_remove_soft();
554
555 return 0;
556}
557
558
[273e668]559/*
560 * remove the hardest protocone and declare it a a jet
561 * - protocones list of protocones (initial jet candidates)
562 * - R2 cone radius (squared)
563 * - ptmin minimal pT allowed for jets
564 * return 0 on success, 1 on error
565 *
566 * The list of remaining particles (and the uncollinear-hard ones)
567 * is updated.
568 */
569int Csplit_merge::add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin){
570
571 int i;
572 Cmomentum *c;
573 Cmomentum *v;
574 double eta, phi;
575 double dx, dy;
576 double R;
577 Cjet jet, jet_candidate;
578 bool found_jet = false;
579
580 if (protocones->size()==0)
581 return 1;
582
583 pt_min2 = ptmin*ptmin;
584 R = sqrt(R2);
585
586 // browse protocones
587 // for each of them, build the list of particles in them
588 for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
589 // initialise variables
590 c = &(*p_it);
591
592 // note: cones have been tested => their (eta,phi) coordinates are computed
593 eta = c->eta;
594 phi = c->phi;
595
596 // NOTE: this is common to this method and add_protocones, so it
597 // could be moved into a 'build_jet_from_protocone' method
598 //
599 // browse particles to create cone contents
600 jet_candidate.v = Cmomentum();
601 jet_candidate.pt_tilde=0;
602 jet_candidate.contents.clear();
603 for (i=0;i<n_left;i++){
604 v = &(p_remain[i]);
605 // for security, avoid including particles with infinite rapidity)
606 // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING
607 //if (fabs(v->pz)!=v->E){
608 dx = eta - v->eta;
609 dy = fabs(phi - v->phi);
610 if (dy>M_PI)
611 dy -= twopi;
612 if (dx*dx+dy*dy<R2){
613 jet_candidate.contents.push_back(v->parent_index);
614 jet_candidate.v+= *v;
615 jet_candidate.pt_tilde+= pt[v->parent_index];
616 v->index=0;
617 }
618 }
619 jet_candidate.n=jet_candidate.contents.size();
620
621 // set the momentum in protocones
622 // (it was only known through eta and phi up to now)
623 *c = jet_candidate.v;
624 c->eta = eta; // restore exact original coords
625 c->phi = phi; // to avoid rounding error inconsistencies
626
627 // set the jet range
628 jet_candidate.range=Ceta_phi_range(eta,phi,R);
629
630 // check that the protojet has large enough pt
631 if (jet_candidate.v.perp2()<pt_min2)
632 continue;
633
634 // assign the split-merge (or progressive-removal) squared scale variable
635 if (_user_scale) {
636 // sm_var2 is the signed square of the user scale returned
637 // for the jet candidate
638 jet_candidate.sm_var2 = (*_user_scale)(jet_candidate);
639 jet_candidate.sm_var2 *= abs(jet_candidate.sm_var2);
640 } else {
641 jet_candidate.sm_var2 = get_sm_var2(jet_candidate.v, jet_candidate.pt_tilde);
642 }
643
644 // now check if it is possibly the hardest
645 if ((! found_jet) ||
646 (_user_scale ? _user_scale->is_larger(jet_candidate, jet)
647 : ptcomparison(jet_candidate, jet))){
648 jet = jet_candidate;
649 found_jet = true;
650 }
651 }
652
653 // make sure at least one of the jets has passed the selection
654 if (!found_jet) return 1;
655
656 // add the jet to the list of jets
657 jets.push_back(jet);
658 jets[jets.size()-1].v.build_etaphi();
659
660#ifdef DEBUG_SPLIT_MERGE
661 cout << "PR-Jet " << jets.size() << " [size " << next_jet.contents.size() << "]:";
662#endif
663
664 // update the list of what particles are left
665 int p_remain_index = 0;
666 int contents_index = 0;
667 //sort(next_jet.contents.begin(),next_jet.contents.end());
668 for (int index=0;index<n_left;index++){
669 if ((contents_index<(int) jet.contents.size()) &&
670 (p_remain[index].parent_index == jet.contents[contents_index])){
671 // this particle belongs to the newly found jet
672 // set pass in initial list
673 particles[p_remain[index].parent_index].index = n_pass;
674#ifdef DEBUG_SPLIT_MERGE
675 cout << " " << jet.contents[contents_index];
676#endif
677 contents_index++;
678 } else {
679 // this particle still has to be clustered
680 p_remain[p_remain_index] = p_remain[index];
681 p_remain[p_remain_index].parent_index = p_remain[index].parent_index;
682 p_remain[p_remain_index].index=1;
683 p_remain_index++;
684 }
685 }
686 p_remain.resize(n_left-jet.contents.size());
687 n_left = p_remain.size();
688 jets[jets.size()-1].pass = particles[jet.contents[0]].index;
689
690 // increase the pass index
691 n_pass++;
692
693#ifdef DEBUG_SPLIT_MERGE
694 cout << endl;
695#endif
696
697 // male sure the list of uncol_hard particles (used for the next
698 // stable cone finding) is updated [could probably be more
699 // efficient]
700 merge_collinear_and_remove_soft();
701
702 return 0;
703}
704
[d7d2da3]705/*
706 * really do the splitting and merging
707 * At the end, the vector jets is filled with the jets found.
708 * the 'contents' field of each jets contains the indices
709 * of the particles included in that jet.
710 * - overlap_tshold threshold for splitting/merging transition
711 * - ptmin minimal pT allowed for jets
712 * return the number of jets is returned
713 ******************************************************************/
714int Csplit_merge::perform(double overlap_tshold, double ptmin){
715 // iterators for the 2 jets
716 cjet_iterator j1;
717 cjet_iterator j2;
718
719 pt_min2 = ptmin*ptmin;
720
721 if (candidates->size()==0)
722 return 0;
723
724 if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
725 ostringstream message;
726 message << "Illegal value for overlap_tshold, f = " << overlap_tshold;
727 message << " (legal values are 0<f<1)";
728 throw Csiscone_error(message.str());
729 }
730
731 // overlap (the contents of this variable depends on the choice for
732 // the split--merge variable.)
733 // Note that the square of the ovelap is used
734 double overlap2;
735
736 // avoid to compute tshold*tshold at each overlap
737 double overlap_tshold2 = overlap_tshold*overlap_tshold;
738
739 do{
740 if (candidates->size()>0){
741 // browse for the first jet
742 j1 = candidates->begin();
743
744 // if hardest jet does not pass threshold then nothing else will
745 // either so one stops the split merge.
746 if (j1->sm_var2<SM_var2_hardest_cut_off) {break;}
747
748 // browse for the second jet
749 j2 = j1;
750 j2++;
751 int j2_relindex = 1; // used only in ifdef, but costs little so keep it outside
752
753 while (j2 != candidates->end()){
754#ifdef DEBUG_SPLIT_MERGE
755 show();
756#endif
757 // check overlapping
758 if (get_overlap(*j1, *j2, &overlap2)){
759 // check if overlapping energy passes threshold
760 // Note that this depends on the ordering variable
761#ifdef DEBUG_SPLIT_MERGE
762 cout << "overlap between cdt 1 and cdt " << j2_relindex+1 << " with overlap "
763 << sqrt(overlap2/j2->sm_var2) << endl<<endl;
764#endif
765 if (overlap2<overlap_tshold2*j2->sm_var2){
766 // split jets
767 split(j1, j2);
768
769 // update iterators
770 j2 = j1 = candidates->begin();
771 j2_relindex = 0;
772 } else {
773 // merge jets
774 merge(j1, j2);
775
776 // update iterators
777 j2 = j1 = candidates->begin();
778 j2_relindex = 0;
779 }
780 }
781 // watch out: split/merge might have caused new jets with pt <
782 // ptmin to disappear, so the total number of jets may
783 // have changed by more than expected and j2 might already by
784 // the end of the candidates list...
785 j2_relindex++;
786 if (j2 != candidates->end()) j2++;
787 } // end of loop on the second jet
788
789 if (j1 != candidates->end()) {
790 // all "second jet" passed without overlapping
791 // (otherwise we won't leave the j2 loop)
792 // => store jet 1 as real jet
793 jets.push_back(*j1);
794 jets[jets.size()-1].v.build_etaphi();
795 // a bug where the contents has non-zero size has been cropping
796 // up in many contexts -- so catch it!
797 assert(j1->contents.size() > 0);
798 jets[jets.size()-1].pass = particles[j1->contents[0]].index;
799#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
800 cand_refs.erase(j1->v.ref);
801#endif
802 candidates->erase(j1);
803
804 //// test that the hardest jet pass the potential cut-off
805 //if ((candidates->size()!=0) &&
806 // (candidates->begin()->sm_var2<SM_var2_hardest_cut_off)){
807 // candidates->clear();
808 //}
809 }
810 }
811 } while (candidates->size()>0);
812
813 // sort jets by pT
814 sort(jets.begin(), jets.end(), jets_pt_less);
815#ifdef DEBUG_SPLIT_MERGE
816 show();
817#endif
818
819 return jets.size();
820}
821
822
823
824// save the event on disk
825// - flux stream used to save jet contents
826//--------------------------------------------
827int Csplit_merge::save_contents(FILE *flux){
828 jet_iterator it_j;
829 Cjet *j1;
830 int i1, i2;
831
832 fprintf(flux, "# %d jets found\n", (int) jets.size());
833 fprintf(flux, "# columns are: eta, phi, pt and number of particles for each jet\n");
834 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
835 j1 = &(*it_j);
836 j1->v.build_etaphi();
837 fprintf(flux, "%f\t%f\t%e\t%d\n",
838 j1->v.eta, j1->v.phi, j1->v.perp(), j1->n);
839 }
840
841 fprintf(flux, "# jet contents\n");
842 fprintf(flux, "# columns are: eta, phi, pt, particle index and jet number\n");
843 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
844 j1 = &(*it_j);
845 for (i2=0;i2<j1->n;i2++)
846 fprintf(flux, "%f\t%f\t%e\t%d\t%d\n",
847 particles[j1->contents[i2]].eta, particles[j1->contents[i2]].phi,
848 particles[j1->contents[i2]].perp(), j1->contents[i2], i1);
849 }
850
851 return 0;
852}
853
854
855// show current jets/candidate status
856//------------------------------------
857int Csplit_merge::show(){
858 jet_iterator it_j;
859 cjet_iterator it_c;
860 Cjet *j;
861 const Cjet *c;
862 int i1, i2;
863
864 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
865 j = &(*it_j);
866 fprintf(stdout, "jet %2d: %e\t%e\t%e\t%e\t", i1+1,
867 j->v.px, j->v.py, j->v.pz, j->v.E);
868 for (i2=0;i2<j->n;i2++)
869 fprintf(stdout, "%d ", j->contents[i2]);
870 fprintf(stdout, "\n");
871 }
872
873 for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
874 c = &(*it_c);
875 fprintf(stdout, "cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1,
876 c->v.px, c->v.py, c->v.pz, c->v.E, sqrt(c->sm_var2));
877 for (i2=0;i2<c->n;i2++)
878 fprintf(stdout, "%d ", c->contents[i2]);
879 fprintf(stdout, "\n");
880 }
881
882 fprintf(stdout, "\n");
883 return 0;
884}
885
886
887// get the overlap between 2 jets
888// - j1 first jet
889// - j2 second jet
890// - overlap2 returned overlap^2 (determined by the choice of SM variable)
891// return true if overlapping, false if disjoint
892//---------------------------------------------------------------------
893bool Csplit_merge::get_overlap(const Cjet &j1, const Cjet &j2, double *overlap2){
894 // check if ranges overlap
895 if (!is_range_overlap(j1.range,j2.range))
896 return false;
897
898 int i1,i2;
899 bool is_overlap;
900
901 // initialise
902 i1=i2=idx_size=0;
903 is_overlap = false;
904 Cmomentum v;
905 double pt_tilde=0.0;
906
907 // compute overlap
908 // at the same time, we store union in indices
909 do{
910 if (j1.contents[i1]<j2.contents[i2]){
911 indices[idx_size] = j1.contents[i1];
912 i1++;
913 } else if (j1.contents[i1]>j2.contents[i2]){
914 indices[idx_size] = j2.contents[i2];
915 i2++;
916 } else { // (j1.contents[i1]==j2.contents[i2])
917 v += particles[j1.contents[i1]];
918 pt_tilde += pt[j1.contents[i1]];
919 indices[idx_size] = j1.contents[i1];
920 i1++;
921 i2++;
922 is_overlap = true;
923 }
924 idx_size++;
925 } while ((i1<j1.n) && (i2<j2.n));
926
927 // finish computing union
928 // (only needed if overlap !)
929 if (is_overlap){
930 while (i1<j1.n){
931 indices[idx_size] = j1.contents[i1];
932 i1++;
933 idx_size++;
934 }
935 while (i2<j2.n){
936 indices[idx_size] = j2.contents[i2];
937 i2++;
938 idx_size++;
939 }
940 }
941
942 // assign the overlapping var as return variable
943 (*overlap2) = get_sm_var2(v, pt_tilde);
944
945 return is_overlap;
946}
947
948
949
950// split the two given jet.
951// during this procedure, the jets j1 & j2 are replaced
952// by 2 new jets. Common particles are associted to the
953// closest initial jet.
954// - it_j1 iterator of the first jet in 'candidates'
955// - it_j2 iterator of the second jet in 'candidates'
956// - j1 first jet (Cjet instance)
957// - j2 second jet (Cjet instance)
958// return true on success, false on error
959////////////////////////////////////////////////////////////////
960bool Csplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
961 int i1, i2;
962 Cjet jet1, jet2;
963 double eta1, phi1, pt1_weight, eta2, phi2, pt2_weight;
964 double dx1, dy1, dx2, dy2;
965 Cmomentum tmp;
966 Cmomentum *v;
967
968 // shorthand to avoid having "->" all over the place
969 const Cjet & j1 = * it_j1;
970 const Cjet & j2 = * it_j2;
971
972 i1=i2=0;
973 jet2.v = jet1.v = Cmomentum();
974 jet2.pt_tilde = jet1.pt_tilde = 0.0;
975
976 // compute centroids
977 // When use_pt_weighted_splitting is activated, the
978 // "geometrical" distance is weighted by the inverse
979 // of the pt of the protojet
980 // This is stored in pt{1,2}_weight
981 tmp = j1.v;
982 tmp.build_etaphi();
983 eta1 = tmp.eta;
984 phi1 = tmp.phi;
985 pt1_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0;
986
987 tmp = j2.v;
988 tmp.build_etaphi();
989 eta2 = tmp.eta;
990 phi2 = tmp.phi;
991 pt2_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0;
992
993 jet1.v = jet2.v = Cmomentum();
994
995 // compute jet splitting
996 do{
997 if (j1.contents[i1]<j2.contents[i2]){
998 // particle i1 belong only to jet 1
999 v = &(particles[j1.contents[i1]]);
1000 jet1.contents.push_back(j1.contents[i1]);
1001 jet1.v += *v;
1002 jet1.pt_tilde += pt[j1.contents[i1]];
1003 i1++;
1004 jet1.range.add_particle(v->eta,v->phi);
1005 } else if (j1.contents[i1]>j2.contents[i2]){
1006 // particle i2 belong only to jet 2
1007 v = &(particles[j2.contents[i2]]);
1008 jet2.contents.push_back(j2.contents[i2]);
1009 jet2.v += *v;
1010 jet2.pt_tilde += pt[j2.contents[i2]];
1011 i2++;
1012 jet2.range.add_particle(v->eta,v->phi);
1013 } else { // (j1.contents[i1]==j2.contents[i2])
1014 // common particle, decide which is the closest centre
1015 v = &(particles[j1.contents[i1]]);
1016
1017 // distance w.r.t. centroid 1
1018 dx1 = eta1 - v->eta;
1019 dy1 = fabs(phi1 - v->phi);
1020 if (dy1>M_PI)
1021 dy1 -= twopi;
1022
1023 // distance w.r.t. centroid 2
1024 dx2 = eta2 - v->eta;
1025 dy2 = fabs(phi2 - v->phi);
1026 if (dy2>M_PI)
1027 dy2 -= twopi;
1028
1029 //? what when == ?
1030 // When use_pt_weighted_splitting is activated, the
1031 // "geometrical" distance is weighted by the inverse
1032 // of the pt of the protojet
1033 double d1sq = (dx1*dx1+dy1*dy1)*pt1_weight;
1034 double d2sq = (dx2*dx2+dy2*dy2)*pt2_weight;
1035 // do bookkeeping on most ambiguous split
1036 if (fabs(d1sq-d2sq) < most_ambiguous_split)
1037 most_ambiguous_split = fabs(d1sq-d2sq);
1038
1039 if (d1sq<d2sq){
1040 // particle i1 belong only to jet 1
1041 jet1.contents.push_back(j1.contents[i1]);
1042 jet1.v += *v;
1043 jet1.pt_tilde += pt[j1.contents[i1]];
1044 jet1.range.add_particle(v->eta,v->phi);
1045 } else {
1046 // particle i2 belong only to jet 2
1047 jet2.contents.push_back(j2.contents[i2]);
1048 jet2.v += *v;
1049 jet2.pt_tilde += pt[j2.contents[i2]];
1050 jet2.range.add_particle(v->eta,v->phi);
1051 }
1052
1053 i1++;
1054 i2++;
1055 }
1056 } while ((i1<j1.n) && (i2<j2.n));
1057
1058 while (i1<j1.n){
1059 v = &(particles[j1.contents[i1]]);
1060 jet1.contents.push_back(j1.contents[i1]);
1061 jet1.v += *v;
1062 jet1.pt_tilde += pt[j1.contents[i1]];
1063 i1++;
1064 jet1.range.add_particle(v->eta,v->phi);
1065 }
1066 while (i2<j2.n){
1067 v = &(particles[j2.contents[i2]]);
1068 jet2.contents.push_back(j2.contents[i2]);
1069 jet2.v += *v;
1070 jet2.pt_tilde += pt[j2.contents[i2]];
1071 i2++;
1072 jet2.range.add_particle(v->eta,v->phi);
1073 }
1074
1075 // finalise jets
1076 jet1.n = jet1.contents.size();
1077 jet2.n = jet2.contents.size();
1078
1079 //jet1.range = j1.range;
1080 //jet2.range = j2.range;
1081
1082 // remove previous jets
1083#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1084 cand_refs.erase(j1.v.ref);
1085 cand_refs.erase(j2.v.ref);
1086#endif
1087 candidates->erase(it_j1);
1088 candidates->erase(it_j2);
1089
1090 // reinsert new ones
1091 insert(jet1);
1092 insert(jet2);
1093
1094 return true;
1095}
1096
1097// merge the two given jet.
1098// during this procedure, the jets j1 & j2 are replaced
1099// by 1 single jets containing both of them.
1100// - it_j1 iterator of the first jet in 'candidates'
1101// - it_j2 iterator of the second jet in 'candidates'
1102// return true on success, false on error
1103////////////////////////////////////////////////////////////////
1104bool Csplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
1105 Cjet jet;
1106 int i;
1107
1108 // build new jet
1109 // note: particles within j1 & j2 have already been stored in indices
1110 for (i=0;i<idx_size;i++){
1111 jet.contents.push_back(indices[i]);
1112 jet.v += particles[indices[i]];
1113 jet.pt_tilde += pt[indices[i]];
1114 }
1115 jet.n = jet.contents.size();
1116
1117 // deal with ranges
1118 jet.range = range_union(it_j1->range, it_j2->range);
1119
1120 // remove old candidates
1121#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1122 if (merge_identical_protocones){
1123 cand_refs.erase(it_j1->v.ref);
1124 cand_refs.erase(it_j2->v.ref);
1125 }
1126#endif
1127 candidates->erase(it_j1);
1128 candidates->erase(it_j2);
1129
1130 // reinsert new candidate
1131 insert(jet);
1132
1133 return true;
1134}
1135
1136/**
1137 * Check whether or not a jet has to be inserted in the
1138 * list of protojets. If it has, set its sm_variable and
1139 * insert it to the list of protojets.
1140 */
1141bool Csplit_merge::insert(Cjet &jet){
1142
1143 // eventually check that no other candidate are present with the
1144 // same cone contents. We recall that this automatic merging of
1145 // identical protocones can lead to infrared-unsafe situations.
1146#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1147 if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
1148 return false;
1149#endif
1150
1151 // check that the protojet has large enough pt
1152 if (jet.v.perp2()<pt_min2)
1153 return false;
1154
1155 // assign SM variable
1156 jet.sm_var2 = get_sm_var2(jet.v, jet.pt_tilde);
1157
1158 // insert the jet.
1159 candidates->insert(jet);
1160
1161 return true;
1162}
1163
1164/**
1165 * given a 4-momentum and its associated pT, return the
1166 * variable that has to be used for SM
1167 * \param v 4 momentum of the protojet
1168 * \param pt_tilde pt_tilde of the protojet
1169 */
1170double Csplit_merge::get_sm_var2(Cmomentum &v, double &pt_tilde){
1171 switch(ptcomparison.split_merge_scale) {
1172 case SM_pt: return v.perp2();
1173 case SM_mt: return v.perpmass2();
1174 case SM_pttilde: return pt_tilde*pt_tilde;
1175 case SM_Et: return v.Et2();
1176 default:
1177 throw Csiscone_error("Unsupported split-merge scale choice: "
1178 + ptcomparison.SM_scale_name());
1179 }
1180
1181 return 0.0;
1182}
1183
1184}
Note: See TracBrowser for help on using the repository browser.