fastjet/contrib/TauComponents.hh

0001 //  Nsubjettiness Package
0002 //  Questions/Comments?  jthaler@jthaler.net
0003 //
0004 //  Copyright (c) 2011-14
0005 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
0006 //
0007 //  $Id: MeasureFunction.hh 742 2014-08-23 15:43:29Z jthaler $
0008 //----------------------------------------------------------------------
0009 // This file is part of FastJet contrib.
0010 //
0011 // It is free software; you can redistribute it and/or modify it under
0012 // the terms of the GNU General Public License as published by the
0013 // Free Software Foundation; either version 2 of the License, or (at
0014 // your option) any later version.
0015 //
0016 // It is distributed in the hope that it will be useful, but WITHOUT
0017 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
0018 // or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
0019 // License for more details.
0020 //
0021 // You should have received a copy of the GNU General Public License
0022 // along with this code. If not, see <http://www.gnu.org/licenses/>.
0023 //----------------------------------------------------------------------
0024
0025 #ifndef __FASTJET_CONTRIB_TAUCOMPONENTS_HH__
0026 #define __FASTJET_CONTRIB_TAUCOMPONENTS_HH__
0027
0028 #include "fastjet/PseudoJet.hh"
0029 #include "fastjet/ClusterSequence.hh"
0030 #include "fastjet/WrappedStructure.hh"
0031
0032
0033 #include <cmath>
0034 #include <vector>
0035 #include <list>
0036 #include <limits>
0037
0038
0039 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0040
0041 namespace contrib{
0042
0043 // Classes defined in this file.
0044 class TauComponents;
0045 class TauPartition;
0046 class NjettinessExtras;
0047
0048 ///------------------------------------------------------------------------
0049 /// \enum TauMode
0050 /// Specified whether tau value has beam region or denominators
0051 ///------------------------------------------------------------------------
0052 enum TauMode {
0053    UNDEFINED_SHAPE = -1, // Added so that constructor would default to some value
0054    UNNORMALIZED_JET_SHAPE = 0,
0055    NORMALIZED_JET_SHAPE = 1,
0056    UNNORMALIZED_EVENT_SHAPE = 2,
0057    NORMALIZED_EVENT_SHAPE = 3,
0058 };
0059
0060 ///////
0061 //
0062 // TauComponents
0063 //
0064 ///////
0065
0066 ///------------------------------------------------------------------------
0067 /// \class TauComponents
0068 /// \brief Output wrapper for supplemental N-(sub)jettiness information
0069 ///
0070 /// This class creates a wrapper for the various tau/subtau values calculated in Njettiness. This class allows Njettiness access to these variables
0071 /// without ever having to do the calculation itself. It takes in subtau numerators and tau denominator from MeasureFunction
0072 /// and outputs tau numerator, and normalized tau and subtau.
0073 ///------------------------------------------------------------------------
0074 class TauComponents {
0075
0076 public:
0077
0078    /// empty constructor necessary to initialize tau_components in Njettiness
0079    /// later set correctly in Njettiness::getTau function
0080    TauComponents() {}
0081
0082    /// This constructor takes input vector and double and calculates all necessary tau components
0083    TauComponents(TauMode tau_mode,
0084                  const std::vector<double> & jet_pieces_numerator,
0085                  double beam_piece_numerator,
0086                  double denominator,
0087                  const std::vector<PseudoJet> & jets,
0088                  const std::vector<PseudoJet> & axes
0089                  );
0090
0091    /// Test for denominator
0092    bool has_denominator() const;
0093    /// Test for beam region
0094    bool has_beam() const;
0095
0096    /// Return tau value
0097    double tau() const { return _tau; }
0098    /// Return jet regions
0099    const std::vector<double>& jet_pieces() const { return _jet_pieces; }
0100    /// Return beam region
0101    double beam_piece() const { return _beam_piece; }
0102
0103    /// Return jet regions (no denominator)
0104    std::vector<double> jet_pieces_numerator() const { return _jet_pieces_numerator; }
0105    /// Return beam regions (no denominator)
0106    double beam_piece_numerator() const { return _beam_piece_numerator; }
0107    /// Return numerator
0108    double numerator() const { return _numerator; }
0109    /// Return denominator
0110    double denominator() const { return _denominator; }
0111
0112    /// Four-vector of total jet (sum of clustered regions)
0113    PseudoJet total_jet() const { return _total_jet;}
0114    /// Four-vector of jet regions
0115    const std::vector<PseudoJet>& jets() const { return _jets;}
0116    /// Four-vector of axes
0117    const std::vector<PseudoJet>& axes() const { return _axes;}
0118
0119    class StructureType;
0120
0121 protected:
0122
0123    /// Defines whether there is a beam or denominator
0124    TauMode _tau_mode;
0125
0126    std::vector<double> _jet_pieces_numerator;   ///< Constructor input (jet region numerator)
0127    double _beam_piece_numerator;                ///< Constructor input (beam region numerator)
0128    double _denominator;                         ///< Constructor input (denominator)
0129
0130    std::vector<double> _jet_pieces;             ///< Derived value (jet regions)
0131    double _beam_piece;                          ///< Derived value (beam region)
0132    double _numerator;                           ///< Derived value (total numerator)
0133    double _tau;                                 ///< Derived value (final value)
0134
0135    PseudoJet _total_jet;                        ///< Total jet four-vector
0136    std::vector<PseudoJet> _jets;                ///< Jet four-vectors
0137    std::vector<PseudoJet> _axes;                ///< AXes four-vectors
0138
0139 };
0140
0141 ///////
0142 //
0143 // TauPartition
0144 //
0145 ///////
0146
0147 ///------------------------------------------------------------------------
0148 /// \class TauPartition
0149 /// \brief Output wrapper for N-(sub)jettiness partitioning information
0150 ///
0151 /// Class for storing partitioning information.
0152 ///------------------------------------------------------------------------
0153 class TauPartition {
0154
0155 public:
0156    /// empty constructor
0157    TauPartition() {}
0158
0159    /// Make partition of size to hold n_jet partitions
0160    TauPartition(int n_jet) {
0161       _jets_list.resize(n_jet);
0162       _jets_partition.resize(n_jet);
0163    }
0164
0165    /// add a particle to the jet
0166    void push_back_jet(int jet_num, const PseudoJet& part_to_add, int part_index) {
0167       _jets_list[jet_num].push_back(part_index);
0168       _jets_partition[jet_num].push_back(part_to_add);
0169    }
0170
0171    /// add a particle to the beam
0172    void push_back_beam(const PseudoJet& part_to_add, int part_index) {
0173       _beam_list.push_back(part_index);
0174       _beam_partition.push_back(part_to_add);
0175    }
0176
0177    /// return jet regions
0178    PseudoJet jet(int jet_num) const { return join(_jets_partition.at(jet_num)); }
0179    /// return beam region
0180    PseudoJet beam() const { return join(_beam_partition);}
0181
0182    /// return jets
0183    std::vector<PseudoJet> jets() const {
0184       std::vector<PseudoJet> jets;
0185       for (unsigned int i = 0; i < _jets_partition.size(); i++) {
0186          jets.push_back(jet(i));
0187       }
0188       return jets;
0189    }
0190
0191    /// jets in list form
0192    const std::list<int> & jet_list(int jet_num) const { return _jets_list.at(jet_num);}
0193    /// beam in list form
0194    const std::list<int> & beam_list() const { return _beam_list;}
0195    /// all jets in list form
0196    const std::vector<std::list<int> > & jets_list() const { return _jets_list;}
0197
0198 private:
0199
0200    std::vector<std::list<int> > _jets_list;   ///<  jets in list form
0201    std::list<int> _beam_list;                 ///<  beam in list form
0202
0203    std::vector<std::vector<PseudoJet> > _jets_partition;  ///< Partition in jet regions
0204    std::vector<PseudoJet> _beam_partition;                ///< Partition in beam region
0205
0206 };
0207
0208
0209 ///////
0210 //
0211 // NjettinessExtras
0212 //
0213 ///////
0214
0215 ///------------------------------------------------------------------------
0216 /// \class NjettinessExtras
0217 /// \brief ClusterSequence add on for N-jettiness information
0218 ///
0219 /// This class contains the same information as TauComponents, but adds additional ways of linking up
0220 /// the jets found in the ClusterSequence::Extras class.
0221 /// This is done in order to help improve the interface for the main NjettinessPlugin class.
0222 ///------------------------------------------------------------------------
0223 class NjettinessExtras : public ClusterSequence::Extras, public TauComponents {
0224
0225 public:
0226    /// Constructor
0227    NjettinessExtras(TauComponents tau_components,
0228                     std::vector<int> cluster_hist_indices)
0229    : TauComponents(tau_components), _cluster_hist_indices(cluster_hist_indices) {}
0230
0231
0232
0233    /// Ask for tau of the whole event, but by querying a jet
0234    double tau(const fastjet::PseudoJet& /*jet*/) const {return _tau;}
0235
0236    /// Ask for tau of an individual jet
0237    double tau_piece(const fastjet::PseudoJet& jet) const {
0238       if (labelOf(jet) == -1) return std::numeric_limits<double>::quiet_NaN(); // nonsense
0239       return _jet_pieces[labelOf(jet)];
0240    }
0241
0242    /// Find axis associated with jet
0243    fastjet::PseudoJet axis(const fastjet::PseudoJet& jet) const {
0244       return _axes[labelOf(jet)];
0245    }
0246
0247    /// Check if extra information is available.
0248    bool has_njettiness_extras(const fastjet::PseudoJet& jet) const {
0249       return (labelOf(jet) >= 0);
0250    }
0251
0252 private:
0253
0254    /// Store cluster history indices to link up with ClusterSequence
0255    std::vector<int> _cluster_hist_indices;
0256
0257    /// Figure out which jet things belonged to
0258    int labelOf(const fastjet::PseudoJet& jet) const {
0259       int thisJet = -1;
0260       for (unsigned int i = 0; i < _jets.size(); i++) {
0261          if (_cluster_hist_indices[i] == jet.cluster_hist_index()) {
0262             thisJet = i;
0263             break;
0264          }
0265       }
0266       return thisJet;
0267    }
0268
0269 public:
0270
0271    // These are old methods for gaining this information
0272    // The recommended interface is given in TauComponents
0273
0274    /// Tau value
0275    double totalTau() const {return _tau;}
0276    /// Jet regions
0277    std::vector<double> subTaus() const {return _jet_pieces;}
0278
0279    /// Tau value
0280    double totalTau(const fastjet::PseudoJet& /*jet*/) const {
0281       return _tau;
0282    }
0283
0284    /// Jet region
0285    double subTau(const fastjet::PseudoJet& jet) const {
0286       if (labelOf(jet) == -1) return std::numeric_limits<double>::quiet_NaN(); // nonsense
0287       return _jet_pieces[labelOf(jet)];
0288    }
0289
0290    /// beam region
0291    double beamTau() const {
0292       return _beam_piece;
0293    }
0294
0295 };
0296
0297
0298 /// Helper function to find out what njettiness_extras are (from jet)
0299 inline const NjettinessExtras * njettiness_extras(const fastjet::PseudoJet& jet) {
0300    const ClusterSequence * myCS = jet.associated_cluster_sequence();
0301    if (myCS == NULL) return NULL;
0302    const NjettinessExtras* extras = dynamic_cast<const NjettinessExtras*>(myCS->extras());
0303    return extras;
0304 }
0305
0306 /// Helper function to find out what njettiness_extras are (from ClusterSequence)
0307 inline const NjettinessExtras * njettiness_extras(const fastjet::ClusterSequence& myCS) {
0308    const NjettinessExtras* extras = dynamic_cast<const NjettinessExtras*>(myCS.extras());
0309    return extras;
0310 }
0311
0312 ///////
0313 //
0314 // TauComponents::StructureType
0315 //
0316 ///////
0317
0318
0319 ///------------------------------------------------------------------------
0320 /// \class TauComponents::StructureType
0321 /// \brief Wrapped structure for jet-based N-(sub)jettiness information
0322 ///
0323 /// Small wrapped structure to store tau information
0324 /// TODO:  Can these be auto-joined?
0325 ///------------------------------------------------------------------------
0326 class TauComponents::StructureType : public WrappedStructure {
0327
0328 public:
0329    /// Constructor
0330    StructureType(const PseudoJet& j) :
0331       WrappedStructure(j.structure_shared_ptr())
0332    {}
0333
0334    /// tau associated with jet
0335    double tau_piece() const { return _tau_piece; }
0336
0337    /// alternative call, though might be confusing
0338    double tau() const { return _tau_piece; }
0339
0340 private:
0341    friend class TauComponents;
0342    double _tau_piece;  ///< tau value associated with jet
0343 };
0344
0345
0346
0347
0348 } //namespace contrib
0349
0350 FASTJET_END_NAMESPACE
0351
0352 #endif  // __FASTJET_CONTRIB_TAUCOMPONENTS_HH__