//**************************************************************************//
//* 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"
#include "yamc.h"
using namespace std;
class RssmRawFile {
public:
int RdRawXml( string fnam );
void RdRawYam( FILE *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;
double temp0, temp1;
double final_setp;
int maj_syntax, min_syntax, maj_outform, min_outform;
};
// Read raw data file, store contents as class variables
void RssmRawFile::RdRawYam( FILE *F_in )
{
COld old;
string lin, key, val;
CCoord co;
CScan sout;
int iz;
double num;
daq_time_step=0;
angles.clear();
CYamlRead y( F_in );
old.as_on_disk = true;
NOlm::OloAdd(&old);
y.checktype( YAML_STREAM_START_EVENT, "no stream\n" );
y.checktype( YAML_DOCUMENT_START_EVENT, "no document\n" );
y.checktype( YAML_MAPPING_START_EVENT, "main: no map\n" );
y.checkvalue( "Meta" );
y.checktype( YAML_MAPPING_START_EVENT, "Meta: no map\n" );
y.checkvalue( "format" );
val = y.getscalar( "format" );
if( strncmp( val.c_str(), "acq5.0 for yaml1", 16 ) )
throw( string("Unexpected format") );
y.skipvalues( 2 );
y.checktype( YAML_MAPPING_END_EVENT, "Meta: no end of map\n" );
y.checkvalue( "History" );
y.checktype( YAML_SEQUENCE_START_EVENT, "History: no seq\n" );
y.checktype( YAML_MAPPING_START_EVENT, "History: meas: no map\n" );
y.checkvalue( "measured on SPHERES" );
old.lDoc.push_back( "measured on SPHERES" );
y.checktype( YAML_MAPPING_START_EVENT, "History: meas: meas: no map\n" );
while( true ){
y.next();
if( y.event.type == YAML_MAPPING_END_EVENT ){
y.done();
break;
}
if( y.event.type != YAML_SCALAR_EVENT )
throw string( "No scalar key in Hist:meas "
"; last valid entry was \"" + y.last_valid + "\"" );
key = (const char*)y.event.data.scalar.value;
y.last_valid = key;
y.done();
y.next();
if( y.event.type != YAML_SCALAR_EVENT )
throw string( "No scalar value in Hist:meas "
"; last valid entry was \"" + y.last_valid + "\"" );
val = (const char*)y.event.data.scalar.value;
y.last_valid = val;
y.done();
old.lDoc.push_back( " " + key + ": " + val );
if ( key=="incremental" ){
if ( val=="true" )
incremental = true;
else if ( val=="false" )
incremental = false;
else
throw "Invalid value of 'incremental'";
} else if ( key=="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 );
}
}
y.checktype( YAML_MAPPING_END_EVENT, "History: meas: no end of map\n" );
while( y.scalar_from_seq( "history", &lin ) ){
old.lDoc.push_back(lin);
}
y.checkvalue( "Param" );
y.checktype( YAML_SEQUENCE_START_EVENT, "Param: no seq\n" );
while( true ){
y.next();
if( y.event.type == YAML_SEQUENCE_END_EVENT )
break;
if( y.event.type != YAML_MAPPING_START_EVENT )
throw string( "Param: no map; last entry was " + y.last_valid );
y.done();
key = y.getscalar( "Param: name or note" );
if( key=="note" ){
y.skipvalues( 1 );
y.checktype( YAML_MAPPING_END_EVENT, "Param: no end of note" );
continue;
}
if( key!="name" )
throw string( "Param entry starting neither with note, "
"nor with name" );
co.name = y.getscalar( "param.name" );
if ( co.name=="daq_time_step" ){
// useful double value
y.checkvalue( "unit" );
co.unit = y.getscalar( "param(" + co.name + ").unit" );
y.checkvalue( "value" );
val = y.getscalar( "param(" + co.name + ").value" );
if( mystd::any2dbl( val, &num ) )
throw string( "param(" ) + co.name + "): invalid value " + val;
old.RPar.push_back( CParam( co, num ) );
y.checktype( YAML_MAPPING_END_EVENT,
"Param "+co.name+": no eom\n" );
if( co.name=="daq_time_step" )
daq_time_step = num;
} else {
// no useful value
while( true ){
y.next();
if( y.event.type == YAML_MAPPING_END_EVENT )
break;
y.done();
}
}
}
y.done();
y.checkvalue( "Detectors" );
y.checktype( YAML_MAPPING_START_EVENT, "Detectors: no map\n" );
y.checkvalue( "note" );
y.skipvalues( 1 );
y.checkvalue( "angles" );
y.checktype( YAML_SEQUENCE_START_EVENT, "Dets:angles: no seq\n" );
while( true ){
y.next();
if( y.event.type == YAML_SEQUENCE_END_EVENT )
break;
if( y.event.type != YAML_SCALAR_EVENT )
throw string ( "Dets:angles: no scalar" );
val = (const char*)y.event.data.scalar.value;
y.done();
if( mystd::any2dbl( val, &num ) )
throw "Dets:angles: invalid value " + val;
}
y.done();
key = y.getscalar( "Note ?" );
}
// Read raw data file, store contents as class variables
int RssmRawFile::RdRawXml( 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 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 xml file");
RssmRawFile R;
int ret;
FILE *F_in;
if ( format==2 ){
ret = R.RdRawXml( file_f );
if( ret==1 )
fprintf( stderr, "Data file not found, or failed to parse.\n" );
} else if ( format==5 ){
if( !(F_in = fopen(file_f.c_str(), "r")) ) {
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) fclose(F_in);
return;
}
} else {
printf( "invalid format" );
return;
}
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: "
"left %i right %i 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_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 format, int flag )
// format values:
// as above
// 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() );
if ( format==2 ) {
ret = RR.back().RdRawXml( fser+"c"+strg(isub) );
if( ret==1 ){
printf( "ignore the above warning\n" );
RR.pop_back();
break;
}
} else if( format==5 ) {
string fnam = fser+"c"+strg(isub);
FILE *F_in;
if( !(F_in = fopen( fnam.c_str(), "r" ) ) )
break;
try{
RR.back().RdRawYam( F_in );
}
catch(exception& e) {
cerr << e.what() << endl;
if (F_in) fclose(F_in);
return;
}
} else {
printf( "invalid format" );
return;
}
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();
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 ||
RR[isub].daq_time_step!=tstep ||
RR[isub].incremental!=incremental_in ){
printf( "SEVERE essential parameters are inconsistent\n" );
return;
}
}
// 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( "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( "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, (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_back(
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_back(
x, count/tmeas/tstep );
}
}
}
}
}
NOlm::SelNew();
if( flag & 1 ){
NOlm::OloAdd(&olf[1]);
}
if( 1 ){
NOlm::OloAdd(&olf[0]);
}
}