ThePEG/Cuts/FuzzyTheta.h

0001 // -*- C++ -*-
0002 //
0003 // FuzzyTheta.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 // Copyright (C) 2009-2019 Simon Platzer
0006 //
0007 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0008 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0009 //
0010 #ifndef ThePEG_FuzzyTheta_H
0011 #define ThePEG_FuzzyTheta_H
0012 //
0013 // This is the declaration of the FuzzyTheta class.
0014 //
0015
0016 #include "ThePEG/Interface/Interfaced.h"
0017 #include "ThePEG/Repository/EventGenerator.h"
0018 #include <cassert>
0019
0020 namespace ThePEG {
0021
0022   namespace CutTypes {
0023
0024     /**
0025      * Identify an energy-type cut
0026      */
0027     struct Energy {};
0028
0029     /**
0030      * Identify a momentum-type cut
0031      */
0032     struct Momentum {};
0033
0034     /**
0035      * Identify a rapidity-type cut
0036      */
0037     struct Rapidity {};
0038
0039     /**
0040      * Identify an azimuth-type cut
0041      */
0042     struct Azimuth {};
0043
0044     /**
0045      * Identify an polar-angle-type cut
0046      */
0047     struct Polar {};
0048
0049   }
0050
0051 /**
0052  * FuzzyTheta implements fuzzy cut prescriptions
0053  *
0054  * @see \ref FuzzyThetaInterfaces "The interfaces"
0055  * defined for FuzzyTheta.
0056  */
0057 class FuzzyTheta: public Interfaced {
0058
0059 public:
0060
0061   /**
0062    * The default constructor.
0063    */
0064   FuzzyTheta();
0065
0066 public:
0067
0068   /**
0069    * Return the (compact) support of the delta approximation
0070    * considered, given its center value. This default version assumes
0071    * a box approximation. All values are assumed to be in units of the
0072    * width considered.
0073    */
0074   virtual pair<double,double> support(double x) const {
0075     return make_pair(x-0.5,x+0.5);
0076   }
0077
0078   /**
0079    * Return the overlap integral of the delta approximation with the
0080    * given box and center. This default version assumes
0081    * a box approximation. All values are assumed to be in units of the
0082    * width considered.
0083    */
0084   virtual double overlap(double x, const pair<double,double>& box) const {
0085     if ( x - 0.5 >= box.first && x + 0.5 <= box.second )
0086       return 1.;
0087     if ( x - 0.5 > box.second || x + 0.5 < box.first )
0088       return 0.;
0089     return min(box.second,x+0.5) - max(box.first,x-0.5);
0090   }
0091
0092   /**
0093    * Return the overlap, optionally considering absolute lower and
0094    * upper bounds.
0095    */
0096   double overlap(double x,
0097          pair<double,double> box,
0098          const pair<double,double>& support) const {
0099     box.first = max(box.first,support.first);
0100     box.second = min(box.second,support.second);
0101     assert(x >= support.first && x <= support.second);
0102     assert(support.second - support.first >= 1.);
0103     if ( x - 0.5 < support.first )
0104       x = support.first + 0.5;
0105     if ( x + 0.5 > support.second )
0106       x = support.second - 0.5;
0107     return overlap(x,box);
0108   }
0109
0110   /**
0111    * Return the bounds for an energy-type cut
0112    */
0113   pair<double,double> bounds(const CutTypes::Energy&) const {
0114     return make_pair(0.,generator()->maximumCMEnergy()/theEnergyWidth);
0115   }
0116
0117   /**
0118    * Return the width for an energy-type cut
0119    */
0120   Energy width(const CutTypes::Energy&) const {
0121     return theEnergyWidth;
0122   }
0123
0124   /**
0125    * Return the bounds for a momentum-type cut
0126    */
0127   pair<double,double> bounds(const CutTypes::Momentum&) const {
0128     return make_pair(0.,0.5*generator()->maximumCMEnergy()/theEnergyWidth);
0129   }
0130
0131   /**
0132    * Return the width for a momentum-type cut
0133    */
0134   Energy width(const CutTypes::Momentum&) const {
0135     return theEnergyWidth;
0136   }
0137
0138   /**
0139    * Return the bounds for a rapidity-type cut
0140    */
0141   pair<double,double> bounds(const CutTypes::Rapidity&) const {
0142     return make_pair(-Constants::MaxRapidity/theRapidityWidth,
0143              Constants::MaxRapidity/theRapidityWidth);
0144   }
0145
0146   /**
0147    * Return the width for a rapidity-type cut
0148    */
0149   double width(const CutTypes::Rapidity&) const {
0150     return theRapidityWidth;
0151   }
0152
0153   /**
0154    * Return the bounds for a azimuth-type cut
0155    */
0156   pair<double,double> bounds(const CutTypes::Azimuth&) const {
0157     return make_pair(0.0,2.*Constants::pi/theAngularWidth);
0158   }
0159
0160   /**
0161    * Return the width for a azimuth-type cut
0162    */
0163   double width(const CutTypes::Azimuth&) const {
0164     return theAngularWidth;
0165   }
0166
0167   /**
0168    * Return the bounds for a polar-angle-type cut
0169    */
0170   pair<double,double> bounds(const CutTypes::Polar&) const {
0171     return make_pair(0.0,Constants::pi/theAngularWidth);
0172   }
0173
0174   /**
0175    * Return the width for a polar-type cut
0176    */
0177   double width(const CutTypes::Polar&) const {
0178     return theAngularWidth;
0179   }
0180
0181   /**
0182    * Check for value inside the given bounds and update the weight
0183    */
0184   template<class CutType, class Value>
0185   bool isInside(const Value& v, const Value& lower, const Value& upper, double& weight) const {
0186     CutType type;
0187     Value w = width(type);
0188     weight *=
0189       overlap(v/w,pair<double,double>(lower/w,upper/w),bounds(type));
0190     if ( weight == 0.0 )
0191       return false;
0192     return true;
0193   }
0194
0195   /**
0196    * Check for value inside the given bounds and update the weight
0197    */
0198   template<class CutType, class Value>
0199   bool isLessThan(const Value& v, const Value& upper, double& weight) const {
0200     CutType type;
0201     Value w = width(type);
0202     pair<double,double> b = bounds(type);
0203     weight *=
0204       overlap(v/w,pair<double,double>(b.first,upper/w),b);
0205     if ( weight == 0.0 )
0206       return false;
0207     return true;
0208   }
0209
0210   /**
0211    * Check for value inside the given bounds and update the weight
0212    */
0213   template<class CutType, class Value>
0214   bool isLargerThan(const Value& v, const Value& lower, double& weight) const {
0215     CutType type;
0216     Value w = width(type);
0217     pair<double,double> b = bounds(type);
0218     weight *=
0219       overlap(v/w,pair<double,double>(lower/w,b.first),b);
0220     if ( weight == 0.0 )
0221       return false;
0222     return true;
0223   }
0224
0225 public:
0226
0227   /** @name Functions used by the persistent I/O system. */
0228   //@{
0229   /**
0230    * Function used to write out object persistently.
0231    * @param os the persistent output stream written to.
0232    */
0233   void persistentOutput(PersistentOStream & os) const;
0234
0235   /**
0236    * Function used to read in object persistently.
0237    * @param is the persistent input stream read from.
0238    * @param version the version number of the object when written.
0239    */
0240   void persistentInput(PersistentIStream & is, int version);
0241   //@}
0242
0243   /**
0244    * The standard Init function used to initialize the interfaces.
0245    * Called exactly once for each class by the class description system
0246    * before the main function starts or
0247    * when this class is dynamically loaded.
0248    */
0249   static void Init();
0250
0251 protected:
0252
0253   /** @name Clone Methods. */
0254   //@{
0255   /**
0256    * Make a simple clone of this object.
0257    * @return a pointer to the new object.
0258    */
0259   virtual IBPtr clone() const;
0260
0261   /** Make a clone of this object, possibly modifying the cloned object
0262    * to make it sane.
0263    * @return a pointer to the new object.
0264    */
0265   virtual IBPtr fullclone() const;
0266   //@}
0267
0268
0269 // If needed, insert declarations of virtual function defined in the
0270 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0271
0272 private:
0273
0274   /**
0275    * The width to be considered for momenta
0276    */
0277   Energy theEnergyWidth;
0278
0279   /**
0280    * The width to be considered for rapidity quantities
0281    */
0282   double theRapidityWidth;
0283
0284   /**
0285    * The width to be considered for angular quantities
0286    */
0287   double theAngularWidth;
0288
0289 private:
0290
0291   /**
0292    * The static object used to initialize the description of this class.
0293    * Indicates that this is a concrete class with persistent data.
0294    */
0295   static ClassDescription<FuzzyTheta> initFuzzyTheta;
0296
0297   /**
0298    * The assignment operator is private and must never be called.
0299    * In fact, it should not even be implemented.
0300    */
0301   FuzzyTheta & operator=(const FuzzyTheta &) = delete;
0302
0303 };
0304
0305 }
0306
0307 #include "ThePEG/Utilities/ClassTraits.h"
0308
0309 namespace ThePEG {
0310
0311 /** @cond TRAITSPECIALIZATIONS */
0312
0313 /** This template specialization informs ThePEG about the
0314  *  base classes of FuzzyTheta. */
0315 template <>
0316 struct BaseClassTrait<FuzzyTheta,1> {
0317   /** Typedef of the first base class of FuzzyTheta. */
0318   typedef Interfaced NthBase;
0319 };
0320
0321 /** This template specialization informs ThePEG about the name of
0322  *  the FuzzyTheta class and the shared object where it is defined. */
0323 template <>
0324 struct ClassTraits<FuzzyTheta>
0325   : public ClassTraitsBase<FuzzyTheta> {
0326   /** Return a platform-independent class name */
0327   static string className() { return "ThePEG::FuzzyTheta"; }
0328 };
0329
0330 /** @endcond */
0331
0332 }
0333
0334 #endif /* ThePEG_FuzzyTheta_H */