diff --git a/cpptools/src/fjcontrib/buildtools/get_fj_contrib.sh b/cpptools/src/fjcontrib/buildtools/get_fj_contrib.sh index b5ea263..486e399 100755 --- a/cpptools/src/fjcontrib/buildtools/get_fj_contrib.sh +++ b/cpptools/src/fjcontrib/buildtools/get_fj_contrib.sh @@ -15,6 +15,7 @@ if [ -d ${srcdir} ]; then [ ! -e fjcontrib-${fjcontrib_version}.tar.gz ] && wget http://fastjet.hepforge.org/contrib/downloads/fjcontrib-${fjcontrib_version}.tar.gz cd - if [ -e ${wdir}/fjcontrib-${fjcontrib_version}.tar.gz ]; then + # RecursiveTools if [ ! -d ${srcdir}/fjcontrib-${fjcontrib_version}/RecursiveTools ]; then cd ${srcdir} tar zxvf ${wdir}/fjcontrib-${fjcontrib_version}.tar.gz fjcontrib-${fjcontrib_version}/RecursiveTools @@ -25,6 +26,7 @@ if [ -d ${srcdir} ]; then cp -v ${srcdir}/custom/Util.* ${srcdir}/fjcontrib-${fjcontrib_version}/RecursiveTools fi + # LundPlane if [ ! -d ${srcdir}/fjcontrib-${fjcontrib_version}/LundPlane ]; then cd ${srcdir} tar zxvf ${wdir}/fjcontrib-${fjcontrib_version}.tar.gz fjcontrib-${fjcontrib_version}/LundPlane @@ -38,12 +40,21 @@ if [ -d ${srcdir} ]; then # cp -v ${srcdir}/custom/GroomerShopUI.* ${srcdir}/fjcontrib-${fjcontrib_version}/LundPlane fi + # ConstituentSubtractor if [ ! -d ${srcdir}/fjcontrib-${fjcontrib_version}/ConstituentSubtractor ]; then cd ${srcdir} tar zxvf ${wdir}/fjcontrib-${fjcontrib_version}.tar.gz fjcontrib-${fjcontrib_version}/ConstituentSubtractor rm fjcontrib-${fjcontrib_version}/ConstituentSubtractor/example_*.cc fi + # Nsubjettiness + if [ ! -d ${srcdir}/fjcontrib-${fjcontrib_version}/Nsubjettiness ]; then + cd ${srcdir} + tar zxvf ${wdir}/fjcontrib-${fjcontrib_version}.tar.gz fjcontrib-${fjcontrib_version}/Nsubjettiness + rm fjcontrib-${fjcontrib_version}/Nsubjettiness/example_*.cc + patch fjcontrib-${fjcontrib_version}/Nsubjettiness/MeasureDefinition.hh -i ${srcdir}/patches/MeasureDefinition.patch + patch fjcontrib-${fjcontrib_version}/Nsubjettiness/AxesDefinition.hh -i ${srcdir}/patches/AxesDefinition.patch + fi fi fi fi diff --git a/cpptools/src/fjcontrib/interface/Nsubjettiness/AxesDefinition_i.hh b/cpptools/src/fjcontrib/interface/Nsubjettiness/AxesDefinition_i.hh new file mode 100644 index 0000000..6c48245 --- /dev/null +++ b/cpptools/src/fjcontrib/interface/Nsubjettiness/AxesDefinition_i.hh @@ -0,0 +1,314 @@ +#include "MeasureDefinition.hh" +#include "ExtraRecombiners.hh" + +#include "fastjet/PseudoJet.hh" +#include + +#include +#include +#include +#include + +class AxesDefinition { +public: + virtual std::vector get_starting_axes(int n_jets, + const std::vector& inputs, + const MeasureDefinition * measure); + virtual std::string short_description(); + virtual std::string description(); + virtual AxesDefinition* create(); + + std::vector get_refined_axes(int n_jets, + const std::vector& inputs, + const std::vector& seedAxes, + const MeasureDefinition * measure = NULL) const; + std::vector get_axes(int n_jets, + const std::vector& inputs, + const MeasureDefinition * measure = NULL) const; + inline std::vector operator() (int n_jets, + const std::vector& inputs, + const MeasureDefinition * measure = NULL) const; + + enum AxesRefiningEnum { + UNDEFINED_REFINE = -1, // added to create a default value + NO_REFINING = 0, + ONE_PASS = 1, + MULTI_PASS = 100, + }; + + int nPass() const; + bool givesRandomizedResults() const; + bool needsManualAxes() const; + void setNPass(int nPass, + int nAttempts = 1000, + double accuracy = 0.0001, + double noise_range = 1.0 // only needed for MultiPass minimization + ); + virtual ~AxesDefinition(); +}; + +class ExclusiveJetAxes : public AxesDefinition { +public: + ExclusiveJetAxes(fastjet::JetDefinition def); + virtual std::vector get_starting_axes(int n_jets, + const std::vector & inputs, + const MeasureDefinition * ) const; + virtual std::string short_description() const; + virtual std::string description() const; + virtual ExclusiveJetAxes* create() const; +}; + +class ExclusiveCombinatorialJetAxes : public AxesDefinition { +public: + ExclusiveCombinatorialJetAxes(fastjet::JetDefinition def, int nExtra = 0); + virtual std::vector get_starting_axes(int n_jets, + const std::vector & inputs, + const MeasureDefinition *measure) const; + virtual std::string short_description() const; + virtual std::string description() const; + virtual ExclusiveCombinatorialJetAxes* create() const; +}; + +class HardestJetAxes : public AxesDefinition { +public: + HardestJetAxes(fastjet::JetDefinition def); + virtual std::vector get_starting_axes(int n_jets, + const std::vector & inputs, + const MeasureDefinition * ) const; + virtual std::string short_description() const; + virtual std::string description() const; + virtual HardestJetAxes* create() const; +// private: +// fastjet::JetDefinition _def; ///< Jet Definition to use. +// static LimitedWarning _too_few_axes_warning; +}; + +class KT_Axes : public ExclusiveJetAxes { +public: + /// Constructor + KT_Axes() + : ExclusiveJetAxes(fastjet::JetDefinition(fastjet::kt_algorithm, + fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant + fastjet::E_scheme, + fastjet::Best) + ); + virtual std::string short_description() const; + virtual std::string description() const; + virtual KT_Axes* create() const; +}; + +class CA_Axes : public ExclusiveJetAxes { +public: + CA_Axes(); + virtual std::string short_description() const; + virtual std::string description() const; + virtual CA_Axes* create() const; +}; + + +class AntiKT_Axes : public HardestJetAxes { +public: + AntiKT_Axes(double R0); + virtual std::string short_description() const; + virtual std::string description() const; + virtual AntiKT_Axes* create() const; +// protected: +// double _R0; ///< AKT jet radius +}; + +class JetDefinitionWrapper { +public: + JetDefinitionWrapper(JetAlgorithm jet_algorithm_in, double R_in, double xtra_param_in, const JetDefinition::Recombiner *recombiner); + JetDefinitionWrapper(JetAlgorithm jet_algorithm_in, double R_in, const JetDefinition::Recombiner *recombiner, fastjet::Strategy strategy_in); + JetDefinition getJetDef(); +// private: +// JetDefinition jet_def; ///< my jet definition +}; + +class WTA_KT_Axes : public ExclusiveJetAxes { +public: + WTA_KT_Axes(); + virtual std::string short_description() const; + virtual std::string description() const; + virtual WTA_KT_Axes* create() const; +}; + +class WTA_CA_Axes : public ExclusiveJetAxes { +public: + WTA_CA_Axes(); + virtual std::string short_description() const; + virtual std::string description() const; + virtual WTA_CA_Axes* create() const; +}; + +class GenKT_Axes : public ExclusiveJetAxes { +public: + GenKT_Axes(double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual GenKT_Axes* create() const; +// protected: +// double _p; ///< genkT power +// double _R0; ///< jet radius +}; + +class WTA_GenKT_Axes : public ExclusiveJetAxes { +public: + WTA_GenKT_Axes(double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual WTA_GenKT_Axes* create() const; +// protected: +// double _p; ///< genkT power +// double _R0; ///< jet radius +}; + +class GenET_GenKT_Axes : public ExclusiveJetAxes { +public: + GenET_GenKT_Axes(double delta, double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual GenET_GenKT_Axes* create() const; +// protected: +// double _delta; ///< Recombination pT weighting +// double _p; ///< GenkT power +// double _R0; ///< jet radius +}; + +class OnePass_KT_Axes : public KT_Axes { +public: + OnePass_KT_Axes(); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_KT_Axes* create() const; +}; + +class OnePass_CA_Axes : public CA_Axes { +public: + OnePass_CA_Axes(); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_CA_Axes* create() const; +}; + +class OnePass_AntiKT_Axes : public AntiKT_Axes { +public: + OnePass_AntiKT_Axes(double R0); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_AntiKT_Axes* create() const; + +}; + +class OnePass_WTA_KT_Axes : public WTA_KT_Axes { +public: + OnePass_WTA_KT_Axes(); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_WTA_KT_Axes* create() const; +}; + +class OnePass_WTA_CA_Axes : public WTA_CA_Axes { +public: + OnePass_WTA_CA_Axes(); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_WTA_CA_Axes* create() const; +}; + +class OnePass_GenKT_Axes : public GenKT_Axes { +public: + OnePass_GenKT_Axes(double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_GenKT_Axes* create() const; +}; + +class OnePass_WTA_GenKT_Axes : public WTA_GenKT_Axes { +public: + OnePass_WTA_GenKT_Axes(double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_WTA_GenKT_Axes* create() const; +}; + +class OnePass_GenET_GenKT_Axes : public GenET_GenKT_Axes { +public: + OnePass_GenET_GenKT_Axes(double delta, double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_GenET_GenKT_Axes* create() const; +}; + +class Manual_Axes : public AxesDefinition { +public: + Manual_Axes(); + virtual std::vector get_starting_axes(int, + const std::vector&, + const MeasureDefinition *) const; + virtual std::string short_description() const; + virtual std::string description() const; + virtual Manual_Axes* create() const; +}; + +class OnePass_Manual_Axes : public Manual_Axes { +public: + OnePass_Manual_Axes(); + virtual std::string short_description() const; + virtual std::string description() const; + virtual OnePass_Manual_Axes* create() const; +}; + +class MultiPass_Axes : public KT_Axes { +public: + MultiPass_Axes(unsigned int Npass); + virtual std::string short_description() const; + virtual std::string description() const; + virtual MultiPass_Axes* create() const; +}; + +class MultiPass_Manual_Axes : public Manual_Axes { +public: + MultiPass_Manual_Axes(unsigned int Npass); + virtual std::string short_description() const; + virtual std::string description() const; + virtual MultiPass_Manual_Axes* create() const; +}; + +class Comb_GenKT_Axes : public ExclusiveCombinatorialJetAxes { +public: + Comb_GenKT_Axes(int nExtra, double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual Comb_GenKT_Axes* create() const; +// private: +// double _nExtra; ///< Number of extra axes +// double _p; ///< GenkT power +// double _R0; ///< jet radius +}; + +class Comb_WTA_GenKT_Axes : public ExclusiveCombinatorialJetAxes { +public: + Comb_WTA_GenKT_Axes(int nExtra, double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual Comb_WTA_GenKT_Axes* create() const; +// private: +// double _nExtra; ///< Number of extra axes +// double _p; ///< GenkT power +// double _R0; ///< jet radius +}; + +class Comb_GenET_GenKT_Axes : public ExclusiveCombinatorialJetAxes { +public: + Comb_GenET_GenKT_Axes(int nExtra, double delta, double p, double R0 = fastjet::JetDefinition::max_allowable_R); + virtual std::string short_description() const; + virtual std::string description() const; + virtual Comb_GenET_GenKT_Axes* create() const; +// private: +// double _nExtra; ///< Number of extra axes +// double _delta; ///< Recombination pT weighting exponent +// double _p; ///< GenkT power +// double _R0; ///< jet radius +}; + \ No newline at end of file diff --git a/cpptools/src/fjcontrib/interface/Nsubjettiness/ExtraRecombiners_i.hh b/cpptools/src/fjcontrib/interface/Nsubjettiness/ExtraRecombiners_i.hh new file mode 100644 index 0000000..04680ff --- /dev/null +++ b/cpptools/src/fjcontrib/interface/Nsubjettiness/ExtraRecombiners_i.hh @@ -0,0 +1,28 @@ +#include "fastjet/PseudoJet.hh" +#include "fastjet/JetDefinition.hh" + +#include +#include +#include +#include +#include +#include +#include + +class GeneralEtSchemeRecombiner : public fastjet::JetDefinition::Recombiner { +public: + GeneralEtSchemeRecombiner(double delta); + virtual std::string description() const; + virtual void recombine(const fastjet::PseudoJet & pa, + const fastjet::PseudoJet & pb, + fastjet::PseudoJet & pab) const; +}; + +class WinnerTakeAllRecombiner : public fastjet::JetDefinition::Recombiner { +public: + WinnerTakeAllRecombiner(double alpha = 1.0); + virtual std::string description() const; + virtual void recombine(const fastjet::PseudoJet & pa, + const fastjet::PseudoJet & pb, + fastjet::PseudoJet & pab) const; +}; diff --git a/cpptools/src/fjcontrib/interface/Nsubjettiness/MeasureDefinition_i.hh b/cpptools/src/fjcontrib/interface/Nsubjettiness/MeasureDefinition_i.hh new file mode 100644 index 0000000..2ddd29e --- /dev/null +++ b/cpptools/src/fjcontrib/interface/Nsubjettiness/MeasureDefinition_i.hh @@ -0,0 +1,178 @@ +#include +#include +#include +#include + +#include "TauComponents.hh" + +class MeasureDefinition { +public: + virtual std::string description() const; + virtual MeasureDefinition* create() const; + virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_distance_squared(const fastjet::PseudoJet& particle) const; + virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_numerator(const fastjet::PseudoJet& particle) const; + virtual double denominator(const fastjet::PseudoJet& particle) const; + virtual std::vector get_one_pass_axes(int n_jets, + const std::vector& inputs, + const std::vector& seedAxes, + int nAttempts = 1000, // cap number of iterations + double accuracy = 0.0001 // cap distance of closest approach + ) const; + double result(const std::vector& particles, const std::vector& axes) const; + inline double operator() (const std::vector& particles, const std::vector& axes) const; + TauComponents component_result(const std::vector& particles, const std::vector& axes) const; + TauPartition get_partition(const std::vector& particles, const std::vector& axes) const; + TauComponents component_result_from_partition(const TauPartition& partition, const std::vector& axes) const; + virtual ~MeasureDefinition(); + TauMode _tau_mode; + bool _useAxisScaling; + MeasureDefinition(); + void setTauMode(TauMode tau_mode); + void setAxisScaling(bool useAxisScaling); + bool has_denominator() const; + bool has_beam() const; + fastjet::PseudoJet lightFrom(const fastjet::PseudoJet& input) const; + static inline double sq(double x); +}; + + +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 LightLikeAxis { +public: + LightLikeAxis(); + LightLikeAxis(double my_rap, double my_phi, double my_weight, double my_mom); + double rap() const; + double phi() const; + double weight() const; + double mom() const; + void set_rap(double my_set_rap); + void set_phi(double my_set_phi); + void set_weight(double my_set_weight); + void set_mom(double my_set_mom); + void reset(double my_rap, double my_phi, double my_weight, double my_mom); + fastjet::PseudoJet ConvertToPseudoJet(); + double DistanceSq(const fastjet::PseudoJet& input) const; + double Distance(const fastjet::PseudoJet& input) const; + double DistanceSq(const LightLikeAxis& input) const; + double Distance(const LightLikeAxis& input) const; +}; + +class DefaultMeasure : public MeasureDefinition { + public: + virtual std::string description() const; + virtual DefaultMeasure* create() const; + virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_distance_squared(const fastjet::PseudoJet& /*particle*/) const; + virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_numerator(const fastjet::PseudoJet& particle) const; + virtual double denominator(const fastjet::PseudoJet& particle) const; + virtual std::vector get_one_pass_axes(int n_jets, + const std::vector& inputs, + const std::vector& seedAxes, + int nAttempts, // cap number of iterations + double accuracy // cap distance of closest approach + ) const; + DefaultMeasure(double beta, double R0, double Rcutoff, DefaultMeasureType measure_type = pt_R); + void setDefaultMeasureType(DefaultMeasureType measure_type); + double energy(const PseudoJet& jet) const; + double angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const; + std::string measure_type_name() const; + template std::vector UpdateAxesFast(const std::vector & old_axes, + const std::vector & inputJets, + double accuracy) const; + std::vector UpdateAxes(const std::vector & old_axes, + const std::vector & inputJets, + double accuracy) const; +}; + + +class NormalizedCutoffMeasure : public DefaultMeasure { +public: + NormalizedCutoffMeasure(double beta, double R0, double Rcutoff, DefaultMeasureType measure_type = pt_R); + virtual std::string description() const; + virtual NormalizedCutoffMeasure* create() const; +}; + +class NormalizedMeasure : public NormalizedCutoffMeasure { +public: + NormalizedMeasure(double beta, double R0, DefaultMeasureType measure_type = pt_R); + virtual std::string description() const; + virtual NormalizedMeasure* create() const; +}; + +class UnnormalizedCutoffMeasure : public DefaultMeasure { +public: + UnnormalizedCutoffMeasure(double beta, double Rcutoff, DefaultMeasureType measure_type = pt_R); + virtual std::string description() const; + virtual UnnormalizedCutoffMeasure* create() const; + +}; + + +class UnnormalizedMeasure : public UnnormalizedCutoffMeasure { +public: + UnnormalizedMeasure(double beta, DefaultMeasureType measure_type = pt_R); + virtual std::string description() const; + virtual UnnormalizedMeasure* create() const; +}; + +class ConicalMeasure : public MeasureDefinition { +public: + ConicalMeasure(double beta, double Rcutoff); + virtual std::string description() const; + virtual ConicalMeasure* create() const; + virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_distance_squared(const fastjet::PseudoJet& /*particle*/) const; + virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_numerator(const fastjet::PseudoJet& particle) const; + virtual double denominator(const fastjet::PseudoJet& /*particle*/) const; +}; + + +class OriginalGeometricMeasure : public MeasureDefinition { +public: + OriginalGeometricMeasure(double Rcutoff); + virtual std::string description() const; + virtual OriginalGeometricMeasure* create() const; + virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_numerator(const fastjet::PseudoJet& particle) const; + virtual double denominator(const fastjet::PseudoJet& /*particle*/) const; +}; + + +class ModifiedGeometricMeasure : public MeasureDefinition { +public: + ModifiedGeometricMeasure(double Rcutoff); + virtual std::string description() const; + virtual ModifiedGeometricMeasure* create() const; + virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_numerator(const fastjet::PseudoJet& particle) const; + virtual double denominator(const fastjet::PseudoJet& /*particle*/) const; +}; + +class ConicalGeometricMeasure : public MeasureDefinition { +public: + ConicalGeometricMeasure(double jet_beta, double beam_gamma, double Rcutoff); + virtual std::string description() const; + virtual ConicalGeometricMeasure* create() const; + virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_distance_squared(const fastjet::PseudoJet& /*particle*/) const; + virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const; + virtual double beam_numerator(const fastjet::PseudoJet& particle) const; + virtual double denominator(const fastjet::PseudoJet& /*particle*/) const; +}; + +class XConeMeasure : public ConicalGeometricMeasure { +public: + XConeMeasure(double jet_beta, double R); + virtual std::string description() const; + virtual XConeMeasure* create() const; +}; diff --git a/cpptools/src/fjcontrib/interface/Nsubjettiness/NjettinessPlugin_i.hh b/cpptools/src/fjcontrib/interface/Nsubjettiness/NjettinessPlugin_i.hh new file mode 100644 index 0000000..3054d56 --- /dev/null +++ b/cpptools/src/fjcontrib/interface/Nsubjettiness/NjettinessPlugin_i.hh @@ -0,0 +1,56 @@ +#include + +#include "Njettiness.hh" +#include "MeasureDefinition.hh" +#include "AxesDefinition.hh" +#include "TauComponents.hh" + +#include "fastjet/ClusterSequence.hh" +#include "fastjet/JetDefinition.hh" + +#include +#include + +class NjettinessPlugin : public JetDefinition::Plugin { +public: + /// Constructor with same arguments as Nsubjettiness (N, AxesDefinition, MeasureDefinition) + NjettinessPlugin(int N, + const AxesDefinition & axes_def, + const MeasureDefinition & measure_def); + virtual std::string description () const; + virtual double R() const; + virtual void run_clustering(ClusterSequence&) const; + void setAxes(const std::vector & myAxes); + virtual ~NjettinessPlugin(); + +// private: +// Njettiness _njettinessFinder; ///< The core Njettiness that does the heavy lifting +// int _N; ///< Number of exclusive jets to find. +// /// Warning if the user tries to use v1.0.3 constructor. +// static LimitedWarning _old_constructor_warning; + +public: + NjettinessPlugin(int N, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode); + NjettinessPlugin(int N, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1); + NjettinessPlugin(int N, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1, + double para2); + NjettinessPlugin(int N, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1, + double para2, + double para3); + NjettinessPlugin(int N, + Njettiness::AxesMode mode, + double beta, + double R0, + double Rcutoff=std::numeric_limits::max()); +}; diff --git a/cpptools/src/fjcontrib/interface/Nsubjettiness/Njettiness_i.hh b/cpptools/src/fjcontrib/interface/Nsubjettiness/Njettiness_i.hh new file mode 100644 index 0000000..d3948a8 --- /dev/null +++ b/cpptools/src/fjcontrib/interface/Nsubjettiness/Njettiness_i.hh @@ -0,0 +1,80 @@ +#include "MeasureDefinition.hh" +#include "AxesDefinition.hh" +#include "TauComponents.hh" + +#include "fastjet/PseudoJet.hh" +#include "fastjet/SharedPtr.hh" +#include + +#include +#include +#include + +class Njettiness { +public: + Njettiness(const AxesDefinition & axes_def, const MeasureDefinition & measure_def); + ~Njettiness(); + void setAxes(const std::vector & myAxes); + TauComponents getTauComponents(unsigned n_jets, const std::vector & inputJets) const; + double getTau(unsigned n_jets, const std::vector & inputJets) const; + TauComponents currentTauComponents() const; + std::vector currentAxes() const; + std::vector seedAxes() const; + std::vector currentJets() const; + fastjet::PseudoJet currentBeam() const; + TauPartition currentPartition() const; + +// private: +// /// AxesDefinition to use. Implemented as SharedPtrs to avoid memory management headaches +// SharedPtr _axes_def; +// /// MeasureDefinition to use. Implemented as SharedPtrs to avoid memory management headaches +// SharedPtr _measure_def; +// // Information about the current information +// // Defined as mutables, so user should be aware that these change when getTau is called. +// // TODO: These are not thread safe and should be fixed somehow +// mutable TauComponents _current_tau_components; //automatically set to have components of 0; these values will be set by the getTau function call +// mutable std::vector _currentAxes; //axes found after minimization +// mutable std::vector _seedAxes; // axes used prior to minimization (if applicable) +// mutable TauPartition _currentPartition; //partitioning information +// /// Warning if the user tries to use v1.0.3 measure style. +// static LimitedWarning _old_measure_warning; +// /// Warning if the user tries to use v1.0.3 axes style. +// static LimitedWarning _old_axes_warning; + +public: + + enum AxesMode { + kt_axes, // exclusive kt axes + ca_axes, // exclusive ca axes + antikt_0p2_axes, // inclusive hardest axes with antikt-0.2 + wta_kt_axes, // Winner Take All axes with kt + wta_ca_axes, // Winner Take All axes with CA + onepass_kt_axes, // one-pass minimization from kt starting point + onepass_ca_axes, // one-pass minimization from ca starting point + onepass_antikt_0p2_axes, // one-pass minimization from antikt-0.2 starting point + onepass_wta_kt_axes, //one-pass minimization of WTA axes with kt + onepass_wta_ca_axes, //one-pass minimization of WTA axes with ca + min_axes, // axes that minimize N-subjettiness (100 passes by default) + manual_axes, // set your own axes with setAxes() + onepass_manual_axes // one-pass minimization from manual starting point + }; + + enum MeasureMode { + normalized_measure, //default normalized measure + unnormalized_measure, //default unnormalized measure + geometric_measure, //geometric measure + normalized_cutoff_measure, //default normalized measure with explicit Rcutoff + unnormalized_cutoff_measure, //default unnormalized measure with explicit Rcutoff + geometric_cutoff_measure //geometric measure with explicit Rcutoff + }; + + Njettiness(AxesMode axes_mode, const MeasureDefinition & measure_def); + Njettiness(AxesMode axes_mode, + MeasureMode measure_mode, + int num_para, + double para1 = std::numeric_limits::quiet_NaN(), + double para2 = std::numeric_limits::quiet_NaN(), + double para3 = std::numeric_limits::quiet_NaN()); + AxesDefinition* createAxesDef(AxesMode axes_mode) const; + MeasureDefinition* createMeasureDef(MeasureMode measure_mode, int num_para, double para1, double para2, double para3) const; +}; diff --git a/cpptools/src/fjcontrib/interface/Nsubjettiness/Nsubjettiness_i.hh b/cpptools/src/fjcontrib/interface/Nsubjettiness/Nsubjettiness_i.hh new file mode 100644 index 0000000..b71f2ea --- /dev/null +++ b/cpptools/src/fjcontrib/interface/Nsubjettiness/Nsubjettiness_i.hh @@ -0,0 +1,96 @@ +#include + +#include "Njettiness.hh" + +#include "fastjet/FunctionOfPseudoJet.hh" +#include +#include + +class Nsubjettiness; +class NsubjettinessRatio; + +class Nsubjettiness : public FunctionOfPseudoJet { +public: + Nsubjettiness(int N, + const AxesDefinition& axes_def, + const MeasureDefinition& measure_def); + double result(const fastjet::PseudoJet& jet) const; + TauComponents component_result(const fastjet::PseudoJet& jet) const; + void setAxes(const std::vector & myAxes); + std::vector seedAxes() const; + std::vector currentAxes() const; + std::vector currentSubjets() const; + TauComponents currentTauComponents() const; + TauPartition currentPartition() const; + +// private: +// /// Core Njettiness code that is called +// Njettiness _njettinessFinder; // TODO: should muck with this so result can be const without this mutable +// /// Number of subjets to find +// int _N; +// /// Warning if the user tries to use v1.0.3 constructor. +// static LimitedWarning _old_constructor_warning; + +public: + Nsubjettiness(int N, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode); + Nsubjettiness(int N, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1); + Nsubjettiness(int N, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1, + double para2); + Nsubjettiness(int N, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1, + double para2, + double para3); + Nsubjettiness(int N, + Njettiness::AxesMode axes_mode, + double beta, + double R0, + double Rcutoff=std::numeric_limits::max()); +}; + + +class NsubjettinessRatio : public FunctionOfPseudoJet { +public: + NsubjettinessRatio(int N, + int M, + const AxesDefinition & axes_def, + const MeasureDefinition & measure_def); + double result(const PseudoJet& jet) const; + +// private: +// Nsubjettiness _nsub_numerator; ///< Function for numerator +// Nsubjettiness _nsub_denominator; ///< Function for denominator + +public: + NsubjettinessRatio(int N, + int M, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode); + NsubjettinessRatio(int N, + int M, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1); + NsubjettinessRatio(int N, + int M, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1, + double para2); + NsubjettinessRatio(int N, + int M, + Njettiness::AxesMode axes_mode, + Njettiness::MeasureMode measure_mode, + double para1, + double para2, + double para3); +}; diff --git a/cpptools/src/fjcontrib/interface/Nsubjettiness/TauComponents_i.hh b/cpptools/src/fjcontrib/interface/Nsubjettiness/TauComponents_i.hh new file mode 100644 index 0000000..e16297a --- /dev/null +++ b/cpptools/src/fjcontrib/interface/Nsubjettiness/TauComponents_i.hh @@ -0,0 +1,84 @@ +#include +#include +#include +#include + +enum TauMode { + UNDEFINED_SHAPE = -1, // Added so that constructor would default to some value + UNNORMALIZED_JET_SHAPE = 0, + NORMALIZED_JET_SHAPE = 1, + UNNORMALIZED_EVENT_SHAPE = 2, + NORMALIZED_EVENT_SHAPE = 3, +}; + +class TauComponents { +public: + TauComponents(); + TauComponents(TauMode tau_mode, + const std::vector & jet_pieces_numerator, + double beam_piece_numerator, + double denominator, + const std::vector & jets, + const std::vector & axes + ); + bool has_denominator() const; + bool has_beam() const; + double tau() const; + const std::vector& jet_pieces() const; + double beam_piece() const; + std::vector jet_pieces_numerator() const; + double beam_piece_numerator() const; + double numerator() const; + double denominator() const; + fastjet::PseudoJet total_jet() const; + const std::vector& jets() const; + const std::vector& axes() const; + + class StructureType : public fastjet::WrappedStructure { + public: + StructureType(const fastjet::PseudoJet& j) : + WrappedStructure(j.structure_shared_ptr()); + double tau_piece() const; + double tau() const; + }; +}; + +class TauPartition { +public: + TauPartition(); + TauPartition(int n_jet); + void push_back_jet(int jet_num, const fastjet::PseudoJet& part_to_add, int part_index); + void push_back_beam(const fastjet::PseudoJet& part_to_add, int part_index); + fastjet::PseudoJet jet(int jet_num) const; + fastjet::PseudoJet beam() const; + std::vector jets() const; + const std::list & jet_list(int jet_num) const; + const std::list & beam_list() const; + const std::vector > & jets_list() const; +}; + + +class NjettinessExtras : public fastjet::ClusterSequence::Extras, public TauComponents { + +public: + /// Constructor + NjettinessExtras(TauComponents tau_components, + std::vector cluster_hist_indices); + double tau(const fastjet::PseudoJet& /*jet*/) const; + double tau_piece(const fastjet::PseudoJet& jet) const; + fastjet::PseudoJet axis(const fastjet::PseudoJet& jet) const; + bool has_njettiness_extras(const fastjet::PseudoJet& jet) const; + + double totalTau() const; + std::vector subTaus() const; + double totalTau(const fastjet::PseudoJet& /*jet*/) const; + double subTau(const fastjet::PseudoJet& jet) const; + double beamTau() const; +}; + + +/// Helper function to find out what njettiness_extras are (from jet) +inline const NjettinessExtras * njettiness_extras(const fastjet::PseudoJet& jet); + +/// Helper function to find out what njettiness_extras are (from ClusterSequence) +inline const NjettinessExtras * njettiness_extras(const fastjet::ClusterSequence& myCS); diff --git a/cpptools/src/fjcontrib/interface/Nsubjettiness/XConePlugin_i.hh b/cpptools/src/fjcontrib/interface/Nsubjettiness/XConePlugin_i.hh new file mode 100644 index 0000000..18a00f0 --- /dev/null +++ b/cpptools/src/fjcontrib/interface/Nsubjettiness/XConePlugin_i.hh @@ -0,0 +1,47 @@ +#include + +#include "NjettinessPlugin.hh" + +#include "fastjet/ClusterSequence.hh" +#include "fastjet/JetDefinition.hh" + +#include +#include + +class XConePlugin : public NjettinessPlugin { +public: + XConePlugin(int N, double R0, double beta = 2.0); + virtual std::string description () const; + virtual double R() const; + virtual ~XConePlugin(); + +// private: +// static double calc_delta(double beta); +// static double calc_power(double beta); +// double _N; ///< Number of desired jets +// double _R0; ///< Jet radius +// double _beta; ///< Angular exponent (beta = 2.0 is dafault, beta = 1.0 is recoil-free) + +// public: + +}; + +class PseudoXConePlugin : public NjettinessPlugin { +public: + PseudoXConePlugin(int N, double R0, double beta = 2.0); + virtual std::string description () const; + virtual double R() const; + // run_clustering is done by NjettinessPlugin + virtual ~PseudoXConePlugin(); + +// private: +// static double calc_delta(double beta); +// /// Static call used within the constructor to set the recommended p value +// static double calc_power(double beta); +// double _N; ///< Number of desired jets +// double _R0; ///< Jet radius +// double _beta; ///< Angular exponent (beta = 2.0 is dafault, beta = 1.0 is recoil-free) + +// public: + +}; diff --git a/cpptools/src/fjcontrib/interface/fjcontrib.i b/cpptools/src/fjcontrib/interface/fjcontrib.i index b51e5c6..bd02c44 100644 --- a/cpptools/src/fjcontrib/interface/fjcontrib.i +++ b/cpptools/src/fjcontrib/interface/fjcontrib.i @@ -6,3 +6,4 @@ %include lundplane.i %include recursivetools.i %include constituentsubtractor.i +%include nsubjettiness.i diff --git a/cpptools/src/fjcontrib/interface/nsubjettiness.i b/cpptools/src/fjcontrib/interface/nsubjettiness.i new file mode 100644 index 0000000..98e3581 --- /dev/null +++ b/cpptools/src/fjcontrib/interface/nsubjettiness.i @@ -0,0 +1,35 @@ +%module nsubjettiness + +%{ + #include + #include + #include "fastjet/ClusterSequence.hh" + #include "fastjet/WrappedStructure.hh" + + #include + #include + #include + #include + #include + #include + #include + #include +%} + +%include "std_string.i" +%include "std_vector.i" + +namespace fastjet{ + namespace contrib{ + + %include "Nsubjettiness/TauComponents_i.hh" + %include "Nsubjettiness/MeasureDefinition_i.hh" + %include "Nsubjettiness/ExtraRecombiners_i.hh" + %include "Nsubjettiness/AxesDefinition_i.hh" + %include "Nsubjettiness/Njettiness_i.hh" + %include "Nsubjettiness/NjettinessPlugin_i.hh" + %include "Nsubjettiness/Nsubjettiness_i.hh" + %include "Nsubjettiness/XConePlugin_i.hh" + + } // namespace contrib +} // namespace fastjet diff --git a/cpptools/src/fjcontrib/patches/AxesDefinition.patch b/cpptools/src/fjcontrib/patches/AxesDefinition.patch new file mode 100644 index 0000000..3369e3a --- /dev/null +++ b/cpptools/src/fjcontrib/patches/AxesDefinition.patch @@ -0,0 +1,34 @@ +--- /Users/ploskon/tmp/fjcontrib-1.042/Nsubjettiness/AxesDefinition.hh 2019-07-18 05:58:16.000000000 -0700 ++++ /Users/ploskon/devel/heppy/cpptools/src/fjcontrib/fjcontrib-1.042/Nsubjettiness/AxesDefinition.hh 2021-02-28 23:45:07.000000000 -0800 +@@ -93,16 +93,17 @@ + /// not be used for iterative refining (since that is the job of MeasureDefinition). + virtual std::vector get_starting_axes(int n_jets, + const std::vector& inputs, +- const MeasureDefinition * measure) const = 0; ++ const MeasureDefinition * measure) const ++ { std::vector v; return v; } // MP + + /// Short description of AxesDefinitions (and any parameters) +- virtual std::string short_description() const = 0; ++ virtual std::string short_description() const {std::string s = "not implemented"; return s;} // MP + + /// Long description of AxesDefinitions (and any parameters) +- virtual std::string description() const = 0; ++ virtual std::string description() const {std::string s = "not implemented"; return s;} // MP + + /// This has to be defined in all derived classes, and allows these to be copied around. +- virtual AxesDefinition* create() const = 0; ++ virtual AxesDefinition* create() const {return (AxesDefinition*) 0x0;}; // MP + + public: + +@@ -186,7 +187,8 @@ + /// Destructor + virtual ~AxesDefinition() {}; + +-protected: ++// protected: ++ public: // MP + + /// Default constructor contains no information. Number of passes has to be set + /// manually by derived classes using setNPass function. diff --git a/cpptools/src/fjcontrib/patches/MeasureDefinition.patch b/cpptools/src/fjcontrib/patches/MeasureDefinition.patch new file mode 100644 index 0000000..af9d3f2 --- /dev/null +++ b/cpptools/src/fjcontrib/patches/MeasureDefinition.patch @@ -0,0 +1,61 @@ +--- /Users/ploskon/tmp/fjcontrib-1.042/Nsubjettiness/MeasureDefinition.hh 2019-07-18 05:58:16.000000000 -0700 ++++ /Users/ploskon/devel/heppy/cpptools/src/fjcontrib/fjcontrib-1.042/Nsubjettiness/MeasureDefinition.hh 2021-02-28 22:57:02.000000000 -0800 +@@ -80,12 +80,12 @@ + class MeasureDefinition { + + public: +- ++ + /// Description of measure and parameters +- virtual std::string description() const = 0; ++ virtual std::string description() const {std::string s = "not implemented"; return s;} // MP + + /// In derived classes, this should return a copy of the corresponding derived class +- virtual MeasureDefinition* create() const = 0; ++ virtual MeasureDefinition* create() const {return (MeasureDefinition*)0x0;} // MP + + //The following five functions define the measure by which tau_N is calculated, + //and are overloaded by the various measures below +@@ -103,12 +103,12 @@ + } + + /// The jet measure used in N-(sub)jettiness +- virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const = 0; ++ virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {return 0.0;} // MP + /// The beam measure used in N-(sub)jettiness +- virtual double beam_numerator(const fastjet::PseudoJet& particle) const = 0; ++ virtual double beam_numerator(const fastjet::PseudoJet& particle) const {return 0.0;} // MP + + /// A possible normalization factor +- virtual double denominator(const fastjet::PseudoJet& particle) const = 0; ++ virtual double denominator(const fastjet::PseudoJet& particle) const {return 0.0;} // MP + + /// 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 +@@ -145,7 +145,8 @@ + /// virtual destructor + virtual ~MeasureDefinition(){} + +-protected: ++// protected: ++ public: + + /// Flag set by derived classes to choose whether or not to use beam/denominator + TauMode _tau_mode; +@@ -260,7 +261,7 @@ + double accuracy // cap distance of closest approach + ) const; + +-protected: ++ protected: + double _beta; ///< Angular exponent + double _R0; ///< Normalization factor + double _Rcutoff; ///< Cutoff radius +@@ -268,6 +269,7 @@ + DefaultMeasureType _measure_type; ///< Type of measure used (i.e. pp style or ee style) + + ++ public: // MP + /// 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)