ExternalData.cc

Go to the documentation of this file.

///
/// @file  ExternalData.cc
/// @brief Container class for values of WC and FF vectors
///

//**** This file is a part of the HAMMER library
//**** Copyright (C) 2016 - 2020 The HAMMER Collaboration
//**** HAMMER is licensed under version 3 of the GPL; see COPYING for details
//**** Please note the MCnet academic guidelines; see GUIDELINES for details

// -*- C++ -*-
#include <iostream>
#include <algorithm>
#include <boost/algorithm/string.hpp>

#include "Hammer/Math/Tensor.hh"
#include "Hammer/ProvidersRepo.hh"
#include "Hammer/ExternalData.hh"
#include "Hammer/Tools/HammerSerial.hh"
#include "Hammer/AmplitudeBase.hh"
#include "Hammer/Config/HammerModules.hh"
#include "Hammer/FormFactorBase.hh"
#include "Hammer/Tools/ParticleUtils.hh"
#include "Hammer/Tools/Logging.hh"
#include "Hammer/Tools/Utils.hh"
#include "Hammer/Math/Histogram.hh"
#include "Hammer/Math/MultiDim/ScalarContainer.hh"

using namespace std;

namespace Hammer {

    namespace MD = MultiDimensional;

    ExternalData::ExternalData(const ProvidersRepo* provs) : _providers{provs} {
        _processWilsonCoefficients.reset(new processWilsonCoefficientsType{});
        _processFFEigenVectors.reset(new processFFEigenVectorsType{});
    }

    ExternalData::~ExternalData() noexcept {
        _processSpecializedWilsonCoefficients.clear();
        if (_processWilsonCoefficients.get() != nullptr) {
            _processWilsonCoefficients.reset();
        }
        if (_processFFEigenVectors.get() != nullptr) {
            _processFFEigenVectors.reset();
        }
        if (_externalVectors.get() != nullptr) {
            _externalVectors.reset();
        }
        _schemes.clear();
    }

    void ExternalData::defineSettings() {
    }

    const Tensor& ExternalData::getExternalVectors(string schemeName, LabelsList labels) const {
        // sort(labels.begin(), labels.end(), std::greater<IndexLabel>());
        if(_externalVectors.get() == nullptr) {
            _externalVectors.reset(new externalVectorsType{});
        }
        auto& ext = (*_externalVectors)[schemeName];
        auto ite = ext.find(labels);
        if (ite != ext.end()) {
            return ite->second;
        }
        if(labels.size() == 0) {
            auto res = ext.insert({labels, Tensor{"External", MD::makeScalar(1.)}});
            return res.first->second;
        }
        vector<pair<MD::SharedTensorData, bool>> elements;
        elements.reserve(labels.size());
        size_t index = (schemeName == "Denominator") ? 1 : 0;
        for (auto elem : labels) {
            bool hc = elem < 0;
            IndexLabel lab = (elem > 0) ? elem : static_cast<IndexLabel>(-elem);
            if(_processWilsonCoefficients.get() != nullptr) {
                auto it = _processWilsonCoefficients->find(lab);
                if (it != _processWilsonCoefficients->end()) {
                    elements.push_back({it->second[index], hc});
                    continue;
                }
            }
            if (_processFFEigenVectors.get() != nullptr) {
                auto it2 = _processFFEigenVectors->find(lab);
                if (it2 != _processFFEigenVectors->end()) {
                    auto it3 = it2->second.find(schemeName);
                    if (it3 != it2->second.end()) {
                        elements.push_back({it3->second, hc});
                        continue;
                    }
                }
            }
        }
        auto res = ext.emplace(labels, Tensor{"External", move(elements)});
        return res.first->second;
    }

    MD::SharedTensorData ExternalData::getWilsonCoefficients(PdgId parent, const std::vector<PdgId>& daughters,
                                const std::vector<PdgId>& granddaughters, WTerm what) const {
        AmplitudeBase* amp = _providers->getAmplitude(parent, daughters, granddaughters);
        return getWilsonCoefficients(amp, what);
    }

