SISCone  2.0.6
split_merge.cpp
1 
2 // File: split_merge.cpp //
3 // Description: source file for splitting/merging (contains the CJet class) //
4 // This file is part of the SISCone project. //
5 // For more details, see http://projects.hepforge.org/siscone //
6 // //
7 // Copyright (c) 2006 Gavin Salam and Gregory Soyez //
8 // //
9 // This program is free software; you can redistribute it and/or modify //
10 // it under the terms of the GNU General Public License as published by //
11 // the Free Software Foundation; either version 2 of the License, or //
12 // (at your option) any later version. //
13 // //
14 // This program is distributed in the hope that it will be useful, //
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
17 // GNU General Public License for more details. //
18 // //
19 // You should have received a copy of the GNU General Public License //
20 // along with this program; if not, write to the Free Software //
21 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
22 // //
23 // $Revision:: 322 $//
24 // $Date:: 2011-11-15 10:12:36 +0100 (Tue, 15 Nov 2011) $//
26 
27 #include "split_merge.h"
28 #include "siscone_error.h"
29 #include "momentum.h"
30 #include <math.h>
31 #include <limits> // for max
32 #include <iostream>
33 #include <algorithm>
34 #include <sstream>
35 #include <cassert>
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 }
57 
58 // default dtor
59 //--------------
61 
62 }
63 
64 // ordering of jets in pt (e.g. used in final jets ordering)
65 //-----------------------------------------------------------
66 bool jets_pt_less(const Cjet &j1, const Cjet &j2){
67  return j1.v.perp2() > j2.v.perp2();
68 }
69 
70 
71 /********************************************************
72  * Csplit_merge_ptcomparison implementation *
73  * This deals with the ordering of the jets candidates *
74  ********************************************************/
75 
76 // odering of two jets
77 // The variable of the ordering is pt or mt
78 // depending on 'split_merge_scale' choice
79 //
80 // with EPSILON_SPLITMERGE defined, this attempts to identify
81 // delicate cases where two jets have identical momenta except for
82 // softish particles -- the difference of pt's may not be correctly
83 // identified normally and so lead to problems for the fate of the
84 // softish particle.
85 //
86 // NB: there is a potential issue in momentum-conserving events,
87 // whereby the harder of two jets becomes ill-defined when a soft
88 // particle is emitted --- this may have a knock-on effect on
89 // subsequent split-merge steps in cases with sufficiently large R
90 // (but we don't know what the limit is...)
91 //------------------------------------------------------------------
92 bool Csplit_merge_ptcomparison::operator ()(const Cjet &jet1, const Cjet &jet2) const{
93  double q1, q2;
94 
95  // compute the value for comparison for both jets
96  // This depends on the choice of variable (mt is the default)
97  q1 = jet1.sm_var2;
98  q2 = jet2.sm_var2;
99 
100  bool res = q1 > q2;
101 
102  // if we enable the refined version of the comparison (see defines.h),
103  // we compute the difference more precisely when the two jets are very
104  // close in the ordering variable.
105 #ifdef EPSILON_SPLITMERGE
106  if ( (fabs(q1-q2) < EPSILON_SPLITMERGE*max(q1,q2)) &&
107  (jet1.v.ref != jet2.v.ref) ) {
108  // get the momentum of the difference
109  Cmomentum difference;
110  double pt_tilde_difference;
111  get_difference(jet1,jet2,&difference,&pt_tilde_difference);
112 
113  // use the following relation: pt1^2 - pt2^2 = (pt1+pt2)*(pt1-pt2)
114  double qdiff;
115  Cmomentum sum = jet1.v ;
116  sum += jet2.v;
117  double pt_tilde_sum = jet1.pt_tilde + jet2.pt_tilde;
118 
119  // depending on the choice of ordering variable, set the result
120  switch (split_merge_scale){
121  case SM_mt:
122  qdiff = sum.E*difference.E - sum.pz*difference.pz;
123  break;
124  case SM_pt:
125  qdiff = sum.px*difference.px + sum.py*difference.py;
126  break;
127  case SM_pttilde:
128  qdiff = pt_tilde_sum*pt_tilde_difference;
129  break;
130  case SM_Et:
131  // diff = E^2 (dpt^2 pz^2- pt^2 dpz^2)
132  // + dE^2 (pt^2+pz^2) pt2^2
133  // where, unless explicitely specified the variables
134  // refer to the first jet or differences jet1-jet2.
135  qdiff = jet1.v.E*jet1.v.E*
136  ((sum.px*difference.px + sum.py*difference.py)*jet1.v.pz*jet1.v.pz
137  -jet1.v.perp2()*sum.pz*difference.pz)
138  +sum.E*difference.E*(jet1.v.perp2()+jet1.v.pz*jet1.v.pz)*jet2.v.perp2();
139  break;
140  default:
141  throw Csiscone_error("Unsupported split-merge scale choice: "
142  + SM_scale_name());
143  }
144  res = qdiff > 0;
145  }
146 #endif // EPSILON_SPLITMERGE
147 
148  return res;
149 }
150 
151 
154 std::string split_merge_scale_name(Esplit_merge_scale sms) {
155  switch(sms) {
156  case SM_pt:
157  return "pt (IR unsafe)";
158  case SM_Et:
159  return "Et (boost dep.)";
160  case SM_mt:
161  return "mt (IR safe except for pairs of identical decayed heavy particles)";
162  case SM_pttilde:
163  return "pttilde (scalar sum of pt's)";
164  default:
165  return "[SM scale without a name]";
166  }
167 }
168 
169 
170 // get the difference between 2 jets
171 // - j1 first jet
172 // - j2 second jet
173 // - v jet1-jet2
174 // - pt_tilde jet1-jet2 pt_tilde
175 // return true if overlapping, false if disjoint
176 //-----------------------------------------------
177 void Csplit_merge_ptcomparison::get_difference(const Cjet &j1, const Cjet &j2, Cmomentum *v, double *pt_tilde) const {
178  int i1,i2;
179 
180  // initialise
181  i1=i2=0;
182  *v = Cmomentum();
183  *pt_tilde = 0.0;
184 
185  // compute overlap
186  // at the same time, we store union in indices
187  do{
188  if (j1.contents[i1]==j2.contents[i2]) {
189  i1++;
190  i2++;
191  } else if (j1.contents[i1]<j2.contents[i2]){
192  (*v) += (*particles)[j1.contents[i1]];
193  (*pt_tilde) += (*pt)[j1.contents[i1]];
194  i1++;
195  } else if (j1.contents[i1]>j2.contents[i2]){
196  (*v) -= (*particles)[j2.contents[i2]];
197  (*pt_tilde) -= (*pt)[j2.contents[i2]];
198  i2++;
199  } else {
200  throw Csiscone_error("get_non_overlap reached part it should never have seen...");
201  }
202  } while ((i1<j1.n) && (i2<j2.n));
203 
204  // deal with particles at the end of the list...
205  while (i1 < j1.n) {
206  (*v) += (*particles)[j1.contents[i1]];
207  (*pt_tilde) += (*pt)[j1.contents[i1]];
208  i1++;
209  }
210  while (i2 < j2.n) {
211  (*v) -= (*particles)[j2.contents[i2]];
212  (*pt_tilde) -= (*pt)[j2.contents[i2]];
213  i2++;
214  }
215 }
216 
217 
218 /********************************************************
219  * class Csplit_merge implementation *
220  * Class used to split and merge jets. *
221  ********************************************************/
222 // default ctor
223 //--------------
225  merge_identical_protocones = false;
226 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
227 #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
228  merge_identical_protocones = true;
229 #endif
230 #endif
231  indices = NULL;
232 
233  // ensure that ptcomparison points to our set of particles (though params not correct)
234  ptcomparison.particles = &particles;
235  ptcomparison.pt = &pt;
236  candidates.reset(new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison));
237 
238  // no hardest cut (col-unsafe)
239  SM_var2_hardest_cut_off = -1.0;
240 
241  // no pt cutoff for the particles to put in p_uncol_hard
242  stable_cone_soft_pt2_cutoff = -1.0;
243 
244  // no pt-weighted splitting
245  use_pt_weighted_splitting = false;
246 }
247 
248 
249 // default dtor
250 //--------------
252  full_clear();
253 }
254 
255 
256 // initialisation function
257 // - _particles list of particles
258 // - protocones list of protocones (initial jet candidates)
259 // - R2 cone radius (squared)
260 // - ptmin minimal pT allowed for jets
261 //-------------------------------------------------------------
262 int Csplit_merge::init(vector<Cmomentum> & /*_particles*/, vector<Cmomentum> *protocones, double R2, double ptmin){
263  // browse protocones
264  return add_protocones(protocones, R2, ptmin);
265 }
266 
267 
268 // initialisation function for particle list
269 // - _particles list of particles
270 //-------------------------------------------------------------
271 int Csplit_merge::init_particles(vector<Cmomentum> &_particles){
272  full_clear();
273 
274  // compute the list of particles
275  // here, we do not need to throw away particles
276  // with infinite rapidity (colinear with the beam)
277  particles = _particles;
278  n = particles.size();
279 
280  // build the vector of particles' pt
281  pt.resize(n);
282  for (int i=0;i<n;i++)
283  pt[i] = particles[i].perp();
284 
285  // ensure that ptcomparison points to our set of particles (though params not correct)
286  ptcomparison.particles = &particles;
287  ptcomparison.pt = &pt;
288 
289  // set up the list of particles left.
290  init_pleft();
291 
292  indices = new int[n];
293 
294  return 0;
295 }
296 
297 
298 // build initial list of remaining particles
299 //------------------------------------------
301  // at this level, we only rule out particles with
302  // infinite rapidity
303  // in the parent particle list, index contain the run
304  // at which particles are puts in jets:
305  // - -1 means infinity rapidity
306  // - 0 means not included
307  // - i mean included at run i
308  int i,j;
309 
310  // copy particles removing the ones with infinite rapidity
311  // determine min,max eta
312  j=0;
313  double eta_min=0.0;
314  double eta_max=0.0;
315  p_remain.clear();
316  for (i=0;i<n;i++){
317  // set ref for checkxor
318  particles[i].ref.randomize();
319 
320  // check if rapidity is not infinite or ill-defined
321  if (fabs(particles[i].pz) < (particles[i].E)){
322  p_remain.push_back(particles[i]);
323  // set up parent index for tracability
324  p_remain[j].parent_index = i;
325  // left particles are marked with a 1
326  // IMPORTANT NOTE: the meaning of index in p_remain is
327  // somehow different as in the initial particle list.
328  // here, within one pass, we use the index to track whether
329  // a particle is included in the current pass (index set to 0
330  // in add_protocones) or still remain (index still 1)
331  p_remain[j].index = 1;
332 
333  j++;
334  // set up parent-particle index
335  particles[i].index = 0;
336 
337  eta_min = min(eta_min, particles[i].eta);
338  eta_max = max(eta_max, particles[i].eta);
339  } else {
340  particles[i].index = -1;
341  }
342  }
343  n_left = p_remain.size();
344  n_pass = 0;
345 
346  Ceta_phi_range epr;
347  epr.eta_min = eta_min-0.01;
348  epr.eta_max = eta_max+0.01;
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<Cjet,Csplit_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<Cmomentum> p_sorted;
406  bool collinear;
407  double dphi;
408 
409  p_uncol_hard.clear();
410 
411  // we first sort the particles according to their rapidity
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_eta_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_pt2_cutoff
423  if (p_sorted[i].perp2()<stable_cone_soft_pt2_cutoff) {
424  i++;
425  continue;
426  }
427 
428  // check if there is(are) particle(s) with the 'same' eta
429  collinear = false;
430  j=i+1;
431  while ((j<n_left) && (fabs(p_sorted[j].eta-p_sorted[i].eta)<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  p_sorted[j] += p_sorted[i];
437  // set collinearity test to true
438  collinear = true;
439  }
440  j++;
441  }
442  // if no collinearity detected, add the particle to our list
443  if (!collinear)
444  p_uncol_hard.push_back(p_sorted[i]);
445  i++;
446  }
447 
448  return 0;
449 }
450 
451 
452 // add a list of protocones
453 // - protocones list of protocones (initial jet candidates)
454 // - R2 cone radius (squared)
455 // - ptmin minimal pT allowed for jets
456 //-------------------------------------------------------------
457 int Csplit_merge::add_protocones(vector<Cmomentum> *protocones, double R2, double ptmin){
458  int i;
459  Cmomentum *c;
460  Cmomentum *v;
461  double eta, phi;
462  double dx, dy;
463  double R;
464  Cjet jet;
465 
466  if (protocones->size()==0)
467  return 1;
468 
469  pt_min2 = ptmin*ptmin;
470  R = sqrt(R2);
471 
472  // browse protocones
473  // for each of them, build the list of particles in them
474  for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
475  // initialise variables
476  c = &(*p_it);
477 
478  // note: cones have been tested => their (eta,phi) coordinates are computed
479  eta = c->eta;
480  phi = c->phi;
481 
482  // browse particles to create cone contents
483  // note that jet is always initialised with default values at this level
484  jet.v = Cmomentum();
485  jet.pt_tilde=0;
486  jet.contents.clear();
487  for (i=0;i<n_left;i++){
488  v = &(p_remain[i]);
489  // for security, avoid including particles with infinite rapidity)
490  // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING
491  //if (fabs(v->pz)!=v->E){
492  dx = eta - v->eta;
493  dy = fabs(phi - v->phi);
494  if (dy>M_PI)
495  dy -= twopi;
496  if (dx*dx+dy*dy<R2){
497  jet.contents.push_back(v->parent_index);
498  jet.v+= *v;
499  jet.pt_tilde+= pt[v->parent_index];
500  v->index=0;
501  }
502  }
503  jet.n=jet.contents.size();
504 
505  // set the momentum in protocones
506  // (it was only known through eta and phi up to now)
507  *c = jet.v;
508  c->eta = eta; // restore exact original coords
509  c->phi = phi; // to avoid rounding error inconsistencies
510 
511  // set the jet range
512  jet.range=Ceta_phi_range(eta,phi,R);
513 
514 #ifdef DEBUG_SPLIT_MERGE
515  cout << "adding jet: ";
516  for (int i2=0;i2<jet.n;i2++)
517  cout << jet.contents[i2] << " ";
518  cout << endl;
519 #endif
520 
521  // add it to the list of jets
522  insert(jet);
523  }
524 
525  // update list of included particles
526  n_pass++;
527 
528 #ifdef DEBUG_SPLIT_MERGE
529  cout << "remaining particles: ";
530 #endif
531  int j=0;
532  for (i=0;i<n_left;i++){
533  if (p_remain[i].index){
534  // copy particle
535  p_remain[j]=p_remain[i];
536  p_remain[j].parent_index = p_remain[i].parent_index;
537  p_remain[j].index=1;
538  // set run in initial list
539  particles[p_remain[j].parent_index].index = n_pass;
540 #ifdef DEBUG_SPLIT_MERGE
541  cout << p_remain[j].parent_index << " ";
542 #endif
543  j++;
544  }
545  }
546 #ifdef DEBUG_SPLIT_MERGE
547  cout << endl;
548 #endif
549  n_left = j;
550  p_remain.resize(j);
551 
552  merge_collinear_and_remove_soft();
553 
554  return 0;
555 }
556 
557 
558 /*
559  * really do the splitting and merging
560  * At the end, the vector jets is filled with the jets found.
561  * the 'contents' field of each jets contains the indices
562  * of the particles included in that jet.
563  * - overlap_tshold threshold for splitting/merging transition
564  * - ptmin minimal pT allowed for jets
565  * return the number of jets is returned
566  ******************************************************************/
567 int Csplit_merge::perform(double overlap_tshold, double ptmin){
568  // iterators for the 2 jets
569  cjet_iterator j1;
570  cjet_iterator j2;
571 
572  pt_min2 = ptmin*ptmin;
573 
574  if (candidates->size()==0)
575  return 0;
576 
577  if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
578  ostringstream message;
579  message << "Illegal value for overlap_tshold, f = " << overlap_tshold;
580  message << " (legal values are 0<f<1)";
581  throw Csiscone_error(message.str());
582  }
583 
584  // overlap (the contents of this variable depends on the choice for
585  // the split--merge variable.)
586  // Note that the square of the ovelap is used
587  double overlap2;
588 
589  // avoid to compute tshold*tshold at each overlap
590  double overlap_tshold2 = overlap_tshold*overlap_tshold;
591 
592  do{
593  if (candidates->size()>0){
594  // browse for the first jet
595  j1 = candidates->begin();
596 
597  // if hardest jet does not pass threshold then nothing else will
598  // either so one stops the split merge.
599  if (j1->sm_var2<SM_var2_hardest_cut_off) {break;}
600 
601  // browse for the second jet
602  j2 = j1;
603  j2++;
604  int j2_relindex = 1; // used only in ifdef, but costs little so keep it outside
605 
606  while (j2 != candidates->end()){
607 #ifdef DEBUG_SPLIT_MERGE
608  show();
609 #endif
610  // check overlapping
611  if (get_overlap(*j1, *j2, &overlap2)){
612  // check if overlapping energy passes threshold
613  // Note that this depends on the ordering variable
614 #ifdef DEBUG_SPLIT_MERGE
615  cout << "overlap between cdt 1 and cdt " << j2_relindex+1 << " with overlap "
616  << sqrt(overlap2/j2->sm_var2) << endl<<endl;
617 #endif
618  if (overlap2<overlap_tshold2*j2->sm_var2){
619  // split jets
620  split(j1, j2);
621 
622  // update iterators
623  j2 = j1 = candidates->begin();
624  j2_relindex = 0;
625  } else {
626  // merge jets
627  merge(j1, j2);
628 
629  // update iterators
630  j2 = j1 = candidates->begin();
631  j2_relindex = 0;
632  }
633  }
634  // watch out: split/merge might have caused new jets with pt <
635  // ptmin to disappear, so the total number of jets may
636  // have changed by more than expected and j2 might already by
637  // the end of the candidates list...
638  j2_relindex++;
639  if (j2 != candidates->end()) j2++;
640  } // end of loop on the second jet
641 
642  if (j1 != candidates->end()) {
643  // all "second jet" passed without overlapping
644  // (otherwise we won't leave the j2 loop)
645  // => store jet 1 as real jet
646  jets.push_back(*j1);
647  jets[jets.size()-1].v.build_etaphi();
648  // a bug where the contents has non-zero size has been cropping
649  // up in many contexts -- so catch it!
650  assert(j1->contents.size() > 0);
651  jets[jets.size()-1].pass = particles[j1->contents[0]].index;
652 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
653  cand_refs.erase(j1->v.ref);
654 #endif
655  candidates->erase(j1);
656 
658  //if ((candidates->size()!=0) &&
659  // (candidates->begin()->sm_var2<SM_var2_hardest_cut_off)){
660  // candidates->clear();
661  //}
662  }
663  }
664  } while (candidates->size()>0);
665 
666  // sort jets by pT
667  sort(jets.begin(), jets.end(), jets_pt_less);
668 #ifdef DEBUG_SPLIT_MERGE
669  show();
670 #endif
671 
672  return jets.size();
673 }
674 
675 
676 
677 // save the event on disk
678 // - flux stream used to save jet contents
679 //--------------------------------------------
681  jet_iterator it_j;
682  Cjet *j1;
683  int i1, i2;
684 
685  fprintf(flux, "# %d jets found\n", (int) jets.size());
686  fprintf(flux, "# columns are: eta, phi, pt and number of particles for each jet\n");
687  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
688  j1 = &(*it_j);
689  j1->v.build_etaphi();
690  fprintf(flux, "%f\t%f\t%e\t%d\n",
691  j1->v.eta, j1->v.phi, j1->v.perp(), j1->n);
692  }
693 
694  fprintf(flux, "# jet contents\n");
695  fprintf(flux, "# columns are: eta, phi, pt, particle index and jet number\n");
696  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
697  j1 = &(*it_j);
698  for (i2=0;i2<j1->n;i2++)
699  fprintf(flux, "%f\t%f\t%e\t%d\t%d\n",
700  particles[j1->contents[i2]].eta, particles[j1->contents[i2]].phi,
701  particles[j1->contents[i2]].perp(), j1->contents[i2], i1);
702  }
703 
704  return 0;
705 }
706 
707 
708 // show current jets/candidate status
709 //------------------------------------
711  jet_iterator it_j;
712  cjet_iterator it_c;
713  Cjet *j;
714  const Cjet *c;
715  int i1, i2;
716 
717  for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
718  j = &(*it_j);
719  fprintf(stdout, "jet %2d: %e\t%e\t%e\t%e\t", i1+1,
720  j->v.px, j->v.py, j->v.pz, j->v.E);
721  for (i2=0;i2<j->n;i2++)
722  fprintf(stdout, "%d ", j->contents[i2]);
723  fprintf(stdout, "\n");
724  }
725 
726  for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
727  c = &(*it_c);
728  fprintf(stdout, "cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1,
729  c->v.px, c->v.py, c->v.pz, c->v.E, sqrt(c->sm_var2));
730  for (i2=0;i2<c->n;i2++)
731  fprintf(stdout, "%d ", c->contents[i2]);
732  fprintf(stdout, "\n");
733  }
734 
735  fprintf(stdout, "\n");
736  return 0;
737 }
738 
739 
740 // get the overlap between 2 jets
741 // - j1 first jet
742 // - j2 second jet
743 // - overlap2 returned overlap^2 (determined by the choice of SM variable)
744 // return true if overlapping, false if disjoint
745 //---------------------------------------------------------------------
746 bool Csplit_merge::get_overlap(const Cjet &j1, const Cjet &j2, double *overlap2){
747  // check if ranges overlap
748  if (!is_range_overlap(j1.range,j2.range))
749  return false;
750 
751  int i1,i2;
752  bool is_overlap;
753 
754  // initialise
755  i1=i2=idx_size=0;
756  is_overlap = false;
757  Cmomentum v;
758  double pt_tilde=0.0;
759 
760  // compute overlap
761  // at the same time, we store union in indices
762  do{
763  if (j1.contents[i1]<j2.contents[i2]){
764  indices[idx_size] = j1.contents[i1];
765  i1++;
766  } else if (j1.contents[i1]>j2.contents[i2]){
767  indices[idx_size] = j2.contents[i2];
768  i2++;
769  } else { // (j1.contents[i1]==j2.contents[i2])
770  v += particles[j1.contents[i1]];
771  pt_tilde += pt[j1.contents[i1]];
772  indices[idx_size] = j1.contents[i1];
773  i1++;
774  i2++;
775  is_overlap = true;
776  }
777  idx_size++;
778  } while ((i1<j1.n) && (i2<j2.n));
779 
780  // finish computing union
781  // (only needed if overlap !)
782  if (is_overlap){
783  while (i1<j1.n){
784  indices[idx_size] = j1.contents[i1];
785  i1++;
786  idx_size++;
787  }
788  while (i2<j2.n){
789  indices[idx_size] = j2.contents[i2];
790  i2++;
791  idx_size++;
792  }
793  }
794 
795  // assign the overlapping var as return variable
796  (*overlap2) = get_sm_var2(v, pt_tilde);
797 
798  return is_overlap;
799 }
800 
801 
802 
803 // split the two given jet.
804 // during this procedure, the jets j1 & j2 are replaced
805 // by 2 new jets. Common particles are associted to the
806 // closest initial jet.
807 // - it_j1 iterator of the first jet in 'candidates'
808 // - it_j2 iterator of the second jet in 'candidates'
809 // - j1 first jet (Cjet instance)
810 // - j2 second jet (Cjet instance)
811 // return true on success, false on error
813 bool Csplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
814  int i1, i2;
815  Cjet jet1, jet2;
816  double eta1, phi1, pt1_weight, eta2, phi2, pt2_weight;
817  double dx1, dy1, dx2, dy2;
818  Cmomentum tmp;
819  Cmomentum *v;
820 
821  // shorthand to avoid having "->" all over the place
822  const Cjet & j1 = * it_j1;
823  const Cjet & j2 = * it_j2;
824 
825  i1=i2=0;
826  jet2.v = jet1.v = Cmomentum();
827  jet2.pt_tilde = jet1.pt_tilde = 0.0;
828 
829  // compute centroids
830  // When use_pt_weighted_splitting is activated, the
831  // "geometrical" distance is weighted by the inverse
832  // of the pt of the protojet
833  // This is stored in pt{1,2}_weight
834  tmp = j1.v;
835  tmp.build_etaphi();
836  eta1 = tmp.eta;
837  phi1 = tmp.phi;
838  pt1_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0;
839 
840  tmp = j2.v;
841  tmp.build_etaphi();
842  eta2 = tmp.eta;
843  phi2 = tmp.phi;
844  pt2_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0;
845 
846  jet1.v = jet2.v = Cmomentum();
847 
848  // compute jet splitting
849  do{
850  if (j1.contents[i1]<j2.contents[i2]){
851  // particle i1 belong only to jet 1
852  v = &(particles[j1.contents[i1]]);
853  jet1.contents.push_back(j1.contents[i1]);
854  jet1.v += *v;
855  jet1.pt_tilde += pt[j1.contents[i1]];
856  i1++;
857  jet1.range.add_particle(v->eta,v->phi);
858  } else if (j1.contents[i1]>j2.contents[i2]){
859  // particle i2 belong only to jet 2
860  v = &(particles[j2.contents[i2]]);
861  jet2.contents.push_back(j2.contents[i2]);
862  jet2.v += *v;
863  jet2.pt_tilde += pt[j2.contents[i2]];
864  i2++;
865  jet2.range.add_particle(v->eta,v->phi);
866  } else { // (j1.contents[i1]==j2.contents[i2])
867  // common particle, decide which is the closest centre
868  v = &(particles[j1.contents[i1]]);
869 
870  // distance w.r.t. centroid 1
871  dx1 = eta1 - v->eta;
872  dy1 = fabs(phi1 - v->phi);
873  if (dy1>M_PI)
874  dy1 -= twopi;
875 
876  // distance w.r.t. centroid 2
877  dx2 = eta2 - v->eta;
878  dy2 = fabs(phi2 - v->phi);
879  if (dy2>M_PI)
880  dy2 -= twopi;
881 
882  //? what when == ?
883  // When use_pt_weighted_splitting is activated, the
884  // "geometrical" distance is weighted by the inverse
885  // of the pt of the protojet
886  double d1sq = (dx1*dx1+dy1*dy1)*pt1_weight;
887  double d2sq = (dx2*dx2+dy2*dy2)*pt2_weight;
888  // do bookkeeping on most ambiguous split
889  if (fabs(d1sq-d2sq) < most_ambiguous_split)
890  most_ambiguous_split = fabs(d1sq-d2sq);
891 
892  if (d1sq<d2sq){
893  // particle i1 belong only to jet 1
894  jet1.contents.push_back(j1.contents[i1]);
895  jet1.v += *v;
896  jet1.pt_tilde += pt[j1.contents[i1]];
897  jet1.range.add_particle(v->eta,v->phi);
898  } else {
899  // particle i2 belong only to jet 2
900  jet2.contents.push_back(j2.contents[i2]);
901  jet2.v += *v;
902  jet2.pt_tilde += pt[j2.contents[i2]];
903  jet2.range.add_particle(v->eta,v->phi);
904  }
905 
906  i1++;
907  i2++;
908  }
909  } while ((i1<j1.n) && (i2<j2.n));
910 
911  while (i1<j1.n){
912  v = &(particles[j1.contents[i1]]);
913  jet1.contents.push_back(j1.contents[i1]);
914  jet1.v += *v;
915  jet1.pt_tilde += pt[j1.contents[i1]];
916  i1++;
917  jet1.range.add_particle(v->eta,v->phi);
918  }
919  while (i2<j2.n){
920  v = &(particles[j2.contents[i2]]);
921  jet2.contents.push_back(j2.contents[i2]);
922  jet2.v += *v;
923  jet2.pt_tilde += pt[j2.contents[i2]];
924  i2++;
925  jet2.range.add_particle(v->eta,v->phi);
926  }
927 
928  // finalise jets
929  jet1.n = jet1.contents.size();
930  jet2.n = jet2.contents.size();
931 
932  //jet1.range = j1.range;
933  //jet2.range = j2.range;
934 
935  // remove previous jets
936 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
937  cand_refs.erase(j1.v.ref);
938  cand_refs.erase(j2.v.ref);
939 #endif
940  candidates->erase(it_j1);
941  candidates->erase(it_j2);
942 
943  // reinsert new ones
944  insert(jet1);
945  insert(jet2);
946 
947  return true;
948 }
949 
950 // merge the two given jet.
951 // during this procedure, the jets j1 & j2 are replaced
952 // by 1 single jets containing both of them.
953 // - it_j1 iterator of the first jet in 'candidates'
954 // - it_j2 iterator of the second jet in 'candidates'
955 // return true on success, false on error
957 bool Csplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
958  Cjet jet;
959  int i;
960 
961  // build new jet
962  // note: particles within j1 & j2 have already been stored in indices
963  for (i=0;i<idx_size;i++){
964  jet.contents.push_back(indices[i]);
965  jet.v += particles[indices[i]];
966  jet.pt_tilde += pt[indices[i]];
967  }
968  jet.n = jet.contents.size();
969 
970  // deal with ranges
971  jet.range = range_union(it_j1->range, it_j2->range);
972 
973  // remove old candidates
974 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
975  if (merge_identical_protocones){
976  cand_refs.erase(it_j1->v.ref);
977  cand_refs.erase(it_j2->v.ref);
978  }
979 #endif
980  candidates->erase(it_j1);
981  candidates->erase(it_j2);
982 
983  // reinsert new candidate
984  insert(jet);
985 
986  return true;
987 }
988 
994 bool Csplit_merge::insert(Cjet &jet){
995 
996  // eventually check that no other candidate are present with the
997  // same cone contents. We recall that this automatic merging of
998  // identical protocones can lead to infrared-unsafe situations.
999 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1000  if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
1001  return false;
1002 #endif
1003 
1004  // check that the protojet has large enough pt
1005  if (jet.v.perp2()<pt_min2)
1006  return false;
1007 
1008  // assign SM variable
1009  jet.sm_var2 = get_sm_var2(jet.v, jet.pt_tilde);
1010 
1011  // insert the jet.
1012  candidates->insert(jet);
1013 
1014  return true;
1015 }
1016 
1023 double Csplit_merge::get_sm_var2(Cmomentum &v, double &pt_tilde){
1024  switch(ptcomparison.split_merge_scale) {
1025  case SM_pt: return v.perp2();
1026  case SM_mt: return v.perpmass2();
1027  case SM_pttilde: return pt_tilde*pt_tilde;
1028  case SM_Et: return v.Et2();
1029  default:
1030  throw Csiscone_error("Unsupported split-merge scale choice: "
1031  + ptcomparison.SM_scale_name());
1032  }
1033 
1034  return 0.0;
1035 }
1036 
1037 }
The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated on Mon May 6 2013 04:18:17 for SISCone by  Doxygen 1.8.1.2