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();

  virtual ~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 */