    MD::SharedTensorData ExternalData::getWilsonCoefficients(AmplitudeBase* ampl, WTerm what) const {
        if(_processWilsonCoefficients.get() != nullptr) {
            auto it = _processWilsonCoefficients->find(ampl->getWCInfo().second);
            if (it != _processWilsonCoefficients->end()) {
                switch (what) {
                case WTerm::NUMERATOR:
                    return it->second[0];
                case WTerm::DENOMINATOR:
                    return it->second[1];
                case WTerm::COMMON:
                    throw Error("Invalid option");
                }
            }
        }
        return nullptr;
    }

    MD::SharedTensorData ExternalData::getSpecializedWilsonCoefficients(AmplitudeBase* ampl) const {
        auto it = _processSpecializedWilsonCoefficients.find(ampl->getWCInfo().second);
        if (it != _processSpecializedWilsonCoefficients.end()) {
            return it->second;
        }
        return nullptr;
    }

    void ExternalData::specializeWCInWeights(const string& prefixName, const map<string, complex<double>>& settings){
        auto ampl = _providers->getWCProvider(prefixName);
        if (ampl != nullptr) {
            ampl->updateWCSettings(settings, WTerm::NUMERATOR);
            auto vec = ampl->getWCVectorFromSettings(WTerm::NUMERATOR);
            specializeWCInWeights(prefixName, vec);
        } else {
             MSG_ERROR("Wilson coefficients " + prefixName + " not found.");
        }
    }

    void ExternalData::specializeWCInWeights(const string& prefixName, const vector<complex<double>>& values){
        auto ampl = _providers->getWCProvider(prefixName);
        if (ampl != nullptr) {
            _isWCSpecializedDict[prefixName] = true;
            auto label = _providers->getWCLabel(prefixName);
            auto itSWC = _processSpecializedWilsonCoefficients.find(label);
            ASSERT(itSWC != _processSpecializedWilsonCoefficients.end());
            ampl->updateWCSettings(values, WTerm::NUMERATOR);
            ampl->updateWCTensor(values, itSWC->second);
            MSG_INFO(prefixName + " Wilson coefficients will be specialized in weights. Histogram specialization will be superceded!");
        } else {
             MSG_ERROR("Wilson coefficients " + prefixName + " not found.");
        }

    }

    void ExternalData::resetSpecializeWCInWeights(const string& prefixName){
        auto ampl = _providers->getWCProvider(prefixName);
        if (ampl != nullptr) {
            _isWCSpecializedDict[prefixName] = false;
        } else {
             MSG_ERROR("Wilson coefficients " + prefixName + " not found.");
        }
    }

    bool ExternalData::isWCSpecialized(AmplitudeBase* ampl) const {
        auto it = _isWCSpecializedDict.find(ampl->getWCInfo().first);
        if(it != _isWCSpecializedDict.end()){
            return it->second;
        }
        return false;
    }

    void ExternalData::setWilsonCoefficients(const string& prefixName, const vector<complex<double>>& values,
                                             WTerm what) {
        auto ampl = _providers->getWCProvider(prefixName);
        if (ampl != nullptr) {
            auto label = _providers->getWCLabel(prefixName);
            ASSERT(_processWilsonCoefficients.get() != nullptr);
            auto itWC = _processWilsonCoefficients->find(label);
            ASSERT(itWC != _processWilsonCoefficients->end());
            ampl->updateWCSettings(values, what);
            switch(what) {
            case WTerm::NUMERATOR:
                ampl->updateWCTensor(values, itWC->second[0]);
                break;
            case WTerm::DENOMINATOR:
                ampl->updateWCTensor(values, itWC->second[1]);
                break;
            case WTerm::COMMON:
                ampl->updateWCTensor(values, itWC->second[0]);
                ampl->updateWCTensor(values, itWC->second[1]);
                break;
            }
        }
    }

    void ExternalData::setWilsonCoefficients(const string& prefixName, const map<string, complex<double>>& values, WTerm what) {
        auto ampl = _providers->getWCProvider(prefixName);
        if (ampl != nullptr) {
            auto label = _providers->getWCLabel(prefixName);
            ASSERT(_processWilsonCoefficients.get() != nullptr);
            auto itWC = _processWilsonCoefficients->find(label);
            ASSERT(itWC != _processWilsonCoefficients->end());
            ampl->updateWCSettings(values, what);
            auto vec = ampl->getWCVectorFromSettings(what);
            switch (what) {
            case WTerm::NUMERATOR:
                ampl->updateWCTensor(vec, itWC->second[0]);
                break;
            case WTerm::DENOMINATOR:
                ampl->updateWCTensor(vec, itWC->second[1]);
                break;
            case WTerm::COMMON:
                ampl->updateWCTensor(vec, itWC->second[0]);
                ampl->updateWCTensor(vec, itWC->second[1]);
                break;
            }
        }
    }

