Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/plugins/SISCone/src/split_merge.cc@ 333

Last change on this file since 333 was 11, checked in by severine ovyn, 16 years ago

Fastjet added; CDFCones directory has been changed

File size: 29.7 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: 1.1 $//
24// $Date: 2008-11-06 14:32:13 $//
25///////////////////////////////////////////////////////////////////////////////
26
27#include "../interface/split_merge.h"
28#include "../interface/siscone_error.h"
29#include "../interface/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
245
246// default dtor
247//--------------
248Csplit_merge::~Csplit_merge(){
249 full_clear();
250}
251
252
253// initialisation function
254// - _particles list of particles
255// - protocones list of protocones (initial jet candidates)
256// - R2 cone radius (squared)
257// - ptmin minimal pT allowed for jets
258//-------------------------------------------------------------
259int Csplit_merge::init(vector<Cmomentum> &_particles, vector<Cmomentum> *protocones, double R2, double ptmin){
260 // browse protocones
261 return add_protocones(protocones, R2, ptmin);
262}
263
264
265// initialisation function for particle list
266// - _particles list of particles
267//-------------------------------------------------------------
268int Csplit_merge::init_particles(vector<Cmomentum> &_particles){
269 full_clear();
270
271 // compute the list of particles
272 // here, we do not need to throw away particles
273 // with infinite rapidity (colinear with the beam)
274 particles = _particles;
275 n = particles.size();
276
277 // build the vector of particles' pt and determine min,max eta
278 double eta_min=0.0; /// for the Ceta_phi_range static member!
279 double eta_max=0.0; /// for the Ceta_phi_range static member!
280 pt.resize(n);
281 for (int i=0;i<n;i++){
282 pt[i] = particles[i].perp();
283 eta_min = min(eta_min, particles[i].eta);
284 eta_max = max(eta_max, particles[i].eta);
285 }
286 Ceta_phi_range epr;
287 epr.eta_min = eta_min-0.01;
288 epr.eta_max = eta_max+0.01;
289
290 // ensure that ptcomparison points to our set of particles (though params not correct)
291 ptcomparison.particles = &particles;
292 ptcomparison.pt = &pt;
293
294 // set up the list of particles left.
295 init_pleft();
296
297 indices = new int[n];
298
299 return 0;
300}
301
302
303// build initial list of remaining particles
304//------------------------------------------
305int Csplit_merge::init_pleft(){
306 // at this level, we only rule out particles with
307 // infinite rapidity
308 // in the parent particle list, index contain the run
309 // at which particles are puts in jets:
310 // - -1 means infinity rapidity
311 // - 0 means not included
312 // - i mean included at run i
313 int i,j;
314
315 // copy particles removing the ones with infinite rapidity
316 j=0;
317 p_remain.clear();
318 for (i=0;i<n;i++){
319 // set ref for checkxor
320 particles[i].ref.randomize();
321
322 // check if rapidity is not infinite or ill-defined
323 if (fabs(particles[i].pz) < (particles[i].E)){
324 p_remain.push_back(particles[i]);
325 // set up parent index for tracability
326 p_remain[j].parent_index = i;
327 // left particles are marked with a 1
328 // IMPORTANT NOTE: the meaning of index in p_remain is
329 // somehow different as in the initial particle list.
330 // here, within one pass, we use the index to track whether
331 // a particle is included in the current pass (index set to 0
332 // in add_protocones) or still remain (index still 1)
333 p_remain[j].index = 1;
334
335 j++;
336 // set up parent-particle index
337 particles[i].index = 0;
338 } else {
339 particles[i].index = -1;
340 }
341 }
342 n_left = p_remain.size();
343 n_pass = 0;
344
345 merge_collinear_and_remove_soft();
346
347 return 0;
348}
349
350
351// partial clearance
352// we want to keep particle list and indices
353// for future usage, so do not clear it !
354// this is done in full_clear
355//----------------------------------------
356int Csplit_merge::partial_clear(){
357 // release jets
358
359 // set up the auto_ptr for the multiset with the _current_ state of
360 // ptcomparison (which may have changed since we constructed the
361 // class)
362 candidates.reset(new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison));
363
364 // start off with huge number
365 most_ambiguous_split = numeric_limits<double>::max();
366
367 jets.clear();
368#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
369 if (merge_identical_protocones)
370 cand_refs.clear();
371#endif
372
373 p_remain.clear();
374
375 return 0;
376}
377
378
379// full clearance
380//----------------
381int Csplit_merge::full_clear(){
382 partial_clear();
383
384 // clear previously allocated memory
385 if (indices != NULL){
386 delete[] indices;
387 }
388 particles.clear();
389
390 return 0;
391}
392
393
394// build the list 'p_uncol_hard' from p_remain by clustering collinear particles
395// note that thins in only used for stable-cone detection
396// so the parent_index field is unnecessary
397//-------------------------------------------------------------------------
398int Csplit_merge::merge_collinear_and_remove_soft(){
399 int i,j;
400 vector<Cmomentum> p_sorted;
401 bool collinear;
402 double dphi;
403
404 p_uncol_hard.clear();
405
406 // we first sort the particles according to their rapidity
407 for (i=0;i<n_left;i++)
408 p_sorted.push_back(p_remain[i]);
409 sort(p_sorted.begin(), p_sorted.end(), momentum_eta_less);
410
411 // then we cluster particles looping over the particles in the following way
412 // if (a particle i has same eta-phi a one after (j))
413 // then add momentum i to j
414 // else add i to the p_uncol_hard list
415 i = 0;
416 while (i<n_left){
417 // check if the particle passes the stable_cone_soft_pt2_cutoff
418 if (p_sorted[i].perp2()<stable_cone_soft_pt2_cutoff) {
419 i++;
420 continue;
421 }
422
423 // check if there is(are) particle(s) with the 'same' eta
424 collinear = false;
425 j=i+1;
426 while ((j<n_left) && (fabs(p_sorted[j].eta-p_sorted[i].eta)<EPSILON_COLLINEAR) && (!collinear)){
427 dphi = fabs(p_sorted[j].phi-p_sorted[i].phi);
428 if (dphi>M_PI) dphi = twopi-dphi;
429 if (dphi<EPSILON_COLLINEAR){
430 // i is collinear with j; add the momentum (i) to the momentum (j)
431 p_sorted[j] += p_sorted[i];
432 // set collinearity test to true
433 collinear = true;
434 }
435 j++;
436 }
437 // if no collinearity detected, add the particle to our list
438 if (!collinear)
439 p_uncol_hard.push_back(p_sorted[i]);
440 i++;
441 }
442
443 return 0;
444}
445
446
447// add a list of protocones
448// - protocones list of protocones (initial jet candidates)
449// - R2 cone radius (squared)
450// - ptmin minimal pT allowed for jets
451//-------------------------------------------------------------
452int Csplit_merge::add_protocones(vector<Cmomentum> *protocones, double R2, double ptmin){
453 int i;
454 Cmomentum *c;
455 Cmomentum *v;
456 double eta, phi;
457 double dx, dy;
458 double R;
459 Cjet jet;
460
461 if (protocones->size()==0)
462 return 1;
463
464 pt_min2 = ptmin*ptmin;
465 R = sqrt(R2);
466
467 // browse protocones
468 // for each of them, build the list of particles in them
469 for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
470 // initialise variables
471 c = &(*p_it);
472
473 // note: cones have been tested => their (eta,phi) coordinates are computed
474 eta = c->eta;
475 phi = c->phi;
476
477 // browse particles to create cone contents
478 // note that jet is always initialised with default values at this level
479 jet.v = Cmomentum();
480 jet.pt_tilde=0;
481 jet.contents.clear();
482 for (i=0;i<n_left;i++){
483 v = &(p_remain[i]);
484 // for security, avoid including particles with infinite rapidity)
485 // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING
486 //if (fabs(v->pz)!=v->E){
487 dx = eta - v->eta;
488 dy = fabs(phi - v->phi);
489 if (dy>M_PI)
490 dy -= twopi;
491 if (dx*dx+dy*dy<R2){
492 jet.contents.push_back(v->parent_index);
493 jet.v+= *v;
494 jet.pt_tilde+= pt[v->parent_index];
495 v->index=0;
496 }
497 }
498 jet.n=jet.contents.size();
499
500 // set the momentum in protocones
501 // (it was only known through eta and phi up to now)
502 *c = jet.v;
503 c->eta = eta; // restore exact original coords
504 c->phi = phi; // to avoid rounding error inconsistencies
505
506 // set the jet range
507 jet.range=Ceta_phi_range(eta,phi,R);
508
509#ifdef DEBUG_SPLIT_MERGE
510 cout << "adding jet: ";
511 for (int i2=0;i2<jet.n;i2++)
512 cout << jet.contents[i2] << " ";
513 cout << endl;
514#endif
515
516 // add it to the list of jets
517 insert(jet);
518 }
519
520 // update list of included particles
521 n_pass++;
522
523#ifdef DEBUG_SPLIT_MERGE
524 cout << "remaining particles: ";
525#endif
526 int j=0;
527 for (i=0;i<n_left;i++){
528 if (p_remain[i].index){
529 // copy particle
530 p_remain[j]=p_remain[i];
531 p_remain[j].parent_index = p_remain[i].parent_index;
532 p_remain[j].index=1;
533 // set run in initial list
534 particles[p_remain[j].parent_index].index = n_pass;
535#ifdef DEBUG_SPLIT_MERGE
536 cout << p_remain[j].parent_index << " ";
537#endif
538 j++;
539 }
540 }
541#ifdef DEBUG_SPLIT_MERGE
542 cout << endl;
543#endif
544 n_left = j;
545 p_remain.resize(j);
546
547 merge_collinear_and_remove_soft();
548
549 return 0;
550}
551
552
553/*
554 * really do the splitting and merging
555 * At the end, the vector jets is filled with the jets found.
556 * the 'contents' field of each jets contains the indices
557 * of the particles included in that jet.
558 * - overlap_tshold threshold for splitting/merging transition
559 * - ptmin minimal pT allowed for jets
560 * return the number of jets is returned
561 ******************************************************************/
562int Csplit_merge::perform(double overlap_tshold, double ptmin){
563 // iterators for the 2 jets
564 cjet_iterator j1;
565 cjet_iterator j2;
566
567 pt_min2 = ptmin*ptmin;
568
569 if (candidates->size()==0)
570 return 0;
571
572 if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
573 ostringstream message;
574 message << "Illegal value for overlap_tshold, f = " << overlap_tshold;
575 message << " (legal values are 0<f<1)";
576 throw Csiscone_error(message.str());
577 }
578
579 // overlap (the contents of this variable depends on the choice for
580 // the split--merge variable.)
581 // Note that the square of the ovelap is used
582 double overlap2;
583
584 // avoid to compute tshold*tshold at each overlap
585 double overlap_tshold2 = overlap_tshold*overlap_tshold;
586
587 do{
588 if (candidates->size()>0){
589 // browse for the first jet
590 j1 = candidates->begin();
591
592 // if hardest jet does not pass threshold then nothing else will
593 // either so one stops the split merge.
594 if (j1->sm_var2<SM_var2_hardest_cut_off) {break;}
595
596 // browse for the second jet
597 j2 = j1;
598 j2++;
599 int j2_relindex = 1; // used only in ifdef, but costs little so keep it outside
600
601 while (j2 != candidates->end()){
602#ifdef DEBUG_SPLIT_MERGE
603 show();
604#endif
605 // check overlapping
606 if (get_overlap(*j1, *j2, &overlap2)){
607 // check if overlapping energy passes threshold
608 // Note that this depends on the ordering variable
609#ifdef DEBUG_SPLIT_MERGE
610 cout << "overlap between cdt 1 and cdt " << j2_relindex+1 << " with overlap "
611 << sqrt(overlap2/j2->sm_var2) << endl<<endl;
612#endif
613 if (overlap2<overlap_tshold2*j2->sm_var2){
614 // split jets
615 split(j1, j2);
616
617 // update iterators
618 j2 = j1 = candidates->begin();
619 j2_relindex = 0;
620 } else {
621 // merge jets
622 merge(j1, j2);
623
624 // update iterators
625 j2 = j1 = candidates->begin();
626 j2_relindex = 0;
627 }
628 }
629 // watch out: split/merge might have caused new jets with pt <
630 // ptmin to disappear, so the total number of jets may
631 // have changed by more than expected and j2 might already by
632 // the end of the candidates list...
633 j2_relindex++;
634 if (j2 != candidates->end()) j2++;
635 } // end of loop on the second jet
636
637 if (j1 != candidates->end()) {
638 // all "second jet" passed without overlapping
639 // (otherwise we won't leave the j2 loop)
640 // => store jet 1 as real jet
641 jets.push_back(*j1);
642 jets[jets.size()-1].v.build_etaphi();
643 // a bug where the contents has non-zero size has been cropping
644 // up in many contexts -- so catch it!
645 assert(j1->contents.size() > 0);
646 jets[jets.size()-1].pass = particles[j1->contents[0]].index;
647#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
648 cand_refs.erase(j1->v.ref);
649#endif
650 candidates->erase(j1);
651
652 //// test that the hardest jet pass the potential cut-off
653 //if ((candidates->size()!=0) &&
654 // (candidates->begin()->sm_var2<SM_var2_hardest_cut_off)){
655 // candidates->clear();
656 //}
657 }
658 }
659 } while (candidates->size()>0);
660
661 // sort jets by pT
662 sort(jets.begin(), jets.end(), jets_pt_less);
663#ifdef DEBUG_SPLIT_MERGE
664 show();
665#endif
666
667 return jets.size();
668}
669
670
671
672// save the event on disk
673// - flux stream used to save jet contents
674//--------------------------------------------
675int Csplit_merge::save_contents(FILE *flux){
676 jet_iterator it_j;
677 Cjet *j1;
678 int i1, i2;
679
680 fprintf(flux, "# %d jets found\n", (int) jets.size());
681 fprintf(flux, "# columns are: eta, phi, pt and number of particles for each jet\n");
682 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
683 j1 = &(*it_j);
684 j1->v.build_etaphi();
685 fprintf(flux, "%lf\t%lf\t%le\t%d\n",
686 j1->v.eta, j1->v.phi, j1->v.perp(), j1->n);
687 }
688
689 fprintf(flux, "# jet contents\n");
690 fprintf(flux, "# columns are: eta, phi, pt, particle index and jet number\n");
691 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
692 j1 = &(*it_j);
693 for (i2=0;i2<j1->n;i2++)
694 fprintf(flux, "%lf\t%lf\t%le\t%d\t%d\n",
695 particles[j1->contents[i2]].eta, particles[j1->contents[i2]].phi,
696 particles[j1->contents[i2]].perp(), j1->contents[i2], i1);
697 }
698
699 return 0;
700}
701
702
703// show current jets/candidate status
704//------------------------------------
705int Csplit_merge::show(){
706 jet_iterator it_j;
707 cjet_iterator it_c;
708 Cjet *j;
709 const Cjet *c;
710 int i1, i2;
711
712 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
713 j = &(*it_j);
714 fprintf(stdout, "jet %2d: %le\t%le\t%le\t%le\t", i1+1,
715 j->v.px, j->v.py, j->v.pz, j->v.E);
716 for (i2=0;i2<j->n;i2++)
717 fprintf(stdout, "%d ", j->contents[i2]);
718 fprintf(stdout, "\n");
719 }
720
721 for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
722 c = &(*it_c);
723 fprintf(stdout, "cdt %2d: %le\t%le\t%le\t%le\t%le\t", i1+1,
724 c->v.px, c->v.py, c->v.pz, c->v.E, sqrt(c->sm_var2));
725 for (i2=0;i2<c->n;i2++)
726 fprintf(stdout, "%d ", c->contents[i2]);
727 fprintf(stdout, "\n");
728 }
729
730 fprintf(stdout, "\n");
731 return 0;
732}
733
734
735// get the overlap between 2 jets
736// - j1 first jet
737// - j2 second jet
738// - overlap2 returned overlap^2 (determined by the choice of SM variable)
739// return true if overlapping, false if disjoint
740//---------------------------------------------------------------------
741bool Csplit_merge::get_overlap(const Cjet &j1, const Cjet &j2, double *overlap2){
742 // check if ranges overlap
743 if (!is_range_overlap(j1.range,j2.range))
744 return false;
745
746 int i1,i2;
747 bool is_overlap;
748
749 // initialise
750 i1=i2=idx_size=0;
751 is_overlap = false;
752 Cmomentum v;
753 double pt_tilde=0.0;
754
755 // compute overlap
756 // at the same time, we store union in indices
757 do{
758 if (j1.contents[i1]<j2.contents[i2]){
759 indices[idx_size] = j1.contents[i1];
760 i1++;
761 } else if (j1.contents[i1]>j2.contents[i2]){
762 indices[idx_size] = j2.contents[i2];
763 i2++;
764 } else { // (j1.contents[i1]==j2.contents[i2])
765 v += particles[j1.contents[i1]];
766 pt_tilde += pt[j1.contents[i1]];
767 indices[idx_size] = j1.contents[i1];
768 i1++;
769 i2++;
770 is_overlap = true;
771 }
772 idx_size++;
773 } while ((i1<j1.n) && (i2<j2.n));
774
775 // finish computing union
776 // (only needed if overlap !)
777 if (is_overlap){
778 while (i1<j1.n){
779 indices[idx_size] = j1.contents[i1];
780 i1++;
781 idx_size++;
782 }
783 while (i2<j2.n){
784 indices[idx_size] = j2.contents[i2];
785 i2++;
786 idx_size++;
787 }
788 }
789
790 // assign the overlapping var as return variable
791 (*overlap2) = get_sm_var2(v, pt_tilde);
792
793 return is_overlap;
794}
795
796
797
798// split the two given jet.
799// during this procedure, the jets j1 & j2 are replaced
800// by 2 new jets. Common particles are associted to the
801// closest initial jet.
802// - it_j1 iterator of the first jet in 'candidates'
803// - it_j2 iterator of the second jet in 'candidates'
804// - j1 first jet (Cjet instance)
805// - j2 second jet (Cjet instance)
806// return true on success, false on error
807////////////////////////////////////////////////////////////////
808bool Csplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
809 int i1, i2;
810 Cjet jet1, jet2;
811 double eta1, phi1, eta2, phi2;
812 double dx1, dy1, dx2, dy2;
813 Cmomentum tmp;
814 Cmomentum *v;
815
816 // shorthand to avoid having "->" all over the place
817 const Cjet & j1 = * it_j1;
818 const Cjet & j2 = * it_j2;
819
820 i1=i2=0;
821 jet2.v = jet1.v = Cmomentum();
822 jet2.pt_tilde = jet1.pt_tilde = 0.0;
823
824 // compute centroids
825 tmp = j1.v;
826 tmp.build_etaphi();
827 eta1 = tmp.eta;
828 phi1 = tmp.phi;
829
830 tmp = j2.v;
831 tmp.build_etaphi();
832 eta2 = tmp.eta;
833 phi2 = tmp.phi;
834
835 jet1.v = jet2.v = Cmomentum();
836
837 // compute jet splitting
838 do{
839 if (j1.contents[i1]<j2.contents[i2]){
840 // particle i1 belong only to jet 1
841 v = &(particles[j1.contents[i1]]);
842 jet1.contents.push_back(j1.contents[i1]);
843 jet1.v += *v;
844 jet1.pt_tilde += pt[j1.contents[i1]];
845 i1++;
846 jet1.range.add_particle(v->eta,v->phi);
847 } else if (j1.contents[i1]>j2.contents[i2]){
848 // particle i2 belong only to jet 2
849 v = &(particles[j2.contents[i2]]);
850 jet2.contents.push_back(j2.contents[i2]);
851 jet2.v += *v;
852 jet2.pt_tilde += pt[j2.contents[i2]];
853 i2++;
854 jet2.range.add_particle(v->eta,v->phi);
855 } else { // (j1.contents[i1]==j2.contents[i2])
856 // common particle, decide which is the closest centre
857 v = &(particles[j1.contents[i1]]);
858
859 // distance w.r.t. centroid 1
860 dx1 = eta1 - v->eta;
861 dy1 = fabs(phi1 - v->phi);
862 if (dy1>M_PI)
863 dy1 -= twopi;
864
865 // distance w.r.t. centroid 2
866 dx2 = eta2 - v->eta;
867 dy2 = fabs(phi2 - v->phi);
868 if (dy2>M_PI)
869 dy2 -= twopi;
870
871 //? what when == ?
872 double d1sq = dx1*dx1+dy1*dy1;
873 double d2sq = dx2*dx2+dy2*dy2;
874 // do bookkeeping on most ambiguous split
875 if (fabs(d1sq-d2sq) < most_ambiguous_split)
876 most_ambiguous_split = fabs(d1sq-d2sq);
877
878 if (d1sq<d2sq){
879 // particle i1 belong only to jet 1
880 jet1.contents.push_back(j1.contents[i1]);
881 jet1.v += *v;
882 jet1.pt_tilde += pt[j1.contents[i1]];
883 jet1.range.add_particle(v->eta,v->phi);
884 } else {
885 // particle i2 belong only to jet 2
886 jet2.contents.push_back(j2.contents[i2]);
887 jet2.v += *v;
888 jet2.pt_tilde += pt[j2.contents[i2]];
889 jet2.range.add_particle(v->eta,v->phi);
890 }
891
892 i1++;
893 i2++;
894 }
895 } while ((i1<j1.n) && (i2<j2.n));
896
897 while (i1<j1.n){
898 v = &(particles[j1.contents[i1]]);
899 jet1.contents.push_back(j1.contents[i1]);
900 jet1.v += *v;
901 jet1.pt_tilde += pt[j1.contents[i1]];
902 i1++;
903 jet1.range.add_particle(v->eta,v->phi);
904 }
905 while (i2<j2.n){
906 v = &(particles[j2.contents[i2]]);
907 jet2.contents.push_back(j2.contents[i2]);
908 jet2.v += *v;
909 jet2.pt_tilde += pt[j2.contents[i2]];
910 i2++;
911 jet2.range.add_particle(v->eta,v->phi);
912 }
913
914 // finalise jets
915 jet1.n = jet1.contents.size();
916 jet2.n = jet2.contents.size();
917
918 //jet1.range = j1.range;
919 //jet2.range = j2.range;
920
921 // remove previous jets
922#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
923 cand_refs.erase(j1.v.ref);
924 cand_refs.erase(j2.v.ref);
925#endif
926 candidates->erase(it_j1);
927 candidates->erase(it_j2);
928
929 // reinsert new ones
930 insert(jet1);
931 insert(jet2);
932
933 return true;
934}
935
936// merge the two given jet.
937// during this procedure, the jets j1 & j2 are replaced
938// by 1 single jets containing both of them.
939// - it_j1 iterator of the first jet in 'candidates'
940// - it_j2 iterator of the second jet in 'candidates'
941// return true on success, false on error
942////////////////////////////////////////////////////////////////
943bool Csplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
944 Cjet jet;
945 int i;
946
947 // build new jet
948 // note: particles within j1 & j2 have already been stored in indices
949 for (i=0;i<idx_size;i++){
950 jet.contents.push_back(indices[i]);
951 jet.v += particles[indices[i]];
952 jet.pt_tilde += pt[indices[i]];
953 }
954 jet.n = jet.contents.size();
955
956 // deal with ranges
957 jet.range = range_union(it_j1->range, it_j2->range);
958
959 // remove old candidates
960#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
961 if (merge_identical_protocones){
962 cand_refs.erase(it_j1->v.ref);
963 cand_refs.erase(it_j2->v.ref);
964 }
965#endif
966 candidates->erase(it_j1);
967 candidates->erase(it_j2);
968
969 // reinsert new candidate
970 insert(jet);
971
972 return true;
973}
974
975/**
976 * Check whether or not a jet has to be inserted in the
977 * list of protojets. If it has, set its sm_variable and
978 * insert it to the list of protojets.
979 */
980bool Csplit_merge::insert(Cjet &jet){
981
982 // eventually check that no other candidate are present with the
983 // same cone contents. We recall that this automatic merging of
984 // identical protocones can lead to infrared-unsafe situations.
985#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
986 if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
987 return false;
988#endif
989
990 // check that the protojet has large enough pt
991 if (jet.v.perp2()<pt_min2)
992 return false;
993
994 // assign SM variable
995 jet.sm_var2 = get_sm_var2(jet.v, jet.pt_tilde);
996
997 // insert the jet.
998 candidates->insert(jet);
999
1000 return true;
1001}
1002
1003/**
1004 * given a 4-momentum and its associated pT, return the
1005 * variable that has to be used for SM
1006 * \param v 4 momentum of the protojet
1007 * \param pt_tilde pt_tilde of the protojet
1008 */
1009double Csplit_merge::get_sm_var2(Cmomentum &v, double &pt_tilde){
1010 switch(ptcomparison.split_merge_scale) {
1011 case SM_pt: return v.perp2();
1012 case SM_mt: return v.perpmass2();
1013 case SM_pttilde: return pt_tilde*pt_tilde;
1014 case SM_Et: return v.Et2();
1015 default:
1016 throw Csiscone_error("Unsupported split-merge scale choice: "
1017 + ptcomparison.SM_scale_name());
1018 }
1019
1020 return 0.0;
1021}
1022
1023}
Note: See TracBrowser for help on using the repository browser.