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 <string>
#include <vector>
#include <algorithm>

#include <libxml/xmlmemory.h>
#include <libxml/parser.h>

#include "olm.h"
#include "asi.h"
#include "gar.h"
#include "rssm.h"

using namespace std;

class RssmRawFile {
 public:
    int ReadRaw( string fnam );
    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;
    double temp0, temp1;
    double final_setp;
    int    maj_syntax, min_syntax, maj_outform, min_outform;
};

// Read raw data file, store contents as class variables
int RssmRawFile::ReadRaw( string fnam )
{
    xmlChar *key, *val, *txt, *syntax, *outform, *format, *axis,
        *contents, *categ, *txt1, *txt2;
    xmlDocPtr doc;
    xmlNodePtr cur, sub;
    char date_string[24];
    struct tm date_broken;

    if( !(doc = xmlParseFile( fnam.c_str() ) ) ){
        xmlFreeDoc( doc );
        return 1;
    }

    if( !(cur = xmlDocGetRootElement( doc ) ) ){
        fprintf( stderr, "Data file has no root element.\n" );
        xmlFreeDoc( doc );
        return -1;
    }
    if( xmlStrcmp( cur->name, (const xmlChar*) "data" ) ){
        fprintf( stderr, "Data file of the wrong type, root node != data\n" );
        xmlFreeDoc( doc );
        return -1;
    }
    if( (syntax = xmlGetProp( cur, (const xmlChar*)"syntax" )) ){
        if( !xmlStrcmp( syntax, (const xmlChar*)"wuml 4.1" ) ){
            maj_syntax = 4;
            min_syntax = 1;
        } else {
            fprintf( stderr,
                     "Recognition of new syntax not yet implemented\n" );
            xmlFreeDoc( doc );
            return -1;
        }
        if( !(outform = xmlGetProp( cur, (const xmlChar*)"generator" )) ){
            fprintf( stderr, "Missing 'generator' information\n" );
            xmlFreeDoc( doc );
            return -1;
        }
        if( !xmlStrcmp( outform, (const xmlChar*)"acq 4.1" ) ){
            maj_outform = 4;
            min_outform = 1;
        } else {
            fprintf( stderr,
                     "Recognition of new output not yet implemented\n" );
            xmlFreeDoc( doc );
            return -1;
        }
    } else {
        if( !(format = xmlGetProp( cur, (const xmlChar*)"format" )) ){
            fprintf( stderr, "Old file has no format specification\n" );
            xmlFreeDoc( doc );
            return -1;
        }
        if( !xmlStrcmp( format, (const xmlChar*)"2007-04-12" ) ){
            maj_syntax = 2;
            min_syntax = 0;
            maj_outform = 2;
            min_outform = 0;
        } else if( !xmlStrcmp( format, (const xmlChar*)"2008-03-03" ) ){
            maj_syntax = 3;
            min_syntax = 0;
            maj_outform = 3;
            min_outform = 0;
        } else {
            fprintf( stderr, "Old file has unexpected format\n" );
            xmlFreeDoc( doc );
            return -1;
        }
    }

    // *** scan parameters *** 

    daq_time_step=0;
    incremental = true;
    measured_until = -999;
    temp0 = -1;
    temp1 = -1;
    final_setp = -1;

    // Scan parameters in header

    sub = cur->xmlChildrenNode;
    while( sub != NULL ){
        if( (!xmlStrcmp( sub->name, (const xmlChar*)"par" ) ) ){
            key = xmlGetProp( sub, (const xmlChar*)"key" );
            val = xmlGetProp( sub, (const xmlChar*)"value" );
            if        ( !xmlStrcmp( key, (const xmlChar*)"temp.temp0" ) ){
                sscanf( (const char*)val, "%lg", &temp0 );
            } else if ( !xmlStrcmp( key, (const xmlChar*)"temp.temp1" ) ){
                sscanf( (const char*)val, "%lg", &temp1 );
            } else if ( !xmlStrcmp( key, (const xmlChar*)"final_setp" ) ){
                sscanf( (const char*)val, "%lg", &final_setp );
            } else if ( !xmlStrcmp( key, (const xmlChar*)"daq_time_step" ) ){
                sscanf( (const char*)val, "%lg", &daq_time_step );
            } else if ( !xmlStrcmp( key, (const xmlChar*)"incremental" ) ){
                sscanf( (const char*)val, "%i", &incremental );
            } else if ( !xmlStrcmp( key, (const xmlChar*)"incremental save" ) 
                ){ // PROVISIONAL: only for data before R51-6
                sscanf( (const char*)val, "%i", &incremental );
            } else if ( !xmlStrcmp( key, (const xmlChar*)"measured_until" ) ){
                sscanf( (const char*)val, "%24c", date_string );
                strptime( date_string, "%F %A %H:%M:%S", &date_broken );
                measured_until = mktime( &date_broken );
                // strftime( tstrg, 30, "%F %b %H:%M:%S",
                //           localtime( &measured_until ) );
                // printf( "measured until %s\n", tstrg );
            // PROVISIONAL: later data include measured_from
            }
            xmlFree( key );
            xmlFree( val );
        }
	sub = sub->next;
    }
    if( !daq_time_step ){
        fprintf( stderr, "Data file has no parameter \"daq_time_step\".\n" );
        xmlFreeDoc( doc );
        return -1;
    }

    // Scan arrays

    sub = cur->xmlChildrenNode;
    int iaxisE=0, iaxisP=4, iaxis;
    if( maj_outform<=2 ){
        angles.resize(19);
        for( int j=0; j<19; ++j )
            angles[j] = j;
    }
    while( sub != NULL ){
        if( (!xmlStrcmp(sub->name, (const xmlChar*)"array" ) ) ){
            axis = xmlGetProp( sub, (const xmlChar*)"xaxis" );
            contents = xmlGetProp( sub, (const xmlChar*)"contents" );
            categ = xmlGetProp( sub, (const xmlChar*)"categorization" );
            if( !axis && !contents ){
                printf( "array has neither old axis nor new contents "
                        "attribute\n" );
                return -1;
            } else if( axis && contents ){
                printf( "array has conflicting attributes: "
                        "old axis and new contents\n" );
                return -1;
            } else if( !xmlStrcmp( axis, (const xmlChar*)"angle" ) ||
                       !xmlStrcmp( contents, 
                                   (const xmlChar*)"detector-info" ) ){
                if( axis ) xmlFree( axis );
                if( contents ) xmlFree( contents );
                txt = xmlNodeListGetString( doc, sub->xmlChildrenNode, 1 );
                vector<double> aux;
                if( mystd::str2vec( (const char*)txt, &aux, 0, false) ){
                    printf( "str2vec( array text ) failed\n" );
                }
                xmlFree( txt );
                angles.clear();
                for( int i=0; i<aux.size(); i+=(axis? 2 : 3 ) )
                    angles.push_back( aux[i] );
                sub = sub->next;
                continue;
            } else if( !xmlStrcmp( axis, (const xmlChar*)"energy" ) ||
                       !xmlStrcmp( contents, (const xmlChar*)"histogram" ) &&
                       !xmlStrcmp( categ, (const xmlChar*)"energy" ) ){
                txt1 = xmlGetProp( sub, (const xmlChar*)"halfspace" );
                if( txt1==NULL )
                    txt1 = 
                        xmlGetProp( sub, (const xmlChar*)"Doppler_halfspace" );
                txt2 = xmlGetProp( sub, (const xmlChar*)"chopper_phase" );
                if( !txt1 || !txt2 ){
                    printf( "array with erg axis missing some attribute\n" );
                    return -1;
                }
                iaxis = iaxisE++;
                if( !(iaxis%2==0 && !xmlStrcmp(txt1,(const xmlChar*)"left") ||
                      iaxis%2==1 && !xmlStrcmp(txt1,(const xmlChar*)"right"))){
                    printf( "unexpected halfspace in energy array\n" );
                    // return -1;
                }
                if( !(iaxis/2==0 && !xmlStrcmp(txt2,(const xmlChar*)"refl") ||
                      iaxis/2==1 && !xmlStrcmp(txt2,(const xmlChar*)"open")) ){
                    printf( "unexpected chopper_phase in energy array\n" );
                    return -1;
                }
                xmlFree( txt1 );
                xmlFree( txt2 );
            } else if( !xmlStrcmp( axis, (const xmlChar*)"chop_phase" ) ||
                       !xmlStrcmp( contents, (const xmlChar*)"histogram" ) &&
                       !xmlStrcmp( categ, (const xmlChar*)"chopper_phase" ) ){
                txt1 = xmlGetProp( sub, (const xmlChar*)"energy_range" );
                if( !txt1 ){
                    printf( "array with cho axis missing energy_range\n" );
                    return -1;
                }
                iaxis = iaxisP++;
                if( !(iaxis==4 && !xmlStrcmp(txt1,(const xmlChar*)"elast") ||
                      iaxis==5 && !xmlStrcmp(txt1,(const xmlChar*)"inela")) ){
                    printf( "unexpected energy_range in chopper array\n" );
                    return -1;
                }
                xmlFree( txt1 );
            } else {
                printf( "unexpected array axis\n" );
                return -1;
            }
            xmlFree( axis );
            txt = xmlNodeListGetString( doc, sub->xmlChildrenNode, 1 );
            if( mystd::str2vec( (const char*)txt,
                                &(rawdata[iaxis]), 0, false) ){
                printf( "str2vec( array text ) failed\n" );
            }
            xmlFree( txt );
        }
	sub = sub->next;
    }
	
    xmlFreeDoc( doc );
    return 0;
}    