    void ExternalData::setWilsonCoefficientsLocal(const std::string& prefixName,
                                                  const std::vector<std::complex<double>>& values) {
        if(_processWilsonCoefficients.get() == nullptr) {
            initWilsonCoefficients();
        }
        auto ampl = _providers->getWCProvider(prefixName);
        if (ampl != nullptr) {
            auto label = _providers->getWCLabel(prefixName);
            ASSERT(_processWilsonCoefficients.get() != nullptr);
            auto itWC = _processWilsonCoefficients->find(label);
            ASSERT(itWC != _processWilsonCoefficients->end());
            ampl->updateWCTensor(values, itWC->second[0]);
        }
    }

    void ExternalData::setWilsonCoefficientsLocal(const std::string& prefixName,
                                                  const std::map<std::string, std::complex<double>>& values) {
        if(_processWilsonCoefficients.get() == nullptr) {
            initWilsonCoefficients();
        }
        auto ampl = _providers->getWCProvider(prefixName);
        if (ampl != nullptr) {
            auto label = _providers->getWCLabel(prefixName);
            ASSERT(_processWilsonCoefficients.get() != nullptr);
            auto itWC = _processWilsonCoefficients->find(label);
            ASSERT(itWC != _processWilsonCoefficients->end());
            auto vec = ampl->getWCVectorFromDict(values);
            vec[0] = 1.;
            ampl->updateWCTensor(vec, itWC->second[0]);
        }
    }

    void ExternalData::resetWilsonCoefficients(const string& prefixName, WTerm what) {
        auto ampl = _providers->getWCProvider(prefixName);
        if (ampl != nullptr) {
            auto label = _providers->getWCLabel(prefixName);
            ASSERT(_processWilsonCoefficients.get() != nullptr);
            auto itWC = _processWilsonCoefficients->find(label);
            ASSERT(itWC != _processWilsonCoefficients->end());
            auto vec = ampl->getWCVectorFromSettings(what);
            auto n = vec.size();
            vec.clear();
            vec.resize(n);
            vec[0]=1.0;
            ampl->updateWCSettings(vec, what);
            switch (what) {
            case WTerm::NUMERATOR:
                ampl->updateWCTensor(vec, itWC->second[0]);
                break;
            case WTerm::DENOMINATOR:
                ampl->updateWCTensor(vec, itWC->second[1]);
                break;
            case WTerm::COMMON:
                ampl->updateWCTensor(vec, itWC->second[0]);
                ampl->updateWCTensor(vec, itWC->second[1]);
                break;
            }
        }
    }

    MD::SharedTensorData ExternalData::getTempWilsonCoefficients(const string& prefixName,
                                                   const vector<complex<double>>& values) const {
        auto ampl = _providers->getWCProvider(prefixName);
        MD::SharedTensorData t{new MD::ScalarContainer{}};
        if (ampl != nullptr) {
            ampl->updateWCTensor(values, t);
        }
        return t;
    }

    MD::SharedTensorData ExternalData::getTempWilsonCoefficients(const string& prefixName,
                                                   const map<string, complex<double>>& values) const {
        auto ampl = _providers->getWCProvider(prefixName);
        MD::SharedTensorData t{new MD::ScalarContainer{}};
        if (ampl != nullptr) {
            auto vec = ampl->getWCVectorFromDict(values);
            ampl->updateWCTensor(vec, t);
        }
        return t;
    }

    MD::SharedTensorData ExternalData::getFFEigenVectors(FormFactorBase* ff, const string& schemeName) const {
        ASSERT(_processFFEigenVectors.get() != nullptr);
        auto it = _processFFEigenVectors->find(ff->getFFErrInfo().second);
        if (it != _processFFEigenVectors->end()) {
            auto it2 = it->second.find(schemeName);
            if(it2 != it->second.end()) {
                return it2->second;
            }
        }
        return nullptr;
    }


