Hammer.cc

Go to the documentation of this file.

///
/// @file  Hammer.cc
/// @brief Main Hammer class
///

//**** 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/Hammer.hh"
#include "Hammer/DictionaryManager.hh"
#include "Hammer/Histos.hh"
#include "Hammer/SettingsHandler.hh"
#include "Hammer/Event.hh"

#include "Hammer/AmplitudeBase.hh"
#include "Hammer/FormFactorBase.hh"
#include "Hammer/Math/Histogram.hh"
#include "Hammer/ProvidersRepo.hh"
#include "Hammer/ExternalData.hh"
#include "Hammer/ProcRates.hh"
#include "Hammer/ProcManager.hh"
#include "Hammer/SchemeDefinitions.hh"
#include "Hammer/ProcessDefinitions.hh"
#include "Hammer/PurePhaseSpaceDefs.hh"
#include "Hammer/Config/HammerModules.hh"

#include "Hammer/Tools/HammerRoot.hh"
#include "Hammer/Tools/HammerSerial.hh"

#include "Hammer/Math/Utils.hh"
#include "Hammer/Math/Units.hh"
#include "Hammer/Tools/ParticleUtils.hh"
#include "Hammer/Tools/Logging.hh"
#include "Hammer/Tools/Utils.hh"
#include "Hammer/Tools/Pdg.hh"


using namespace std;

namespace Hammer {

    Hammer::Hammer() : _containers{new DictionaryManager{}},
                       _settings{new SettingsHandler{}},
                       _histograms{new Histos{_containers.get()}},
                       _builder{new flatbuffers::FlatBufferBuilder{16*1024*1024}},
                       _event{new Event{_histograms.get(), _containers.get()}} {
        Log::resetWarningCounters();
        this->setSettingsHandler(*_settings);
        _histograms->setSettingsHandler(*_settings);
        _event->setSettingsHandler(*_settings);
        _containers->setSettingsHandler(*_settings);
        // ProcManager{};
        _histograms->init();
    }

    Hammer::~Hammer() noexcept {
        _builder->Clear();
        _event->clear();
    }

    void Hammer::initEvent(double weight) {
        _event->clear();
        if(_histograms->isValidHistogram("Total Sum of Weights", 1)){
            _event->setHistogramBin("Total Sum of Weights", IndexList{0});
        }
        if(!isZero(weight - 1.0)) {
            _event->setEventBaseWeight(weight);
        }
    }

    size_t Hammer::addProcess(Process& p) {
        return _event->addProcess(p);
    }

    void Hammer::removeProcess(size_t id) {
        _event->removeProcess(id);
    }

    void Hammer::setEventHistogramBin(const string& name, const IndexList& bins) {
        _event->setHistogramBin(name, bins);
    }

    void Hammer::fillEventHistogram(const string& name, const vector<double>& values) {
        if(_histograms->canFill(name, values)){
            setEventHistogramBin(name, _histograms->getBinIndices(name, values));
        }
    }

    void Hammer::setEventBaseWeight(double weight) {
        if(!isZero(weight - 1.0)) {
            _event->setEventBaseWeight(weight);
        }
    }

    bool Hammer::loadEventWeights(IOBuffer& buffer, bool merge) {
        if(buffer.kind == RecordType::EVENT) {
            const Serial::FBEvent* buf = flatbuffers::GetRoot<Serial::FBEvent>(buffer.start);
            if (buf != nullptr) {
                return _event->read(buf, merge);
            }
        }
        MSG_ERROR("Invalid event record on read!");
        return false;
    }

    IOBuffer Hammer::saveEventWeights() const {
        _builder->Clear();
        _event->write(&(*_builder));
        IOBuffer result;
        result.kind = RecordType::EVENT;
        result.start = _builder->GetBufferPointer();
        result.length = _builder->GetSize();
        return result;
    }

    bool Hammer::loadRunHeader(IOBuffer& buffer, bool merge) {
        if(buffer.kind == RecordType::HEADER) {
            const Serial::FBHeader* buf = flatbuffers::GetRoot<Serial::FBHeader>(buffer.start);
            if (buf != nullptr) {
                bool result = _settings->read(buf, merge);
                result &= _containers->read(buf, merge);
                return result;
            }
        }
        MSG_ERROR("Invalid header record on read!");
        return false;
    }

