SISCone  3.0.5
split_merge.cpp
1 // 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:: 390 $//
24 // $Date:: 2016-03-03 11:06:52 +0100 (Thu, 03 Mar 2016) $//
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 
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 //--------------
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 //-----------------------------------------------------------
68 bool 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 //------------------------------------------------------------------
94 bool 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 
156 std::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 //-----------------------------------------------
179 void 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 //--------------
227  merge_identical_protocones = false;
228 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
229 #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
230  merge_identical_protocones = 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)
237  ptcomparison.particles = &particles;
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
245  stable_cone_soft_pt2_cutoff = -1.0;
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 //-------------------------------------------------------------
265 int 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 //-------------------------------------------------------------
274 int 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)
289  ptcomparison.particles = &particles;
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 
353  merge_collinear_and_remove_soft();
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
377  if (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 //----------------
390  partial_clear();
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 //-------------------------------------------------------------
460 int 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 
555  merge_collinear_and_remove_soft();
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  */
571 int 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]
702  merge_collinear_and_remove_soft();
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  ******************************************************************/
716 int 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 //---------------------------------------------------------------------
895 bool 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
962 bool 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
1106 bool 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
1124  if (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 
1143 bool 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 
1172 double Csplit_merge::get_sm_var2(Cmomentum &v, double &pt_tilde){
1173  switch(ptcomparison.split_merge_scale) {
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: "
1180  + ptcomparison.SM_scale_name());
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
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
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...
int init(std::vector< Cmomentum > &_particles, std::vector< Cmomentum > *protocones, double R2, double ptmin=0.0)
initialisation function
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
int partial_clear()
partial clearance
int show()
show jets/candidates status
int init_pleft()
build initial list of left particles
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.
int add_protocones(std::vector< Cmomentum > *protocones, double R2, double ptmin=0.0)
add a list of protocones
~Csplit_merge()
default dtor
int save_contents(FILE *flux)
save final jets
Csplit_merge()
default ctor
int full_clear()
full clearance
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