    void ExternalData::setFFEigenVectors(const FFPrefixGroup& process,
                                         const vector<double>& values) {
        auto ff = _providers->getFFErrProvider(process);
        if (ff != nullptr) {
            auto label = _providers->getFFErrLabel(process);
            ASSERT(_processFFEigenVectors.get() != nullptr);
            auto itFFE1 = _processFFEigenVectors->find(label);
            ASSERT(itFFE1 != _processFFEigenVectors->end());
            ff->updateFFErrSettings(values);
            auto schemes = _providers->schemeNamesFromPrefixAndGroup(process);
            vector<double> tmpvalues;
            tmpvalues.reserve(values.size() + 1);
            tmpvalues.push_back(1.);
            tmpvalues.insert(tmpvalues.end(), values.begin(),values.end());
            for(auto& elem: schemes) {
                auto itFFE2 = itFFE1->second.find(elem);
                if(itFFE2 != itFFE1->second.end()) {
                    ff->updateFFErrTensor(tmpvalues, itFFE2->second);
                }
            }
        }
    }

    void ExternalData::setFFEigenVectors(const FFPrefixGroup& process, const map<string, double>& values) {
        auto ff = _providers->getFFErrProvider(process);
        if (ff != nullptr) {
            auto label = _providers->getFFErrLabel(process);
            ASSERT(_processFFEigenVectors.get() != nullptr);
            auto itFFE1 = _processFFEigenVectors->find(label);
            ASSERT(itFFE1 != _processFFEigenVectors->end());
            ff->updateFFErrSettings(values);
            auto vec = ff->getErrVectorFromSettings();
            auto schemes = _providers->schemeNamesFromPrefixAndGroup(process);
            for (auto& elem : schemes) {
                auto itFFE2 = itFFE1->second.find(elem);
                if (itFFE2 != itFFE1->second.end()) {
                    ff->updateFFErrTensor(vec, itFFE2->second);
                }
            }
        }
    }

    void ExternalData::setFFEigenVectorsLocal(const FFPrefixGroup& process, const vector<double>& values) {
        if (_processFFEigenVectors.get() == nullptr) {
            initFormFactorErrors();
        }
        auto ff = _providers->getFFErrProvider(process);
        if (ff != nullptr) {
            auto label = _providers->getFFErrLabel(process);
            ASSERT(_processFFEigenVectors.get() != nullptr);
            auto itFFE1 = _processFFEigenVectors->find(label);
            ASSERT(itFFE1 != _processFFEigenVectors->end());
            auto schemes = _providers->schemeNamesFromPrefixAndGroup(process);
            vector<double> tmpvalues;
            tmpvalues.reserve(values.size() + 1);
            tmpvalues.push_back(1.);
            tmpvalues.insert(tmpvalues.end(), values.begin(), values.end());
            for (auto& elem : schemes) {
                auto itFFE2 = itFFE1->second.find(elem);
                if (itFFE2 != itFFE1->second.end()) {
                    ff->updateFFErrTensor(tmpvalues, itFFE2->second);
                }
            }
        }
    }

    void ExternalData::setFFEigenVectorsLocal(const FFPrefixGroup& process, const map<string, double>& values) {
        if (_processFFEigenVectors.get() == nullptr) {
            initFormFactorErrors();
        }
        auto ff = _providers->getFFErrProvider(process);
        if (ff != nullptr) {
            auto label = _providers->getFFErrLabel(process);
            ASSERT(_processFFEigenVectors.get() != nullptr);
            auto itFFE1 = _processFFEigenVectors->find(label);
            ASSERT(itFFE1 != _processFFEigenVectors->end());
            auto vec = ff->getErrVectorFromDict(values);
            auto schemes = _providers->schemeNamesFromPrefixAndGroup(process);
            for (auto& elem : schemes) {
                auto itFFE2 = itFFE1->second.find(elem);
                if (itFFE2 != itFFE1->second.end()) {
                    ff->updateFFErrTensor(vec, itFFE2->second);
                }
            }
        }
    }


