DalitzIntegrator.h

Go to the documentation of this file.

#ifndef AMPGEN_DALITZINTEGRATOR_H
#define AMPGEN_DALITZINTEGRATOR_H
 
#include <stddef.h>
#include <chrono>
#include <cmath>
#include <complex>
#include <functional>
#include <string>
#include <utility>
 
#include "AmpGen/simd/utils.h"
#include "AmpGen/NumericalIntegration.h"
#include "AmpGen/ProfileClock.h"
 
class TGraph; 
class TH1D;
class TH2D;
 
namespace AmpGen
{
  class Event;
  class Projection2D;
  class Projection;
 
  class DalitzIntegrator
  {
    public:
 
      typedef std::pair<real_v, real_v> sqCo;
      DalitzIntegrator( const double& s0, const double& s1, const double& s2, const double& s3);
 
      template <typename FCN> double integrateDP( FCN&& fcn, const double& s) const;
      template <typename FCN> double integrateDP( FCN&& fcn ) const
      {
        return integrateDP( fcn, m_s0 );
      }
      real_v getMAB(const sqCo& coords)   const;
      real_v J(const sqCo& coords) const;
      real_v getMAB(const sqCo& coords  , const double& s) const;
      real_v J(const sqCo& coords, const double& s) const;
      double sqDp1(const Event& evt) const;
      double sqDp2(const Event& evt) const;
      real_v safe_sqrt(const real_v& x ) const { 
        #if ENABLE_AVX 
        return select( x > 0. , sqrt(x) , real_v(0.) ); 
        #else 
          return x > 0 ? std::sqrt(x) : 0;
        #endif
      }
      void setEvent(const sqCo& x, real_v* event, const double& s) const;
      
      void debug() const; 
      void setEvent(const sqCo& x, real_v* event) const;
       
      
      void set(const double& s0, const double& s1, const double& s2, const double& s3);
      void setMin();
      void setMother(const double& s);
 
      TH1D* makePlot( const std::function<double(const double*)>& fcn, const Projection& projection,
          const std::string& name, const size_t& nSamples = 1000000 );
 
      TH2D* makePlot( const std::function<double(const double*)>& fcn, const Projection2D& projection,
          const std::string& name, const size_t& nSamples = 1000000 );
      
      sqCo getCoordinates( const Event& evt ) const;
 
 
      TGraph* makeBoundaryGraph( const Projection2D& ) const; 
    private:
      double    m_min;
      double    m_max;
      double    m_s0;
      double    m_s1;
      double    m_s2;
      double    m_s3;
   
  };
  template <typename FCN>
  double DalitzIntegrator::integrateDP( FCN&& fcn, const double& s) const
  {
    #if INSTRUCTION_SET != 0 && INSTRUCTION_SET != INSTRUCTION_SET_AVX2d
      #pragma message("WARNING: DalitzIntegrator only supports scalar or AVX2(d) instruction sets")
    #else 
    real_v event[12] = {0.};
    ProfileClock pc1; 
    auto i1 = integrate<2>(  [&]( const std::array<real_v,2>& x ) {
      setEvent( *reinterpret_cast<const sqCo*>(&x) , event, s);
      return J( *reinterpret_cast<const sqCo*>(&x) , s) * real( fcn(event) ); }, std::array<double,2>{0., 0.}, std::array<double, 2>{1.,1.} ) /s;
    return i1; 
    #endif 
    return 0; 
  }
 
 
} // namespace AmpGen
 
#endif