Changeset 35cdc46 in git for external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc

Timestamp:

Sep 3, 2014, 3:18:54 PM (10 years ago)

Author:

Pavel Demin <demin@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

be2222c

Parents:

5b5a56b

Message:

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

File:

• external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc (modified) (3 diffs)

external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: MeasureFunction.cc 670 2014-06-06 01:24:42Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -36 +37 @@
 
 // Return all of the necessary TauComponents for specific input particles and axes
-TauComponents MeasureFunction::result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) {
+TauComponents MeasureFunction::result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
+   
+   // first find partition
+   // this sets jetPartitionStorage and beamPartitionStorage
+   PseudoJet beamPartitionStorage;
+   std::vector<fastjet::PseudoJet> jetPartitionStorage = get_partition(particles,axes,&beamPartitionStorage);
+   
+   // then return result calculated from partition
+   return result_from_partition(jetPartitionStorage,axes,&beamPartitionStorage);
+}
 
-   std::vector<double> jetPieces(axes.size(), 0.0);
-   double beamPiece = 0.0;
+std::vector<fastjet::PseudoJet> MeasureFunction::get_partition(const std::vector<fastjet::PseudoJet>& particles,
+                                                               const std::vector<fastjet::PseudoJet>& axes,
+                                                               PseudoJet * beamPartitionStorage) const {
    
-   // Calculates the unnormalized sub-tau values, i.e. a std::vector of the contributions to tau_N of each Voronoi region (or region within R_0)
+   std::vector<std::vector<PseudoJet> > jetPartition(axes.size());
+   std::vector<PseudoJet> beamPartition;
+   
+   // Figures out the partiting of the input particles into the various jet pieces
+   // Based on which axis the parition is closest to
    for (unsigned i = 0; i < particles.size(); i++) {
       
@@ -60 +75 @@
       
       if (j_min == -1) {
-         if (_has_beam) beamPiece += beam_numerator(particles[i]);
+         if (_has_beam) beamPartition.push_back(particles[i]);
          else assert(_has_beam);  // this should never happen.
       } else {
-         jetPieces[j_min] += jet_numerator(particles[i],axes[j_min]);
+         jetPartition[j_min].push_back(particles[i]);
       }
    }
    
-   // Calculates normalization for tau and subTau if _has_denominator is true, otherwise returns 1.0 (i.e. no normalization)
-   double tauDen = 0.0;
-   if (_has_denominator) {
-      for (unsigned i = 0; i < particles.size(); i++) {
-         tauDen += denominator(particles[i]);
-      }
-   } else {
-      tauDen = 1.0; // if no denominator, then 1.0 for no normalization factor
+   // Store beam partition
+   if (beamPartitionStorage) {
+      *beamPartitionStorage = join(beamPartition);
+   }
+
+   // Store jet partitions
+   std::vector<PseudoJet> jetPartitionStorage(axes.size(),PseudoJet(0,0,0,0));
+   for (unsigned j = 0; j < axes.size(); j++) {
+      jetPartitionStorage[j] = join(jetPartition[j]);
    }
    
+   return jetPartitionStorage;
+}
+
+// does partition, but only stores index of PseudoJets
+std::vector<std::list<int> > MeasureFunction::get_partition_list(const std::vector<fastjet::PseudoJet>& particles,
+                                                                 const std::vector<fastjet::PseudoJet>& axes) const {
+
+   std::vector<std::list<int> > jetPartition(axes.size());
+   
+   // Figures out the partiting of the input particles into the various jet pieces
+   // Based on which axis the parition is closest to
+   for (unsigned i = 0; i < particles.size(); i++) {
+      
+      // find minimum distance; start with beam (-1) for reference
+      int j_min = -1;
+      double minRsq;
+      if (_has_beam) minRsq = beam_distance_squared(particles[i]);
+      else minRsq = std::numeric_limits<double>::max(); // make it large value
+      
+      // check to see which axis the particle is closest to
+      for (unsigned j = 0; j < axes.size(); j++) {
+         double tempRsq = jet_distance_squared(particles[i],axes[j]); // delta R distance
+         if (tempRsq < minRsq) {
+            minRsq = tempRsq;
+            j_min = j;
+         }
+      }
+      
+      if (j_min == -1) {
+         assert(_has_beam); // consistency check
+      } else {
+         jetPartition[j_min].push_back(i);
+      }
+   }
+   
+   return jetPartition;
+}
+   
+
+// Uses existing partition and calculates result
+// TODO:  Can we cache this for speed up when doing area subtraction?
+TauComponents MeasureFunction::result_from_partition(const std::vector<fastjet::PseudoJet>& jet_partition,
+                                                     const std::vector<fastjet::PseudoJet>& axes,
+                                                     PseudoJet * beamPartitionStorage) const {
+   
+   std::vector<double> jetPieces(axes.size(), 0.0);
+   double beamPiece = 0.0;
+   
+   double tauDen = 0.0;
+   if (!_has_denominator) tauDen = 1.0;  // if no denominator, then 1.0 for no normalization factor
+   
+   // first find jet pieces
+   for (unsigned j = 0; j < axes.size(); j++) {
+      std::vector<PseudoJet> thisPartition = jet_partition[j].constituents();
+      for (unsigned i = 0; i < thisPartition.size(); i++) {
+         jetPieces[j] += jet_numerator(thisPartition[i],axes[j]); //numerator jet piece
+         if (_has_denominator) tauDen += denominator(thisPartition[i]); // denominator
+      }
+   }
+   
+   // then find beam piece
+   if (_has_beam) {
+      assert(beamPartitionStorage); // make sure I have beam information
+      std::vector<PseudoJet> beamPartition = beamPartitionStorage->constituents();
+
+      for (unsigned i = 0; i < beamPartition.size(); i++) {
+         beamPiece += beam_numerator(beamPartition[i]); //numerator beam piece
+         if (_has_denominator) tauDen += denominator(beamPartition[i]); // denominator
+      }
+   }
    return TauComponents(jetPieces, beamPiece, tauDen, _has_denominator, _has_beam);
 }
+
+   
+   
+   
    
 } //namespace contrib