siscone is hosted by Hepforge, IPPP Durham

The SISCone Jet Algorithm

Version 3.0.6



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

The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated on Tue Jun 20 2023 18:08:37 for SISCone by  Doxygen 1.9.4