rssm.cpp

//**************************************************************************************************
//*  FRIDA: fast reliable interactive data analysis
//*  (C) Joachim Wuttke 1990-, v2(C++) 2001-
//*  http://apps.jcns.fz-juelich.de/frida
//**************************************************************************************************

//! \file  rssm.cpp
//! \brief NRSSM: directly read SPHERES raw data

#include "defs.hpp"

#include <cstring>
#include <fstream>
#include <experimental/filesystem>
#include <yaml-cpp/yaml.h>

#include "../readplus/ask.hpp"
#include "../trivia/string_ops.hpp"

#include "mem.hpp"
#include "obj.hpp"
#include "olf.hpp"
#include "rssm.hpp"
#include "slice.hpp"

namespace fs = std::experimental::filesystem;


//! All relevant information extracted from one raw-data file from SPHERES.

class CRawfileSpheres
{
public:
    void RdRawYam(std::ifstream& F_in);
    double daq_time_step;
    time_t measured_from;
    time_t measured_until;
    int measured_at;
    double measured_during;
    bool incremental;
    int ndet;
    vector<double> rawdata[6];
    int maj_syntax, min_syntax, maj_outform, min_outform;
};

//! Read raw data file, store contents as class variables

void CRawfileSpheres::RdRawYam(std::ifstream& F_in)
{
    string key, val;
    CCoord co;
    CSpec sout;
    double num;

    maj_outform = 5;
    min_outform = 1;

    // YAML::Parser parser(F_in);
    const YAML::Node& doc = YAML::Load(F_in);

    // read Meta:
    if (!doc["Meta"])
        throw S("no Meta");
    try {
        YAML::NodeType::value Metatype = doc["Meta"].Type();
        if (Metatype != YAML::NodeType::Map && doc["Meta"].size() > 0)
            throw S("DATA BUG: Meta is not a MAP");
        if (!doc["Meta"]["format"])
            throw S("DATA BUG: no format in Meta");
        string val;
        val = doc["Meta"]["format"].as<string>();
        if (strncmp(val.c_str(), "acq5.2 for yaml1", 16))
            throw "UNEXPECTED DATA: format: " + val + " instead of acq5.2";
    } catch (std::exception& e) {
        throw "failed to read Meta section: " + string(e.what());
    }

    // ignore History.

    // read Shortpar:
    try {
        const YAML::Node& spar = doc["Shortpar"];
        if (!spar)
            throw S("DATA BUG: no Shortpar");

        if (!spar["incremental"])
            throw S("DATA BUG: no Shortpar 'incremental'");
        string val = spar["incremental"].as<string>();
        if (val == "true")
            incremental = true;
        else if (val == "false")
            incremental = false;
        else
            throw "DATA BUG: invalid value [" + val + "] of 'incremental'";

        if (!spar["subs_until"])
            throw S("DATA BUG: no Shortpar 'subs_until'");
        val = spar["subs_until"].as<string>();
        char date_string[24];
        sscanf(val.c_str(), "%24c", date_string);
        struct tm date_broken;
        strptime(date_string, "%F %A %H:%M:%S", &date_broken);
        measured_until = mktime(&date_broken);

        if (!spar["ndet"])
            throw S("DATA BUG: no Shortpar 'ndet'");
        ndet = spar["ndet"].as<int>();

        if (!spar["daq_time_step"])
            throw S("DATA BUG: no Shortpar 'daq_time_step'");
        val = spar["daq_time_step"].as<string>();
        string numval, unit;
        triv::string_extract_word(val, &numval, &unit);
        if (unit != "s")
            throw "DATA BUG: daq_time_step has unexpected unit '" + unit + "'";
        triv::any2dbl(numval, &daq_time_step);

    } catch (std::exception& e) {
        throw "DATA BUG: failed to read Shortpar section: " + string(e.what());
    }

    // read EnergyHistograms:
    try {
        if (!doc["EnergyHistograms"])
            throw S("DATA BUG: no EnergyHistograms");
        if (doc["EnergyHistograms"].Type() != YAML::NodeType::Sequence)
            throw S("DATA BUG: EnergyHistograms is not a SEQUENCE");
        int nEne = doc["EnergyHistograms"].size();
        for (int iEne = 0; iEne < nEne; iEne++) {
            if (doc["EnergyHistograms"][iEne].Type() != YAML::NodeType::Sequence)
                throw "DATA BUG: EnergyHistogram[" + S(iEne) + "] is not a SEQUENCE";
            val = doc["EnergyHistograms"][iEne][0].as<string>();
            if (!triv::any2dbl(val, &num))
                throw "E-Hist: invalid x-value " + val;
            for (int i = 0; i < 4; ++i)
                rawdata[i].push_back(num);
            for (int i = 0; i < 4; ++i) { // counts
                // inline Node& Node::operator=(const T& rhs);

                const YAML::Node arr = doc["EnergyHistograms"][iEne][i + 1];
                if (arr.Type() != YAML::NodeType::Sequence)
                    throw "DATA BUG: EnergyHistogram[" + S(iEne) + "][" + S(i + 1)
                        + "] is not a SEQUENCE";
                for (int j = 0; j < ndet; ++j) {
                    val = (arr)[j].as<string>();
                    //(*arr)[j] >> val;
                    if (!triv::any2dbl(val, &num))
                        throw "E-Hist: invalid count " + val;
                    rawdata[i].push_back(num);
                }
            }
            for (int i = 0; i < 4; ++i) { // tsteps
                const YAML::Node arr = (doc["EnergyHistograms"][iEne][i + 1 + 4]);
                if (arr.Type() != YAML::NodeType::Sequence)
                    throw "DATA BUG: EnergyHistogram[" + S(iEne) + "][" + S(i + 1 + 4)
                        + "] is not a SEQUENCE";
                for (int j = 0; j < ndet; ++j) {

                    val = (arr)[j].as<string>();
                    if (!triv::any2dbl(val, &num))
                        throw "E-Hist: invalid tstep " + val;
                    rawdata[i].push_back(num);
                }
            }
        }
    } catch (std::exception& e) {
        throw "BAD BUG: failed to read EnergyHistograms section: " + string(e.what());
    }

    // read ChopperHistograms:
    try {
        if (!doc["ChopperHistograms"])
            throw S("no ChopperHistograms");
        if (doc["ChopperHistograms"].Type() != YAML::NodeType::Sequence)
            throw S("DATA BUG: ChopperHistograms is not a SEQUENCE");
        int nCho = doc["ChopperHistograms"].size();
        for (int iCho = 0; iCho < nCho; iCho++) {
            if (doc["ChopperHistograms"][iCho].Type() != YAML::NodeType::Sequence)
                throw "DATA BUG: ChopperHistograms " + S(iCho) + " is not a SEQUENCE";
            val = doc["ChopperHistograms"][iCho][0].as<string>();
            if (!triv::any2dbl(val, &num))
                throw "DATA BUG: C-Hist: invalid x-value " + val;
            for (int i = 4; i < 6; ++i)
                rawdata[i].push_back(num);
            for (int i = 4; i < 6; ++i) { // 2 cnts lines
                const YAML::Node arr = (doc["ChopperHistograms"][iCho][i - 4 + 1]);
                if (arr.Type() != YAML::NodeType::Sequence)
                    throw "DATA BUG: ChopperHistogram[" + S(iCho) + "][" + S(i - 4 + 1)
                        + "] is not a SEQUENCE";
                for (int j = 0; j < ndet; ++j) {
                    val = (arr)[j].as<string>();
                    if (!triv::any2dbl(val, &num))
                        throw "E-Hist: invalid count " + val;
                    rawdata[i].push_back(num);
                }
            }
            for (int i = 4; i < 6; ++i) { // 2 time lines
                const YAML::Node arr = (doc["ChopperHistograms"][iCho][i - 4 + 1 + 2]);
                if (arr.Type() != YAML::NodeType::Sequence)
                    throw "DATA BUG: ChopperHistogram[" + S(iCho) + "][" + S(i - 4 + 1 + 2)
                        + "] is not a SEQUENCE";
                for (int j = 0; j < ndet; ++j) {
                    val = (arr)[j].as<string>();
                    if (!triv::any2dbl(val, &num))
                        throw "E-Hist: invalid count " + val;
                    rawdata[i].push_back(num);
                }
            }
        }
    } catch (std::exception& e) {
        throw "DATA BUG: failed to read ChopperHistograms section: " + string(e.what());
    }
}

