Changeset 273e668 in git for external/fastjet/plugins/SISCone/split_merge.cc

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:

• external/fastjet/plugins/SISCone/split_merge.cc (modified) (6 diffs)

external/fastjet/plugins/SISCone/split_merge.cc

@@ -21 +21 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 322                                                          $//
-// $Date:: 2011-11-15 10:12:36 +0100 (Tue, 15 Nov 2011)                     $//
+// $Revision:: 370                                                          $//
+// $Date:: 2014-09-04 17:03:15 +0200 (Thu, 04 Sep 2014)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -28 +28 @@
 #include "siscone_error.h"
 #include "momentum.h"
-#include <math.h>
 #include <limits>   // for max
 #include <iostream>
@@ -34 +33 @@
 #include <sstream>
 #include <cassert>
+#include <cmath>
 
 namespace siscone{
@@ -229 +229 @@
 #endif
 #endif
+  _user_scale = NULL;
   indices = NULL;
 
@@ -237 +238 @@
 
   // no hardest cut (col-unsafe)
-  SM_var2_hardest_cut_off = -1.0;
+  SM_var2_hardest_cut_off = -numeric_limits<double>::max();
 
   // no pt cutoff for the particles to put in p_uncol_hard
@@ -555 +556 @@
 }
 
+
+/*
+ * remove the hardest protocone and declare it a a jet 
+ *  - protocones  list of protocones (initial jet candidates)
+ *  - R2          cone radius (squared)
+ *  - ptmin       minimal pT allowed for jets
+ * return 0 on success, 1 on error
+ *
+ * The list of remaining particles (and the uncollinear-hard ones)
+ * is updated.
+ */
+int Csplit_merge::add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin){
+
+  int i;
+  Cmomentum *c;
+  Cmomentum *v;
+  double eta, phi;
+  double dx, dy;
+  double R;
+  Cjet jet, jet_candidate;
+  bool found_jet = false;
+
+  if (protocones->size()==0)
+    return 1;
+
+  pt_min2 = ptmin*ptmin;
+  R = sqrt(R2);
+
+  // browse protocones
+  // for each of them, build the list of particles in them
+  for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
+    // initialise variables
+    c = &(*p_it);
+
+    // note: cones have been tested => their (eta,phi) coordinates are computed
+    eta = c->eta;
+    phi = c->phi;
+
+    // NOTE: this is common to this method and add_protocones, so it
+    // could be moved into a 'build_jet_from_protocone' method
+    //
+    // browse particles to create cone contents
+    jet_candidate.v = Cmomentum();
+    jet_candidate.pt_tilde=0;
+    jet_candidate.contents.clear();
+    for (i=0;i<n_left;i++){
+      v = &(p_remain[i]);
+      // for security, avoid including particles with infinite rapidity)
+      // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING
+      //if (fabs(v->pz)!=v->E){
+      dx = eta - v->eta;
+      dy = fabs(phi - v->phi);
+      if (dy>M_PI) 
+        dy -= twopi;
+      if (dx*dx+dy*dy<R2){
+        jet_candidate.contents.push_back(v->parent_index);
+        jet_candidate.v+= *v;
+        jet_candidate.pt_tilde+= pt[v->parent_index];
+        v->index=0;
+      }
+    }
+    jet_candidate.n=jet_candidate.contents.size();
+
+    // set the momentum in protocones 
+    // (it was only known through eta and phi up to now)
+    *c = jet_candidate.v;
+    c->eta = eta; // restore exact original coords
+    c->phi = phi; // to avoid rounding error inconsistencies
+
+    // set the jet range
+    jet_candidate.range=Ceta_phi_range(eta,phi,R);
+
+    // check that the protojet has large enough pt
+    if (jet_candidate.v.perp2()<pt_min2)
+      continue;
+
+    // assign the split-merge (or progressive-removal) squared scale variable
+    if (_user_scale) {
+      // sm_var2 is the signed square of the user scale returned
+      // for the jet candidate
+      jet_candidate.sm_var2 = (*_user_scale)(jet_candidate);
+      jet_candidate.sm_var2 *= abs(jet_candidate.sm_var2);
+    } else {
+      jet_candidate.sm_var2 = get_sm_var2(jet_candidate.v, jet_candidate.pt_tilde);
+    }
+
+    // now check if it is possibly the hardest
+    if ((! found_jet) ||
+        (_user_scale ? _user_scale->is_larger(jet_candidate, jet)
+                     : ptcomparison(jet_candidate, jet))){
+      jet = jet_candidate;
+      found_jet = true;
+    }
+  }
+
+  // make sure at least one of the jets has passed the selection
+  if (!found_jet) return 1;  
+  
+  // add the jet to the list of jets
+  jets.push_back(jet);
+  jets[jets.size()-1].v.build_etaphi();
+
+#ifdef DEBUG_SPLIT_MERGE
+  cout << "PR-Jet " << jets.size() << " [size " << next_jet.contents.size() << "]:";
+#endif
+    
+  // update the list of what particles are left
+  int p_remain_index = 0;
+  int contents_index = 0;
+  //sort(next_jet.contents.begin(),next_jet.contents.end());
+  for (int index=0;index<n_left;index++){
+    if ((contents_index<(int) jet.contents.size()) &&
+        (p_remain[index].parent_index == jet.contents[contents_index])){
+      // this particle belongs to the newly found jet
+      // set pass in initial list
+      particles[p_remain[index].parent_index].index = n_pass;
+#ifdef DEBUG_SPLIT_MERGE
+      cout << " " << jet.contents[contents_index];
+#endif
+      contents_index++;
+    } else {
+      // this particle still has to be clustered
+      p_remain[p_remain_index] = p_remain[index];
+      p_remain[p_remain_index].parent_index = p_remain[index].parent_index;
+      p_remain[p_remain_index].index=1;
+      p_remain_index++;
+    }
+  }
+  p_remain.resize(n_left-jet.contents.size());
+  n_left = p_remain.size();
+  jets[jets.size()-1].pass = particles[jet.contents[0]].index;
+
+  // increase the pass index
+  n_pass++;
+
+#ifdef DEBUG_SPLIT_MERGE
+  cout << endl;
+#endif
+
+  // male sure the list of uncol_hard particles (used for the next
+  // stable cone finding) is updated [could probably be more
+  // efficient]
+  merge_collinear_and_remove_soft();
+  
+  return 0;
+}
 
 /*