ProcRates.cc

Go to the documentation of this file.

///
/// @file  ProcRates.cc
/// @brief Container class for process rate tensors
///

//**** 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/Tools/HammerSerial.hh"
#include "Hammer/ProcRates.hh"
#include "Hammer/ExternalData.hh"
#include "Hammer/Tools/ParticleUtils.hh"
#include "Hammer/Tools/Logging.hh"
#include "Hammer/Tools/Utils.hh"
#include "Hammer/Math/MultiDim/ScalarContainer.hh"

using namespace std;

namespace Hammer {

    namespace MD = MultiDimensional;

    ProcRates::ProcRates(const ExternalData* ext) : _external{ext} {
    }

    ProcRates::~ProcRates() noexcept {
        _processRates.clear();
        _external = nullptr;
    }

    void ProcRates::defineSettings() {
    }

    void ProcRates::init() {
        _processRates.clear();
    }

    SchemeDict<Tensor>* ProcRates::getProcessRates(HashId id) {
        auto it = _processRates.find(id);
        if(it == _processRates.end()) {
            auto res = _processRates.insert({id, SchemeDict<Tensor>{}});
            if(res.second) {
                return &(res.first->second);
            }
            else {
                return nullptr;
            }
        }
        return &(it->second);
    }

    double ProcRates::getVertexRate(const HashId& id, const string& scheme) const {
        auto it = _processRates.find(id);
        if (it != _processRates.end()){
            auto it2 = it->second.find(scheme);
            if(it2 != it->second.end()){
                if(it2->second.labels().size() > 0) {
                    auto t = _external->getExternalVectors(scheme, it2->second.labels());
                    t.dot(it2->second);
                    if(t.rank() != 0) {
                        MSG_ERROR("Invalid rank at end of evaluation");
                    }
                    return t.element().real();
                }
                else {
                    if(it2->second.rank() != 0) {
                        MSG_ERROR("Invalid rank at end of evaluation");
                    }
                    return it2->second.element().real();
                }
            }
            MSG_ERROR("Rate not found for specified FF scheme. Oh, the humanity.");
            return 0.;
        }
        MSG_ERROR("Rate not found for specified vertex ID. Oh, the humanity.");
        return 0.;
    }

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

    bool ProcRates::read(const Serial::FBRates* msgreader, bool merge) {
        auto rates = msgreader->rates();
        if(!merge) {
            _processRates.clear();
        }
        bool result = true;
        for(unsigned int i = 0; i< rates->size() && result; ++i) {
            auto res = _processRates.insert({rates->Get(i)->id(), SchemeDict<Tensor>{}});
            if(merge && !res.second) {
                res.first = _processRates.find(rates->Get(i)->id());
                if(res.first == _processRates.end()) {
                    return false;
                }
            }
            auto vals = rates->Get(i)->ratevals();
            auto chs = rates->Get(i)->ratenames();
            for(unsigned int j = 0; j < vals->size() && result; ++j) {
                auto it = res.first->second.find(chs->Get(j)->c_str());
                if(merge && it != res.first->second.end() && strcmp(chs->Get(j)->c_str(), "Denominator") != 0) {
                    MSG_ERROR("Try to merge two rates with factor schemes with same name '" +
                              string(chs->Get(j)->c_str()) + "'!");
                    result = false;

                }
                else {
                    auto res2 = res.first->second.insert({chs->Get(j)->c_str(), Tensor{}});
                    res2.first->second.read(vals->Get(j));
                }
            }
        }
        return result;
    }

    void ProcRates::write(flatbuffers::FlatBufferBuilder* msgwriter) const {
        vector<flatbuffers::Offset<Serial::FBRate>> rates;
        rates.reserve(_processRates.size());
        for (auto& elem : _processRates) {
            vector<flatbuffers::Offset<flatbuffers::String>> chs;
            vector<flatbuffers::Offset<Serial::FBTensor>> rts;
            for(auto& elem2: elem.second) {
                auto serCh = msgwriter->CreateString(elem2.first);
                flatbuffers::Offset<Serial::FBTensor> val;
                elem2.second.write(&(*msgwriter),&val);
                chs.push_back(serCh);
                rts.push_back(val);
            }
            auto serialChs = msgwriter->CreateVector(chs);
            auto serialrts = msgwriter->CreateVector(rts);
            Serial::FBRateBuilder serialFF{*msgwriter};
            serialFF.add_id(elem.first);
            serialFF.add_ratevals(serialrts);
            serialFF.add_ratenames(serialChs);
            auto resFF = serialFF.Finish();
            rates.push_back(resFF);
        }
        auto serialrates = msgwriter->CreateVector(rates);
        Serial::FBRatesBuilder serialout{*msgwriter};
        serialout.add_rates(serialrates);
        auto headoffset = serialout.Finish();
        msgwriter->Finish(headoffset);
    }


} // namespace Hammer