//! Read one raw data file.

void NRSSM::read_spec(int flag)

// flag values (to be OR'ed):
//    1 save also open
//    2 save also chop histo
//    4 save spectra also separately for both directions
//    8 do not save regular spectrum
//   16 no normalization

{
    string file_f = wask("Read SPHERES data from file");
    CRawfileSpheres R;
    std::ifstream F_in;

    F_in.open(file_f.c_str(), std::fstream::in);
    if (!(F_in.is_open()))
        throw "Cannot read " + file_f;
    R.RdRawYam(F_in);

    double tstep = R.daq_time_step;
    vector<double>* raw = R.rawdata;
    int ndet = R.ndet;

    // *** set file headers ***

    POld olf[8];
    int iolf;
    for (iolf = 0; iolf < 8; ++iolf) {
        olf[iolf] = POld(new COld);
        olf[iolf]->xco = iolf < 4 || iolf > 5 ? CCoord("E", "ueV") : CCoord("phase", "deg");
        olf[iolf]->yco = CCoord("cts", "sec-1");
        olf[iolf]->ZCo.push_back(CCoord("det", ""));
    }
    string name = fs::path(file_f).filename().stem();
    olf[0]->name = name + "left";
    olf[1]->name = name + "right";
    olf[2]->name = name + "open_left";
    olf[3]->name = name + "open_right";
    olf[4]->name = name + "chop_elast";
    olf[5]->name = name + "chop_inela";
    olf[6]->name = name;
    olf[7]->name = name + "open";
    string doc = "ry" + (flag == 0 ? "" : S(flag + 1)) + " " + file_f;
    olf[0]->log_action(doc + " # refl, left halfspace");
    olf[1]->log_action(doc + " # refl, right halfspace");
    olf[2]->log_action(doc + " # open, left halfspace");
    olf[3]->log_action(doc + " # open, right halfspace");
    olf[4]->log_action(doc + " # elastic");
    olf[5]->log_action(doc + " # inelastic");
    olf[6]->log_action(doc + " # refl");
    olf[7]->log_action(doc + " # open");

    for (iolf = 0; iolf < 8; ++iolf) {
        for (int j = 0; j < ndet; ++j) {
            PSpec s(new CSpec);
            s->z.resize(1);
            s->z[0] = PObjInt(new CObjInt(j));
            olf[iolf]->V.push_back(move(s));
        }
    }

    if (R.maj_outform <= 2) { // old data format
        double x, tmeas, count;
        // transfer data from raw arrays:
        for (iolf = 0; iolf < 6; ++iolf) {
            for (int iraw = 0; iraw < raw[iolf].size(); iraw += (ndet + 2)) {
                x = raw[iolf][iraw];
                tmeas = raw[iolf][iraw + ndet + 1];
                if (tmeas < 1)
                    continue;
                for (int j = 0; j < ndet; ++j) {
                    count = raw[iolf][iraw + 1 + j];
                    olf[iolf]->VS(j)->push_xy(x, count / tmeas / tstep);
                }
            }
        }
        // sum over halfspaces:
        for (int ichop = 0; ichop < 2; ++ichop) {
            int iolf0 = ichop * 2;
            int iolf1 = ichop * 2 + 1;
            int nline = raw[iolf0].size();
            if (raw[iolf1].size() != nline)
                throw S("inconsistent raw array sizes");
            for (int iraw = 0; iraw < nline; iraw += (ndet + 2)) {
                x = raw[iolf0][iraw];
                tmeas = raw[iolf0][iraw + ndet + 1] + raw[iolf1][iraw + ndet + 1];
                if (tmeas < 1)
                    continue;
                for (int j = 0; j < ndet; ++j) {
                    count = raw[iolf0][iraw + 1 + j] + raw[iolf1][iraw + 1 + j];
                    olf[6 + ichop]->VS(j)->push_xy(x, count / tmeas / tstep);
                }
            }
        }

    } else { // as of 2008-03-03
        double x, tmeas, count;
        // transfer data from raw arrays:
        for (iolf = 0; iolf < 6; ++iolf) {
            for (int iraw = 0; iraw < raw[iolf].size(); iraw += (1 + ndet * 2)) {
                x = raw[iolf][iraw];
                for (int j = 0; j < ndet; ++j) {
                    count = raw[iolf][iraw + 1 + j];
                    if (!(flag & 16)) {
                        tmeas = raw[iolf][iraw + 1 + ndet + j];
                        if (tmeas < 1)
                            continue;
                        count /= (tmeas * tstep);
                    }
                    olf[iolf]->VS(j)->push_xy(x, count);
                }
            }
        }
        // sum over halfspaces:
        for (int ichop = 0; ichop < 2; ++ichop) {
            int iolf0 = ichop * 2;
            int iolf1 = ichop * 2 + 1;
            int nline = raw[iolf0].size();
            if (raw[iolf1].size() != nline)
                throw S("inconsistent raw array sizes");
            for (int iraw = 0; iraw < nline; iraw += (1 + ndet * 2)) {
                x = raw[iolf0][iraw];
                for (int j = 0; j < ndet; ++j) {
                    count = raw[iolf0][iraw + 1 + j] + raw[iolf1][iraw + 1 + j];
                    if (!(flag & 16)) {
                        tmeas = raw[iolf0][iraw + 1 + ndet + j] + raw[iolf1][iraw + 1 + ndet + j];
                        if (tmeas < 1)
                            continue;
                        count /= (tmeas * tstep);
                    }
                    olf[6 + ichop]->VS(j)->push_xy(x, count);
                }
            }
        }
    }

    if (flag & 2) {
        SMem::instance()->mem_store(move(olf[4]));
        SMem::instance()->mem_store(move(olf[5]));
    }
    if (flag & 2 && flag & 4) {
        SMem::instance()->mem_store(move(olf[2]));
        SMem::instance()->mem_store(move(olf[3]));
    }
    if (flag & 4) {
        SMem::instance()->mem_store(move(olf[0]));
        SMem::instance()->mem_store(move(olf[1]));
    }
    if (flag & 1) {
        SMem::instance()->mem_store(move(olf[7]));
    }
    if (!(flag & 8)) {
        SMem::instance()->mem_store(move(olf[6]));
    }
}