    IOBuffer Hammer::saveRunHeader() const {
        _builder->Clear();
        _containers->write(&(*_builder));
        IOBuffer result;
        result.kind = RecordType::HEADER;
        result.start = _builder->GetBufferPointer();
        result.length = _builder->GetSize();
        return result;
    }

    string Hammer::loadHistogramDefinition(IOBuffer& buffer, bool merge) {
        if(buffer.kind == RecordType::HISTOGRAM_DEFINITION) {
            const Serial::FBHistoDefinition* buf = flatbuffers::GetRoot<Serial::FBHistoDefinition>(buffer.start);
            if (buf != nullptr) {
                return _histograms->readDefinition(buf, merge);
            }
        }
        MSG_ERROR("Invalid histogram record on read!");
        return "";
    }

    HistoInfo Hammer::loadHistogram(IOBuffer& buffer, bool merge) {
        if(buffer.kind == RecordType::HISTOGRAM) {
            const Serial::FBHistogram* buf = flatbuffers::GetRoot<Serial::FBHistogram>(buffer.start);
            if (buf != nullptr) {
                return _histograms->readHistogram(buf, merge);
            }
        }
        MSG_ERROR("Invalid histogram record on read!");
        return HistoInfo{"", "", EventUIDGroup{}};
    }

    IOBuffers Hammer::saveHistogram(const string& name) const {
        unique_ptr<Serial::DetachedBuffers> resultVec{new Serial::DetachedBuffers{}};
        _histograms->checkExists(name);
        _builder->Clear();
        if(_histograms->writeDefinition(&(*_builder), name)) {
            resultVec->add(_builder->Release(), RecordType::HISTOGRAM_DEFINITION);
            for (auto& schemeName : _containers->schemeDefs().getFFSchemeNames()) {
                for (auto& eventIDs : _histograms->getEventIDRepsForHisto(name, schemeName)) {
                    _builder->Clear();
                    _histograms->writeHistogram(&(*_builder), name, schemeName, eventIDs);
                    resultVec->add(_builder->Release(), RecordType::HISTOGRAM);
                }
            }
        }
        return IOBuffers{move(resultVec)};
    }

    IOBuffers Hammer::saveHistogram(const string& name, const string& scheme) const {
        unique_ptr<Serial::DetachedBuffers> resultVec{new Serial::DetachedBuffers{}};
        _histograms->checkExists(name);
        _builder->Clear();
        if(_histograms->writeDefinition(&(*_builder), name)) {
            resultVec->add(_builder->Release(), RecordType::HISTOGRAM_DEFINITION);
            for(auto& eventIDs: _histograms->getEventIDRepsForHisto(name, scheme)) {
                _builder->Clear();
                _histograms->writeHistogram(&(*_builder), name, scheme, eventIDs);
                resultVec->add(_builder->Release(), RecordType::HISTOGRAM);
            }
            if(resultVec->size() == 1) {
                MSG_ERROR("No scheme '" + scheme + "' present for histogram '" + name + "'. Only definition will be saved.");
            }
        }
        return IOBuffers{move(resultVec)};
    }

    IOBuffers Hammer::saveHistogram(const string& name, const EventUIDGroup& eventIDs) const {
        unique_ptr<Serial::DetachedBuffers> resultVec{new Serial::DetachedBuffers{}};
        _histograms->checkExists(name);
        _builder->Clear();
        if (_histograms->writeDefinition(&(*_builder), name)) {
            resultVec->add(_builder->Release(), RecordType::HISTOGRAM_DEFINITION);
            for (auto& schemeName : _containers->schemeDefs().getFFSchemeNames()) {
                _builder->Clear();
                _histograms->writeHistogram(&(*_builder), name, schemeName, *eventIDs.begin());
                resultVec->add(_builder->Release(), RecordType::HISTOGRAM);
            }
            if(resultVec->size() == 1) {
                string groupId = "[ ";
                for(auto& elem: eventIDs) {
                    groupId += "[ ";
                    for(auto& elem2: elem) {
                        groupId += to_string(elem2) + " ";
                    }
                    groupId += "] ";
                }
                groupId += "] ";
                MSG_ERROR("No event ID " + groupId + "present for histogram '" + name + "'. Only definition will be saved.");
            }
        }
        return IOBuffers{move(resultVec)};
    }

