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 // WARNING: this is not the main SISCone trunk but //
6 // an adaptation to spherical coordinates //
7 // For more details, see http://projects.hepforge.org/siscone //
8 // //
9 // Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez //
10 // //
11 // This program is free software; you can redistribute it and/or modify //
12 // it under the terms of the GNU General Public License as published by //
13 // the Free Software Foundation; either version 2 of the License, or //
14 // (at your option) any later version. //
15 // //
16 // This program is distributed in the hope that it will be useful, //
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
19 // GNU General Public License for more details. //
20 // //
21 // You should have received a copy of the GNU General Public License //
22 // along with this program; if not, write to the Free Software //
23 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
24 // //
25 // $Revision:: 390 $//
26 // $Date:: 2016-03-03 11:06:52 +0100 (Thu, 03 Mar 2016) $//
28 
29 #include <siscone/siscone_error.h>
30 #include "split_merge.h"
31 #include "momentum.h"
32 #include <limits> // for max
33 #include <iostream>
34 #include <algorithm>
35 #include <sstream>
36 #include <cassert>
37 #include <cmath>
38 
39 namespace siscone_spherical{
40 
41 using namespace std;
42 
43 /********************************************************
44  * class CSphjet implementation *
45  * real Jet information. *
46  * This class contains information for one single jet. *
47  * That is, first, its momentum carrying information *
48  * about its centre and pT, and second, its particle *
49  * contents *
50  ********************************************************/
51 // default ctor
52 //--------------
54  n = 0;
55  v = CSphmomentum();
56  E_tilde = 0.0;
57  sm_var2 = 0.0;
58  pass = CJET_INEXISTENT_PASS; // initialised to a value that should
59  // notappear in the end (after clustering)
60 }
61 
62 // default dtor
63 //--------------
65 
66 }
67 
68 // ordering of jets in E (e.g. used in final jets ordering)
69 //----------------------------------------------------------
70 bool jets_E_less(const CSphjet &j1, const CSphjet &j2){
71  return j1.v.E > j2.v.E;
72 }
73 
74 
75 /********************************************************
76  * CSphsplit_merge_ptcomparison implementation *
77  * This deals with the ordering of the jets candidates *
78  ********************************************************/
79 
80 // odering of two jets
81 // The variable of the ordering is pt or mt
82 // depending on 'split_merge_scale' choice
83 //
84 // with EPSILON_SPLITMERGE defined, this attempts to identify
85 // delicate cases where two jets have identical momenta except for
86 // softish particles -- the difference of pt's may not be correctly
87 // identified normally and so lead to problems for the fate of the
88 // softish particle.
89 //
90 // NB: there is a potential issue in momentum-conserving events,
91 // whereby the harder of two jets becomes ill-defined when a soft
92 // particle is emitted --- this may have a knock-on effect on
93 // subsequent split-merge steps in cases with sufficiently large R
94 // (but we don't know what the limit is...)
95 //------------------------------------------------------------------
96 bool CSphsplit_merge_ptcomparison::operator ()(const CSphjet &jet1, const CSphjet &jet2) const{
97  double q1, q2;
98 
99  // compute the value for comparison for both jets
100  // This depends on the choice of variable (mt is the default)
101  q1 = jet1.sm_var2;
102  q2 = jet2.sm_var2;
103 
104  bool res = q1 > q2;
105 
106  // if we enable the refined version of the comparison (see defines.h),
107  // we compute the difference more precisely when the two jets are very
108  // close in the ordering variable.
109 #ifdef EPSILON_SPLITMERGE
110  if ( (fabs(q1-q2) < EPSILON_SPLITMERGE*max(q1,q2)) &&
111  (jet1.v.ref != jet2.v.ref) ) {
112 #ifdef DEBUG_SPLIT_MERGE
113  cout << "Using high-precision ordering tests" << endl;
114 #endif
115  // get the momentum of the difference
116  CSphmomentum difference;
117  double E_tilde_difference;
118  get_difference(jet1,jet2,&difference,&E_tilde_difference);
119 
120  // use the following relation: pt1^2 - pt2^2 = (pt1+pt2)*(pt1-pt2)
121  double qdiff;
122  CSphmomentum sum = jet1.v ;
123  sum += jet2.v;
124  double E_tilde_sum = jet1.E_tilde + jet2.E_tilde;
125 
126  // depending on the choice of ordering variable, set the result
127  switch (split_merge_scale){
128  case SM_Etilde:
129  qdiff = E_tilde_sum*E_tilde_difference;
130  break;
131  case SM_E:
132  qdiff = sum.E*difference.E;
133  break;
134  default:
135  throw siscone::Csiscone_error("Unsupported split-merge scale choice: "
136  + SM_scale_name());
137  }
138  res = qdiff > 0;
139  }
140 #endif // EPSILON_SPLITMERGE
141 
142  return res;
143 }
144 
145 
148 std::string split_merge_scale_name(Esplit_merge_scale sms) {
149  switch(sms) {
150  case SM_E:
151  return "E (IR unsafe for pairs of identical decayed heavy particles)";
152  case SM_Etilde:
153  return "Etilde (sum of E.[1+sin^2(theta_{i,jet})])";
154  default:
155  return "[SM scale without a name]";
156  }
157 }
158 
159 
160 // get the difference between 2 jets
161 // - j1 first jet
162 // - j2 second jet
163 // - v jet1-jet2
164 // - pt_tilde jet1-jet2 pt_tilde
165 // return true if overlapping, false if disjoint
166 //-----------------------------------------------
168  CSphmomentum *v, double *E_tilde) const {
169  int i1,i2;
170 
171  // initialise
172  i1=i2=0;
173  *v = CSphmomentum();
174  *E_tilde = 0.0;
175 
176  CSph3vector jet1_axis = j1.v;
177  //jet1_axis /= j1.v._norm;
178  jet1_axis /= j1.v.E;
179  CSph3vector jet2_axis = j2.v;
180  //jet2_axis /= j2.v._norm;
181  jet2_axis /= j2.v.E;
182 
183  // compute overlap
184  // at the same time, we store union in indices
185  // note tat for Etilde, we'll add the direct energy contributino at the end
186  do{
187  if (j1.contents[i1]==j2.contents[i2]) {
188  const CSphmomentum & p = (*particles)[j1.contents[i1]];
189  (*E_tilde) += p.E*((norm2_cross_product3(p,jet1_axis)-norm2_cross_product3(p,jet2_axis))/(*particles_norm2)[j1.contents[i1]]);
190  i1++;
191  i2++;
192  } else if (j1.contents[i1]<j2.contents[i2]){
193  const CSphmomentum & p = (*particles)[j1.contents[i1]];
194  (*v) += p;
195  (*E_tilde) += p.E*norm2_cross_product3(p,jet1_axis)/(*particles_norm2)[j1.contents[i1]];
196  i1++;
197  } else if (j1.contents[i1]>j2.contents[i2]){
198  const CSphmomentum &p = (*particles)[j2.contents[i2]];
199  (*v) -= p;
200  (*E_tilde) -= p.E*norm2_cross_product3(p,jet2_axis)/(*particles_norm2)[j2.contents[i2]];
201  i2++;
202  } else {
203  throw siscone::Csiscone_error("get_non_overlap reached part it should never have seen...");
204  }
205  } while ((i1<j1.n) && (i2<j2.n));
206 
207  // deal with particles at the end of the list...
208  while (i1 < j1.n) {
209  const CSphmomentum &p = (*particles)[j1.contents[i1]];
210  (*v) += p;
211  (*E_tilde) += p.E*norm2_cross_product3(p,jet1_axis)/(*particles_norm2)[j1.contents[i1]];
212  i1++;
213  }
214  while (i2 < j2.n) {
215  const CSphmomentum &p = (*particles)[j2.contents[i2]];
216  (*v) -= p;
217  (*E_tilde) -= p.E*norm2_cross_product3(p,jet2_axis)/(*particles_norm2)[j2.contents[i2]];
218  i2++;
219  }
220 
221  // add the direct energy contribution to Etilde
222  (*E_tilde) += v->E;
223 }
224 
225 
226 /********************************************************
227  * class CSphsplit_merge implementation *
228  * Class used to split and merge jets. *
229  ********************************************************/
230 // default ctor
231 //--------------
233  merge_identical_protocones = false;
234 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
235 #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
236  merge_identical_protocones = true;
237 #endif
238 #endif
239  _user_scale = NULL;
240  indices = NULL;
241 
242  // ensure that ptcomparison points to our set of particles (though params not correct)
243  ptcomparison.particles = &particles;
244  ptcomparison.particles_norm2 = &particles_norm2;
245  candidates.reset(new multiset<CSphjet,CSphsplit_merge_ptcomparison>(ptcomparison));
246 
247  // no hardest cut (col-unsafe)
248  SM_var2_hardest_cut_off = -numeric_limits<double>::max();
249 
250  // no energy cutoff for the particles to put in p_uncol_hard
251  stable_cone_soft_E2_cutoff = -1.0;
252 
253  // no pt-weighted splitting
254  use_E_weighted_splitting = false;
255 }
256 
257 
258 // default dtor
259 //--------------
261  full_clear();
262 }
263 
264 
265 // initialisation function
266 // - _particles list of particles
267 // - protocones list of protocones (initial jet candidates)
268 // - R2 cone radius (squared)
269 // - Emin minimal energy allowed for jets
270 //-------------------------------------------------------------
271 int CSphsplit_merge::init(vector<CSphmomentum> & /*_particles*/, vector<CSphmomentum> *protocones, double R2, double Emin){
272  // browse protocones
273  return add_protocones(protocones, R2, Emin);
274 }
275 
276 
277 // initialisation function for particle list
278 // - _particles list of particles
279 //-------------------------------------------------------------
280 int CSphsplit_merge::init_particles(vector<CSphmomentum> &_particles){
281  full_clear();
282 
283  // compute the list of particles
284  // here, we do not need to throw away particles
285  // with infinite rapidity (colinear with the beam)
286  particles = _particles;
287  n = particles.size();
288 
289  // store the particle norm^2
290  particles_norm2.resize(n);
291  for (int i=0;i<n;i++){
292  particles_norm2[i] = particles[i].norm2();
293  }
294 
295  // ensure that ptcomparison points to our set of particles (though params not correct)
296  ptcomparison.particles = &particles;
297  ptcomparison.particles_norm2 = &particles_norm2;
298 
299  // set up the list of particles left.
300  init_pleft();
301 
302  indices = new int[n];
303 
304  return 0;
305 }
306 
307 
308 // build initial list of remaining particles
309 //------------------------------------------
311  // at this level, we only rule out particles with
312  // infinite rapidity
313  // in the parent particle list, index contain the run
314  // at which particles are puts in jets:
315  // - -1 means infinity rapidity
316  // - 0 means not included
317  // - i mean included at run i
318  int i,j;
319 
320  // copy particles removing the ones with infinite rapidity
321  j=0;
322  p_remain.clear();
323  for (i=0;i<n;i++){
324  // set ref for checkxor
325  particles[i].ref.randomize();
326 
327  //REMOVED: check if rapidity is not infinite or ill-defined
328  //if (fabs(particles[i].pz) < (particles[i].E)){
329  p_remain.push_back(particles[i]);
330  // set up parent index for tracability
331  p_remain[j].parent_index = i;
332  // left particles are marked with a 1
333  // IMPORTANT NOTE: the meaning of index in p_remain is
334  // somehow different as in the initial particle list.
335  // here, within one pass, we use the index to track whether
336  // a particle is included in the current pass (index set to 0
337  // in add_protocones) or still remain (index still 1)
338  p_remain[j].index = 1;
339 
340  j++;
341  // set up parent-particle index
342  particles[i].index = 0;
343  //} else {
344  // particles[i].index = -1;
345  //}
346  }
347  n_left = p_remain.size();
348  n_pass = 0;
349 
350  merge_collinear_and_remove_soft();
351 
352  return 0;
353 }
354 
355 
356 // partial clearance
357 // we want to keep particle list and indices
358 // for future usage, so do not clear it !
359 // this is done in full_clear
360 //----------------------------------------
362  // release jets
363 
364  // set up the auto_ptr for the multiset with the _current_ state of
365  // ptcomparison (which may have changed since we constructed the
366  // class)
367  candidates.reset(new multiset<CSphjet,CSphsplit_merge_ptcomparison>(ptcomparison));
368 
369  // start off with huge number
370  most_ambiguous_split = numeric_limits<double>::max();
371 
372  jets.clear();
373 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
374  if (merge_identical_protocones)
375  cand_refs.clear();
376 #endif
377 
378  p_remain.clear();
379 
380  return 0;
381 }
382 
383 
384 // full clearance
385 //----------------
387  partial_clear();
388 
389  // clear previously allocated memory
390  if (indices != NULL){
391  delete[] indices;
392  }
393  particles.clear();
394 
395  return 0;
396 }
397 
398 
399 // build the list 'p_uncol_hard' from p_remain by clustering collinear particles
400 // note that thins in only used for stable-cone detection
401 // so the parent_index field is unnecessary
402 //-------------------------------------------------------------------------
404  int i,j;
405  vector<CSphmomentum> p_sorted;
406  bool collinear;
407  double dphi;
408 
409  p_uncol_hard.clear();
410 
411  // we first sort the particles according to their theta angle
412  for (i=0;i<n_left;i++)
413  p_sorted.push_back(p_remain[i]);
414  sort(p_sorted.begin(), p_sorted.end(), momentum_theta_less);
415 
416  // then we cluster particles looping over the particles in the following way
417  // if (a particle i has same eta-phi a one after (j))
418  // then add momentum i to j
419  // else add i to the p_uncol_hard list
420  i = 0;
421  while (i<n_left){
422  // check if the particle passes the stable_cone_soft_E2_cutoff
423  if (p_sorted[i].E*p_sorted[i].E<stable_cone_soft_E2_cutoff) {
424  i++;
425  continue;
426  }
427 
428  // check if there is(are) particle(s) with the 'same' theta
429  collinear = false;
430  j=i+1;
431  while ((j<n_left) && (fabs(p_sorted[j]._theta-p_sorted[i]._theta)<EPSILON_COLLINEAR) && (!collinear)){
432  dphi = fabs(p_sorted[j]._phi-p_sorted[i]._phi);
433  if (dphi>M_PI) dphi = twopi-dphi;
434  if (dphi<EPSILON_COLLINEAR){
435  // i is collinear with j; add the momentum (i) to the momentum (j)
436 #ifdef DEBUG_SPLIT_MERGE
437  cout << "# collinear merging at point " << p_sorted[i]._theta << ", " << p_sorted[j]._phi << endl;
438 #endif
439  p_sorted[j] += p_sorted[i];
440  //p_sorted[j].build_thetaphi();
441  p_sorted[j].build_norm();
442  // set collinearity test to true
443  collinear = true;
444  }
445  j++;
446  }
447  // if no collinearity detected, add the particle to our list
448  if (!collinear)
449  p_uncol_hard.push_back(p_sorted[i]);
450  i++;
451  }
452 
453  return 0;
454 }
455 
456 
457 // add a list of protocones
458 // - protocones list of protocones (initial jet candidates)
459 // - R2 cone radius (squared)
460 // - Emin minimal energy allowed for jets
461 //-------------------------------------------------------------
462 int CSphsplit_merge::add_protocones(vector<CSphmomentum> *protocones, double R2, double Emin){
463  int i;
464  CSphmomentum *c;
465  CSphmomentum *v;
466  double tan2R;
467  CSphjet jet;
468 
469  if (protocones->size()==0)
470  return 1;
471 
472  E_min = Emin;
473  double R = sqrt(R2);
474  tan2R = tan(R);
475  tan2R *= tan2R;
476 
477 #ifdef DEBUG_SPLIT_MERGE
478  cout << "particle list: ";
479  for (int i2=0;i2<n_left;i2++)
480  cout << p_remain[i2].parent_index << " "
481  << p_remain[i2].px << " " << p_remain[i2].py << " "
482  << p_remain[i2].pz << " " << p_remain[i2].E << endl;
483  cout << endl;
484 #endif
485 
486  // browse protocones
487  // for each of them, build the list of particles in them
488  for (vector<CSphmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
489  // initialise variables
490  c = &(*p_it);
491 
492  // browse particles to create cone contents
493  // note that jet is always initialised with default values at this level
494  jet.v = CSphmomentum();
495  jet.contents.clear();
496  for (i=0;i<n_left;i++){
497  v = &(p_remain[i]);
498  if (is_closer(v, c, tan2R)){
499  jet.contents.push_back(v->parent_index);
500  jet.v+= *v;
501  v->index=0;
502  }
503  }
504  jet.n=jet.contents.size();
505 
506  // compute Etilde for that jet.
507  // we can't do that before as it requires knowledge of the jet axis
508  // which has just been computed.
509  compute_Etilde(jet);
510 
511  // set the momentum in protocones
512  // (it was only known through its spatial coordinates up to now)
513  *c = jet.v;
514  c->build_thetaphi();
515 
516  // set the jet range
517  jet.range=CSphtheta_phi_range(c->_theta,c->_phi,R);
518 
519 #ifdef DEBUG_SPLIT_MERGE
520  cout << "adding protojet: ";
521 
522  unsigned int phirange=jet.range.phi_range;
523  for (unsigned int i2=0;i2<32;i2++) fprintf(stdout, "%d", (phirange&(1<<i2)) >> i2 );
524  fprintf(stdout, "\t");
525  unsigned int thetarange=jet.range.theta_range;
526  for (unsigned int i2=0;i2<32;i2++) fprintf(stdout, "%d", (thetarange&(1<<i2)) >> i2);
527  fprintf(stdout, "\t");
528 
529  for (int i2=0;i2<jet.n;i2++)
530  cout << jet.contents[i2] << " ";
531  cout << endl;
532 #endif
533 
534  // add it to the list of jets
535  insert(jet);
536  }
537 
538  // update list of included particles
539  n_pass++;
540 
541 #ifdef DEBUG_SPLIT_MERGE
542  cout << "remaining particles: ";
543 #endif
544  int j=0;
545  for (i=0;i<n_left;i++){
546  if (p_remain[i].index){
547  // copy particle
548  p_remain[j]=p_remain[i];
549  p_remain[j].parent_index = p_remain[i].parent_index;
550  p_remain[j].index=1;
551  // set run in initial list
552  particles[p_remain[j].parent_index].index = n_pass;
553 #ifdef DEBUG_SPLIT_MERGE
554  cout << p_remain[j].parent_index << " ";
555 #endif
556  j++;
557  }
558  }
559 #ifdef DEBUG_SPLIT_MERGE
560  cout << endl;
561 #endif
562  n_left = j;
563  p_remain.resize(j);
564 
565  merge_collinear_and_remove_soft();
566 
567  return 0;
568 }
569 
570 /*
571  * remove the hardest protocone and declare it a a jet
572  * - protocones list of protocones (initial jet candidates)
573  * - R2 cone radius (squared)
574 // - Emin minimal energy allowed for jets
575  * return 0 on success, 1 on error
576  *
577  * The list of remaining particles (and the uncollinear-hard ones)
578  * is updated.
579  */
580 int CSphsplit_merge::add_hardest_protocone_to_jets(std::vector<CSphmomentum> *protocones, double R2, double Emin){
581 
582  int i;
583  CSphmomentum *c;
584  CSphmomentum *v;
585  double R, tan2R;
586  CSphjet jet, jet_candidate;
587  bool found_jet = false;
588 
589  if (protocones->size()==0)
590  return 1;
591 
592  E_min = Emin;
593  R = sqrt(R2);
594  tan2R = tan(R);
595  tan2R *= tan2R;
596 
597  // browse protocones
598  // for each of them, build the list of particles in them
599  for (vector<CSphmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
600  // initialise variables
601  c = &(*p_it);
602 
603  // browse particles to create cone contents
604  // note that jet is always initialised with default values at this level
605  jet_candidate.v = CSphmomentum();
606  jet_candidate.contents.clear();
607  for (i=0;i<n_left;i++){
608  v = &(p_remain[i]);
609  if (is_closer(v, c, tan2R)){
610  jet_candidate.contents.push_back(v->parent_index);
611  jet_candidate.v+= *v;
612  v->index=0;
613  }
614  }
615  jet_candidate.n=jet_candidate.contents.size();
616 
617  // compute Etilde for that jet.
618  // we can't do that before as it requires knowledge of the jet axis
619  // which has just been computed.
620  compute_Etilde(jet_candidate);
621 
622  // set the momentum in protocones
623  // (it was only known through its spatial coordinates up to now)
624  *c = jet_candidate.v;
625  c->build_thetaphi();
626 
627  // set the jet range
628  jet_candidate.range=CSphtheta_phi_range(c->_theta,c->_phi,R);
629 
630  // check that the protojet has large enough pt
631  if (jet_candidate.v.E<E_min)
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.E_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_thetaphi();
659  jets[jets.size()-1].v.build_norm();
660 
661 #ifdef DEBUG_SPLIT_MERGE
662  cout << "PR-Jet " << jets.size() << " [size " << jet.contents.size() << "]:";
663 #endif
664 
665  // update the list of what particles are left
666  int p_remain_index = 0;
667  int contents_index = 0;
668  //sort(next_jet.contents.begin(),next_jet.contents.end());
669  for (int index=0;index<n_left;index++){
670  if ((contents_index<(int) jet.contents.size()) &&
671  (p_remain[index].parent_index == jet.contents[contents_index])){
672  // this particle belongs to the newly found jet
673  // set pass in initial list
674  particles[p_remain[index].parent_index].index = n_pass;
675 #ifdef DEBUG_SPLIT_MERGE
676  cout << " " << jet.contents[contents_index];
677 #endif
678  contents_index++;
679  } else {
680  // this particle still has to be clustered
681  p_remain[p_remain_index] = p_remain[index];
682  p_remain[p_remain_index].parent_index = p_remain[index].parent_index;
683  p_remain[p_remain_index].index=1;
684  p_remain_index++;
685  }
686  }
687  p_remain.resize(n_left-jet.contents.size());
688  n_left = p_remain.size();
689  jets[jets.size()-1].pass = particles[jet.contents[0]].index;
690 
691  // update list of included particles
692  n_pass++;
693 
694 #ifdef DEBUG_SPLIT_MERGE
695  cout << endl;
696 #endif
697 
698  // male sure the list of uncol_hard particles (used for the next
699  // stable cone finding) is updated [could probably be more
700  // efficient]
701  merge_collinear_and_remove_soft();
702 
703  return 0;
704 }
705 
706 /*
707  * really do the splitting and merging
708  * At the end, the vector jets is filled with the jets found.
709  * the 'contents' field of each jets contains the indices
710  * of the particles included in that jet.
711  * - overlap_tshold threshold for splitting/merging transition
712  * - Emin minimal energy allowed for jets
713  * return the number of jets is returned
714  ******************************************************************/
715 int CSphsplit_merge::perform(double overlap_tshold, double Emin){
716  // iterators for the 2 jets
717  cjet_iterator j1;
718  cjet_iterator j2;
719 
720  E_min = Emin;
721 
722  if (candidates->size()==0)
723  return 0;
724 
725  if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
726  ostringstream message;
727  message << "Illegal value for overlap_tshold, f = " << overlap_tshold;
728  message << " (legal values are 0<f<1)";
729  throw siscone::Csiscone_error(message.str());
730  }
731 
732  // overlap (the contents of this variable depends on the choice for
733  // the split--merge variable.)
734  // Note that the square of the ovelap is used
735  double overlap2;
736 
737  // avoid to compute tshold*tshold at each overlap
738  double overlap_tshold2 = overlap_tshold*overlap_tshold;
739 
740  do{
741  if (candidates->size()>0){
742  // browse for the first jet
743  j1 = candidates->begin();
744 
745  // if hardest jet does not pass threshold then nothing else will
746  // either so one stops the split merge.
747  //if (j1->sm_var2<SM_var2_hardest_cut_off) {break;}
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  if (j2_relindex==1) show();
758  cout << "check overlap between cdt 1 and cdt " << j2_relindex+1 << " with overlap " << endl;
759 #endif
760  // check overlapping
761  if (get_overlap(*j1, *j2, &overlap2)){
762  // check if overlapping energy passes threshold
763  // Note that this depends on the ordering variable
764 #ifdef DEBUG_SPLIT_MERGE
765  cout << "overlap between cdt 1 and cdt " << j2_relindex+1 << " with overlap "
766  << sqrt(overlap2)/j2->v.E << endl<<endl;
767 #endif
768  // We use the energy for the overlap computation
769  if (overlap2<overlap_tshold2*sqr(j2->v.E)){
770 #ifdef DEBUG_SPLIT_MERGE
771  cout << " --> split" << endl<<endl;
772 #endif
773  // split jets
774  split(j1, j2);
775 
776  // update iterators
777  j2 = j1 = candidates->begin();
778  j2_relindex = 0;
779  } else {
780 #ifdef DEBUG_SPLIT_MERGE
781  cout << " --> merge" << endl<<endl;
782 #endif
783  // merge jets
784  merge(j1, j2);
785 
786  // update iterators
787  j2 = j1 = candidates->begin();
788  j2_relindex = 0;
789  }
790  }
791  // watch out: split/merge might have caused new jets with E <
792  // Emin to disappear, so the total number of jets may
793  // have changed by more than expected and j2 might already by
794  // the end of the candidates list...
795  j2_relindex++;
796  if (j2 != candidates->end()) j2++;
797  } // end of loop on the second jet
798 
799  if (j1 != candidates->end()) {
800  // all "second jet" passed without overlapping
801  // (otherwise we won't leave the j2 loop)
802  // => store jet 1 as real jet
803  jets.push_back(*j1);
804  jets[jets.size()-1].v.build_thetaphi();
805  jets[jets.size()-1].v.build_norm();
806  // a bug where the contents has non-zero size has been cropping
807  // up in many contexts -- so catch it!
808  assert(j1->contents.size() > 0);
809  jets[jets.size()-1].pass = particles[j1->contents[0]].index;
810 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
811  cand_refs.erase(j1->v.ref);
812 #endif
813  candidates->erase(j1);
814  }
815  }
816  } while (candidates->size()>0);
817 
818  // sort jets by Energy
819  sort(jets.begin(), jets.end(), jets_E_less);
820 #ifdef DEBUG_SPLIT_MERGE
821  show();
822 #endif
823 
824  return jets.size();
825 }
826 
827 
828 
829 // save the event on disk
830 // - flux stream used to save jet contents
831 //--------------------------------------------
833  jet_iterator it_j;
834  CSphjet *j1;
835  int i1, i2;
836 
837  fprintf(flux, "# %d jets found\n", (int) jets.size());
838  fprintf(flux, "# columns are: px, py, pz, E and number of particles for each jet\n");
839  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
840  j1 = &(*it_j);
841  fprintf(flux, "%e\t%e\t%e\t%e\t%d\n",
842  j1->v.px, j1->v.py, j1->v.pz, j1->v.E, j1->n);
843  }
844 
845  fprintf(flux, "# jet contents\n");
846  fprintf(flux, "# columns are: px, py, pz, E, particle index and jet number\n");
847  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
848  j1 = &(*it_j);
849  for (i2=0;i2<j1->n;i2++)
850  fprintf(flux, "%e\t%e\t%e\t%e\t%d\t%d\n",
851  particles[j1->contents[i2]].px, particles[j1->contents[i2]].py,
852  particles[j1->contents[i2]].pz, particles[j1->contents[i2]].E,
853  j1->contents[i2], i1);
854  }
855 
856  return 0;
857 }
858 
859 
860 // show current jets/candidate status
861 //------------------------------------
863  jet_iterator it_j;
864  cjet_iterator it_c;
865  CSphjet *j;
866  const CSphjet *c;
867  int i1, i2;
868 
869  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
870  j = &(*it_j);
871  fprintf(stdout, "jet %2d: %e\t%e\t%e\t%e\t", i1+1,
872  j->v.px, j->v.py, j->v.pz, j->v.E);
873 
874  unsigned int phirange=j->range.phi_range;
875  for (i2=0;i2<32;i2++) fprintf(stdout, "%d", (phirange&(1<<i2)) >> i2 );
876  fprintf(stdout, "\t");
877  unsigned int thetarange=j->range.theta_range;
878  for (i2=0;i2<32;i2++) fprintf(stdout, "%d", (thetarange&(1<<i2)) >> i2);
879  fprintf(stdout, "\t");
880 
881  for (i2=0;i2<j->n;i2++)
882  fprintf(stdout, "%d ", j->contents[i2]);
883  fprintf(stdout, "\n");
884  }
885 
886  for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
887  c = &(*it_c);
888  fprintf(stdout, "cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1,
889  c->v.px, c->v.py, c->v.pz, c->v.E, sqrt(c->sm_var2));
890 
891  unsigned int phirange=c->range.phi_range;
892  for (i2=0;i2<32;i2++) fprintf(stdout, "%d", (phirange&(1<<i2)) >> i2 );
893  fprintf(stdout, "\t");
894  unsigned int thetarange=c->range.theta_range;
895  for (i2=0;i2<32;i2++) fprintf(stdout, "%d", (thetarange&(1<<i2)) >> i2);
896  fprintf(stdout, "\t");
897 
898  for (i2=0;i2<c->n;i2++)
899  fprintf(stdout, "%d ", c->contents[i2]);
900  fprintf(stdout, "\n");
901  }
902 
903  fprintf(stdout, "\n");
904  return 0;
905 }
906 
907 
908 // get the overlap between 2 jets
909 // - j1 first jet
910 // - j2 second jet
911 // - overlap2 returned overlap^2 (determined by the choice of SM variable)
912 // return true if overlapping, false if disjoint
913 //---------------------------------------------------------------------
914 bool CSphsplit_merge::get_overlap(const CSphjet &j1, const CSphjet &j2, double *overlap2){
915  // check if ranges overlap
916  if (!is_range_overlap(j1.range,j2.range))
917  return false;
918 
919  int i1,i2;
920  bool is_overlap;
921 
922  // initialise
923  i1=i2=idx_size=0;
924  is_overlap = false;
925  CSphmomentum v;
926 
927  // compute overlap
928  // at the same time, we store union in indices
929  do{
930  if (j1.contents[i1]<j2.contents[i2]){
931  indices[idx_size] = j1.contents[i1];
932  i1++;
933  } else if (j1.contents[i1]>j2.contents[i2]){
934  indices[idx_size] = j2.contents[i2];
935  i2++;
936  } else { // (j1.contents[i1]==j2.contents[i2])
937  v += particles[j2.contents[i2]];
938  indices[idx_size] = j2.contents[i2];
939  i1++;
940  i2++;
941  is_overlap = true;
942  }
943  idx_size++;
944  } while ((i1<j1.n) && (i2<j2.n));
945 
946  // finish computing union
947  // (only needed if overlap !)
948  if (is_overlap){
949  while (i1<j1.n){
950  indices[idx_size] = j1.contents[i1];
951  i1++;
952  idx_size++;
953  }
954  while (i2<j2.n){
955  indices[idx_size] = j2.contents[i2];
956  i2++;
957  idx_size++;
958  }
959  }
960 
961  // assign the overlapping var as return variable
962  (*overlap2) = sqr(v.E); //get_sm_var2(v, E_tilde);
963 
964  return is_overlap;
965 }
966 
967 
968 
969 // split the two given jet.
970 // during this procedure, the jets j1 & j2 are replaced
971 // by 2 new jets. Common particles are associted to the
972 // closest initial jet.
973 // - it_j1 iterator of the first jet in 'candidates'
974 // - it_j2 iterator of the second jet in 'candidates'
975 // - j1 first jet (CSphjet instance)
976 // - j2 second jet (CSphjet instance)
977 // return true on success, false on error
979 bool CSphsplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
980  int i1, i2;
981  CSphjet jet1, jet2;
982  double E1_weight, E2_weight;
983  CSphmomentum tmp;
984  CSphmomentum *v;
985 
986  // shorthand to avoid having "->" all over the place
987  const CSphjet & j1 = * it_j1;
988  const CSphjet & j2 = * it_j2;
989 
990  i1=i2=0;
991  jet2.v = jet1.v = CSphmomentum();
992 
993  // compute centroids
994  // When use_E_weighted_splitting is activated, the
995  // "geometrical" distance is weighted by the inverse
996  // of the E of the protojet
997  // This is stored in E{1,2}_weight
998  E1_weight = (use_E_weighted_splitting) ? 1.0/j1.v.E/j1.v.E : 1.0;
999  E2_weight = (use_E_weighted_splitting) ? 1.0/j2.v.E/j2.v.E : 1.0;
1000 
1001  // compute jet splitting
1002  do{
1003  if (j1.contents[i1]<j2.contents[i2]){
1004  // particle i1 belong only to jet 1
1005  v = &(particles[j1.contents[i1]]);
1006  jet1.contents.push_back(j1.contents[i1]);
1007  jet1.v += *v;
1008  //jet1.pt_tilde += pt[j1.contents[i1]];
1009  i1++;
1010  jet1.range.add_particle(v->_theta,v->_phi);
1011  } else if (j1.contents[i1]>j2.contents[i2]){
1012  // particle i2 belong only to jet 2
1013  v = &(particles[j2.contents[i2]]);
1014  jet2.contents.push_back(j2.contents[i2]);
1015  jet2.v += *v;
1016  //jet2.pt_tilde += pt[j2.contents[i2]];
1017  i2++;
1018  jet2.range.add_particle(v->_theta,v->_phi);
1019  } else { // (j1.contents[i1]==j2.contents[i2])
1020  // common particle, decide which is the closest centre
1021  v = &(particles[j1.contents[i1]]);
1022 
1023  //TODO: improve this brutal use of atan2 and sqrt !!!!
1024 
1025  //? what when == ?
1026  // When use_E_weighted_splitting is activated, the
1027  // "geometrical" distance is weighted by the inverse
1028  // of the E of the protojet
1029  double d1 = get_distance(&(j1.v), v)*E1_weight;
1030  double d2 = get_distance(&(j2.v), v)*E2_weight;
1031  // do bookkeeping on most ambiguous split
1032  if (fabs(d1-d2) < most_ambiguous_split)
1033  most_ambiguous_split = fabs(d1-d2);
1034 
1035  if (d1<d2){
1036  // particle i1 belong only to jet 1
1037  jet1.contents.push_back(j1.contents[i1]);
1038  jet1.v += *v;
1039  //jet1.pt_tilde += pt[j1.contents[i1]];
1040  jet1.range.add_particle(v->_theta,v->_phi);
1041  } else {
1042  // particle i2 belong only to jet 2
1043  jet2.contents.push_back(j2.contents[i2]);
1044  jet2.v += *v;
1045  //jet2.pt_tilde += pt[j2.contents[i2]];
1046  jet2.range.add_particle(v->_theta,v->_phi);
1047  }
1048 
1049  i1++;
1050  i2++;
1051  }
1052  } while ((i1<j1.n) && (i2<j2.n));
1053 
1054  while (i1<j1.n){
1055  v = &(particles[j1.contents[i1]]);
1056  jet1.contents.push_back(j1.contents[i1]);
1057  jet1.v += *v;
1058  //jet1.pt_tilde += pt[j1.contents[i1]];
1059  i1++;
1060  jet1.range.add_particle(v->_theta,v->_phi);
1061  }
1062  while (i2<j2.n){
1063  v = &(particles[j2.contents[i2]]);
1064  jet2.contents.push_back(j2.contents[i2]);
1065  jet2.v += *v;
1066  //jet2.pt_tilde += pt[j2.contents[i2]];
1067  i2++;
1068  jet2.range.add_particle(v->_theta,v->_phi);
1069  }
1070 
1071  // finalise jets
1072  jet1.n = jet1.contents.size();
1073  jet2.n = jet2.contents.size();
1074 
1075  // now the jet axis is known, we can compute Etilde
1076  compute_Etilde(jet1);
1077  compute_Etilde(jet2);
1078 
1079  // remove previous jets
1080 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1081  cand_refs.erase(j1.v.ref);
1082  cand_refs.erase(j2.v.ref);
1083 #endif
1084  candidates->erase(it_j1);
1085  candidates->erase(it_j2);
1086 
1087  // reinsert new ones
1088  insert(jet1);
1089  insert(jet2);
1090 
1091  return true;
1092 }
1093 
1094 // merge the two given jet.
1095 // during this procedure, the jets j1 & j2 are replaced
1096 // by 1 single jets containing both of them.
1097 // - it_j1 iterator of the first jet in 'candidates'
1098 // - it_j2 iterator of the second jet in 'candidates'
1099 // return true on success, false on error
1101 bool CSphsplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
1102  CSphjet jet;
1103  int i;
1104 
1105  // build new jet
1106  // note: particles within j1 & j2 have already been stored in indices
1107  for (i=0;i<idx_size;i++){
1108  jet.contents.push_back(indices[i]);
1109  jet.v += particles[indices[i]];
1110  //jet.pt_tilde += pt[indices[i]];
1111  }
1112  jet.n = jet.contents.size();
1113 
1114  compute_Etilde(jet);
1115 
1116  // deal with ranges
1117  jet.range = range_union(it_j1->range, it_j2->range);
1118 
1119  // remove old candidates
1120 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1121  if (merge_identical_protocones){
1122  cand_refs.erase(it_j1->v.ref);
1123  cand_refs.erase(it_j2->v.ref);
1124  }
1125 #endif
1126  candidates->erase(it_j1);
1127  candidates->erase(it_j2);
1128 
1129  // reinsert new candidate
1130  insert(jet);
1131 
1132  return true;
1133 }
1134 
1140 bool CSphsplit_merge::insert(CSphjet &jet){
1141 
1142  // eventually check that no other candidate are present with the
1143  // same cone contents. We recall that this automatic merging of
1144  // identical protocones can lead to infrared-unsafe situations.
1145 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1146  if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
1147  return false;
1148 #endif
1149 
1150  // check that the protojet has large enough energy
1151  if (jet.v.E<E_min)
1152  return false;
1153 
1154  // assign SM variable
1155  jet.sm_var2 = get_sm_var2(jet.v, jet.E_tilde);
1156 
1157  // insert the jet.
1158  candidates->insert(jet);
1159 
1160  return true;
1161 }
1162 
1169 double CSphsplit_merge::get_sm_var2(CSphmomentum &v, double &E_tilde){
1170  switch(ptcomparison.split_merge_scale) {
1171  case SM_E: return v.E*v.E;
1172  case SM_Etilde: return E_tilde*E_tilde;
1173  default:
1174  throw siscone::Csiscone_error("Unsupported split-merge scale choice: "
1175  + ptcomparison.SM_scale_name());
1176  }
1177 
1178  //return 0.0;
1179 }
1180 
1181 
1182 
1184 void CSphsplit_merge::compute_Etilde(CSphjet &jet){
1185  jet.v.build_norm();
1186  jet.E_tilde=0.0;
1187  CSph3vector jet_axis = jet.v;
1188  //if (jet.v._norm==0){
1189  // jet_axis = CSph3vector(0.0,0.0,0.0);
1190  //} else {
1191  jet_axis/=jet.v.E;
1192  //}
1193  //cout << "~~~ Axis: " << jet.v.px << " " << jet.v.py << " " << jet.v.pz << " " << jet.v._norm << endl;
1194  //cout << "~~~ Axis: " << jet_axis.px << " " << jet_axis.py << " " << jet_axis.pz << endl;
1195  for (vector<int>::iterator cont_it=jet.contents.begin(); cont_it!=jet.contents.end(); cont_it++){
1196  const CSphmomentum &p = particles[*cont_it];
1197  jet.E_tilde+=p.E*(1.0+norm2_cross_product3(p,jet_axis)/particles_norm2[*cont_it]);
1198  }
1199 }
1200 
1201 }
class corresponding to errors that will be thrown by siscone
Definition: siscone_error.h:38
base class for managing the spatial part of Cmomentum (defined after)
Definition: momentum.h:54
double _theta
particle theta angle (available ONLY after a call to build_thetaphi)
Definition: momentum.h:136
double _phi
particle phi angle (available ONLY after a call to build_thetaphi)
Definition: momentum.h:137
siscone::Creference ref
reference number for the vector
Definition: momentum.h:142
void build_thetaphi()
just a useful tool to store theta and phi locally (in _theta and _phi) in case you need repeated acce...
Definition: momentum.cpp:153
real Jet information.
Definition: split_merge.h:56
double sm_var2
ordering variable used for ordering and overlap in the split–merge.
Definition: split_merge.h:80
double E_tilde
sum of E_i [ 1 +|p_i x p_J|^2/(|p_i|^2 E_J^2)]
Definition: split_merge.h:65
CSphmomentum v
jet momentum
Definition: split_merge.h:64
std::vector< int > contents
particle contents (list of indices)
Definition: split_merge.h:67
CSphtheta_phi_range range
covered range in eta-phi
Definition: split_merge.h:83
int n
number of particles inside
Definition: split_merge.h:66
base class for dynamic coordinates management
Definition: momentum.h:158
int index
internal particle number
Definition: momentum.h:211
int parent_index
particle number in the parent list
Definition: momentum.h:210
bool operator()(const CSphjet &jet1, const CSphjet &jet2) const
comparison of 2 CSphjet
Definition: split_merge.cpp:96
void get_difference(const CSphjet &j1, const CSphjet &j2, CSphmomentum *v, double *E_tilde) const
get the difference between 2 jets, calculated such that rounding errors will not affect the result ev...
int show()
show jets/candidates status
int perform(double overlap_tshold, double Emin=0.0)
really do the splitting and merging At the end, the vector jets is filled with the jets found.
int partial_clear()
partial clearance
int add_hardest_protocone_to_jets(std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
remove the hardest protocone and declare it a a jet
int init_pleft()
build initial list of left particles
int init(std::vector< CSphmomentum > &_particles, std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
initialisation function
int save_contents(FILE *flux)
save final jets
int merge_collinear_and_remove_soft()
build the list 'p_uncol_hard' from p_remain by clustering collinear particles note that thins in only...
int init_particles(std::vector< CSphmomentum > &_particles)
initialisation function for particle list
int add_protocones(std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
add a list of protocones
class for holding a covering range in eta-phi
Definition: geom_2d.h:52
unsigned int theta_range
theta range as a binary coding of covered cells
Definition: geom_2d.h:75
unsigned int phi_range
phi range as a binary coding of covered cells
Definition: geom_2d.h:78
#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