//! Read a series of raw data files.

void NRSSM::read_series(int flag)

// flag values (to be OR'ed):
//    1 save also open

{
    char tstrg[30];
    string fser, fnam;
    time_t t;
    int isub;
    vector<CRawfileSpheres> RR;

    // Read consolidated files:
    fser = wask("Read SPHERES data from series");
    for (isub = 0;; ++isub) {
        fnam = fser + "a" + S(isub);
        std::ifstream F_in;
        F_in.open(fnam.c_str(), std::fstream::in);
        if (!(F_in.is_open()))
            break;
        // cout << "successfully opened "<<fnam<<"\n";
        CRawfileSpheres rf;
        rf.RdRawYam(F_in);
        RR.push_back(rf);
    }

    int nsub = RR.size();
    if (nsub == 0)
        throw "could not open file " + fnam;
    else if (nsub == 1)
        throw S("series contains just one spectrum -> use a simpler method");
    printf("successfully read %i files\n", nsub);

    // Correct time:
    double mean_time = (RR[nsub - 1].measured_until - RR[0].measured_until) / nsub;
    printf("mean measuring time %g seconds per file\n", mean_time);
    t = (time_t)(RR[0].measured_until - mean_time);
    struct tm tlocal;
    strftime(tstrg, 30, "%F %b %H:%M:%S", localtime_r(&t, &tlocal));
    printf("measurement times relative to start time %s\n", tstrg);
    for (isub = 0; isub < nsub; ++isub)
        RR[isub].measured_at =
            RR[isub].measured_until - RR[0].measured_until + (int)(mean_time / 2);

    // Get parameters, check consistency:
    int ndet = RR[0].ndet;
    printf("used %i detectors\n", ndet);
    double tstep = RR[0].daq_time_step;
    bool incremental_in = RR[0].incremental;

    for (isub = 1; isub < nsub; ++isub) {
        if (RR[isub].ndet != ndet)
            throw "inconsistent no of det: file " + S(isub) + " has " + S(RR[isub].ndet)
                + " instead of " + S(ndet);
        if (RR[isub].daq_time_step != tstep)
            throw S("inconsistent time step");
        if (RR[isub].incremental != incremental_in)
            throw S("inconsistent increment status");
    }

    // Set file headers:
    POld olf[2];
    int iolf;
    for (int iolf = 0; iolf < 2; ++iolf)
        olf[iolf] = POld(new COld);
    int j, nE;

    for (iolf = 0; iolf < 2; ++iolf) {
        olf[iolf]->log_action("acquire data from series" + fser);
        olf[iolf]->xco = CCoord("E", "ueV");
        olf[iolf]->yco = CCoord("cts", "sec-1");
        olf[iolf]->ZCo.push_back(CCoord("#sub", ""));
        olf[iolf]->ZCo.push_back(CCoord("det", ""));
    }
    string name = fs::path(fser).filename().stem();
    olf[0]->name = name;
    olf[1]->name = name + "open";
    olf[0]->log_action("refl, average over both halfspaces");
    olf[1]->log_action("open, average over both halfspaces");


    for (isub = 0; isub < nsub; ++isub) {
        // ? = (double)RR[isub].measured_at/3600;
        for (j = 0; j < ndet; ++j) {
            PSpec s(new CSpec);
            s->z.resize(2);
            s->z[0] = PObjInt(new CObjInt(isub));
            s->z[1] = PObjInt(new CObjInt(j));
            for (iolf = 0; iolf < 2; ++iolf)
                olf[iolf]->V.push_back(move(s));
        }
    }

    if (RR[0].maj_outform >= 3) {
        printf("new format\n");
        nE = RR[0].rawdata[0].size() / (1 + ndet * 2);
        // Reverse incrementation?
        if (incremental_in) {
            printf("reversing incrementation\n");
            for (isub = nsub - 1; isub >= 1; --isub) {
                for (int iarr = 0; iarr < 4; ++iarr) {
                    for (int iE = 0; iE < nE; ++iE) {
                        for (j = 1; j < 1 + ndet * 2; ++j) {
                            RR[isub].rawdata[iarr][iE * (1 + ndet * 2) + j] -=
                                RR[isub - 1].rawdata[iarr][iE * (1 + ndet * 2) + j];
                        }
                    }
                }
            }
        }

        // Transfer data from raw arrays, sum over halfspaces:
        vector<double>* raw;
        double x, tmeas, count;
        int iraw;
        for (isub = 0; isub < nsub; ++isub) {
            raw = RR[isub].rawdata;
            for (int ichop = 0; ichop < 2; ++ichop) {
                int iolf0 = ichop * 2;
                int iolf1 = ichop * 2 + 1;
                if (raw[iolf1].size() / (1 + ndet * 2) != nE) {
                    printf(
                        "inconsistent raw array sizes:\n  chop %i raw1 %i nE %i\n", ichop,
                        (int)raw[iolf1].size(), nE);
                    return;
                }
                for (int iE = 0; iE < nE; ++iE) {
                    iraw = iE * (1 + ndet * 2);
                    x = raw[iolf0][iraw];
                    for (j = 0; j < ndet; ++j) {
                        tmeas = raw[iolf0][iraw + ndet + 1 + j] + raw[iolf1][iraw + ndet + 1 + j];
                        if (tmeas < 1)
                            continue;
                        count = raw[iolf0][iraw + 1 + j] + raw[iolf1][iraw + 1 + j];
                        olf[ichop]->VS(isub * ndet + j)->push_xy(x, count / tmeas / tstep);
                    }
                }
            }
        }

    } else {
        printf("old format\n");
        nE = RR[0].rawdata[0].size() / (ndet + 2);
        // Reverse incrementation?
        if (incremental_in) {
            printf("reversing incrementation\n");
            for (isub = nsub - 1; isub >= 1; --isub) {
                for (int iarr = 0; iarr < 4; ++iarr) {
                    for (int iE = 0; iE < nE; ++iE) {
                        for (j = 1; j < ndet + 2; ++j) {
                            RR[isub].rawdata[iarr][iE * (ndet + 2) + j] -=
                                RR[isub - 1].rawdata[iarr][iE * (ndet + 2) + j];
                        }
                    }
                }
            }
        }

        // Transfer data from raw arrays, sum over halfspaces:
        vector<double>* raw;
        double x, tmeas, count;
        int iraw;
        for (isub = 0; isub < nsub; ++isub) {
            raw = RR[isub].rawdata;
            for (int ichop = 0; ichop < 2; ++ichop) {
                int iolf0 = ichop * 2;
                int iolf1 = ichop * 2 + 1;
                if (raw[iolf1].size() / (ndet + 2) != nE) {
                    printf("inconsistent raw array sizes\n");
                    return;
                }
                for (int iE = 0; iE < nE; ++iE) {
                    iraw = iE * (ndet + 2);
                    x = raw[iolf0][iraw];
                    tmeas = raw[iolf0][iraw + ndet + 1] + raw[iolf1][iraw + ndet + 1];
                    if (tmeas < 1)
                        continue;
                    for (j = 0; j < ndet; ++j) {
                        count = raw[iolf0][iraw + 1 + j] + raw[iolf1][iraw + 1 + j];
                        olf[ichop]->VS(isub * ndet + j)->push_xy(x, count / tmeas / tstep);
                    }
                }
            }
        }
    }

    if (flag & 1) {
        SMem::instance()->mem_store(move(olf[1]));
    }
    if (1) {
        SMem::instance()->mem_store(move(olf[0]));
    }
}