SISCone  2.0.6
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:: 322 $//
26 // $Date:: 2011-11-15 10:12:36 +0100 (Tue, 15 Nov 2011) $//
28 
29 #include <siscone/siscone_error.h>
30 #include "split_merge.h"
31 #include "momentum.h"
32 #include <math.h>
33 #include <limits> // for max
34 #include <iostream>
35 #include <algorithm>
36 #include <sstream>
37 #include <cassert>
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 }
59 
60 // default dtor
61 //--------------
63 
64 }
65 
66 // ordering of jets in E (e.g. used in final jets ordering)
67 //----------------------------------------------------------
68 bool jets_E_less(const CSphjet &j1, const CSphjet &j2){
69  return j1.v.E > j2.v.E;
70 }
71 
72 
73 /********************************************************
74  * CSphsplit_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 CSphsplit_merge_ptcomparison::operator ()(const CSphjet &jet1, const CSphjet &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 #ifdef DEBUG_SPLIT_MERGE
111  cout << "Using high-precision ordering tests" << endl;
112 #endif
113  // get the momentum of the difference
114  CSphmomentum difference;
115  double E_tilde_difference;
116  get_difference(jet1,jet2,&difference,&E_tilde_difference);
117 
118  // use the following relation: pt1^2 - pt2^2 = (pt1+pt2)*(pt1-pt2)
119  double qdiff;
120  CSphmomentum sum = jet1.v ;
121  sum += jet2.v;
122  double E_tilde_sum = jet1.E_tilde + jet2.E_tilde;
123 
124  // depending on the choice of ordering variable, set the result
125  switch (split_merge_scale){
126  case SM_Etilde:
127  qdiff = E_tilde_sum*E_tilde_difference;
128  break;
129  case SM_E:
130  qdiff = sum.E*difference.E;
131  break;
132  default:
133  throw siscone::Csiscone_error("Unsupported split-merge scale choice: "
134  + SM_scale_name());
135  }
136  res = qdiff > 0;
137  }
138 #endif // EPSILON_SPLITMERGE
139 
140  return res;
141 }
142 
143 
146 std::string split_merge_scale_name(Esplit_merge_scale sms) {
147  switch(sms) {
148  case SM_E:
149  return "E (IR unsafe for pairs of identical decayed heavy particles)";
150  case SM_Etilde:
151  return "Etilde (sum of E.[1+sin^2(theta_{i,jet})])";
152  default:
153  return "[SM scale without a name]";
154  }
155 }
156 
157 
158 // get the difference between 2 jets
159 // - j1 first jet
160 // - j2 second jet
161 // - v jet1-jet2
162 // - pt_tilde jet1-jet2 pt_tilde
163 // return true if overlapping, false if disjoint
164 //-----------------------------------------------
166  CSphmomentum *v, double *E_tilde) const {
167  int i1,i2;
168 
169  // initialise
170  i1=i2=0;
171  *v = CSphmomentum();
172  *E_tilde = 0.0;
173 
174  CSph3vector jet1_axis = j1.v;
175  //jet1_axis /= j1.v._norm;
176  jet1_axis /= j1.v.E;
177  CSph3vector jet2_axis = j2.v;
178  //jet2_axis /= j2.v._norm;
179  jet2_axis /= j2.v.E;
180 
181  // compute overlap
182  // at the same time, we store union in indices
183  // note tat for Etilde, we'll add the direct energy contributino at the end
184  do{
185  if (j1.contents[i1]==j2.contents[i2]) {
186  const CSphmomentum & p = (*particles)[j1.contents[i1]];
187  (*E_tilde) += p.E*((norm2_cross_product3(p,jet1_axis)-norm2_cross_product3(p,jet2_axis))/(*particles_norm2)[j1.contents[i1]]);
188  i1++;
189  i2++;
190  } else if (j1.contents[i1]<j2.contents[i2]){
191  const CSphmomentum & p = (*particles)[j1.contents[i1]];
192  (*v) += p;
193  (*E_tilde) += p.E*norm2_cross_product3(p,jet1_axis)/(*particles_norm2)[j1.contents[i1]];
194  i1++;
195  } else if (j1.contents[i1]>j2.contents[i2]){
196  const CSphmomentum &p = (*particles)[j2.contents[i2]];
197  (*v) -= p;
198  (*E_tilde) -= p.E*norm2_cross_product3(p,jet2_axis)/(*particles_norm2)[j2.contents[i2]];
199  i2++;
200  } else {
201  throw siscone::Csiscone_error("get_non_overlap reached part it should never have seen...");
202  }
203  } while ((i1<j1.n) && (i2<j2.n));
204 
205  // deal with particles at the end of the list...
206  while (i1 < j1.n) {
207  const CSphmomentum &p = (*particles)[j1.contents[i1]];
208  (*v) += p;
209  (*E_tilde) += p.E*norm2_cross_product3(p,jet1_axis)/(*particles_norm2)[j1.contents[i1]];
210  i1++;
211  }
212  while (i2 < j2.n) {
213  const CSphmomentum &p = (*particles)[j2.contents[i2]];
214  (*v) -= p;
215  (*E_tilde) -= p.E*norm2_cross_product3(p,jet2_axis)/(*particles_norm2)[j2.contents[i2]];
216  i2++;
217  }
218 
219  // add the direct energy contribution to Etilde
220  (*E_tilde) += v->E;
221 }
222 
223 
224 /********************************************************
225  * class CSphsplit_merge implementation *
226  * Class used to split and merge jets. *
227  ********************************************************/
228 // default ctor
229 //--------------
231  merge_identical_protocones = false;
232 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
233 #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
234  merge_identical_protocones = true;
235 #endif
236 #endif
237  indices = NULL;
238 
239  // ensure that ptcomparison points to our set of particles (though params not correct)
240  ptcomparison.particles = &particles;
241  ptcomparison.particles_norm2 = &particles_norm2;
242  candidates.reset(new multiset<CSphjet,CSphsplit_merge_ptcomparison>(ptcomparison));
243 
244  // no hardest cut (col-unsafe)
245  SM_var2_hardest_cut_off = -1.0;
246 
247  // no energy cutoff for the particles to put in p_uncol_hard
248  stable_cone_soft_E2_cutoff = -1.0;
249 
250  // no pt-weighted splitting
251  use_E_weighted_splitting = false;
252 }
253 
254 
255 // default dtor
256 //--------------
258  full_clear();
259 }
260 
261 
262 // initialisation function
263 // - _particles list of particles
264 // - protocones list of protocones (initial jet candidates)
265 // - R2 cone radius (squared)
266 // - Emin minimal energy allowed for jets
267 //-------------------------------------------------------------
268 int CSphsplit_merge::init(vector<CSphmomentum> & /*_particles*/, vector<CSphmomentum> *protocones, double R2, double Emin){
269  // browse protocones
270  return add_protocones(protocones, R2, Emin);
271 }
272 
273 
274 // initialisation function for particle list
275 // - _particles list of particles
276 //-------------------------------------------------------------
277 int CSphsplit_merge::init_particles(vector<CSphmomentum> &_particles){
278  full_clear();
279 
280  // compute the list of particles
281  // here, we do not need to throw away particles
282  // with infinite rapidity (colinear with the beam)
283  particles = _particles;
284  n = particles.size();
285 
286  // store the particle norm^2
287  particles_norm2.resize(n);
288  for (int i=0;i<n;i++){
289  particles_norm2[i] = particles[i].norm2();
290  }
291 
292  // ensure that ptcomparison points to our set of particles (though params not correct)
293  ptcomparison.particles = &particles;
294  ptcomparison.particles_norm2 = &particles_norm2;
295 
296  // set up the list of particles left.
297  init_pleft();
298 
299  indices = new int[n];
300 
301  return 0;
302 }
303 
304 
305 // build initial list of remaining particles
306 //------------------------------------------
308  // at this level, we only rule out particles with
309  // infinite rapidity
310  // in the parent particle list, index contain the run
311  // at which particles are puts in jets:
312  // - -1 means infinity rapidity
313  // - 0 means not included
314  // - i mean included at run i
315  int i,j;
316 
317  // copy particles removing the ones with infinite rapidity
318  j=0;
319  p_remain.clear();
320  for (i=0;i<n;i++){
321  // set ref for checkxor
322  particles[i].ref.randomize();
323 
324  //REMOVED: check if rapidity is not infinite or ill-defined
325  //if (fabs(particles[i].pz) < (particles[i].E)){
326  p_remain.push_back(particles[i]);
327  // set up parent index for tracability
328  p_remain[j].parent_index = i;
329  // left particles are marked with a 1
330  // IMPORTANT NOTE: the meaning of index in p_remain is
331  // somehow different as in the initial particle list.
332  // here, within one pass, we use the index to track whether
333  // a particle is included in the current pass (index set to 0
334  // in add_protocones) or still remain (index still 1)
335  p_remain[j].index = 1;
336 
337  j++;
338  // set up parent-particle index
339  particles[i].index = 0;
340  //} else {
341  // particles[i].index = -1;
342  //}
343  }
344  n_left = p_remain.size();
345  n_pass = 0;
346 
347  merge_collinear_and_remove_soft();
348 
349  return 0;
350 }
351 
352 
353 // partial clearance
354 // we want to keep particle list and indices
355 // for future usage, so do not clear it !
356 // this is done in full_clear
357 //----------------------------------------
359  // release jets
360 
361  // set up the auto_ptr for the multiset with the _current_ state of
362  // ptcomparison (which may have changed since we constructed the
363  // class)
364  candidates.reset(new multiset<CSphjet,CSphsplit_merge_ptcomparison>(ptcomparison));
365 
366  // start off with huge number
367  most_ambiguous_split = numeric_limits<double>::max();
368 
369  jets.clear();
370 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
371  if (merge_identical_protocones)
372  cand_refs.clear();
373 #endif
374 
375  p_remain.clear();
376 
377  return 0;
378 }
379 
380 
381 // full clearance
382 //----------------
384  partial_clear();
385 
386  // clear previously allocated memory
387  if (indices != NULL){
388  delete[] indices;
389  }
390  particles.clear();
391 
392  return 0;
393 }
394 
395 
396 // build the list 'p_uncol_hard' from p_remain by clustering collinear particles
397 // note that thins in only used for stable-cone detection
398 // so the parent_index field is unnecessary
399 //-------------------------------------------------------------------------
401  int i,j;
402  vector<CSphmomentum> p_sorted;
403  bool collinear;
404  double dphi;
405 
406  p_uncol_hard.clear();
407 
408  // we first sort the particles according to their theta angle
409  for (i=0;i<n_left;i++)
410  p_sorted.push_back(p_remain[i]);
411  sort(p_sorted.begin(), p_sorted.end(), momentum_theta_less);
412 
413  // then we cluster particles looping over the particles in the following way
414  // if (a particle i has same eta-phi a one after (j))
415  // then add momentum i to j
416  // else add i to the p_uncol_hard list
417  i = 0;
418  while (i<n_left){
419  // check if the particle passes the stable_cone_soft_E2_cutoff
420  if (p_sorted[i].E*p_sorted[i].E<stable_cone_soft_E2_cutoff) {
421  i++;
422  continue;
423  }
424 
425  // check if there is(are) particle(s) with the 'same' theta
426  collinear = false;
427  j=i+1;
428  while ((j<n_left) && (fabs(p_sorted[j]._theta-p_sorted[i]._theta)<EPSILON_COLLINEAR) && (!collinear)){
429  dphi = fabs(p_sorted[j]._phi-p_sorted[i]._phi);
430  if (dphi>M_PI) dphi = twopi-dphi;
431  if (dphi<EPSILON_COLLINEAR){
432  // i is collinear with j; add the momentum (i) to the momentum (j)
433 #ifdef DEBUG_SPLIT_MERGE
434  cout << "# collinear merging at point " << p_sorted[i]._theta << ", " << p_sorted[j]._phi << endl;
435 #endif
436  p_sorted[j] += p_sorted[i];
437  //p_sorted[j].build_thetaphi();
438  p_sorted[j].build_norm();
439  // set collinearity test to true
440  collinear = true;
441  }
442  j++;
443  }
444  // if no collinearity detected, add the particle to our list
445  if (!collinear)
446  p_uncol_hard.push_back(p_sorted[i]);
447  i++;
448  }
449 
450  return 0;
451 }
452 
453 
454 // add a list of protocones
455 // - protocones list of protocones (initial jet candidates)
456 // - R2 cone radius (squared)
457 // - Emin minimal energy allowed for jets
458 //-------------------------------------------------------------
459 int CSphsplit_merge::add_protocones(vector<CSphmomentum> *protocones, double R2, double Emin){
460  int i;
461  CSphmomentum *c;
462  CSphmomentum *v;
463  double tan2R;
464  CSphjet jet;
465 
466  if (protocones->size()==0)
467  return 1;
468 
469  E_min = Emin;
470  double R = sqrt(R2);
471  tan2R = tan(R);
472  tan2R *= tan2R;
473 
474 #ifdef DEBUG_SPLIT_MERGE
475  cout << "particle list: ";
476  for (int i2=0;i2<n_left;i2++)
477  cout << p_remain[i2].parent_index << " "
478  << p_remain[i2].px << " " << p_remain[i2].py << " "
479  << p_remain[i2].pz << " " << p_remain[i2].E << endl;
480  cout << endl;
481 #endif
482 
483  // browse protocones
484  // for each of them, build the list of particles in them
485  for (vector<CSphmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
486  // initialise variables
487  c = &(*p_it);
488 
489  // browse particles to create cone contents
490  // note that jet is always initialised with default values at this level
491  jet.v = CSphmomentum();
492  jet.contents.clear();
493  for (i=0;i<n_left;i++){
494  v = &(p_remain[i]);
495  if (is_closer(v, c, tan2R)){
496  jet.contents.push_back(v->parent_index);
497  jet.v+= *v;
498  v->index=0;
499  }
500  }
501  jet.n=jet.contents.size();
502 
503  // compute Etilde for that jet.
504  // we can't do that before as it requires knowledge of the jet axis
505  // which has just been computed.
506  compute_Etilde(jet);
507 
508  // set the momentum in protocones
509  // (it was only known through its spatial coordinates up to now)
510  *c = jet.v;
511  c->build_thetaphi();
512 
513  // set the jet range
514  jet.range=CSphtheta_phi_range(c->_theta,c->_phi,R);
515 
516 #ifdef DEBUG_SPLIT_MERGE
517  cout << "adding protojet: ";
518  for (int i2=0;i2<jet.n;i2++)
519  cout << jet.contents[i2] << " ";
520  cout << endl;
521 #endif
522 
523  // add it to the list of jets
524  insert(jet);
525  }
526 
527  // update list of included particles
528  n_pass++;
529 
530 #ifdef DEBUG_SPLIT_MERGE
531  cout << "remaining particles: ";
532 #endif
533  int j=0;
534  for (i=0;i<n_left;i++){
535  if (p_remain[i].index){
536  // copy particle
537  p_remain[j]=p_remain[i];
538  p_remain[j].parent_index = p_remain[i].parent_index;
539  p_remain[j].index=1;
540  // set run in initial list
541  particles[p_remain[j].parent_index].index = n_pass;
542 #ifdef DEBUG_SPLIT_MERGE
543  cout << p_remain[j].parent_index << " ";
544 #endif
545  j++;
546  }
547  }
548 #ifdef DEBUG_SPLIT_MERGE
549  cout << endl;
550 #endif
551  n_left = j;
552  p_remain.resize(j);
553 
554  merge_collinear_and_remove_soft();
555 
556  return 0;
557 }
558 
559 
560 /*
561  * really do the splitting and merging
562  * At the end, the vector jets is filled with the jets found.
563  * the 'contents' field of each jets contains the indices
564  * of the particles included in that jet.
565  * - overlap_tshold threshold for splitting/merging transition
566  * - Emin minimal energy allowed for jets
567  * return the number of jets is returned
568  ******************************************************************/
569 int CSphsplit_merge::perform(double overlap_tshold, double Emin){
570  // iterators for the 2 jets
571  cjet_iterator j1;
572  cjet_iterator j2;
573 
574  E_min = Emin;
575 
576  if (candidates->size()==0)
577  return 0;
578 
579  if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
580  ostringstream message;
581  message << "Illegal value for overlap_tshold, f = " << overlap_tshold;
582  message << " (legal values are 0<f<1)";
583  throw siscone::Csiscone_error(message.str());
584  }
585 
586  // overlap (the contents of this variable depends on the choice for
587  // the split--merge variable.)
588  // Note that the square of the ovelap is used
589  double overlap2;
590 
591  // avoid to compute tshold*tshold at each overlap
592  double overlap_tshold2 = overlap_tshold*overlap_tshold;
593 
594  do{
595  if (candidates->size()>0){
596  // browse for the first jet
597  j1 = candidates->begin();
598 
599  // if hardest jet does not pass threshold then nothing else will
600  // either so one stops the split merge.
601  //if (j1->sm_var2<SM_var2_hardest_cut_off) {break;}
602  if (j1->sm_var2<SM_var2_hardest_cut_off) {break;}
603 
604  // browse for the second jet
605  j2 = j1;
606  j2++;
607  int j2_relindex = 1; // used only in ifdef, but costs little so keep it outside
608 
609  while (j2 != candidates->end()){
610 #ifdef DEBUG_SPLIT_MERGE
611  show();
612 #endif
613  // check overlapping
614  if (get_overlap(*j1, *j2, &overlap2)){
615  // check if overlapping energy passes threshold
616  // Note that this depends on the ordering variable
617 #ifdef DEBUG_SPLIT_MERGE
618  cout << "overlap between cdt 1 and cdt " << j2_relindex+1 << " with overlap "
619  << sqrt(overlap2)/j2->v.E << endl<<endl;
620 #endif
621  // We use the energy for the overlap computation
622  if (overlap2<overlap_tshold2*sqr(j2->v.E)){
623  // split jets
624  split(j1, j2);
625 
626  // update iterators
627  j2 = j1 = candidates->begin();
628  j2_relindex = 0;
629  } else {
630  // merge jets
631  merge(j1, j2);
632 
633  // update iterators
634  j2 = j1 = candidates->begin();
635  j2_relindex = 0;
636  }
637  }
638  // watch out: split/merge might have caused new jets with E <
639  // Emin to disappear, so the total number of jets may
640  // have changed by more than expected and j2 might already by
641  // the end of the candidates list...
642  j2_relindex++;
643  if (j2 != candidates->end()) j2++;
644  } // end of loop on the second jet
645 
646  if (j1 != candidates->end()) {
647  // all "second jet" passed without overlapping
648  // (otherwise we won't leave the j2 loop)
649  // => store jet 1 as real jet
650  jets.push_back(*j1);
651  jets[jets.size()-1].v.build_thetaphi();
652  jets[jets.size()-1].v.build_norm();
653  // a bug where the contents has non-zero size has been cropping
654  // up in many contexts -- so catch it!
655  assert(j1->contents.size() > 0);
656  jets[jets.size()-1].pass = particles[j1->contents[0]].index;
657 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
658  cand_refs.erase(j1->v.ref);
659 #endif
660  candidates->erase(j1);
661  }
662  }
663  } while (candidates->size()>0);
664 
665  // sort jets by Energy
666  sort(jets.begin(), jets.end(), jets_E_less);
667 #ifdef DEBUG_SPLIT_MERGE
668  show();
669 #endif
670 
671  return jets.size();
672 }
673 
674 
675 
676 // save the event on disk
677 // - flux stream used to save jet contents
678 //--------------------------------------------
680  jet_iterator it_j;
681  CSphjet *j1;
682  int i1, i2;
683 
684  fprintf(flux, "# %d jets found\n", (int) jets.size());
685  fprintf(flux, "# columns are: px, py, pz, E and number of particles for each jet\n");
686  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
687  j1 = &(*it_j);
688  fprintf(flux, "%e\t%e\t%e\t%e\t%d\n",
689  j1->v.px, j1->v.py, j1->v.pz, j1->v.E, j1->n);
690  }
691 
692  fprintf(flux, "# jet contents\n");
693  fprintf(flux, "# columns are: px, py, pz, E, particle index and jet number\n");
694  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
695  j1 = &(*it_j);
696  for (i2=0;i2<j1->n;i2++)
697  fprintf(flux, "%e\t%e\t%e\t%e\t%d\t%d\n",
698  particles[j1->contents[i2]].px, particles[j1->contents[i2]].py,
699  particles[j1->contents[i2]].pz, particles[j1->contents[i2]].E,
700  j1->contents[i2], i1);
701  }
702 
703  return 0;
704 }
705 
706 
707 // show current jets/candidate status
708 //------------------------------------
710  jet_iterator it_j;
711  cjet_iterator it_c;
712  CSphjet *j;
713  const CSphjet *c;
714  int i1, i2;
715 
716  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
717  j = &(*it_j);
718  fprintf(stdout, "jet %2d: %e\t%e\t%e\t%e\t", i1+1,
719  j->v.px, j->v.py, j->v.pz, j->v.E);
720  for (i2=0;i2<j->n;i2++)
721  fprintf(stdout, "%d ", j->contents[i2]);
722  fprintf(stdout, "\n");
723  }
724 
725  for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
726  c = &(*it_c);
727  fprintf(stdout, "cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1,
728  c->v.px, c->v.py, c->v.pz, c->v.E, sqrt(c->sm_var2));
729  for (i2=0;i2<c->n;i2++)
730  fprintf(stdout, "%d ", c->contents[i2]);
731  fprintf(stdout, "\n");
732  }
733 
734  fprintf(stdout, "\n");
735  return 0;
736 }
737 
738 
739 // get the overlap between 2 jets
740 // - j1 first jet
741 // - j2 second jet
742 // - overlap2 returned overlap^2 (determined by the choice of SM variable)
743 // return true if overlapping, false if disjoint
744 //---------------------------------------------------------------------
745 bool CSphsplit_merge::get_overlap(const CSphjet &j1, const CSphjet &j2, double *overlap2){
746  // check if ranges overlap
747  if (!is_range_overlap(j1.range,j2.range))
748  return false;
749 
750  int i1,i2;
751  bool is_overlap;
752 
753  // initialise
754  i1=i2=idx_size=0;
755  is_overlap = false;
756  CSphmomentum v;
757 
758  // compute overlap
759  // at the same time, we store union in indices
760  do{
761  if (j1.contents[i1]<j2.contents[i2]){
762  indices[idx_size] = j1.contents[i1];
763  i1++;
764  } else if (j1.contents[i1]>j2.contents[i2]){
765  indices[idx_size] = j2.contents[i2];
766  i2++;
767  } else { // (j1.contents[i1]==j2.contents[i2])
768  v += particles[j2.contents[i2]];
769  indices[idx_size] = j2.contents[i2];
770  i1++;
771  i2++;
772  is_overlap = true;
773  }
774  idx_size++;
775  } while ((i1<j1.n) && (i2<j2.n));
776 
777  // finish computing union
778  // (only needed if overlap !)
779  if (is_overlap){
780  while (i1<j1.n){
781  indices[idx_size] = j1.contents[i1];
782  i1++;
783  idx_size++;
784  }
785  while (i2<j2.n){
786  indices[idx_size] = j2.contents[i2];
787  i2++;
788  idx_size++;
789  }
790  }
791 
792  // assign the overlapping var as return variable
793  (*overlap2) = sqr(v.E); //get_sm_var2(v, E_tilde);
794 
795  return is_overlap;
796 }
797 
798 
799 
800 // split the two given jet.
801 // during this procedure, the jets j1 & j2 are replaced
802 // by 2 new jets. Common particles are associted to the
803 // closest initial jet.
804 // - it_j1 iterator of the first jet in 'candidates'
805 // - it_j2 iterator of the second jet in 'candidates'
806 // - j1 first jet (CSphjet instance)
807 // - j2 second jet (CSphjet instance)
808 // return true on success, false on error
810 bool CSphsplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
811  int i1, i2;
812  CSphjet jet1, jet2;
813  double E1_weight, E2_weight;
814  CSphmomentum tmp;
815  CSphmomentum *v;
816 
817  // shorthand to avoid having "->" all over the place
818  const CSphjet & j1 = * it_j1;
819  const CSphjet & j2 = * it_j2;
820 
821  i1=i2=0;
822  jet2.v = jet1.v = CSphmomentum();
823 
824  // compute centroids
825  // When use_E_weighted_splitting is activated, the
826  // "geometrical" distance is weighted by the inverse
827  // of the E of the protojet
828  // This is stored in E{1,2}_weight
829  E1_weight = (use_E_weighted_splitting) ? 1.0/j1.v.E/j1.v.E : 1.0;
830  E2_weight = (use_E_weighted_splitting) ? 1.0/j2.v.E/j2.v.E : 1.0;
831 
832  // compute jet splitting
833  do{
834  if (j1.contents[i1]<j2.contents[i2]){
835  // particle i1 belong only to jet 1
836  v = &(particles[j1.contents[i1]]);
837  jet1.contents.push_back(j1.contents[i1]);
838  jet1.v += *v;
839  //jet1.pt_tilde += pt[j1.contents[i1]];
840  i1++;
841  jet1.range.add_particle(v->_theta,v->_phi);
842  } else if (j1.contents[i1]>j2.contents[i2]){
843  // particle i2 belong only to jet 2
844  v = &(particles[j2.contents[i2]]);
845  jet2.contents.push_back(j2.contents[i2]);
846  jet2.v += *v;
847  //jet2.pt_tilde += pt[j2.contents[i2]];
848  i2++;
849  jet2.range.add_particle(v->_theta,v->_phi);
850  } else { // (j1.contents[i1]==j2.contents[i2])
851  // common particle, decide which is the closest centre
852  v = &(particles[j1.contents[i1]]);
853 
854  //TODO: improve this brutal use of atan2 and sqrt !!!!
855 
856  //? what when == ?
857  // When use_E_weighted_splitting is activated, the
858  // "geometrical" distance is weighted by the inverse
859  // of the E of the protojet
860  double d1 = get_distance(&(j1.v), v)*E1_weight;
861  double d2 = get_distance(&(j2.v), v)*E2_weight;
862  // do bookkeeping on most ambiguous split
863  if (fabs(d1-d2) < most_ambiguous_split)
864  most_ambiguous_split = fabs(d1-d2);
865 
866  if (d1<d2){
867  // particle i1 belong only to jet 1
868  jet1.contents.push_back(j1.contents[i1]);
869  jet1.v += *v;
870  //jet1.pt_tilde += pt[j1.contents[i1]];
871  jet1.range.add_particle(v->_theta,v->_phi);
872  } else {
873  // particle i2 belong only to jet 2
874  jet2.contents.push_back(j2.contents[i2]);
875  jet2.v += *v;
876  //jet2.pt_tilde += pt[j2.contents[i2]];
877  jet2.range.add_particle(v->_theta,v->_phi);
878  }
879 
880  i1++;
881  i2++;
882  }
883  } while ((i1<j1.n) && (i2<j2.n));
884 
885  while (i1<j1.n){
886  v = &(particles[j1.contents[i1]]);
887  jet1.contents.push_back(j1.contents[i1]);
888  jet1.v += *v;
889  //jet1.pt_tilde += pt[j1.contents[i1]];
890  i1++;
891  jet1.range.add_particle(v->_theta,v->_phi);
892  }
893  while (i2<j2.n){
894  v = &(particles[j2.contents[i2]]);
895  jet2.contents.push_back(j2.contents[i2]);
896  jet2.v += *v;
897  //jet2.pt_tilde += pt[j2.contents[i2]];
898  i2++;
899  jet2.range.add_particle(v->_theta,v->_phi);
900  }
901 
902  // finalise jets
903  jet1.n = jet1.contents.size();
904  jet2.n = jet2.contents.size();
905 
906  // now the jet axis is known, we can compute Etilde
907  compute_Etilde(jet1);
908  compute_Etilde(jet2);
909 
910  // remove previous jets
911 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
912  cand_refs.erase(j1.v.ref);
913  cand_refs.erase(j2.v.ref);
914 #endif
915  candidates->erase(it_j1);
916  candidates->erase(it_j2);
917 
918  // reinsert new ones
919  insert(jet1);
920  insert(jet2);
921 
922  return true;
923 }
924 
925 // merge the two given jet.
926 // during this procedure, the jets j1 & j2 are replaced
927 // by 1 single jets containing both of them.
928 // - it_j1 iterator of the first jet in 'candidates'
929 // - it_j2 iterator of the second jet in 'candidates'
930 // return true on success, false on error
932 bool CSphsplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
933  CSphjet jet;
934  int i;
935 
936  // build new jet
937  // note: particles within j1 & j2 have already been stored in indices
938  for (i=0;i<idx_size;i++){
939  jet.contents.push_back(indices[i]);
940  jet.v += particles[indices[i]];
941  //jet.pt_tilde += pt[indices[i]];
942  }
943  jet.n = jet.contents.size();
944 
945  compute_Etilde(jet);
946 
947  // deal with ranges
948  jet.range = range_union(it_j1->range, it_j2->range);
949 
950  // remove old candidates
951 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
952  if (merge_identical_protocones){
953  cand_refs.erase(it_j1->v.ref);
954  cand_refs.erase(it_j2->v.ref);
955  }
956 #endif
957  candidates->erase(it_j1);
958  candidates->erase(it_j2);
959 
960  // reinsert new candidate
961  insert(jet);
962 
963  return true;
964 }
965 
971 bool CSphsplit_merge::insert(CSphjet &jet){
972 
973  // eventually check that no other candidate are present with the
974  // same cone contents. We recall that this automatic merging of
975  // identical protocones can lead to infrared-unsafe situations.
976 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
977  if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
978  return false;
979 #endif
980 
981  // check that the protojet has large enough energy
982  if (jet.v.E<E_min)
983  return false;
984 
985  // assign SM variable
986  jet.sm_var2 = get_sm_var2(jet.v, jet.E_tilde);
987 
988  // insert the jet.
989  candidates->insert(jet);
990 
991  return true;
992 }
993 
1000 double CSphsplit_merge::get_sm_var2(CSphmomentum &v, double &E_tilde){
1001  switch(ptcomparison.split_merge_scale) {
1002  case SM_E: return v.E*v.E;
1003  case SM_Etilde: return E_tilde*E_tilde;
1004  default:
1005  throw siscone::Csiscone_error("Unsupported split-merge scale choice: "
1006  + ptcomparison.SM_scale_name());
1007  }
1008 
1009  return 0.0;
1010 }
1011 
1012 
1013 
1015 void CSphsplit_merge::compute_Etilde(CSphjet &jet){
1016  jet.v.build_norm();
1017  jet.E_tilde=0.0;
1018  CSph3vector jet_axis = jet.v;
1019  //if (jet.v._norm==0){
1020  // jet_axis = CSph3vector(0.0,0.0,0.0);
1021  //} else {
1022  jet_axis/=jet.v.E;
1023  //}
1024  //cout << "~~~ Axis: " << jet.v.px << " " << jet.v.py << " " << jet.v.pz << " " << jet.v._norm << endl;
1025  //cout << "~~~ Axis: " << jet_axis.px << " " << jet_axis.py << " " << jet_axis.pz << endl;
1026  for (vector<int>::iterator cont_it=jet.contents.begin(); cont_it!=jet.contents.end(); cont_it++){
1027  const CSphmomentum &p = particles[*cont_it];
1028  jet.E_tilde+=p.E*(1.0+norm2_cross_product3(p,jet_axis)/particles_norm2[*cont_it]);
1029  }
1030 }
1031 
1032 }
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 n
number of particles inside
Definition: split_merge.h:63
double pz
z-momentum
Definition: momentum.h:133
bool operator()(const CSphjet &jet1, const CSphjet &jet2) const
comparison of 2 CSphjet
Definition: split_merge.cpp:94
double sm_var2
ordering variable used for ordering and overlap in the split–merge.
Definition: split_merge.h:77
int merge_collinear_and_remove_soft()
build the list 'p_uncol_hard' from p_remain by clustering collinear particles note that thins in only...
class corresponding to errors that will be thrown by siscone
Definition: siscone_error.h:38
double _phi
particle phi angle (available ONLY after a call to build_thetaphi)
Definition: momentum.h:137
#define EPSILON_SPLITMERGE
The following define enables you to allow for identical protocones to be merged automatically after e...
Definition: defines.h:111
siscone::Creference ref
reference number for the vector
Definition: momentum.h:142
int init_particles(std::vector< CSphmomentum > &_particles)
initialisation function for particle list
const double twopi
definition of 2*M_PI which is useful a bit everyhere!
Definition: defines.h:114
int save_contents(FILE *flux)
save final jets
CSphtheta_phi_range range
covered range in eta-phi
Definition: split_merge.h:80
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
double py
y-momentum
Definition: momentum.h:132
int add_protocones(std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
add a list of protocones
CSphmomentum v
jet momentum
Definition: split_merge.h:61
class for holding a covering range in eta-phi
Definition: geom_2d.h:52
double E_tilde
sum of E_i [ 1 +|p_i x p_J|^2/(|p_i|^2 E_J^2)]
Definition: split_merge.h:62
int parent_index
particle number in the parent list
Definition: momentum.h:210
int init(std::vector< CSphmomentum > &_particles, std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
initialisation function
std::vector< int > contents
particle contents (list of indices)
Definition: split_merge.h:64
double _theta
particle theta angle (available ONLY after a call to build_thetaphi)
Definition: momentum.h:136
int show()
show jets/candidates status
int index
internal particle number
Definition: momentum.h:211
int partial_clear()
partial clearance
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 init_pleft()
build initial list of left particles
base class for dynamic coordinates management
Definition: momentum.h:158
#define EPSILON_COLLINEAR
The following parameter controls collinear safety.
Definition: defines.h:75
base class for managing the spatial part of Cmomentum (defined after)
Definition: momentum.h:54
double px
x-momentum
Definition: momentum.h:131
real Jet information.
Definition: split_merge.h:53
The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated on Mon Oct 19 2015 20:25:39 for SISCone by  Doxygen 1.8.9.1