rssm.cpp

//**************************************************************************//
//* FRIDA: fast reliable interactive data analysis                         *//
//* rssm.cpp: read SPHERES raw data                                        *//
//* (C) Joachim Wuttke 1990-, v2(C++) 2001-                                *//
//* http://frida.sourceforge.net                                           *//
//**************************************************************************//

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <iostream>
#include <fstream>

#include "mystd.h"
#include "olm.h"
#include <ask_simple_value.h>
#include "rssm.h"
#include <yaml-cpp.h>
#include <string.h>

using namespace std;

class RssmRawFile {
 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;
    vector<double> rawdata[6];
    vector<double> angles;
    int    maj_syntax, min_syntax, maj_outform, min_outform;
};

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

//void RssmRawFile::RdRawYam( FILE *F_in )
void RssmRawFile::RdRawYam( ifstream& F_in )
{
    string lin, key, val, bla, blub;
    CCoord co;
    CSpec sout;
    int ndet;
    double num;
    vector<double> numvec;

    maj_outform = 5;
    min_outform = 1;

    daq_time_step=0;
    angles.clear();

    YAML::Parser parser(F_in);
    YAML::Node doc;
    parser.GetNextDocument(doc);

    //start to read Meta
    if(!doc.FindValue("Meta"))
        throw string("no Meta" );
    try {
        YAML::CONTENT_TYPE Metatype = doc["Meta"].GetType();
        if ( Metatype != YAML::CT_MAP && doc["Meta"].size() > 0)
            throw string("Meta is not a MAP type" );
        if(!doc["Meta"].FindValue("format"))
            throw string("no format in Meta" );
        doc["Meta"]["format"] >> val;
        if( strncmp( val.c_str(), "acq5.1 for yaml1", 16 ) )
            throw( string("Unexpected format: ") + val + " instead of acq5.1" );
    } catch(YAML::ParserException& e) {
        std::cout << "no Meta" << e.what() << "\n";
    }
    
    //do nothing in History
    
    //start to read Shortpar
    if(!doc.FindValue("Shortpar"))
        throw string("no Shortpar" );
    for(YAML::Iterator it=doc["Shortpar"].begin();
    it!=doc["Shortpar"].end();++it) {
        it.first() >> key;
        it.second() >> val;
        if      ( key=="incremental" ){
            if      ( val=="true" )
                incremental = true;
            else if ( val=="false" )
                incremental = false;
            else
                throw string( "Invalid value [" ) + val +
                "] of 'incremental'";
        } else if ( key=="measured_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( !mystd::any2int( val, &ndet ) )
                throw string( "invalid ndet" );
        } else if ( key=="daq_time_step" ){
            mystd::string_extract_word( val, &bla, &blub );
            mystd::any2dbl( bla, &daq_time_step );
        }
    }
    if( !daq_time_step )
        throw string( "daq time step not found" );
    
    //do nothing in Longpar
    
    //start to read Detectors
    if(!doc.FindValue("Detectors"))
        throw string("no Detectors" );
    if ( doc["Detectors"].GetType() != YAML::CT_SEQUENCE &&
    doc["Detectors"].size() > 0 )
        throw string("Detectors is not a SEQUENCE type" );
    for (unsigned int iDet = 0; iDet < (doc["Detectors"].size());
    iDet++) {
        if ( doc["Detectors"][iDet].GetType() != YAML::CT_MAP )
            throw string("Detectors " + strg(iDet) +
             " is not a MAP type" );
        if( const YAML::Node *pName = 
        doc["Detectors"][iDet].FindValue("angle") ){
            *pName >> val;
            if( !mystd::any2dbl( val, &num ) ){
                num = 0;
                printf( "warning: invalid angle %s\n", 
                val.c_str() );
            }
        }else
            throw string( "missing angle" );
        // path_offset ?
        angles.push_back( num );
    }
    
    //start to read EnergyHistograms
    if(!doc.FindValue("EnergyHistograms"))
        throw string("no EnergyHistograms" );
    if ( doc["EnergyHistograms"].GetType() != YAML::CT_SEQUENCE &&
         doc["EnergyHistograms"].size() > 0 )
        throw string("EnergyHistograms is not a SEQUENCE type" );
    for (unsigned int iEne = 0; 
         iEne < (doc["EnergyHistograms"].size()); iEne++) {
        if ( doc["EnergyHistograms"][iEne].GetType() != YAML::CT_SEQUENCE )
            throw "EnergyHistograms " + strg(iEne) + 
                  " is not a SEQUENCE type"; 
        doc["EnergyHistograms"][iEne][0] >> val;
        if( !mystd::any2dbl( val, &num ) )
            throw string( "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 ){
            doc["EnergyHistograms"][iEne][i+1] >> val;// cnts line
            if( !mystd::str2vec( val, &numvec ) )
                throw string( "cnts line extraction failed" );
            if( numvec.size()!=ndet )
                throw string( "hist: ndet incompatibility" );
            for( int ii=0; ii<numvec.size(); ++ii )
                rawdata[i].push_back( numvec[ii] );
        }
        for( int i=0; i<4; ++i ){
            doc["EnergyHistograms"][iEne][i+1+4] >> val;// time line
            if( !mystd::str2vec( val, &numvec ) )
                throw string( "time line extraction failed" );
            for( int ii=0; ii<numvec.size(); ++ii )
                rawdata[i].push_back( numvec[ii] );
        }
    }
    
    //start to read ChopperHistograms
    if(!doc.FindValue("ChopperHistograms"))
        throw string("no ChopperHistograms" );
    if ( doc["ChopperHistograms"].GetType() != YAML::CT_SEQUENCE &&
         doc["ChopperHistograms"].size() > 0 )
        throw string("ChopperHistograms is not a SEQUENCE type" );
    for (unsigned int iCho = 0; 
         iCho < (doc["ChopperHistograms"].size()); iCho++) {
        if ( doc["ChopperHistograms"][iCho].GetType() != YAML::CT_SEQUENCE )
            throw string("ChopperHistograms " + strg(iCho) +
                   " is not a SEQUENCE type" );         
        doc["ChopperHistograms"][iCho][0] >> val;
        if( !mystd::any2dbl( val, &num ) )
            throw string( "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
            doc["ChopperHistograms"][iCho][i-4+1] >> val;//cnts line
            if( !mystd::str2vec( val, &numvec ) )
                throw string( "cnts line extraction failed" );
            if( numvec.size()!=ndet )
                throw string( "hist: ndet incompatibility" );
            for( int ii=0; ii<numvec.size(); ++ii )
                rawdata[i].push_back( numvec[ii] );
        }
        for( int i=4; i<6; ++i ){            // 2 time lines
            doc["ChopperHistograms"][iCho][i-4+1+2] >> val;//time line
            if( !mystd::str2vec( val, &numvec ) )
                throw string( "time line extraction failed" );
            for( int ii=0; ii<numvec.size(); ++ii )
                rawdata[i].push_back( numvec[ii] );
        }
    }
}    