    IOBuffers Hammer::saveHistogram(const string& name, const string& scheme, const EventUIDGroup& eventIDs) const {
        unique_ptr<Serial::DetachedBuffers> resultVec{new Serial::DetachedBuffers{}};
        _histograms->checkExists(name);
        _builder->Clear();
        if (_histograms->writeDefinition(&(*_builder), name)) {
            resultVec->add(_builder->Release(), RecordType::HISTOGRAM_DEFINITION);
            _builder->Clear();
            _histograms->writeHistogram(&(*_builder), name, scheme, *eventIDs.begin());
            resultVec->add(_builder->Release(), RecordType::HISTOGRAM);
            if (resultVec->size() == 1) {
                string groupId = "[ ";
                for (auto& elem : eventIDs) {
                    groupId += "[ ";
                    for (auto& elem2 : elem) {
                        groupId += to_string(elem2) + " ";
                    }
                    groupId += "] ";
                }
                groupId += "] ";
                MSG_ERROR("No event ID " + groupId + "/ scheme '" + scheme + "' combination present for histogram '" + name +
                          "'. Only definition will be saved.");
            }
        }
        return IOBuffers{move(resultVec)};
    }

    IOBuffers Hammer::saveHistogram(const HistoInfo& info) const {
        if(info.scheme.size() > 0 && info.eventGroupId.size() > 0) {
            return saveHistogram(info.name, info.scheme, info.eventGroupId);
        }
        else if(info.scheme.size() > 0) {
            return saveHistogram(info.name, info.scheme);
        }
        else if(info.eventGroupId.size() > 0) {
            return saveHistogram(info.name, info.eventGroupId);
        }
        else {
            return saveHistogram(info.name);
        }
    }

    bool Hammer::loadRates(IOBuffer& buffer, bool merge) {
        if(buffer.kind == RecordType::RATE) {
            const Serial::FBRates* buf = flatbuffers::GetRoot<Serial::FBRates>(buffer.start);
            if (buf != nullptr) {
                return _containers->rates().read(buf, merge);
            }
        }
        MSG_ERROR("Invalid rate record on read!");
        return false;
    }

    IOBuffer Hammer::saveRates() const {
        _builder->Clear();
        _containers->rates().write(&(*_builder));
        IOBuffer result;
        result.kind = RecordType::RATE;
        result.start = _builder->GetBufferPointer();
        result.length = _builder->GetSize();
        return result;
    }

    void Hammer::processEvent() {
        if(isOn("CalcProcesses")) {
            _event->calc();
        }
        if(_histograms->size() > 0 && isOn("CalcHistograms")) {
            _event->fillHistograms();
        }
    }

    void Hammer::readCards(const string& fileDecays, const string& fileOptions) {
        _containers->readDecays(fileDecays);
        _settings->readSettings(fileOptions);
    }

    void Hammer::saveOptionCard(const string& fileOptions, bool useDefault) const {
        _settings->saveSettings(fileOptions, useDefault);
    }

    void Hammer::saveHeaderCard(const string& fileDecays) const {
        _containers->saveDecays(fileDecays);
    }

    void Hammer::saveReferences(const std::string& fileRefs) const {
        _settings->saveReferences(fileRefs);
    }

    void Hammer::setOptions(const string& options) {
        _settings->parseSettings(options);
    }

    void Hammer::setHeader(const string& options) {
        _containers->parseDecays(options);
    }

    void Hammer::initRun() {
//        addHistogram("Total Sum of Weights", IndexList{1}, false);
        _containers->init();
        _event->init();
        _histograms->init();
    }

    void Hammer::addTotalSumOfWeights(const bool compress, const bool witherrors){
        addHistogram("Total Sum of Weights", IndexList{1}, false);
        if(compress) {collapseProcessesInHistogram("Total Sum of Weights");}
        keepErrorsInHistogram("Total Sum of Weights", witherrors);
    }

