Fork me on GitHub

Ignore:
Timestamp:
Dec 9, 2014, 1:27:13 PM (10 years ago)
Author:
Michele <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
37deb3b, 9e991f8
Parents:
f6b6ee7 (diff), e7e90df (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'TestFastJet310b1'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • external/fastjet/contribs/Nsubjettiness/MeasureFunction.hh

    rf6b6ee7 r49234af  
    55//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
    66//
     7//  $Id: MeasureFunction.hh 678 2014-06-12 20:43:03Z jthaler $
    78//----------------------------------------------------------------------
    89// This file is part of FastJet contrib.
    … …  
    3132#include <limits>
    3233
     34
    3335FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
    3436
    3537namespace contrib{
    3638
    37 inline double sq(double x) {return x*x;}
    38 
    39 ///////
    40 //
    41 // Measure Function
    42 //
    43 ///////
    44 
     39///////
     40//
     41// TauComponents
     42// (eventually we might want to put this in a separate header file)
     43//
     44///////
     45   
    4546/// \class TauComponents
    4647// This class creates a wrapper for the various tau/subtau values calculated in Njettiness. This class allows Njettiness access to these variables
    4748// without ever having to do the calculation itself. It takes in subtau numerators and tau denominator from MeasureFunction
    4849// and outputs tau numerator, and normalized tau and subtau.
    49 // TODO:  Consider merging with NjettinessExtras.  Add axes information?
    5050class TauComponents {
    51 private:
    52    
    53    // these values are input in the constructor
    54    std::vector<double> _jet_pieces_numerator;
    55    double _beam_piece_numerator;
    56    double _denominator;
    57    bool _has_denominator; //added so that TauComponents knows if denominator is used or not
    58    bool _has_beam; //added so that TauComponents knows if beam regions is used or not
    59    
    60    // these values are derived from above values
    61    std::vector<double> _jet_pieces;
    62    double _beam_piece;
    63    double _numerator;
    64    double _tau;
    65    
    6651   
    6752public:
    … …  
    115100   double beam_piece() const { return _beam_piece; }
    116101   double tau() const { return _tau; }
    117    
    118 };
     102
     103private:
     104   
     105   // these values are input in the constructor
     106   std::vector<double> _jet_pieces_numerator;
     107   double _beam_piece_numerator;
     108   double _denominator;
     109   bool _has_denominator; //added so that TauComponents knows if denominator is used or not
     110   bool _has_beam; //added so that TauComponents knows if beam regions is used or not
     111   
     112   // these values are derived from above values
     113   std::vector<double> _jet_pieces;
     114   double _beam_piece;
     115   double _numerator;
     116   double _tau;
     117   
     118};
     119
     120///////
     121//
     122// Measure Function
     123//
     124///////
     125
    119126
    120127//------------------------------------------------------------------------
    … …  
    125132class MeasureFunction {
    126133   
     134public:
     135   //These functions define the measure by which tau_N is calculated,
     136   //and are overloaded by the various measures below
     137   
     138   // Distanes to axes.  These are called many times, so need to be as fast as possible
     139   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const = 0;
     140   virtual double beam_distance_squared(const fastjet::PseudoJet& particle) const = 0;
     141   
     142   // The actual measures used in N-(sub)jettiness
     143   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const  = 0;
     144   virtual double beam_numerator(const fastjet::PseudoJet& particle) const = 0;
     145   
     146   // a possible normalization factor
     147   virtual double denominator(const fastjet::PseudoJet& particle) const = 0;
     148   
     149   //------
     150   // The functions below call the above functions and are not virtual
     151   //------
     152   
     153   // Return all of the necessary TauComponents for specific input particles and axes
     154   // Also have optional pointers to get out information about partitioning
     155   TauComponents result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;
     156
     157   // Just getting tau value if that is all that is needed
     158   double tau(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
     159      return result(particles,axes).tau();
     160   }
     161   
     162   // Create the partitioning and stores internally
     163   std::vector<fastjet::PseudoJet> get_partition(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes, PseudoJet * beamPartitionStorage = NULL) const;
     164
     165   // Essentially same as get_partition, but in the form needed for the jet algorithm
     166   std::vector<std::list<int> > get_partition_list(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;
     167
     168   // calculates the tau result using an existing partition
     169   TauComponents result_from_partition(const std::vector<fastjet::PseudoJet>& jet_partitioning, const std::vector<fastjet::PseudoJet>& axes, PseudoJet * beamPartitionStorage = NULL) const;
     170
     171   // shorthand for squaring
     172   static inline double sq(double x) {return x*x;}
     173
     174   //virtual destructor
     175   virtual ~MeasureFunction(){}
     176   
    127177protected:
    128178   //bool set by derived classes to choose whether or not to use the denominator
    … …  
    133183   MeasureFunction(bool has_denominator = true, bool has_beam = true) : _has_denominator(has_denominator), _has_beam(has_beam) {}
    134184   
    135 public:
    136    virtual ~MeasureFunction(){}
    137    
    138    //These functions define the measure by which tau_N is calculated,
    139    //and are overloaded by the various measures below
    140    
    141    // Distanes to axes.  These are called many times, so need to be as fast as possible
    142    virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) = 0;
    143    virtual double beam_distance_squared(const fastjet::PseudoJet& particle) = 0;
    144    
    145    // The actual measures used in N-(sub)jettiness
    146    virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) = 0;
    147    virtual double beam_numerator(const fastjet::PseudoJet& particle) = 0;
    148    
    149    // a possible normalization factor
    150    virtual double denominator(const fastjet::PseudoJet& particle) = 0;
    151    
    152    
    153    // These functions call the above functions and are not virtual
    154    
    155    // Do I cluster a particle into a jet?
    156    bool do_cluster(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
    157       return (jet_distance_squared(particle,axis) <= beam_distance_squared(particle));
    158    }
    159    
    160    // Return all of the necessary TauComponents for specific input particles and axes
    161    TauComponents result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes);
    162 
    163    double tau(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) {
    164       return result(particles,axes).tau();
    165    }
    166 
    167    
    168 };
    169 
    170 
    171 /// \class DefaultNormalizedMeasure
     185};
     186
     187
     188/// \class DefaultNormalizedMeasureFunction
    172189// This class is the default measure, inheriting from the class above. This class will calculate tau_N
    173190// of a jet according to this measure. This measure is defined as the pT of the particle multiplied by deltaR
    174191// to the power of beta. This class includes the normalization factor determined by R0
    175 class DefaultNormalizedMeasure : public MeasureFunction {
    176 
    177    private:
    178       double _beta;
    179       double _R0;
    180       double _Rcutoff;
    181 
    182    public:
    183 
    184       DefaultNormalizedMeasure(double beta, double R0, double Rcutoff, bool normalized = true)
    185       : MeasureFunction(normalized), _beta(beta), _R0(R0), _Rcutoff(Rcutoff) {}
    186 
    187       virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
    188          return particle.squared_distance(axis);
    189       }
    190    
    191       virtual double beam_distance_squared(const fastjet::PseudoJet& particle) {
    192          return sq(_Rcutoff);
    193       }
    194 
    195       virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
    196          return particle.perp() * std::pow(jet_distance_squared(particle,axis),_beta/2.0);
    197       }
    198    
    199       virtual double beam_numerator(const fastjet::PseudoJet& particle) {
    200          return particle.perp() * std::pow(_Rcutoff,_beta);
    201       }
    202 
    203       virtual double denominator(const fastjet::PseudoJet& particle) {
    204          return particle.perp() * std::pow(_R0,_beta);
    205       }
    206 
     192class DefaultNormalizedMeasureFunction : public MeasureFunction {
     193
     194public:
     195
     196   DefaultNormalizedMeasureFunction(double beta, double R0, double Rcutoff, bool normalized = true)
     197   : MeasureFunction(normalized), _beta(beta), _R0(R0), _Rcutoff(Rcutoff) {}
     198
     199   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
     200      return particle.squared_distance(axis);
     201   }
     202
     203   virtual double beam_distance_squared(const fastjet::PseudoJet& /*particle*/) const {
     204      return sq(_Rcutoff);
     205   }
     206
     207   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const{
     208      return particle.perp() * std::pow(jet_distance_squared(particle,axis),_beta/2.0);
     209   }
     210
     211   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
     212      return particle.perp() * std::pow(_Rcutoff,_beta);
     213   }
     214
     215   virtual double denominator(const fastjet::PseudoJet& particle) const {
     216      return particle.perp() * std::pow(_R0,_beta);
     217   }
     218   
     219private:
     220   double _beta;
     221   double _R0;
     222   double _Rcutoff;
     223
     224   
    207225};
    208226
    209227//------------------------------------------------------------------------
    210 /// \class DefaultUnnormalizedMeasure
     228/// \class DefaultUnnormalizedMeasureFunction
    211229// This class is the unnormalized default measure, inheriting from the class above. The only difference from above
    212230// is that the denominator is defined to be 1.0 by setting _has_denominator to false.
    213 class DefaultUnnormalizedMeasure : public DefaultNormalizedMeasure {
    214 
    215    public:
    216       // Since all methods are identical, UnnormalizedMeasure inherits directly from NormalizedMeasure. R0 is defaulted to NAN since the value of R0 is unecessary for this class.
    217       // the "false" flag sets _has_denominator in MeasureFunction to false so no denominator is used.
    218       DefaultUnnormalizedMeasure(double beta, double Rcutoff) : DefaultNormalizedMeasure(beta, NAN, Rcutoff, false) {}
    219 
    220       
     231class DefaultUnnormalizedMeasureFunction : public DefaultNormalizedMeasureFunction {
     232
     233public:
     234   // Since all methods are identical, UnnormalizedMeasure inherits directly
     235   // from NormalizedMeasure. R0 is a dummy value since the value of R0 is unecessary for this class,
     236   // and the "false" flag sets _has_denominator in MeasureFunction to false so no denominator is used.
     237   DefaultUnnormalizedMeasureFunction(double beta, double Rcutoff)
     238   : DefaultNormalizedMeasureFunction(beta, std::numeric_limits<double>::quiet_NaN(), Rcutoff, false) {}
    221239};
    222240
    223241//------------------------------------------------------------------------
    224 /// \class GeometricMeasure
     242/// \class GeometricMeasureFunction
    225243// This class is the geometic measure, inheriting from the class above. This class will calculate tau_N
    226244// of a jet according to this measure. This measure is defined by the Lorentz dot product between
    227245// the particle and the axis. This class includes normalization of tau_N.
    228 class GeometricMeasure : public MeasureFunction {
    229 
    230    private:
    231       double _jet_beta;
    232       double _beam_beta;
    233       double _Rcutoff;
    234 
    235       // create light-like axis
    236       fastjet::PseudoJet lightFrom(const fastjet::PseudoJet& input) const {
    237          double length = sqrt(pow(input.px(),2) + pow(input.py(),2) + pow(input.pz(),2));
    238          return fastjet::PseudoJet(input.px()/length,input.py()/length,input.pz()/length,1.0);
    239       }
    240 
    241    public:
    242       // Right now, we are hard coded for beam_beta = 1.0, but that will need to change
    243       GeometricMeasure(double jet_beta, double Rcutoff) : _jet_beta(jet_beta), _beam_beta(1.0), _Rcutoff(Rcutoff) {}
    244    
    245       virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
    246          fastjet::PseudoJet lightAxis = lightFrom(axis);
    247          double pseudoRsquared = 2.0*dot_product(lightFrom(axis),particle)/(lightAxis.pt()*particle.pt());
    248          return pseudoRsquared;
    249       }
    250    
    251       virtual double beam_distance_squared(const fastjet::PseudoJet& particle) {
    252          return sq(_Rcutoff);
    253       }
    254 
    255       virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
    256          fastjet::PseudoJet lightAxis = lightFrom(axis);
    257          double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(lightAxis.pt(),_beam_beta - 1.0);
    258          return particle.pt() * weight * std::pow(jet_distance_squared(particle,axis),_jet_beta/2.0);
    259       }
    260    
    261       virtual double beam_numerator(const fastjet::PseudoJet& particle) {
    262          double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(particle.pt()/particle.e(),_beam_beta - 1.0);
    263          return particle.pt() * weight * std::pow(_Rcutoff,_jet_beta);
    264       }
    265 
    266       virtual double denominator(const fastjet::PseudoJet& particle) {
    267          return 1.0;
    268       }
    269 };
    270 
    271 
     246class GeometricMeasureFunction : public MeasureFunction {
     247
     248public:
     249   // Right now, we are hard coded for beam_beta = 1.0, but that will need to change
     250   GeometricMeasureFunction(double jet_beta, double Rcutoff) :
     251     MeasureFunction(false), // doesn't have denominator
     252     _jet_beta(jet_beta), _beam_beta(1.0), _Rcutoff(Rcutoff) {}
     253
     254   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
     255      fastjet::PseudoJet lightAxis = lightFrom(axis);
     256      double pseudoRsquared = 2.0*dot_product(lightFrom(axis),particle)/(lightAxis.pt()*particle.pt());
     257      return pseudoRsquared;
     258   }
     259
     260   virtual double beam_distance_squared(const fastjet::PseudoJet&  /*particle*/) const {
     261      return sq(_Rcutoff);
     262   }
     263
     264   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
     265      fastjet::PseudoJet lightAxis = lightFrom(axis);
     266      double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(lightAxis.pt(),_beam_beta - 1.0);
     267      return particle.pt() * weight * std::pow(jet_distance_squared(particle,axis),_jet_beta/2.0);
     268   }
     269
     270   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
     271      double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(particle.pt()/particle.e(),_beam_beta - 1.0);
     272      return particle.pt() * weight * std::pow(_Rcutoff,_jet_beta);
     273   }
     274
     275   virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
     276      return std::numeric_limits<double>::quiet_NaN();
     277   }
     278   
     279   
     280private:
     281   double _jet_beta;
     282   double _beam_beta;
     283   double _Rcutoff;
     284   
     285   // create light-like axis
     286   fastjet::PseudoJet lightFrom(const fastjet::PseudoJet& input) const {
     287      double length = sqrt(pow(input.px(),2) + pow(input.py(),2) + pow(input.pz(),2));
     288      return fastjet::PseudoJet(input.px()/length,input.py()/length,input.pz()/length,1.0);
     289   }
     290
     291};
     292   
     293   
    272294} //namespace contrib
    273295
Note: See TracChangeset for help on using the changeset viewer.