source: HiSusy/trunk/Delphes/Delphes-3.0.9/external/fastjet/plugins/SISCone/split_merge.cc @ 5

Last change on this file since 5 was 5, checked in by zerwas, 11 years ago

update to Delphes-3.0.9

File size: 30.3 KB
Line 
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:: 859                                                          $//
24// $Date:: 2012-11-28 02:49:23 +0100 (Wed, 28 Nov 2012)                     $//
25///////////////////////////////////////////////////////////////////////////////
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
37namespace siscone{
38
39using namespace std;
40
41/********************************************************
42 * class Cjet implementation                            *
43 * real Jet information.                                *
44 * This class contains information for one single jet.  *
45 * That is, first, its momentum carrying information    *
46 * about its centre and pT, and second, its particle    *
47 * contents                                             *
48 ********************************************************/
49// default ctor
50//--------------
51Cjet::Cjet(){
52  n = 0;
53  v = Cmomentum();
54  pt_tilde = 0.0;
55  sm_var2 = 0.0;
56}
57
58// default dtor
59//--------------
60Cjet::~Cjet(){
61
62}
63
64// ordering of jets in pt (e.g. used in final jets ordering)
65//-----------------------------------------------------------
66bool 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//------------------------------------------------------------------
92bool 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
152/// return a name for the sm scale choice
153/// NB: used internally and also by fastjet
154std::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//-----------------------------------------------
177void 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//--------------
224Csplit_merge::Csplit_merge(){
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//--------------
251Csplit_merge::~Csplit_merge(){
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//-------------------------------------------------------------
262int 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//-------------------------------------------------------------
271int 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//------------------------------------------
300int Csplit_merge::init_pleft(){
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;  /// for the Ceta_phi_range static member!
314  double eta_max=0.0;  /// for the Ceta_phi_range static member!
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//----------------------------------------
361int Csplit_merge::partial_clear(){
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//----------------
386int Csplit_merge::full_clear(){
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//-------------------------------------------------------------------------
403int Csplit_merge::merge_collinear_and_remove_soft(){
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//-------------------------------------------------------------
457int 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 ******************************************************************/
567int 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
657        //// test that the hardest jet pass the potential cut-off
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//--------------------------------------------
680int Csplit_merge::save_contents(FILE *flux){
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//------------------------------------
710int Csplit_merge::show(){
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//---------------------------------------------------------------------
746bool 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
812////////////////////////////////////////////////////////////////
813bool 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
956////////////////////////////////////////////////////////////////
957bool 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
989/**
990 * Check whether or not a jet has to be inserted in the
991 * list of protojets. If it has, set its sm_variable and
992 * insert it to the list of protojets.
993 */
994bool 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
1017/**
1018 * given a 4-momentum and its associated pT, return the
1019 * variable that has to be used for SM
1020 * \param v          4 momentum of the protojet
1021 * \param pt_tilde   pt_tilde of the protojet
1022 */
1023double 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}
Note: See TracBrowser for help on using the repository browser.