    void Hammer::addHistogram(const string& name, const IndexList& binSizes, bool hasUnderOverFlow,
                              const vector<pair<double, double>>& ranges) {
        _histograms->addHistogramDefinition(name, binSizes, hasUnderOverFlow, ranges);
    }

    void Hammer::addHistogram(const string& name, const vector<vector<double>>& binEdges,
                              bool hasUnderOverFlow) {
        _histograms->addHistogramDefinition(name, binEdges, hasUnderOverFlow);
    }

    void Hammer::keepErrorsInHistogram(const string& name, bool value) {
        _histograms->setHistogramKeepErrors(name, value);
    }

    void Hammer::collapseProcessesInHistogram(const string& name) {
        _histograms->setHistogramCompression(name);
    }

    void Hammer::specializeWCInWeights(const string& process, const vector<complex<double>>& values){
        _containers->externalData().specializeWCInWeights(process, values);
    }

    void Hammer::specializeWCInWeights(const string& process, const map<string, complex<double>>& settings){
        _containers->externalData().specializeWCInWeights(process, settings);
    }

    void Hammer::resetSpecializeWCInWeights(const string& process){
        _containers->externalData().resetSpecializeWCInWeights(process);
    }

    void Hammer::specializeWCInHistogram(const string& name, const string& process,
                                 const vector<complex<double>>& values) {
        auto data = _containers->externalData().getTempWilsonCoefficients(process, values);
        _histograms->addHistogramFixedData(name, data);
    }

    void Hammer::specializeWCInHistogram(const string& name, const string& process,
                                 const map<string, complex<double>>& settings) {
        auto data = _containers->externalData().getTempWilsonCoefficients(process, settings);
        _histograms->addHistogramFixedData(name, data);
    }

    void Hammer::specializeFFInHistogram(const string& name, const string& process, const string& group,
                                 const vector<double>& values) {
        auto data = _containers->externalData().getTempFFEigenVectors({process, group}, values);
        _histograms->addHistogramFixedData(name, data);
    }

    void Hammer::specializeFFInHistogram(const string& name, const string& process, const string& group,
                                 const map<string, double>& settings) {
        auto data = _containers->externalData().getTempFFEigenVectors({process, group}, settings);
        _histograms->addHistogramFixedData(name, data);
    }

    void Hammer::resetSpecializationInHistogram(const string& name) {
        _histograms->resetHistogramFixedData(name);
    }

    void Hammer::createProjectedHistogram(const string& oldName, const string& newName,
                                          const set<uint16_t>& collapsedIndexPositions) {
        _histograms->createProjectedHistogram(oldName, newName, collapsedIndexPositions);
    }


    void Hammer::removeHistogram(const string& name) {
        _histograms->removeHistogram(name);
    }

    void Hammer::addFFScheme(const string& schemeName, const map<string, string>& schemes) {
        _containers->schemeDefs().addFFScheme(schemeName, schemes);
    }

    void Hammer::setFFInputScheme(const map<string, string>& schemes) {
        _containers->schemeDefs().setFFInputScheme(schemes);
    }

    void Hammer::removeFFScheme(const string& schemeName) {
        _containers->schemeDefs().removeFFScheme(schemeName);
    }

    vector<string> Hammer::getFFSchemeNames() const {
        return _containers->schemeDefs().getFFSchemeNames();
    }

    void Hammer::includeDecay(const vector<string>& names) {
        _containers->processDefs().addIncludedDecay(names);
    }

    void Hammer::includeDecay(const string& name) {
        _containers->processDefs().addIncludedDecay({name});
    }

    void Hammer::forbidDecay(const vector<string>& names) {
        _containers->processDefs().addForbiddenDecay(names);
    }

    void Hammer::forbidDecay(const string& name) {
        _containers->processDefs().addForbiddenDecay({name});
    }

    void Hammer::addPurePSVertices(const set<string>& vertices, WTerm what) {
        _containers->purePSDefs().addPurePhaseSpaceVertices(vertices, what);
    }

    void Hammer::clearPurePSVertices(WTerm what) {
        _containers->purePSDefs().clearPurePhaseSpaceVertices(what);
    }