void NRSSM::ReadSpec( int format, int flag )

// format values:
//    2 acq2-4 xml
//    5 acq5   yaml
// 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

{
    string file_f, name;
    file_f = wask("Read SPHERES data from file");
    RssmRawFile R;
    int ret;
    ifstream F_in;

    if        ( format==2 ){
        throw string( "XML support has been withdrawn from frida2 on 19may09" );
    } else if ( format==5 ){
        F_in.open (file_f.c_str(),fstream::in);
        if( !(F_in.is_open()) ) {
            cerr << "Cannot read " << file_f << "\n";
            return;
        }
        try{
            R.RdRawYam( F_in );
        }
        catch( string &e) {
            cerr << "Cannot read " << file_f << ":\n";
            cerr << e << endl;
            if (F_in.is_open()) F_in.close();
            return;
        }
        ret = 0;
    } else {
        throw string( "invalid format" );
    }
    if( ret )
        return;

    double tstep = R.daq_time_step;
    vector<double> *raw = R.rawdata;
    int ndet = R.angles.size();

    // *** set file headers ***

    COld olf[8];
    int iolf;
    
    for( iolf=0; iolf<8; ++iolf ){
        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", "") );
    }
    mystd::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 ? "" : strg(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( int j=0; j<ndet; ++j ){
        PSpec S( new CSpec );
        S->z.resize( 1 );
        S->z[0] = R.angles[j];
        for( iolf=0; iolf<8; ++iolf )
            olf[iolf].V.push_back(S);
    }

    if( R.maj_outform<=2 ){
        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 ){
                printf( "inconsistent raw array sizes\n" );
                return;
            }
            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 ){
                    tmeas = raw[iolf][iraw+1+ndet+j];
                    if( tmeas<1 ) continue;
                    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 ){
                printf( "inconsistent raw array sizes: "
                        "left %zu right %zu lines\n",
                        (raw+iolf0)->size(), (raw+iolf1)->size() );
                return;
            }
            for( int iraw=0; iraw<nline; iraw+=(1+ndet*2) ){
                x = raw[iolf0][iraw];
                for ( int j=0; j<ndet; ++j ){
                    tmeas = raw[iolf0][iraw+1+ndet+j] +
                        raw[iolf1][iraw+1+ndet+j];
                    if( tmeas<1 ) continue;
                    count = raw[iolf0][iraw+1+j] + raw[iolf1][iraw+1+j];
                    olf[6+ichop].VS(j)->push_xy( x, count/tmeas/tstep );
                }
            }
        }
    }
    
    NOlm::SelNew();
    if( flag & 2 ){
        NOlm::OloAdd(&olf[4]);
        NOlm::OloAdd(&olf[5]);
    }
    if( flag & 2 && flag & 4 ){
        NOlm::OloAdd(&olf[2]);
        NOlm::OloAdd(&olf[3]);
    }
    if( flag & 4 ){
        NOlm::OloAdd(&olf[0]);
        NOlm::OloAdd(&olf[1]);
    }
    if( flag & 1 ){
        NOlm::OloAdd(&olf[7]);
    }
    if( ! (flag & 8 ) ){
        NOlm::OloAdd(&olf[6]);
    }

}


