include/HepMC3/Feature.h

0001 // -*- C++ -*-
0002 //
0003 // This file is part of HepMC
0004 // Copyright (C) 2014-2023 The HepMC collaboration (see AUTHORS for details)
0005 //
0006 ///
0007 /// @file Feature.h
0008 /// @brief Defines Feature interface for selecting Particles according to extracted Features.
0009 ///
0010
0011 #ifndef HEPMC3_FEATURE_H
0012 #define HEPMC3_FEATURE_H
0013
0014 #include <functional>
0015 #include <memory>
0016 #include <limits>
0017 #include "HepMC3/GenParticle.h"
0018 #include "HepMC3/Filter.h"
0019
0020
0021 namespace HepMC3 {
0022
0023 //////////////////////////////////////////////////////////////////////
0024
0025 /**
0026  *  @brief GenericFeature defines the Feature interface
0027  *  GenericFeature is not intended to be used directly.  The
0028  *  derived Feature class and its specialisations should be used.
0029  *
0030  *  A Feature wraps a function object that can extract a
0031  *  generic Feature_type from a ConstGenParticlePtr.  Usually the
0032  *  Feature_type would be something like int (e.g. status) or
0033  *  double (e.g. pT), but it could in principle be any attribute of a
0034  *  particle so long as there are well defined <, <=, >, >=, == and
0035  *  != operators for that attribute, as well as an abs function.
0036  *
0037  *  Once a Feature is defined, you can obtain Filters that select
0038  *  Particles according to that Feature by e.g.
0039  *  Feature<int> status([](ConstGenParticlePtr p)->int{return p->status();});
0040  *  bool is_stable = (status == 1)(p);
0041  *  Filter is_beam = (status == 4);
0042  *  bool beam = is_beam(p);
0043  *
0044  *  An abs function is also defined, so abs(Feature) works as you'd
0045  *  expect, e.g.
0046  *  Feature<double> rapidity([](ConstGenParticlePtr p)->double{return p->momentum().rap();});
0047  *  Filter rapCut = abs(rapidity) < 2.5;
0048  *
0049  *  Please also see the Selector interface, which defines an
0050  *  abstract interface to Feature that is free of the template params
0051  *  and also includes some standard Features such as
0052  *
0053  *  Selector::STATUS;
0054  *  Selector::PDG_ID;
0055  *  Selector::PT;
0056  *  Selector::RAPIDITY;
0057  */
0058 template<typename Feature_type>
0059 class GenericFeature {
0060 public:
0061     /// @brief evaluator type
0062     using Evaluator_type = std::function<Feature_type(ConstGenParticlePtr)>;
0063     /// @brief shared pointer for evaluator type
0064     using EvaluatorPtr   =  std::shared_ptr<Evaluator_type>;
0065
0066     /// @brief access the underlying feature value
0067     Feature_type operator()(ConstGenParticlePtr input)const {
0068         return (*m_internal)(input);
0069     }
0070
0071     /// @brief greater than operator
0072     /// @return Filter function
0073     Filter operator > (Feature_type value) const {
0074         EvaluatorPtr functor = m_internal;
0075         return [value, functor](ConstGenParticlePtr input)->bool{return (*functor)(input) >  value;};
0076     }
0077     /// @brief less than operator
0078     /// @return Filter function
0079     Filter operator < (Feature_type value) const {
0080         EvaluatorPtr functor = m_internal;
0081         return [value, functor](ConstGenParticlePtr input)->bool{return (*functor)(input) <  value;};
0082     }
0083
0084     /// @brief greater than or equals operator
0085     /// @return Filter function
0086     Filter operator >= (Feature_type value) const {
0087         EvaluatorPtr functor = m_internal;
0088         return [value, functor](ConstGenParticlePtr input)->bool{return (*functor)(input) >=  value;};
0089     }
0090
0091     /// @brief less than or equals operator
0092     /// @return Filter function
0093     Filter operator <= (Feature_type value) const {
0094         EvaluatorPtr functor = m_internal;
0095         return [value, functor](ConstGenParticlePtr input)->bool{return (*functor)(input) <=  value;};
0096     }
0097
0098     /// @brief equality operator
0099     /// @return Filter function
0100     virtual Filter operator == (Feature_type value) const {
0101         EvaluatorPtr functor = m_internal;
0102         return [value, functor](ConstGenParticlePtr input)->bool{return (*functor)(input) == value;};
0103     }
0104
0105     /// @brief inequality operator
0106     /// @return Filter function
0107     virtual Filter operator != (Feature_type value) const {
0108         EvaluatorPtr functor = m_internal;
0109         return [value, functor](ConstGenParticlePtr input)->bool{return (*functor)(input) != value;};
0110     }
0111
0112 protected:
0113     /// Hide the constructor so no one can use GenericFeature directly
0114     GenericFeature(Evaluator_type functor):m_internal(std::make_shared<Evaluator_type>(functor)) {}
0115
0116     /// Hide the copy constructor
0117     GenericFeature(const GenericFeature &copy) : m_internal(copy.m_internal) {}
0118
0119     /** @brief Move constructor */
0120     GenericFeature(GenericFeature && ) = default;
0121     /** @brief = */
0122     GenericFeature& operator=(const GenericFeature&) = default;
0123     /** @brief = */
0124     GenericFeature& operator=(GenericFeature&&) = default;
0125
0126
0127     /// @brief internal copy of func for evaluation
0128     /// on the heap so will persist in resulting Filters even if
0129     /// parent Feature object was destroyed
0130     EvaluatorPtr m_internal;
0131 };
0132
0133 //////////////////////////////////////////////////////////////////////
0134
0135 /** @brief Expose GenericFeature interface to derived Feature class
0136  *
0137  *  This will get used for generic class types that aren't integral or
0138  *  floating point types.
0139  *
0140  *  A Feature wraps a function object that can extract a
0141  *  generic Feature_type from a ConstGenParticlePtr.  Usually the
0142  *  Feature_type would be something like int (e.g. status) or
0143  *  double (e.g. pT), but it could in principle be any attribute of a
0144  *  particle so long as there are well defined <, <=, >, >=, == and
0145  *  != operators for that attribute, as well as an abs function.
0146  *
0147  *  Once a Feature is defined, you can obtain Filters that select
0148  *  Particles according to that Feature by e.g.
0149  *  Feature<int> status([](ConstGenParticlePtr p)->int{return p->status();});
0150  *  bool is_stable = (status == 1)(p);
0151  *  Filter is_beam = (status == 4);
0152  *  bool beam = is_beam(p);
0153  *
0154  *  An abs function is also defined, so abs(Feature) works as you'd
0155  *  expect, e.g.
0156  *  Feature<double> rapidity([](ConstGenParticlePtr p)->double{return p->momentum().rap();});
0157  *  Filter rapCut = abs(rapidity) < 2.5;
0158  *
0159  *  Please also see the Selector interface, which defines an
0160  *  abstract interface to Feature that is free of the template params
0161  *  and also includes some standard Features such as
0162  *
0163  *  Selector::STATUS;
0164  *  Selector::PDG_ID;
0165  *  Selector::PT;
0166  *  Selector::RAPIDITY;
0167
0168  */
0169 template<typename Feature_type, typename Dummy = void>
0170 class Feature : public GenericFeature<Feature_type> {
0171 public:
0172     using typename GenericFeature<Feature_type>::Evaluator_type;
0173     using typename GenericFeature<Feature_type>::EvaluatorPtr;
0174     using GenericFeature<Feature_type>::m_internal;
0175
0176     using GenericFeature<Feature_type>::operator ();
0177     using GenericFeature<Feature_type>::operator >;
0178     using GenericFeature<Feature_type>::operator >=;
0179     using GenericFeature<Feature_type>::operator <;
0180     using GenericFeature<Feature_type>::operator <=;
0181     using GenericFeature<Feature_type>::operator ==;
0182     using GenericFeature<Feature_type>::operator !=;
0183
0184     /// @brief  Feature
0185     Feature(Evaluator_type functor) : GenericFeature<Feature_type>(functor) {}
0186     /// @brief  Copy
0187     Feature(const Feature &copy) : GenericFeature<Feature_type>(copy) {}
0188     /** @brief Move constructor */
0189     Feature(Feature && ) = default;
0190     /** @brief = */
0191     Feature& operator=(const Feature&) = default;
0192     /** @brief = */
0193     Feature& operator=(Feature&&) = default;
0194
0195     /// @brief  Abs function
0196     Feature<Feature_type> abs() const {
0197         EvaluatorPtr functor = m_internal;
0198         Evaluator_type absfunctor = [functor](ConstGenParticlePtr p)->Feature_type{return ::abs((*functor)(p));};
0199         return Feature<Feature_type>(absfunctor);
0200     }
0201 };
0202
0203 //////////////////////////////////////////////////////////////////////
0204
0205 /** @brief Specialisation of Feature for integral types
0206  *
0207  *  It is a valid operator to compare an int to a float, but the
0208  *  generic version of these operators in the base class will
0209  *  first cast input float to an int, then compare that.  In some cases
0210  *  the comparison will be incorrect because of rounding the float.
0211  *  e.g. int x=5; float y=5.5; bool result = x<y; would be wrong
0212  *  because y first gets converted to int 5.
0213  *
0214  *  To solve this, we provide specialised comparison operators for
0215  *  integral type and double.  Note that the opposite specialisation
0216  *  in which the Feature_type is floating_point is not necessary
0217  */
0218 template<typename Feature_type>
0219 class Feature<Feature_type, typename std::enable_if<std::is_integral<Feature_type>::value, void>::type> : public GenericFeature<Feature_type> {
0220 public:
0221     using GenericFeature<Feature_type>::operator ();
0222     using GenericFeature<Feature_type>::operator >;
0223     using GenericFeature<Feature_type>::operator >=;
0224     using GenericFeature<Feature_type>::operator <;
0225     using GenericFeature<Feature_type>::operator <=;
0226     using GenericFeature<Feature_type>::operator ==;
0227     using GenericFeature<Feature_type>::operator !=;
0228
0229     using typename GenericFeature<Feature_type>::Evaluator_type;
0230     using typename GenericFeature<Feature_type>::EvaluatorPtr;
0231
0232     using GenericFeature<Feature_type>::m_internal;
0233
0234     /// @brief  Feature
0235     Feature(Evaluator_type functor) : GenericFeature<Feature_type>(functor) {}
0236     /// @brief  Feature
0237     Feature(const Feature &copy) : GenericFeature<Feature_type>(copy) {}
0238     /** @brief Move constructor */
0239     Feature(Feature && ) = default;
0240     /** @brief = */
0241     Feature& operator=(const Feature&) = default;
0242     /** @brief = */
0243     Feature& operator=(Feature&&) = default;
0244
0245
0246     /// @brief abs function
0247     Feature<Feature_type> abs() const {
0248         EvaluatorPtr functor = m_internal;
0249         Evaluator_type absfunctor = [functor](ConstGenParticlePtr p)->Feature_type{return ::abs((*functor)(p));};
0250         return Feature<Feature_type>(absfunctor);
0251     }
0252
0253     /// @brief greater operator
0254     Filter operator > (double value) const {
0255         EvaluatorPtr functor = m_internal;
0256         return [value, functor](ConstGenParticlePtr input)->bool{return (*functor)(input) >  value;};
0257     }
0258
0259     /// @brief less operator
0260     Filter operator < (double value) const {
0261         EvaluatorPtr functor = m_internal;
0262         return [value, functor](ConstGenParticlePtr input)->bool{return (*functor)(input) <  value;};
0263     }
0264
0265     /// @brief equal operator
0266     Filter operator == (double value) const {
0267         EvaluatorPtr functor = m_internal;
0268         return [value, functor](ConstGenParticlePtr input)->bool{
0269             Feature_type local = (*functor)(input);
0270             return std::abs(local - value) <  std::numeric_limits<double>::epsilon();
0271         };
0272     }
0273
0274     /// @brief greater or equal operator
0275     Filter operator >= (double value) const { return !( (*this) < value );}
0276
0277     /// @brief less or equal operator
0278     Filter operator <= (double value) const { return !( (*this) > value );}
0279
0280     /// @brief not equal operator
0281     Filter operator != (double value) const {
0282         return !( (*this) == value );
0283     }
0284 };
0285
0286 //////////////////////////////////////////////////////////////////////
0287
0288 /** @brief specialisation of Feature for floating point type
0289  *
0290  *  Test of equality of floating point types is not safe.  Here we
0291  *  provide a "reasonable" definition of equality based on the
0292  *  floating point precision.
0293  */
0294
0295 template<typename Feature_type>
0296 class Feature<Feature_type, typename std::enable_if<std::is_floating_point<Feature_type>::value, void>::type> : public GenericFeature<Feature_type> {
0297 public:
0298     using typename GenericFeature<Feature_type>::Evaluator_type;
0299     using typename GenericFeature<Feature_type>::EvaluatorPtr;
0300
0301     using GenericFeature<Feature_type>::operator ();
0302     using GenericFeature<Feature_type>::operator >;
0303     using GenericFeature<Feature_type>::operator >=;
0304     using GenericFeature<Feature_type>::operator <;
0305     using GenericFeature<Feature_type>::operator <=;
0306
0307     using GenericFeature<Feature_type>::m_internal;
0308
0309     /// @brief Feature
0310     Feature(Evaluator_type functor) : GenericFeature<Feature_type>(functor) {}
0311     /// @brief Copy
0312     Feature(const Feature &copy) : GenericFeature<Feature_type>(copy) {}
0313
0314     /** @brief Move constructor */
0315     Feature(Feature && ) = default;
0316     /** @brief = */
0317     Feature& operator=(const Feature&) = default;
0318     /** @brief = */
0319     Feature& operator=(Feature&&) = default;
0320
0321
0322     /// @brief abs function
0323     Feature<Feature_type> abs() const {
0324         EvaluatorPtr functor = m_internal;
0325         Evaluator_type absfunctor = [functor](ConstGenParticlePtr p)->Feature_type{return std::abs((*functor)(p));};
0326         return Feature<Feature_type>(absfunctor);
0327     }
0328
0329     Filter operator == (Feature_type value) const override {
0330         EvaluatorPtr functor = m_internal;
0331         return [value, functor](ConstGenParticlePtr input)->bool{
0332             Feature_type local = (*functor)(input);
0333             return std::less_equal<Feature_type>{}(fabs(local - value) ,  std::numeric_limits<Feature_type>::epsilon());
0334         };
0335     }
0336
0337     Filter operator != (Feature_type value) const override {
0338         return !( (*this) == value );
0339     }
0340 };
0341
0342 //////////////////////////////////////////////////////////////////////
0343
0344 /**
0345  *  @brief Obtain the absolute value of a Feature.
0346  *  This works as you'd expect.  If foo is a valid Feature, then
0347  *  abs(foo) returns a new Feature that corresponds to the absolute
0348  *  value of the foo feature.  You can construct a Filter from that in
0349  *  the usual way with e.g. Filter f = abs(foo) > 10.;
0350  */
0351 template<typename Feature_type>
0352 Feature<Feature_type> abs(const Feature<Feature_type> &input) {
0353     return input.abs();
0354 }
0355
0356 }
0357
0358 #endif