    void Hammer::setUnits(string name) {
        _settings->changeSetting<string>("Hammer", "Units", name);
        Units& units = Units::instance();
        _mcunits = units.getUnitsRescalingToMC(name, "GeV");
    }

    void Hammer::defineSettings() {
        setPath("Hammer");
        addSetting<bool>("CalcHistograms", true);
        addSetting<bool>("CalcProcesses", true);
        addSetting<string>("Units", "GeV");
    }

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

    void Hammer::setWilsonCoefficients(const string& process, const vector<complex<double>>& values, WTerm what) {
        _containers->externalData().setWilsonCoefficients(process, values, what);
    }

    void Hammer::setWilsonCoefficients(const string& process, const map<string, complex<double>>& values, WTerm what) {
        _containers->externalData().setWilsonCoefficients(process, values, what);
    }

    void Hammer::setWilsonCoefficientsLocal(const string& process, const vector<complex<double>>& values) {
        _containers->externalData().setWilsonCoefficientsLocal(process, values);
    }

    void Hammer::setWilsonCoefficientsLocal(const string& process, const map<string, complex<double>>& values) {
        _containers->externalData().setWilsonCoefficientsLocal(process, values);
    }

    void Hammer::resetWilsonCoefficients(const string& process, WTerm what) {
        _containers->externalData().resetWilsonCoefficients(process, what);
    }

    void Hammer::setFFEigenvectors(const string& process, const string& group, const vector<double>& values) {
        FFPrefixGroup tmp{process, group};
        if(_containers->providers().checkFFPrefixAndGroup(tmp)) {
            _containers->externalData().setFFEigenVectors(tmp, values);
        }
    }

    void Hammer::setFFEigenvectors(const string& process, const string& group, const map<string, double>& settings) {
        FFPrefixGroup tmp{process, group};
        if(_containers->providers().checkFFPrefixAndGroup(tmp)) {
            _containers->externalData().setFFEigenVectors(tmp, settings);
        }
    }

    void Hammer::setFFEigenvectorsLocal(const string& process, const string& group, const vector<double>& values) {
        FFPrefixGroup tmp{process, group};
        if (_containers->providers().checkFFPrefixAndGroup(tmp)) {
            _containers->externalData().setFFEigenVectorsLocal(tmp, values);
        }
    }

    void Hammer::setFFEigenvectorsLocal(const string& process, const string& group, const map<string, double>& settings) {
        FFPrefixGroup tmp{process, group};
        if (_containers->providers().checkFFPrefixAndGroup(tmp)) {
            _containers->externalData().setFFEigenVectorsLocal(tmp, settings);
        }
    }

    void Hammer::resetFFEigenvectors(const string& process, const string& group) {
        FFPrefixGroup tmp{process, group};
        if (_containers->providers().checkFFPrefixAndGroup(tmp)) {
            _containers->externalData().resetFFEigenVectors(tmp);
        }
    }

    double Hammer::getWeight(const string& scheme, const vector<size_t>& processes) const {
        auto tempProcesses = processes;
        if(tempProcesses.size() == 0) {
            copy(_event->getEventId().begin(), _event->getEventId().end(), back_inserter(tempProcesses));
        }
        double result = _event->getEventBaseWeight();
        for(auto elem: tempProcesses) {
            try{
                result *= _event->getWeight(scheme, elem);
            }
            catch(RangeError& err) {
                MSG_ERROR("Weight not found: " << err.what() << ". Skipping");
                continue;
            }
        }
        return result;
    }

    double Hammer::getWeight(const string& scheme, const vector<vector<string>>& processes) const {
        const auto procIdSet = _containers->processDefs().decayStringsToProcessIds(processes);
        const vector<HashId> processesId(procIdSet.begin(), procIdSet.end());
        return getWeight(scheme, processesId);
    }

    map<size_t, double> Hammer::getWeights(const string& scheme) const {
        return _event->getWeights(scheme);
    }

    double Hammer::getRate(const HashId& vertexid, const string& scheme) const {
        return _containers->rates().getVertexRate(vertexid, scheme)*_mcunits;
    }

    double Hammer::getRate(const PdgId& parent, const vector<PdgId>& daughters, const string& scheme) const {
        auto tmp = combineDaughters(daughters, {});
        HashId vertexid = processID(parent, tmp);
        //Integral
        return _containers->rates().getVertexRate(vertexid, scheme)*_mcunits;
    }

