MoMEMta.h

/*
 *  MoMEMta: a modular implementation of the Matrix Element Method
 *  Copyright (C) 2016  Universite catholique de Louvain (UCL), Belgium
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */


#pragma once

#include <memory>
#include <vector>

#include <momemta/Module.h>
#include <momemta/ParameterSet.h>
#include <momemta/Particle.h>
#include <momemta/Pool.h>
#include <momemta/Types.h>

class Configuration;
class SharedLibrary;

namespace momemta {
class ComputationGraph;
}

class MoMEMta {
    public:
        enum class IntegrationStatus {
            SUCCESS, 
            ACCURACY_NOT_REACHED, 
            FAILED, 
            ABORTED, 
            DIM_OUT_OF_RANGE, 
            NONE 
        };

        MoMEMta(const Configuration& configuration);
        virtual ~MoMEMta();

        std::vector<std::pair<double, double>> computeWeights(const std::vector<momemta::Particle>& particles,
                                                              const LorentzVector& met=LorentzVector());

        void setEvent(const std::vector<momemta::Particle>& particles, const LorentzVector& met=LorentzVector());
        
        std::vector<double> evaluateIntegrand(const std::vector<double>& psPoints);
        
        IntegrationStatus getIntegrationStatus() const;

        const Pool& getPool() const;

    private:
        class integrands_output_error: public std::runtime_error {
            using std::runtime_error::runtime_error;
        };
        class cuba_configuration_error: public std::runtime_error {
            using std::runtime_error::runtime_error;
        };
        class unphysical_lorentzvector_error: public std::exception {
        public:
            unphysical_lorentzvector_error(const LorentzVector& p4);
            virtual const char* what() const noexcept override;
        private:
            std::string _what;
        };
        class invalid_inputs: public std::runtime_error {
            using std::runtime_error::runtime_error;
        };
        class integrands_nonfinite_error: public std::runtime_error {
            using std::runtime_error::runtime_error;
        };

        void checkIfPhysical(const LorentzVector& p4);

        void initPool(const Configuration& configuration);

        int integrand(const double* psPoints, double* results, const double* weights=nullptr);

        static int CUBAIntegrand(const int *nDim, const double* psPoint, const int *nComp, double *value, void *inputs, const int *nVec, const int *core);
        static int CUBAIntegrandWeighted(const int *nDim, const double* psPoint, const int *nComp, double *value, void *inputs, const int *nVec, const int *core, const double *weight);
        static void cuba_logging(const char*);

        PoolPtr m_pool;
        std::shared_ptr<momemta::ComputationGraph> m_computation_graph;

        using SharedLibraryPtr = std::shared_ptr<SharedLibrary>;
        std::vector<SharedLibraryPtr> m_libraries;

        std::size_t m_n_dimensions;
        std::size_t m_n_components;
        ParameterSet m_cuba_configuration;

        IntegrationStatus integration_status = IntegrationStatus::NONE;

        // Pool inputs
        std::shared_ptr<std::vector<double>> m_ps_points;
        std::shared_ptr<double> m_ps_weight;

        std::unordered_map<std::string, std::shared_ptr<LorentzVector>> m_inputs_p4;
        std::unordered_map<std::string, std::shared_ptr<int64_t>> m_inputs_type;
        std::shared_ptr<LorentzVector> m_met;
        std::vector<Value<double>> m_integrands;
};