void NRSSM::ReadScan( 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

{
    string file_f, name;
    file_f = wask("Read SPHERES data from xml file");
    RssmRawFile R;
    int ret = R.ReadRaw( file_f );
    if( ret==1 )
        fprintf( stderr, "Data file not parsed successfully.\n" );
    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].lDoc.push_back( "acquire data from "+file_f );
        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";
    olf[0].lDoc.push_back( "refl, left halfspace" );
    olf[1].lDoc.push_back( "refl, right halfspace" );
    olf[2].lDoc.push_back( "open, left halfspace" );
    olf[3].lDoc.push_back( "open, right halfspace" );
    olf[4].lDoc.push_back( "elastic region" );
    olf[5].lDoc.push_back( "inelastic region" );
    olf[6].lDoc.push_back( "refl, average over both halfspaces" );
    olf[7].lDoc.push_back( "open, average over both halfspaces" );

    CScan S;
    for( int j=0; j<ndet; ++j ){
        S.Clear();
        S.z.resize( 1 );
        S.z[0] = R.angles[j];
        for( iolf=0; iolf<8; ++iolf )
            olf[iolf].VS.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_back( 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_back( 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_back( 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+=(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_back( 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 flag )

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

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

    // Read consolidated files:
    fser = wask("Read SPHERES data from series");
    for( isub=0; ; ++isub ){
        RR.push_back( RssmRawFile() );
        ret = RR.back().ReadRaw( fser+"c"+strg(isub) );
        if( ret==1 ){
            printf( "ignore the above warning\n" );
            RR.pop_back();
            break;
        }
        if( ret )
            return;
    }
    int nsub = RR.size();
    if( nsub==0 )
        return;
    else if( nsub==1 ){
        printf( "series contains just one scan -> use a simpler method\n" );
        return;
    }
    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();
    double tstep = RR[0].daq_time_step;
    bool incremental_in = RR[0].incremental;
    int nE = RR[0].rawdata[0].size()/(ndet+2);

    for( isub=1; isub<nsub; ++isub ){
        if( RR[isub].angles.size()!=ndet ||
            RR[isub].daq_time_step!=tstep ||
            RR[isub].incremental!=incremental_in ||
            RR[0].rawdata[0].size()/(ndet+2)!=nE ){
            printf( "SEVERE essential parameters are inconsistent\n" );
            return;
        }
    }

    // Set file headers:
    COld olf[2];
    int iolf, j;

    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( "T0", "K") );
        olf[iolf].ZCo.push_back( CCoord( "T1", "K") );
        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 ){
            CScan S;
            S.z.resize( 3 );
            S.z[0] = RR[isub].temp0;
            S.z[1] = RR[isub].temp1;
            S.z[2] = RR[isub].angles[j];
            for( iolf=0; iolf<2; ++iolf )
                olf[iolf].VS.push_back( S );
        }
    }

    if( RR[0].maj_outform>=3 ){
        printf( "series / format080303 not yet implemented\n" );
        return;
    }

    // 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_back( 
                        x, count/tmeas/tstep );
                }
            }
        }
    }
    
    NOlm::SelNew();
    if( flag & 1 ){
        NOlm::OloAdd(&olf[1]);
    }
    if( 1 ){
        NOlm::OloAdd(&olf[0]);
    }

}