Forked from
mlz / Frida
1841 commits behind the upstream repository.
-
Wuttke, Joachim authoredWuttke, Joachim authored
rssm.cpp 19.12 KiB
//**************************************************************************//
//* FRIDA: flexible rapid interactive data analysis *//
//* rssm.cpp: read SPHERES raw data *//
//* (C) Joachim Wuttke 1990-, v2(C++) 2001- *//
//* http://www.messen-und-deuten.de/frida *//
//**************************************************************************//
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include <yaml-cpp/yaml.h>
#include <ask_simple_value.h>
#include "mystd.h"
#include "olf.h"
#include "mem.h"
#include "rssm.h"
using namespace std;
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;
uint 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( FILE *F_in )
void CRawfileSpheres::RdRawYam( ifstream& F_in )
{
string lin, key, val, bla, blub;
CCoord co;
CSpec sout;
double num;
vector<double> numvec;
maj_outform = 5;
min_outform = 1;
daq_time_step=0;
YAML::Parser parser(F_in);
YAML::Node doc;
const YAML::Node *arr;
parser.GetNextDocument(doc);
// read Meta:
if(!doc.FindValue("Meta"))
throw "no Meta";
try {
YAML::NodeType::value Metatype = doc["Meta"].Type();
if ( Metatype != YAML::NodeType::Map && doc["Meta"].size() > 0)
throw "DATA BUG: Meta is not a MAP";
if(!doc["Meta"].FindValue("format"))
throw "DATA BUG: no format in Meta";
doc["Meta"]["format"] >> val;
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.FindValue("Shortpar"))
throw "DATA BUG: 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 "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( !mystd::any2uint( val, &ndet ) )
throw string( "DATA BUG: 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( "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.FindValue("EnergyHistograms") )
throw "DATA BUG: no EnergyHistograms";
if ( doc["EnergyHistograms"].Type() != YAML::NodeType::Sequence )
throw "DATA BUG: EnergyHistograms is not a SEQUENCE";
uint nEne = doc["EnergyHistograms"].size();
for ( uint iEne = 0; iEne<nEne; iEne++ ) {
if ( doc["EnergyHistograms"][iEne].Type() !=
YAML::NodeType::Sequence )
throw "DATA BUG: EnergyHistogram[" + strg(iEne) +
"] is not a SEQUENCE";
doc["EnergyHistograms"][iEne][0] >> val;
if( !mystd::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
arr = &(doc["EnergyHistograms"][iEne][i+1]);
if( arr->Type() != YAML::NodeType::Sequence )
throw "DATA BUG: EnergyHistogram[" + strg(iEne) + "][" +
strg(i+1) + "] is not a SEQUENCE";
for( int j=0; j<ndet; ++j ) {
*arr >> val;
if( !mystd::any2dbl( val, &num ) )
throw "E-Hist: invalid count " + val;
rawdata[i].push_back( num ); }
}
for( int i=0; i<4; ++i ){ // tsteps
arr = &(doc["EnergyHistograms"][iEne][i+1+4]);
if( arr->Type() != YAML::NodeType::Sequence )
throw "DATA BUG: EnergyHistogram[" + strg(iEne) + "][" +
strg(i+1+4) + "] is not a SEQUENCE";
for( int j=0; j<ndet; ++j ) {
*arr >> val;
if( !mystd::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.FindValue("ChopperHistograms"))
throw string("no ChopperHistograms" );
if ( doc["ChopperHistograms"].Type() != YAML::NodeType::Sequence )
throw "DATA BUG: ChopperHistograms is not a SEQUENCE";
uint nCho = doc["ChopperHistograms"].size();
for (uint iCho = 0; iCho<nCho; iCho++) {
if ( doc["ChopperHistograms"][iCho].Type() != YAML::NodeType::Sequence )
throw "DATA BUG: ChopperHistograms " + strg(iCho) +
" is not a SEQUENCE";
doc["ChopperHistograms"][iCho][0] >> val;
if( !mystd::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
doc["ChopperHistograms"][iCho][i-4+1] >> val;//cnts line
if( !mystd::str2vec( val, &numvec ) )
throw "DATA BUG: cnts line extraction failed";
if( numvec.size()!=ndet )
throw "DATA BUG: 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 "DATA BUG: time line extraction failed";
for( int ii=0; ii<numvec.size(); ++ii )
rawdata[i].push_back( numvec[ii] );
}
}
} catch( exception& e ) {
throw "DATA BUG: failed to read ChopperHistograms section: " +
string(e.what());
}
}
void NRSSM::ReadSpec( 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 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;
uint 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", "") );
}
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( iolf=0; iolf<8; ++iolf ){
for( int j=0; j<ndet; ++j ){
PSpec S( new CSpec );
S->z.resize( 1 );
S->z[0] = 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 string( "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 ){
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 )
throw string( "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 ){
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 );
}
}
}
}
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, 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"+strg(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 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].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 string( "inconsistent no of det: file " ) + strg(isub)
+ " has " + strg(RR[isub].ndet) + " 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:
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", "") );
}
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] = 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 );
}
}
}
}
}
if( flag & 1 ){
NOlm::OloAdd( olf[1] );
}
if( 1 ){
NOlm::OloAdd( olf[0] );
}
}