EIC code displayed by LXR master/​include/​fastjet/​contrib/​MeasureDefinition.hh	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-01-18 09:57:14

 //  Nsubjettiness Package
 //  Questions/Comments?  jthaler@jthaler.net
 //
 //  Copyright (c) 2011-14
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
 //  $Id: MeasureDefinition.hh 828 2015-07-20 14:52:06Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
 //
 // It is free software; you can redistribute it and/or modify it under
 // the terms of the GNU General Public License as published by the
 // Free Software Foundation; either version 2 of the License, or (at
 // your option) any later version.
 //
 // It is distributed in the hope that it will be useful, but WITHOUT
 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 // or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
 // License for more details.
 //
 // You should have received a copy of the GNU General Public License
 // along with this code. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
 
 #ifndef __FASTJET_CONTRIB_MEASUREDEFINITION_HH__
 #define __FASTJET_CONTRIB_MEASUREDEFINITION_HH__
 
 #include "fastjet/PseudoJet.hh"
 #include <cmath>
 #include <vector>
 #include <list>
 #include <limits>
 
 #include "TauComponents.hh"
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
 namespace contrib{
 
 
 
 // The following Measures are available (and their relevant arguments):
 // Recommended for usage as jet shapes
 class DefaultMeasure;               // Default measure from which next classes derive (should not be called directly)
 class NormalizedMeasure;            // (beta,R0)
 class UnnormalizedMeasure;          // (beta)
 class NormalizedCutoffMeasure;      // (beta,R0,Rcutoff)
 class UnnormalizedCutoffMeasure;    // (beta,Rcutoff)
 
 // New measures as of v2.2
 // Recommended for usage as event shapes (or for jet finding)
 class ConicalMeasure;               // (beta,R)
 class OriginalGeometricMeasure;     // (R)
 class ModifiedGeometricMeasure;     // (R)
 class ConicalGeometricMeasure;      // (beta, gamma, R)
 class XConeMeasure;                 // (beta, R)
    
 // Formerly GeometricMeasure, now no longer recommended, kept commented out only for cross-check purposes
 //class DeprecatedGeometricMeasure;         // (beta)
 //class DeprecatedGeometricCutoffMeasure;   // (beta,Rcutoff)
 
    
 ///////
 //
 // MeasureDefinition
 //
 ///////
 
 //This is a helper class for the Minimum Axes Finders. It is defined later.
 class LightLikeAxis;                                          
    
 ///------------------------------------------------------------------------
 /// \class MeasureDefinition
 /// \brief Base class for measure definitions
 ///
 /// This is the base class for measure definitions.  Derived classes will calculate
 /// the tau_N of a jet given a specific measure and a set of axes.  The measure is
 /// determined by various jet and beam distances (and possible normalization factors).
 ///------------------------------------------------------------------------
 class MeasureDefinition {
    
 public:
    
    /// Description of measure and parameters
    virtual std::string description() const = 0;
    
    /// In derived classes, this should return a copy of the corresponding derived class
    virtual MeasureDefinition* create() const = 0;
 
    //The following five functions define the measure by which tau_N is calculated,
    //and are overloaded by the various measures below
    
    /// Distanes to jet axis.  This is called many times, so needs to be as fast as possible
    /// Unless overloaded, it just calls jet_numerator
    virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
       return jet_numerator(particle,axis);
    }
 
    /// Distanes to beam.  This is called many times, so needs to be as fast as possible
    /// Unless overloaded, it just calls beam_numerator
    virtual double beam_distance_squared(const fastjet::PseudoJet& particle) const {
       return beam_numerator(particle);
    }
    
    /// The jet measure used in N-(sub)jettiness
    virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const = 0;
    /// The beam measure used in N-(sub)jettiness
    virtual double beam_numerator(const fastjet::PseudoJet& particle) const = 0;
    
    /// A possible normalization factor
    virtual double denominator(const fastjet::PseudoJet& particle) const = 0;
    
    /// Run a one-pass minimization routine.  There is a generic one-pass minimization that works for a wide variety of measures.
    /// This should be overloaded to create a measure-specific minimization scheme
    virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
                                                              const std::vector<fastjet::PseudoJet>& inputs,
                                                              const std::vector<fastjet::PseudoJet>& seedAxes,
                                                              int nAttempts = 1000,         // cap number of iterations
                                                              double accuracy = 0.0001      // cap distance of closest approach
    ) const;
    
 public:
    
    /// Returns the tau value for a give set of particles and axes
    double result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
       return component_result(particles,axes).tau();
    }
    
    /// Short-hand for the result() function
    inline double operator() (const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
       return result(particles,axes);
    }
    
    /// Return all of the TauComponents for specific input particles and axes
    TauComponents component_result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;
    
    /// Create the partitioning according the jet/beam distances and stores them a TauPartition
    TauPartition get_partition(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;
 
    /// Calculate the tau result using an existing partition
    TauComponents component_result_from_partition(const TauPartition& partition, const std::vector<fastjet::PseudoJet>& axes) const;
 
    
 
    /// virtual destructor
    virtual ~MeasureDefinition(){}
    
 protected:
    
    /// Flag set by derived classes to choose whether or not to use beam/denominator
    TauMode _tau_mode;
    
    /// Flag set by derived classes to say whether cheap get_one_pass_axes method can be used (true by default)
    bool _useAxisScaling;
 
    /// This is the only constructor, which requires _tau_mode and _useAxisScaling to be manually set by derived classes.
    MeasureDefinition() : _tau_mode(UNDEFINED_SHAPE), _useAxisScaling(true) {}
 
    
    /// Used by derived classes to set whether or not to use beam/denominator information
    void setTauMode(TauMode tau_mode) {
       _tau_mode = tau_mode;
    }
 
    /// Used by derived classes to say whether one can use cheap get_one_pass_axes
    void setAxisScaling(bool useAxisScaling) {
       _useAxisScaling = useAxisScaling;
    }
 
    /// Uses denominator information?
    bool has_denominator() const { return (_tau_mode == NORMALIZED_JET_SHAPE || _tau_mode == NORMALIZED_EVENT_SHAPE);}
    /// Uses beam information?
    bool has_beam() const {return (_tau_mode == UNNORMALIZED_EVENT_SHAPE || _tau_mode == NORMALIZED_EVENT_SHAPE);}
 
    /// Create light-like axis (used in default one-pass minimization routine)
    fastjet::PseudoJet lightFrom(const fastjet::PseudoJet& input) const {
       double length = sqrt(pow(input.px(),2) + pow(input.py(),2) + pow(input.pz(),2));
       return fastjet::PseudoJet(input.px()/length,input.py()/length,input.pz()/length,1.0);
    }
    
    /// Shorthand for squaring
    static inline double sq(double x) {return x*x;}
 
 };
    
 
 ///////
 //
 // Default Measures
 //
 ///////
    
    
 ///------------------------------------------------------------------------
 /// \enum DefaultMeasureType
 /// \brief Options for default measure
 ///
 /// Can be used to switch between pp and ee measure types in the DefaultMeasure
 ///------------------------------------------------------------------------
 enum DefaultMeasureType {
    pt_R,       ///  use transverse momenta and boost-invariant angles,
    E_theta,    ///  use energies and angles,
    lorentz_dot, ///  use dot product inspired measure
    perp_lorentz_dot /// use conical geometric inspired measures
 };
 
 ///------------------------------------------------------------------------
 /// \class DefaultMeasure
 /// \brief Base class for default N-subjettiness measure definitions
 ///
 /// This class is the default measure as defined in the original N-subjettiness papers.
 /// Based on the conical measure, but with a normalization factor
 /// This measure is defined as the pT of the particle multiplied by deltaR
 /// to the power of beta. This class includes the normalization factor determined by R0
 ///------------------------------------------------------------------------
 class DefaultMeasure : public MeasureDefinition {
    
 public:
 
    /// Description
    virtual std::string description() const;
    /// To allow copying around of these objects
    virtual DefaultMeasure* create() const {return new DefaultMeasure(*this);}
 
    /// fast jet distance
    virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
       return angleSquared(particle, axis);
    }
    
    /// fast beam distance
    virtual double beam_distance_squared(const fastjet::PseudoJet& /*particle*/) const {
       return _RcutoffSq;
    }
    
    /// true jet distance (given by general definitions of energy and angle)
    virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const{
       double jet_dist = angleSquared(particle, axis);
       if (jet_dist > 0.0) {
          return energy(particle) * std::pow(jet_dist,_beta/2.0);
       } else {
          return 0.0;
       }
    }
    
    /// true beam distance
    virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
       return energy(particle) * std::pow(_Rcutoff,_beta);
    }
    
    /// possible denominator for normalization
    virtual double denominator(const fastjet::PseudoJet& particle) const {
       return energy(particle) * std::pow(_R0,_beta);
    }
 
    /// Special minimization routine (from v1.0 of N-subjettiness)
    virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
                                                              const std::vector<fastjet::PseudoJet>& inputs,
                                                              const std::vector<fastjet::PseudoJet>& seedAxes,
                                                              int nAttempts,   // cap number of iterations
                                                              double accuracy  // cap distance of closest approach
                                                              ) const;
    
 protected:
    double _beta;      ///< Angular exponent
    double _R0;        ///< Normalization factor
    double _Rcutoff;   ///< Cutoff radius
    double _RcutoffSq; ///< Cutoff radius squared
    DefaultMeasureType _measure_type;  ///< Type of measure used (i.e. pp style or ee style)
    
    
    /// Constructor is protected so that no one tries to call this directly.
    DefaultMeasure(double beta, double R0, double Rcutoff, DefaultMeasureType measure_type = pt_R)
    : MeasureDefinition(), _beta(beta), _R0(R0), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)), _measure_type(measure_type)
    {
       if (beta <= 0) throw Error("DefaultMeasure:  You must choose beta > 0.");
       if (R0 <= 0) throw Error("DefaultMeasure:  You must choose R0 > 0.");
       if (Rcutoff <= 0) throw Error("DefaultMeasure:  You must choose Rcutoff > 0.");
    }
    
    /// Added set measure method in case it becomes useful later
    void setDefaultMeasureType(DefaultMeasureType measure_type) {
       _measure_type = measure_type;
    }
 
    /// Generalized energy value (determined by _measure_type)
    double energy(const PseudoJet& jet) const;
    /// Generalized angle value (determined by _measure_type)
    double angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const;
 
    /// Name of _measure_type, so description will include the measure type
    std::string measure_type_name() const {
       if (_measure_type == pt_R) return "pt_R";
       else if (_measure_type == E_theta) return "E_theta";
       else if (_measure_type == lorentz_dot) return "lorentz_dot";
       else if (_measure_type == perp_lorentz_dot) return "perp_lorentz_dot";
       else return "Measure Type Undefined";
    }      
 
    /// templated for speed (TODO: probably should remove, since not clear that there is a speed gain)
    template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
                                                               const std::vector <fastjet::PseudoJet> & inputJets,
                                                               double accuracy) const;
    
    /// called by get_one_pass_axes to update axes iteratively
    std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
                                          const std::vector <fastjet::PseudoJet> & inputJets,
                                          double accuracy) const;
 };
    
 
 ///------------------------------------------------------------------------
 /// \class NormalizedCutoffMeasure
 /// \brief Dimensionless default measure, with radius cutoff
 ///
 /// This measure is just a wrapper for DefaultMeasure
 ///------------------------------------------------------------------------
 class NormalizedCutoffMeasure : public DefaultMeasure {
 
 public:
 
    /// Constructor
    NormalizedCutoffMeasure(double beta, double R0, double Rcutoff, DefaultMeasureType measure_type = pt_R) 
    : DefaultMeasure(beta, R0, Rcutoff, measure_type) {
       setTauMode(NORMALIZED_JET_SHAPE);
    }
 
    /// Description
    virtual std::string description() const;
 
    /// For copying purposes
    virtual NormalizedCutoffMeasure* create() const {return new NormalizedCutoffMeasure(*this);}
 
 };
 
 ///------------------------------------------------------------------------
 /// \class NormalizedMeasure
 /// \brief Dimensionless default measure, with no cutoff
 ///
 /// This measure is the same as NormalizedCutoffMeasure, with Rcutoff taken to infinity.
 ///------------------------------------------------------------------------
 class NormalizedMeasure : public NormalizedCutoffMeasure {
 
 public:
 
    /// Constructor
    NormalizedMeasure(double beta, double R0, DefaultMeasureType measure_type = pt_R)
    : NormalizedCutoffMeasure(beta, R0, std::numeric_limits<double>::max(), measure_type) {
       _RcutoffSq = std::numeric_limits<double>::max();
       setTauMode(NORMALIZED_JET_SHAPE);
    }
    
    /// Description
    virtual std::string description() const;
    /// For copying purposes
    virtual NormalizedMeasure* create() const {return new NormalizedMeasure(*this);}
 
 };
    
 
 ///------------------------------------------------------------------------
 /// \class UnnormalizedCutoffMeasure
 /// \brief Dimensionful default measure, with radius cutoff
 ///
 /// This class is the unnormalized conical measure. The only difference from NormalizedCutoffMeasure
 /// is that the denominator is defined to be 1.0 by setting _has_denominator to false.
 /// class UnnormalizedCutoffMeasure : public NormalizedCutoffMeasure {
 ///------------------------------------------------------------------------
 class UnnormalizedCutoffMeasure : public DefaultMeasure {
    
 public:
    
    /// Since all methods are identical, UnnormalizedMeasure inherits directly
    /// from NormalizedMeasure. R0 is a dummy value since the value of R0 is unecessary for this class,
    /// and the "false" flag sets _has_denominator in MeasureDefinition to false so no denominator is used.
    UnnormalizedCutoffMeasure(double beta, double Rcutoff, DefaultMeasureType measure_type = pt_R)
    : DefaultMeasure(beta, std::numeric_limits<double>::quiet_NaN(), Rcutoff, measure_type) {
       setTauMode(UNNORMALIZED_EVENT_SHAPE);
    }
 
    /// Description
    virtual std::string description() const;
    /// For copying purposes
    virtual UnnormalizedCutoffMeasure* create() const {return new UnnormalizedCutoffMeasure(*this);}
 
 };
 
    
 ///------------------------------------------------------------------------
 /// \class UnnormalizedMeasure
 /// \brief Dimensionless default measure, with no cutoff
 ///
 /// This measure is the same as UnnormalizedCutoffMeasure, with Rcutoff taken to infinity.
 ///------------------------------------------------------------------------
 class UnnormalizedMeasure : public UnnormalizedCutoffMeasure {
    
 public:
    /// Since all methods are identical, UnnormalizedMeasure inherits directly
    /// from NormalizedMeasure. R0 is a dummy value since the value of R0 is unecessary for this class,
    /// and the "false" flag sets _has_denominator in MeasureDefinition to false so no denominator is used.
    UnnormalizedMeasure(double beta, DefaultMeasureType measure_type = pt_R)
    : UnnormalizedCutoffMeasure(beta, std::numeric_limits<double>::max(), measure_type) {
       _RcutoffSq = std::numeric_limits<double>::max();
       setTauMode(UNNORMALIZED_JET_SHAPE);
    }
 
    /// Description
    virtual std::string description() const;
    
    /// For copying purposes
    virtual UnnormalizedMeasure* create() const {return new UnnormalizedMeasure(*this);}
    
 };
 
 
 ///------------------------------------------------------------------------
 /// \class ConicalMeasure
 /// \brief Dimensionful event-shape measure, with radius cutoff
 ///
 /// Very similar to UnnormalizedCutoffMeasure, but with different normalization convention
 /// and using the new default one-pass minimization algorithm.
 /// Axes are also made to be light-like to ensure sensible behavior
 /// Intended to be used as an event shape.
 ///------------------------------------------------------------------------
 class ConicalMeasure : public MeasureDefinition {
    
 public:
    
    /// Constructor
    ConicalMeasure(double beta, double Rcutoff)
    :   MeasureDefinition(), _beta(beta), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
       if (beta <= 0) throw Error("ConicalMeasure:  You must choose beta > 0.");
       if (Rcutoff <= 0) throw Error("ConicalMeasure:  You must choose Rcutoff > 0.");
       setTauMode(UNNORMALIZED_EVENT_SHAPE);
    }
    
    /// Description
    virtual std::string description() const;
    /// For copying purposes
    virtual ConicalMeasure* create() const {return new ConicalMeasure(*this);}
    
    /// fast jet distance
    virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
       PseudoJet lightAxis = lightFrom(axis);
       return particle.squared_distance(lightAxis);
    }
    
    /// fast beam distance
    virtual double beam_distance_squared(const fastjet::PseudoJet&  /*particle*/) const {
       return _RcutoffSq;
    }
    
    
    /// true jet distance
    virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
       PseudoJet lightAxis = lightFrom(axis);
       double jet_dist = particle.squared_distance(lightAxis)/_RcutoffSq;
       double jet_perp = particle.perp();
       
       if (_beta == 2.0) {
          return jet_perp * jet_dist;
       } else {
          return jet_perp * pow(jet_dist,_beta/2.0);
       }
    }
    
    /// true beam distance
    virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
       return particle.perp();
    }
    
    /// no denominator used for this measure
    virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
       return std::numeric_limits<double>::quiet_NaN();
    }
    
 protected:
    double _beta;     ///< angular exponent
    double _Rcutoff;  ///< effective jet radius
    double _RcutoffSq;///< effective jet radius squared
 };
    
 
 
 ///------------------------------------------------------------------------
 /// \class OriginalGeometricMeasure
 /// \brief Dimensionful event-shape measure, with dot-product distances
 ///
 /// This class is the original (and hopefully now correctly coded) geometric measure.
 /// This measure is defined by the Lorentz dot product between
 /// the particle and the axis.  This class does not include normalization of tau_N.
 /// New in Nsubjettiness v2.2
 /// NOTE: This is defined differently from the DeprecatedGeometricMeasure which are now commented out.
 ///------------------------------------------------------------------------
 class OriginalGeometricMeasure : public MeasureDefinition {
 
 public:
    /// Constructor
    OriginalGeometricMeasure(double Rcutoff)
    : MeasureDefinition(), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
       if (Rcutoff <= 0) throw Error("OriginalGeometricMeasure:  You must choose Rcutoff > 0.");
       setTauMode(UNNORMALIZED_EVENT_SHAPE);
       setAxisScaling(false);  // No need to rescale axes (for speed up in one-pass minimization)
    }
 
    /// Description
    virtual std::string description() const;
    /// For copying purposes
    virtual OriginalGeometricMeasure* create() const {return new OriginalGeometricMeasure(*this);}
    
    // This class uses the default jet_distance_squared and beam_distance_squared
 
    /// true jet measure
    virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
       return dot_product(lightFrom(axis), particle)/_RcutoffSq;
    }
 
    /// true beam measure
    virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
       fastjet::PseudoJet beam_a(0,0,1,1);
       fastjet::PseudoJet beam_b(0,0,-1,1);
       double min_perp = std::min(dot_product(beam_a, particle),dot_product(beam_b, particle));
       return min_perp;
    }
 
    /// no denominator needed for this measure.
    virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
       return std::numeric_limits<double>::quiet_NaN();
    }
    
 protected:
    double _Rcutoff;   ///< Effective jet radius (rho = R^2)
    double _RcutoffSq; ///< Effective jet radius squared
 
 };
 
 
 ///------------------------------------------------------------------------
 /// \class ModifiedGeometricMeasure
 /// \brief Dimensionful event-shape measure, with dot-product distances, modified beam measure
 ///
 /// This class is the Modified geometric measure.  This jet measure is defined by the Lorentz dot product between
 /// the particle and the axis, as in the Original Geometric Measure. The beam measure is defined differently from
 /// the above OriginalGeometric to allow for more conical jets. New in Nsubjettiness v2.2
 ///------------------------------------------------------------------------
 class ModifiedGeometricMeasure : public MeasureDefinition {
 
 public:
    /// Constructor
    ModifiedGeometricMeasure(double Rcutoff)
    :  MeasureDefinition(), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
       if (Rcutoff <= 0) throw Error("ModifiedGeometricMeasure:  You must choose Rcutoff > 0.");
       setTauMode(UNNORMALIZED_EVENT_SHAPE);
       setAxisScaling(false);  // No need to rescale axes (for speed up in one-pass minimization)
    }
 
    /// Description
    virtual std::string description() const;
    /// For copying purposes
    virtual ModifiedGeometricMeasure* create() const {return new ModifiedGeometricMeasure(*this);}
    
    // This class uses the default jet_distance_squared and beam_distance_squared
 
    /// True jet measure
    virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
       return dot_product(lightFrom(axis), particle)/_RcutoffSq;
    }
 
    /// True beam measure
    virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
       fastjet::PseudoJet lightParticle = lightFrom(particle);
       return 0.5*particle.mperp()*lightParticle.pt();
    }
 
    /// This measure does not require a denominator
    virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
       return std::numeric_limits<double>::quiet_NaN();
    }
    
 protected:
    double _Rcutoff;   ///< Effective jet radius (rho = R^2)
    double _RcutoffSq; ///< Effective jet radius squared
 
 
 };
 
 ///------------------------------------------------------------------------
 /// \class ConicalGeometricMeasure
 /// \brief Dimensionful event-shape measure, basis for XCone jet algorithm
 ///
 /// This class is the Conical Geometric measure.  This measure is defined by the Lorentz dot product between
 /// the particle and the axis normalized by the axis and particle pT, as well as a factor of cosh(y) to vary
 /// the rapidity depepdence of the beam. New in Nsubjettiness v2.2, and the basis for the XCone jet algorithm
 ///------------------------------------------------------------------------
 class ConicalGeometricMeasure : public MeasureDefinition {
 
 public:
    
    /// Constructor
    ConicalGeometricMeasure(double jet_beta, double beam_gamma, double Rcutoff)
    :  MeasureDefinition(),  _jet_beta(jet_beta), _beam_gamma(beam_gamma), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)){
       if (jet_beta <= 0)   throw Error("ConicalGeometricMeasure:  You must choose beta > 0.");
       if (beam_gamma <= 0) throw Error("ConicalGeometricMeasure:  You must choose gamma > 0.");
       if (Rcutoff <= 0)   throw Error("ConicalGeometricMeasure:  You must choose Rcutoff > 0.");
       setTauMode(UNNORMALIZED_EVENT_SHAPE);
    }
 
    /// Description
    virtual std::string description() const;
    /// For copying purposes
    virtual ConicalGeometricMeasure* create() const {return new ConicalGeometricMeasure(*this);}
    
    /// fast jet measure
    virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
       fastjet::PseudoJet lightAxis = lightFrom(axis);
       double pseudoRsquared = 2.0*dot_product(lightFrom(axis),particle)/(lightAxis.pt()*particle.pt());
       return pseudoRsquared;
    }
 
    /// fast beam measure
    virtual double beam_distance_squared(const fastjet::PseudoJet&  /*particle*/) const {
       return _RcutoffSq;
    }
 
    /// true jet measure
    virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
       double jet_dist = jet_distance_squared(particle,axis)/_RcutoffSq;
       if (jet_dist > 0.0) {
          fastjet::PseudoJet lightParticle = lightFrom(particle);
          double weight = (_beam_gamma == 1.0) ? 1.0 : std::pow(0.5*lightParticle.pt(),_beam_gamma - 1.0);
          return particle.pt() * weight * std::pow(jet_dist,_jet_beta/2.0);
       } else {
          return 0.0;
       }
    }
 
    /// true beam measure
    virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
       fastjet::PseudoJet lightParticle = lightFrom(particle);
       double weight = (_beam_gamma == 1.0) ? 1.0 : std::pow(0.5*lightParticle.pt(),_beam_gamma - 1.0);
       return particle.pt() * weight;
    }
 
    /// no denominator needed
    virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
       return std::numeric_limits<double>::quiet_NaN();
    }
       
 protected:
    double _jet_beta;   ///< jet angular exponent
    double _beam_gamma; ///< beam angular exponent (gamma = 1.0 is recommended)
    double _Rcutoff;    ///< effective jet radius
    double _RcutoffSq;  ///< effective jet radius squared
    
 };
 
 ///------------------------------------------------------------------------
 /// \class XConeMeasure
 /// \brief Dimensionful event-shape measure used in XCone jet algorithm
 ///
 /// This class is the XCone Measure.  This is the default measure for use with the
 /// XCone algorithm. It is identical to the conical geometric measure but with gamma = 1.0.
 ///------------------------------------------------------------------------
 class XConeMeasure : public ConicalGeometricMeasure {
 
 public:
    /// Constructor
    XConeMeasure(double jet_beta, double R)
    :   ConicalGeometricMeasure(jet_beta,
                                1.0, // beam_gamma, hard coded at gamma = 1.0 default
                                R    // Rcutoff scale
                                ) { }
 
    /// Description
    virtual std::string description() const;
    /// For copying purposes
    virtual XConeMeasure* create() const {return new XConeMeasure(*this);}
 
 };
 
 ///------------------------------------------------------------------------
 /// \class LightLikeAxis
 /// \brief Helper class to define light-like axes directions
 ///
 /// This is a helper class for the minimization routines.
 /// It creates a convenient way of defining axes in order to better facilitate calculations.
 ///------------------------------------------------------------------------
 class LightLikeAxis {
 
 public:
    /// Bare constructor
    LightLikeAxis() : _rap(0.0), _phi(0.0), _weight(0.0), _mom(0.0) {}
    /// Constructor
    LightLikeAxis(double my_rap, double my_phi, double my_weight, double my_mom) :
    _rap(my_rap), _phi(my_phi), _weight(my_weight), _mom(my_mom) {}
    
    /// Rapidity
    double rap() const {return _rap;}
    /// Azimuth
    double phi() const {return _phi;}
    /// weight factor
    double weight() const {return _weight;}
    /// pt momentum
    double mom() const {return _mom;}
    
    /// set rapidity
    void set_rap(double my_set_rap) {_rap = my_set_rap;}
    /// set azimuth
    void set_phi(double my_set_phi) {_phi = my_set_phi;}
    /// set weight factor
    void set_weight(double my_set_weight) {_weight = my_set_weight;}
    /// set pt momentum
    void set_mom(double my_set_mom) {_mom = my_set_mom;}
    /// set all kinematics
    void reset(double my_rap, double my_phi, double my_weight, double my_mom) {_rap=my_rap; _phi=my_phi; _weight=my_weight; _mom=my_mom;}
 
    /// Return PseudoJet version
    fastjet::PseudoJet ConvertToPseudoJet();
    
    /// Squared distance to PseudoJet
    double DistanceSq(const fastjet::PseudoJet& input) const {
       return DistanceSq(input.rap(),input.phi());
    }
 
    /// Distance to PseudoJet
    double Distance(const fastjet::PseudoJet& input) const {
       return std::sqrt(DistanceSq(input));
    }
 
    /// Squared distance to Lightlikeaxis
    double DistanceSq(const LightLikeAxis& input) const {
       return DistanceSq(input.rap(),input.phi());
    }
 
    /// Distance to Lightlikeaxis
    double Distance(const LightLikeAxis& input) const {
       return std::sqrt(DistanceSq(input));
    }
 
 private:
    double _rap;    ///< rapidity
    double _phi;    ///< azimuth
    double _weight; ///< weight factor
    double _mom;    ///< pt momentum
    
    /// Internal squared distance calculation
    double DistanceSq(double rap2, double phi2) const {
       double rap1 = _rap;
       double phi1 = _phi;
       
       double distRap = rap1-rap2;
       double distPhi = std::fabs(phi1-phi2);
       if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
       return distRap*distRap + distPhi*distPhi;
    }
    
    /// Internal distance calculation
    double Distance(double rap2, double phi2) const {
       return std::sqrt(DistanceSq(rap2,phi2));
    }
       
 };
    
    
 ////------------------------------------------------------------------------
 ///// \class DeprecatedGeometricCutoffMeasure
 //// This class is the old, incorrectly coded geometric measure.
 //// It is kept in case anyone wants to check old code, but should not be used for production purposes.
 //class DeprecatedGeometricCutoffMeasure : public MeasureDefinition {
 //
 //public:
 //
 //   // Please, please don't use this.
 //   DeprecatedGeometricCutoffMeasure(double jet_beta, double Rcutoff)
 //   :   MeasureDefinition(),
 //      _jet_beta(jet_beta),
 //      _beam_beta(1.0), // This is hard coded, since alternative beta_beam values were never checked.
 //      _Rcutoff(Rcutoff),
 //      _RcutoffSq(sq(Rcutoff)) {
 //         setTauMode(UNNORMALIZED_EVENT_SHAPE);
 //         setAxisScaling(false);
 //         if (jet_beta != 2.0) {
 //         throw Error("Geometric minimization is currently only defined for beta = 2.0.");
 //      }
 //   }
 //
 //   virtual std::string description() const;
 //
 //   virtual DeprecatedGeometricCutoffMeasure* create() const {return new DeprecatedGeometricCutoffMeasure(*this);}
 //
 //   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
 //      fastjet::PseudoJet lightAxis = lightFrom(axis);
 //      double pseudoRsquared = 2.0*dot_product(lightFrom(axis),particle)/(lightAxis.pt()*particle.pt());
 //      return pseudoRsquared;
 //   }
 //
 //   virtual double beam_distance_squared(const fastjet::PseudoJet&  /*particle*/) const {
 //      return _RcutoffSq;
 //   }
 //
 //   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
 //      fastjet::PseudoJet lightAxis = lightFrom(axis);
 //      double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(lightAxis.pt(),_beam_beta - 1.0);
 //      return particle.pt() * weight * std::pow(jet_distance_squared(particle,axis),_jet_beta/2.0);
 //   }
 //
 //   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
 //      double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(particle.pt()/particle.e(),_beam_beta - 1.0);
 //      return particle.pt() * weight * std::pow(_Rcutoff,_jet_beta);
 //   }
 //
 //   virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
 //      return std::numeric_limits<double>::quiet_NaN();
 //   }
 //
 //   virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
 //                                                             const std::vector<fastjet::PseudoJet>& inputs,
 //                                                             const std::vector<fastjet::PseudoJet>& seedAxes,
 //                                                             int nAttempts,    // cap number of iterations
 //                                                             double accuracy   // cap distance of closest approach
 //                                                            ) const;
 //
 //protected:
 //   double _jet_beta;
 //   double _beam_beta;
 //   double _Rcutoff;
 //   double _RcutoffSq;
 //
 //};
 //
 //// ------------------------------------------------------------------------
 //// / \class DeprecatedGeometricMeasure
 //// Same as DeprecatedGeometricMeasureCutoffMeasure, but with Rcutoff taken to infinity.
 //// NOTE:  This class should not be used for production purposes.
 //class DeprecatedGeometricMeasure : public DeprecatedGeometricCutoffMeasure {
 //
 //public:
 //   DeprecatedGeometricMeasure(double beta)
 //   : DeprecatedGeometricCutoffMeasure(beta,std::numeric_limits<double>::max()) {
 //      _RcutoffSq = std::numeric_limits<double>::max();
 //      setTauMode(UNNORMALIZED_JET_SHAPE);
 //   }
 //
 //   virtual std::string description() const;
 //
 //   virtual DeprecatedGeometricMeasure* create() const {return new DeprecatedGeometricMeasure(*this);}
 //};
 
 
 } //namespace contrib
 
 FASTJET_END_NAMESPACE
 
 #endif  // __FASTJET_CONTRIB_MEASUREDEFINITION_HH__

This page was automatically generated by the 2.3.7 LXR engine. The LXR team