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 <cstring>
#include <fstream>
#include <vector>

#include <yaml-cpp/yaml.h>

#include <../trivia/string_ops.hpp>
#include <../trivia/file_ops.hpp>
#include <../readplus/ask.hpp>

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


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

class CRawfileSpheres {
 public:
    void RdRawYam( 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( ifstream& F_in )
{
    string key, val, bla, blub;
    CCoord co;
    CSpec sout;
    double num;

    maj_outform = 5;
    min_outform = 1;

    daq_time_step=0;

    //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( exception& e) {
        throw "failed to read Meta section: " + string(e.what());
    }

    // ignore History.

    // read Shortpar:
    try{
        if(!doc["Shortpar"])
            throw S("DATA BUG: no Shortpar");
        for(auto it=doc["Shortpar"].begin(); it!=doc["Shortpar"].end(); ++it) {
            key = it->first.as<string>();
             val=it->second.as<string>();
             if      ( key=="incremental" ){
                if      ( val=="true" )
                    incremental = true;
                else if ( val=="false" )
                    incremental = false;
                else
                    throw "DATA BUG: invalid value [" + val + "] of 'incremental'";
            } else if ( key=="subs_until" ){
                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 );
            } else if ( key=="ndet" ){
                if( !triv::any2int( val, &ndet ) )
                    throw S("DATA BUG: invalid ndet");
            } else if ( key=="daq_time_step" ){
                triv::string_extract_word( val, &bla, &blub );
                triv::any2dbl( bla, &daq_time_step );
            }
        }
        if( !daq_time_step )
            throw S("DATA BUG: daq time step not found");
    } catch( 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( 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( 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, name;
    file_f = wask("Read SPHERES data from file");
    CRawfileSpheres R;
    ifstream F_in;

    F_in.open (file_f.c_str(),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", "") );
    }
    triv::fname_divide( file_f, 0, &name, 0 );
    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]->lDoc.push_back( doc + " # refl, left halfspace" );
    olf[1]->lDoc.push_back( doc + " # refl, right halfspace" );
    olf[2]->lDoc.push_back( doc + " # open, left halfspace" );
    olf[3]->lDoc.push_back( doc + " # open, right halfspace" );
    olf[4]->lDoc.push_back( doc + " # elastic" );
    olf[5]->lDoc.push_back( doc + " # inelastic" );
    olf[6]->lDoc.push_back( doc + " # refl" );
    olf[7]->lDoc.push_back( 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(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 ){
        NOlm::mem_store( olf[4] );
        NOlm::mem_store( olf[5] );
    }
    if( flag & 2 && flag & 4 ){
        NOlm::mem_store( olf[2] );
        NOlm::mem_store( olf[3] );
    }
    if( flag & 4 ){
        NOlm::mem_store( olf[0] );
        NOlm::mem_store( olf[1] );
    }
    if( flag & 1 ){
        NOlm::mem_store( olf[7] );
    }
    if( ! (flag & 8 ) ){
        NOlm::mem_store( 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, name;
    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);
        ifstream F_in;
        F_in.open(fnam.c_str(),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]->lDoc.push_back( "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", "") );
    }
    triv::fname_divide( fser, 0, &name, 0 );
    olf[0]->name = name;
    olf[1]->name = name+"open";
    olf[0]->lDoc.push_back( "refl, average over both halfspaces" );
    olf[1]->lDoc.push_back( "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( 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 ){
        NOlm::mem_store( olf[1] );
    }
    if( 1 ){
        NOlm::mem_store( olf[0] );
    }
}