void NRSSM::ReadSeries( int format, int flag )

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

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

    // Read consolidated files:
    fser = wask("Read SPHERES data from series");
    for( isub=0; ; ++isub ){
        if       ( format==2 ) {
            throw string(
                "XML support has been withdrawn from frida2 on 19may09" );
        } else if( format==5 ) {
            fnam = fser+"a"+strg(isub);
            ifstream F_in;
            F_in.open(fnam.c_str(),fstream::in);
            if( !(F_in.is_open()) )
                break;
            // cout << "successfully opened "<<fnam<<"\n";
            try{ 
                RssmRawFile rf; 
                rf.RdRawYam( F_in );
                RR.push_back( rf );
            }
            catch(exception& e) {
                cerr << e.what() << endl;
                if (F_in.is_open()) F_in.close();
                return;
            }
            ret = 0;
        } else {
            throw string( "invalid format parameter" );
        } 
        if ( ret )
            return;
    }
    int nsub = RR.size();
    if( nsub==0 )
        throw string( "could not open file " ) + fnam;
    else if( nsub==1 )
        throw string( "series contains just one spectrum -> use a simpler method");
    printf( "successfully read %d 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);
    strftime( tstrg, 30, "%F %b %H:%M:%S", localtime( &t ) );
    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].angles.size();
    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].angles.size()!=ndet )
            throw string( "inconsistent no of det: file " ) + strg(isub)
                + " has " + strg(RR[isub].angles.size()) + " instead of "
                + strg(ndet);
        if( RR[isub].daq_time_step!=tstep )
            throw string( "inconsistent time step" );
        if( RR[isub].incremental!=incremental_in )
            throw string( "inconsistent increment status" );
    }

    // Set file headers:
    COld olf[2];
    int iolf, 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", "") );
    }
    mystd::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] = isub;
            S->z[1] = RR[isub].angles[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 %zu nE %i\n",
                            ichop, (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 );
                    }
                }
            }
        }
    }

    NOlm::SelNew();
    if( flag & 1 ){
        NOlm::OloAdd(&olf[1]);
    }
    if( 1 ){
        NOlm::OloAdd(&olf[0]);
    }

}