Changeset 273e668 in git for external/fastjet/plugins/SISCone/split_merge.cc
- Timestamp:
- Oct 15, 2014, 10:55:55 AM (10 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 35b9204, b25d4cf
- Parents:
- f14bd6a
- git-author:
- Pavel Demin <pavel.demin@…> (10/10/14 08:56:40)
- git-committer:
- Pavel Demin <pavel.demin@…> (10/15/14 10:55:55)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
external/fastjet/plugins/SISCone/split_merge.cc
rf14bd6a r273e668 21 21 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // 22 22 // // 23 // $Revision:: 3 22$//24 // $Date:: 201 1-11-15 10:12:36 +0100 (Tue, 15 Nov 2011) $//23 // $Revision:: 370 $// 24 // $Date:: 2014-09-04 17:03:15 +0200 (Thu, 04 Sep 2014) $// 25 25 /////////////////////////////////////////////////////////////////////////////// 26 26 … … 28 28 #include "siscone_error.h" 29 29 #include "momentum.h" 30 #include <math.h>31 30 #include <limits> // for max 32 31 #include <iostream> … … 34 33 #include <sstream> 35 34 #include <cassert> 35 #include <cmath> 36 36 37 37 namespace siscone{ … … 229 229 #endif 230 230 #endif 231 _user_scale = NULL; 231 232 indices = NULL; 232 233 … … 237 238 238 239 // no hardest cut (col-unsafe) 239 SM_var2_hardest_cut_off = - 1.0;240 SM_var2_hardest_cut_off = -numeric_limits<double>::max(); 240 241 241 242 // no pt cutoff for the particles to put in p_uncol_hard … … 555 556 } 556 557 558 559 /* 560 * remove the hardest protocone and declare it a a jet 561 * - protocones list of protocones (initial jet candidates) 562 * - R2 cone radius (squared) 563 * - ptmin minimal pT allowed for jets 564 * return 0 on success, 1 on error 565 * 566 * The list of remaining particles (and the uncollinear-hard ones) 567 * is updated. 568 */ 569 int Csplit_merge::add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin){ 570 571 int i; 572 Cmomentum *c; 573 Cmomentum *v; 574 double eta, phi; 575 double dx, dy; 576 double R; 577 Cjet jet, jet_candidate; 578 bool found_jet = false; 579 580 if (protocones->size()==0) 581 return 1; 582 583 pt_min2 = ptmin*ptmin; 584 R = sqrt(R2); 585 586 // browse protocones 587 // for each of them, build the list of particles in them 588 for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){ 589 // initialise variables 590 c = &(*p_it); 591 592 // note: cones have been tested => their (eta,phi) coordinates are computed 593 eta = c->eta; 594 phi = c->phi; 595 596 // NOTE: this is common to this method and add_protocones, so it 597 // could be moved into a 'build_jet_from_protocone' method 598 // 599 // browse particles to create cone contents 600 jet_candidate.v = Cmomentum(); 601 jet_candidate.pt_tilde=0; 602 jet_candidate.contents.clear(); 603 for (i=0;i<n_left;i++){ 604 v = &(p_remain[i]); 605 // for security, avoid including particles with infinite rapidity) 606 // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING 607 //if (fabs(v->pz)!=v->E){ 608 dx = eta - v->eta; 609 dy = fabs(phi - v->phi); 610 if (dy>M_PI) 611 dy -= twopi; 612 if (dx*dx+dy*dy<R2){ 613 jet_candidate.contents.push_back(v->parent_index); 614 jet_candidate.v+= *v; 615 jet_candidate.pt_tilde+= pt[v->parent_index]; 616 v->index=0; 617 } 618 } 619 jet_candidate.n=jet_candidate.contents.size(); 620 621 // set the momentum in protocones 622 // (it was only known through eta and phi up to now) 623 *c = jet_candidate.v; 624 c->eta = eta; // restore exact original coords 625 c->phi = phi; // to avoid rounding error inconsistencies 626 627 // set the jet range 628 jet_candidate.range=Ceta_phi_range(eta,phi,R); 629 630 // check that the protojet has large enough pt 631 if (jet_candidate.v.perp2()<pt_min2) 632 continue; 633 634 // assign the split-merge (or progressive-removal) squared scale variable 635 if (_user_scale) { 636 // sm_var2 is the signed square of the user scale returned 637 // for the jet candidate 638 jet_candidate.sm_var2 = (*_user_scale)(jet_candidate); 639 jet_candidate.sm_var2 *= abs(jet_candidate.sm_var2); 640 } else { 641 jet_candidate.sm_var2 = get_sm_var2(jet_candidate.v, jet_candidate.pt_tilde); 642 } 643 644 // now check if it is possibly the hardest 645 if ((! found_jet) || 646 (_user_scale ? _user_scale->is_larger(jet_candidate, jet) 647 : ptcomparison(jet_candidate, jet))){ 648 jet = jet_candidate; 649 found_jet = true; 650 } 651 } 652 653 // make sure at least one of the jets has passed the selection 654 if (!found_jet) return 1; 655 656 // add the jet to the list of jets 657 jets.push_back(jet); 658 jets[jets.size()-1].v.build_etaphi(); 659 660 #ifdef DEBUG_SPLIT_MERGE 661 cout << "PR-Jet " << jets.size() << " [size " << next_jet.contents.size() << "]:"; 662 #endif 663 664 // update the list of what particles are left 665 int p_remain_index = 0; 666 int contents_index = 0; 667 //sort(next_jet.contents.begin(),next_jet.contents.end()); 668 for (int index=0;index<n_left;index++){ 669 if ((contents_index<(int) jet.contents.size()) && 670 (p_remain[index].parent_index == jet.contents[contents_index])){ 671 // this particle belongs to the newly found jet 672 // set pass in initial list 673 particles[p_remain[index].parent_index].index = n_pass; 674 #ifdef DEBUG_SPLIT_MERGE 675 cout << " " << jet.contents[contents_index]; 676 #endif 677 contents_index++; 678 } else { 679 // this particle still has to be clustered 680 p_remain[p_remain_index] = p_remain[index]; 681 p_remain[p_remain_index].parent_index = p_remain[index].parent_index; 682 p_remain[p_remain_index].index=1; 683 p_remain_index++; 684 } 685 } 686 p_remain.resize(n_left-jet.contents.size()); 687 n_left = p_remain.size(); 688 jets[jets.size()-1].pass = particles[jet.contents[0]].index; 689 690 // increase the pass index 691 n_pass++; 692 693 #ifdef DEBUG_SPLIT_MERGE 694 cout << endl; 695 #endif 696 697 // male sure the list of uncol_hard particles (used for the next 698 // stable cone finding) is updated [could probably be more 699 // efficient] 700 merge_collinear_and_remove_soft(); 701 702 return 0; 703 } 557 704 558 705 /*
Note:
See TracChangeset
for help on using the changeset viewer.