FuzzyTheta.h

// -*- C++ -*-
//
// FuzzyTheta.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
// Copyright (C) 2009-2019 Simon Platzer
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_FuzzyTheta_H
#define ThePEG_FuzzyTheta_H
//
// This is the declaration of the FuzzyTheta class.
//
 
#include "ThePEG/Interface/Interfaced.h"
#include "ThePEG/Repository/EventGenerator.h"
#include <cassert>
 
namespace ThePEG {
 
  namespace CutTypes {
 
    struct Energy {};
 
    struct Momentum {};
 
    struct Rapidity {};
 
    struct Azimuth {};
 
    struct Polar {};
 
  }
 
class FuzzyTheta: public Interfaced {
 
public:
 
  FuzzyTheta();
 
public:
 
  virtual pair<double,double> support(double x) const {
    return make_pair(x-0.5,x+0.5);
  }
 
  virtual double overlap(double x, const pair<double,double>& box) const {
    if ( x - 0.5 >= box.first && x + 0.5 <= box.second )
      return 1.;
    if ( x - 0.5 > box.second || x + 0.5 < box.first )
      return 0.;
    return min(box.second,x+0.5) - max(box.first,x-0.5);
  }
 
  double overlap(double x,
                 pair<double,double> box,
                 const pair<double,double>& support) const {
    box.first = max(box.first,support.first);
    box.second = min(box.second,support.second);
    assert(x >= support.first && x <= support.second);
    assert(support.second - support.first >= 1.);
    if ( x - 0.5 < support.first )
      x = support.first + 0.5;
    if ( x + 0.5 > support.second )
      x = support.second - 0.5;
    return overlap(x,box);
  }
 
  pair<double,double> bounds(const CutTypes::Energy&) const {
    return make_pair(0.,generator()->maximumCMEnergy()/theEnergyWidth);
  }
 
  Energy width(const CutTypes::Energy&) const {
    return theEnergyWidth;
  }
 
  pair<double,double> bounds(const CutTypes::Momentum&) const {
    return make_pair(0.,0.5*generator()->maximumCMEnergy()/theEnergyWidth);
  }
 
  Energy width(const CutTypes::Momentum&) const {
    return theEnergyWidth;
  }
 
  pair<double,double> bounds(const CutTypes::Rapidity&) const {
    return make_pair(-Constants::MaxRapidity/theRapidityWidth,
                     Constants::MaxRapidity/theRapidityWidth);
  }
 
  double width(const CutTypes::Rapidity&) const {
    return theRapidityWidth;
  }
 
  pair<double,double> bounds(const CutTypes::Azimuth&) const {
    return make_pair(0.0,2.*Constants::pi/theAngularWidth);
  }
 
  double width(const CutTypes::Azimuth&) const {
    return theAngularWidth;
  }
 
  pair<double,double> bounds(const CutTypes::Polar&) const {
    return make_pair(0.0,Constants::pi/theAngularWidth);
  }
 
  double width(const CutTypes::Polar&) const {
    return theAngularWidth;
  }
 
  template<class CutType, class Value>
  bool isInside(const Value& v, const Value& lower, const Value& upper, double& weight) const {
    CutType type;
    Value w = width(type);
    weight *=
      overlap(v/w,pair<double,double>(lower/w,upper/w),bounds(type));
    if ( weight == 0.0 )
      return false;
    return true;
  }
 
  template<class CutType, class Value>
  bool isLessThan(const Value& v, const Value& upper, double& weight) const {
    CutType type;
    Value w = width(type);
    pair<double,double> b = bounds(type);
    weight *=
      overlap(v/w,pair<double,double>(b.first,upper/w),b);
    if ( weight == 0.0 )
      return false;
    return true;
  }
 
  template<class CutType, class Value>
  bool isLargerThan(const Value& v, const Value& lower, double& weight) const {
    CutType type;
    Value w = width(type);
    pair<double,double> b = bounds(type);
    weight *=
      overlap(v/w,pair<double,double>(lower/w,b.first),b);
    if ( weight == 0.0 )
      return false;
    return true;
  }
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
 
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
 
private:
 
  Energy theEnergyWidth;
 
  double theRapidityWidth;
 
  double theAngularWidth;
 
private:
 
  static ClassDescription<FuzzyTheta> initFuzzyTheta;
 
  FuzzyTheta & operator=(const FuzzyTheta &) = delete;
 
};
 
}
 
#include "ThePEG/Utilities/ClassTraits.h"
 
namespace ThePEG {
 
template <>
struct BaseClassTrait<FuzzyTheta,1> {
  typedef Interfaced NthBase;
};
 
template <>
struct ClassTraits<FuzzyTheta>
  : public ClassTraitsBase<FuzzyTheta> {
  static string className() { return "ThePEG::FuzzyTheta"; }
};
 
}
 
#endif /* ThePEG_FuzzyTheta_H */