    void ExternalData::resetFFEigenVectors(const FFPrefixGroup& process) {
        auto ff = _providers->getFFErrProvider(process);
        if (ff != nullptr) {
            auto label = _providers->getFFErrLabel(process);
            ASSERT(_processFFEigenVectors.get() != nullptr);
            auto itFFE1 = _processFFEigenVectors->find(label);
            ASSERT(itFFE1 != _processFFEigenVectors->end());
            auto n = ff->getErrVectorFromSettings().size();
            vector<double> tmpvalues(n-1);
            tmpvalues.reserve(n);
            ff->updateFFErrSettings(tmpvalues);
            tmpvalues.insert(tmpvalues.begin(),1.);
            auto schemes = _providers->schemeNamesFromPrefixAndGroup(process);
            for(auto& elem: schemes) {
                auto itFFE2 = itFFE1->second.find(elem);
                if(itFFE2 != itFFE1->second.end()) {
                    ff->updateFFErrTensor(tmpvalues, itFFE2->second);
                }
            }
        }
    }

    MD::SharedTensorData ExternalData::getTempFFEigenVectors(const FFPrefixGroup& process,
                                               const vector<double>& values) const {
        auto ff = _providers->getFFErrProvider(process);
        MD::SharedTensorData t{new MD::ScalarContainer{}};
        if (ff != nullptr) {
            ff->updateFFErrTensor(values, t);
        }
        return t;
    }

    MD::SharedTensorData ExternalData::getTempFFEigenVectors(const FFPrefixGroup& process,
                                               const map<string, double>& values) const {
        auto ff = _providers->getFFErrProvider(process);
        MD::SharedTensorData t{new MD::ScalarContainer{}};
        if (ff != nullptr) {
            auto vec = ff->getErrVectorFromDict(values);
            ff->updateFFErrTensor(vec, t);
        }
        return t;
    }


    void ExternalData::initWilsonCoefficients() {
        _processWilsonCoefficients.reset(new processWilsonCoefficientsType{});
        _processSpecializedWilsonCoefficients.clear();
        for (auto& elem : _providers->getAllWCProviders()) {
            auto res = _processWilsonCoefficients->emplace(elem.first, array<MD::SharedTensorData, 2>{});
            if (res.second) {
                auto vecN = elem.second->getWCVectorFromSettings(WTerm::NUMERATOR);
                elem.second->updateWCTensor(vecN, res.first->second[0]);
                auto vecD = elem.second->getWCVectorFromSettings(WTerm::DENOMINATOR);
                elem.second->updateWCTensor(vecD, res.first->second[1]);
            }
        }
        for (auto& elem: *_processWilsonCoefficients){
            _processSpecializedWilsonCoefficients.insert({elem.first, elem.second[0]});
        }
    }

    void ExternalData::initFormFactorErrors() {
        _processFFEigenVectors.reset(new processFFEigenVectorsType{});
        for (auto& elem : _providers->getAllFFErrProviders()) {
            auto res = _processFFEigenVectors->emplace(elem.first, SchemeDict<MD::SharedTensorData>{});
            if (res.second) {
                for (auto& elem2 : elem.second) {
                    auto res2 = res.first->second.emplace(elem2.first, MD::SharedTensorData{});
                    if (res2.second) {
                        bool useDefault = (elem2.first == "Denominator");
                        auto vecN = elem2.second->getErrVectorFromSettings(useDefault);
                        elem2.second->updateFFErrTensor(vecN, res2.first->second);
                    }
                }
            }
        }
    }

    void ExternalData::initExternalVectors() {
        _externalVectors.reset(new externalVectorsType{});
        _externalVectors->reserve(_schemes.size() + 1);
        for (auto& elem : _schemes) {
            _externalVectors->insert({elem, UMap<LabelsList, Tensor>{}});
        }
        _externalVectors->insert({"Denominator", UMap<LabelsList, Tensor>{}});
    }

    void ExternalData::init(vector<string> schemeNames) {
        SettingsHandler* settings = getSettingsHandler();
        _schemes = schemeNames;
        initExternalVectors();
        if (settings != nullptr) {
            initWilsonCoefficients();
            initFormFactorErrors();
        }
    }


    Log& ExternalData::getLog() const {
        return Log::getLog("Hammer.ExternalData");
    }


} // namespace Hammer