Fork me on GitHub

Ignore:
Timestamp:
Oct 15, 2014, 10:55:55 AM (10 years ago)
Author:
Pavel Demin <pavel.demin@…>
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)
Message:

upgrade FastJet to version 3.1.0

File:
1 edited

Legend:

Unmodified
Added
Removed
  • external/fastjet/plugins/SISCone/split_merge.cc

    rf14bd6a r273e668  
    2121// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
    2222//                                                                           //
    23 // $Revision:: 322                                                          $//
    24 // $Date:: 2011-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)                     $//
    2525///////////////////////////////////////////////////////////////////////////////
    2626
     
    2828#include "siscone_error.h"
    2929#include "momentum.h"
    30 #include <math.h>
    3130#include <limits>   // for max
    3231#include <iostream>
     
    3433#include <sstream>
    3534#include <cassert>
     35#include <cmath>
    3636
    3737namespace siscone{
     
    229229#endif
    230230#endif
     231  _user_scale = NULL;
    231232  indices = NULL;
    232233
     
    237238
    238239  // no hardest cut (col-unsafe)
    239   SM_var2_hardest_cut_off = -1.0;
     240  SM_var2_hardest_cut_off = -numeric_limits<double>::max();
    240241
    241242  // no pt cutoff for the particles to put in p_uncol_hard
     
    555556}
    556557
     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 */
     569int 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}
    557704
    558705/*
Note: See TracChangeset for help on using the changeset viewer.