    double Hammer::getRate(const string& vertex, const string& scheme) const {
        PID& pdg = PID::instance();
        auto vertexid = pdg.expandToValidVertexUIDs(vertex);
       if (vertexid.size() == 1){
            return _containers->rates().getVertexRate(vertexid[0], scheme)*_mcunits;
        } else {
            MSG_ERROR("Vertex string does not correspond to a (charge) unique process. ");
            return 0.;
        }
    }

    double Hammer::getDenominatorRate(const HashId& vertexid) const {
        return getRate(vertexid, "Denominator");
    }

    double Hammer::getDenominatorRate(const PdgId& parent, const vector<PdgId>& daughters) const {
        return getRate(parent, daughters, "Denominator");
    }

    double Hammer::getDenominatorRate(const string& vertex) const {
        return getRate(vertex, "Denominator");
    }


    IOHistogram Hammer::getHistogram(const string& name, const string& scheme) const {
        return  _histograms->getHistogram(name, scheme);
    }

    EventIdGroupDict<IOHistogram> Hammer::getHistograms(const string& name, const string& scheme) const {
        return  _histograms->getHistograms(name, scheme);
    }

    EventUIDGroup Hammer::getHistogramEventIds(const string& name, const string& scheme) const {
        return  _histograms->getHistogramEventIds(name, scheme);
    }

    vector<vector<double>> Hammer::getHistogramBinEdges(const string& name) const {
        return _histograms->getHistogramEdges(name);
    }

    IndexList Hammer::getHistogramShape(const string& name) const {
        return _histograms->getHistogramShape(name);
    }

    bool Hammer::histogramHasUnderOverFlows(const string& name) const {
        return _histograms->getUnderOverFlows(name);
    }

#ifdef HAVE_ROOT

    unique_ptr<TH1D> Hammer::getHistogram1D(const string& name, const string& scheme) const {
        return _histograms->getHistogram1D(name, scheme);
    }

    unique_ptr<TH2D> Hammer::getHistogram2D(const string& name, const string& scheme) const {
        return _histograms->getHistogram2D(name, scheme);
    }

    unique_ptr<TH3D> Hammer::getHistogram3D(const string& name, const string& scheme) const {
        return _histograms->getHistogram3D(name, scheme);
    }

    EventIdGroupDict<unique_ptr<TH1D>> Hammer::getHistograms1D(const string& name, const string& scheme) const {
        return _histograms->getHistograms1D(name, scheme);
    }

    EventIdGroupDict<unique_ptr<TH2D>> Hammer::getHistograms2D(const string& name, const string& scheme) const {
        return _histograms->getHistograms2D(name, scheme);
    }

    EventIdGroupDict<unique_ptr<TH3D>> Hammer::getHistograms3D(const string& name, const string& scheme) const {
        return _histograms->getHistograms3D(name, scheme);
    }

    void Hammer::setHistogram1D(const string& name, const string& scheme, TH1D& histogram) const {
        _histograms->setHistogram1D(name, scheme, histogram);
    }

    void Hammer::setHistogram2D(const string& name, const string& scheme, TH2D& histogram) const {
        _histograms->setHistogram2D(name, scheme, histogram);
    }

    void Hammer::setHistogram3D(const string& name, const string& scheme, TH3D& histogram) const {
        _histograms->setHistogram3D(name, scheme, histogram);
    }

    void Hammer::setHistograms1D(const string& name, const string& scheme,
                                 EventIdGroupDict<unique_ptr<TH1D>>& histograms) const {
        _histograms->setHistograms1D(name, scheme, histograms);
    }

    void Hammer::setHistograms2D(const string& name, const string& scheme,
                                 EventIdGroupDict<unique_ptr<TH2D>>& histograms) const {
        _histograms->setHistograms2D(name, scheme, histograms);
    }

    void Hammer::setHistograms3D(const string& name, const string& scheme,
                                 EventIdGroupDict<unique_ptr<TH3D>>& histograms) const {
        _histograms->setHistograms3D(name, scheme, histograms);
    }

#endif

} // namespace Hammer