Integrator.h

Go to the documentation of this file.

#ifndef AMPGEN_INTEGRATOR_H
#define AMPGEN_INTEGRATOR_H
 
#include "AmpGen/Types.h"
#include "AmpGen/EventList.h"
#include <array>
#include <complex>
#include "AmpGen/simd/utils.h"
#include "AmpGen/Store.h"
#include "AmpGen/EventListSIMD.h"
#include "AmpGen/EventList.h"
 
namespace AmpGen
{ 
  class Integrator
  {
    #if ENABLE_AVX 
    using EventList_t = EventListSIMD;
    #else 
    using EventList_t = EventList; 
    #endif
    struct QueuedIntegral 
    {
      QueuedIntegral() = default; 
      QueuedIntegral(complex_t* result, const unsigned& i, const unsigned& j) 
        : result(result), i(i), j(j) {}
      complex_t* result = {nullptr};
      unsigned i = {0};
      unsigned j = {0};
    };
    public:
    Integrator() = default; 
 
    template <typename EventList_type, typename T> Integrator( const EventList_type* events, const std::vector<T>& expressions ={}) : m_events(events)
    {
      if( events == nullptr ) {
        WARNING("No events specified, returning");
        return; 
      }
      m_cache.allocate(events, expressions); 
      m_weight.resize(events->nBlocks() );
      real_v norm_acc = 0.;
      for( size_t i = 0 ; i < events->nBlocks(); ++i )
      {
        m_weight[i] = events->weight(i) / events->genPDF(i);
        norm_acc   += m_weight[i]; 
      }
      m_norm = utils::sum_elements(norm_acc);
    }
 
    bool isReady()            const;
    void queueIntegral( complex_t* result, const unsigned& i, const unsigned& j );
    void flush();
    
    template <class return_type> return_type get( const unsigned& index, const unsigned& evt ) const ;
    template <class T> unsigned getCacheIndex( const T& t ) const { return m_cache.find(t.name())[0] ; }
    double norm() const { return m_norm; }
 
    template <class T> void updateCache(const T& expression){ if( isReady() ) m_cache.update(expression); }
    template <class T> const T* events() const { return static_cast<const T*>(m_events); }
 
    const auto& cache() const { return m_cache; }
    private:
    static constexpr size_t             N         = {8}; 
    size_t                              m_counter = {0};  
    std::array<QueuedIntegral, N>       m_integrals;
    const void*                         m_events  = {nullptr};
    std::vector<real_v>                 m_weight; 
    FunctionCache<EventList_t, complex_v, Alignment::SoA>    m_cache;      
    double                              m_norm    = {0};
    void integrateBlock();
  };
  
  class Bilinears 
  {
    private: 
      size_t rows;
      size_t cols;
      std::vector<complex_t> norms;
      std::vector<bool> markAsZero; 
      std::vector<bool> calculate;   
    public:
      Bilinears( const size_t& r = 0, const size_t& c = 0 );
      complex_t get(const size_t& x, const size_t& y) const;
      complex_t get(const size_t& x, const size_t& y, Integrator* integ = nullptr, const size_t& kx=0, const size_t& ky=0){
        if( integ != nullptr ) integ->queueIntegral(&norms[x*cols+y], kx, ky );
        return norms[x*cols+y];
      }
      void      set(const size_t& x, const size_t& y,  const complex_t& f );
      void  setZero(const size_t& x, const size_t& y);
      void resetCalculateFlags();
      complex_t& operator()( const size_t& x, const size_t& y );
      bool   isZero(const size_t& x, const size_t& y);
      bool workToDo(const size_t& x, const size_t& y) const;
      void   resize(const size_t& r, const size_t& c = 1 );
  };
 
} // namespace AmpGen
#endif