dataFiles.cpp

Go to the documentation of this file.

/* dataFiles.cpp: Atom probe data file handling routines
 * Copyright (C) 2014  Daniel Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "atomprobe/io/dataFiles.h"
#include "atomprobe/primitives/ionHit.h"
#include "atomprobe/helper/stringFuncs.h"
#include "atomprobe/helper/endianTest.h"
#include "atomprobe/helper/sampling.h"
#include "helper/helpFuncs.h"

#include "atomprobe/helper/aptAssert.h"

#include <vector>
#include <string>
#include <fstream>
#include <map>
#include <algorithm>
#include <sys/time.h>


namespace AtomProbe{

using std::vector;
using std::string;
using std::ifstream;
using std::ofstream;
 
//Pos file structure
//===========
typedef struct IONHIT
{
float pos[3];
float massToCharge;
} IONHIT;
//=========


enum
{
    //Pos file constants
    POS_ALLOC_FAIL=1,
    POS_OPEN_FAIL,
    POS_SIZE_MODULUS_ERR,
    POS_SIZE_EMPTY_ERR, 
    POS_READ_FAIL,
    POS_NAN_LOAD_ERROR,
    POS_FILE_ENUM_END
};

enum
{
    //tapsim loading error constants
    TAPSIM_FILE_FORMAT_FAIL=1,
    TAPSIM_OPEN_FAIL,
};


const char *ATO_ERR_STRINGS[] = { "",
                "Error opening file",
                "File is empty",
                "Filesize does not match expected format",
                "File version number not <4, as expected",
                "Unable to allocate memory to store data",
                "Unable to detect endian-ness in file"
                };

const char *RECORDREAD_ERR_STRINGS[] = { "",
                    "Unable to determine filesize",
                    "Filesize indicates that file contains a non-integer number of entries",
                    "Unable to open file",
                    "Unable to allocate memory for reading file contents",
                    "Unable to perform read  operation on file",
                    "Read past end of file requested",
                    "Entry in file appears to be invalid",
    };

template<class T>
unsigned int fixedRecordReader(const char *filename, bool (*recordReader)(const char *bufRead, const char *destBuf),
            size_t recordSize, std::vector<T> &outputData)
{
    size_t fileSize;
    if(!getFilesize(filename, fileSize))
        return RECORDREAD_ERR_GET_FILESIZE;
    
    
    //Check we have an even number of entries
    if(fileSize % recordSize )
        return RECORDREAD_ERR_FILESIZE_MODULO;

    std::ifstream f(filename);

    if(!f)
        return RECORDREAD_ERR_FILE_OPEN;

    
    size_t numFileEntries = fileSize/recordSize;
    
    char *buffer;   
    try{
        buffer = new char[recordSize];
        outputData.resize(numFileEntries);
    }
    catch(std::bad_alloc)
    {
        return RECORDREAD_ERR_NOMEM;
    }

    for(size_t ui=0;ui<numFileEntries;ui++)
    {
    
        if(!f.read(buffer,recordSize))
        {
            delete[] buffer;
            return RECORDREAD_BAD_FILEREAD;
        }
    
        if(!(*recordReader)(buffer,(char *)(&outputData[ui])))
        {
            delete[] buffer;
            return RECORDREAD_BAD_RECORD;
        }
    }
    
    delete[] buffer;    
    return 0;
}

template<class T>
unsigned int fixedRecordChunkReader(const char *filename, bool (*recordReader)(const char *bufRead, const char *destBuf),
    size_t recordSize, std::vector<T> &outputData, 
    unsigned int chunkSize,unsigned int chunkOffset, unsigned int &nEntriesLeft)
{
    //Chunk size is in units of records
    //Chunk offset is in units of chunks
    //nEntriesLeft is in units of records

    //Obtain the size of the file
    size_t fileSize;
    if(!getFilesize(filename, fileSize))
        return RECORDREAD_ERR_GET_FILESIZE;
    
    //Check we have an even number of entries
    if(fileSize % recordSize )
        return RECORDREAD_ERR_FILESIZE_MODULO;

    //Open file
    std::ifstream f(filename);
    if(!f)
        return RECORDREAD_ERR_FILE_OPEN;

    
    size_t numFileEntries = fileSize/recordSize;
    unsigned int entryOffset=chunkOffset*chunkSize; //offset in records

    //Caller should not specify read start from beyond file end
    if(entryOffset > numFileEntries)
        return RECORDREAD_ERR_CHUNKOFFSET;


    //Read either the minimum number of entries, or a chunk of entries
    size_t nEntriesToRead = std::min(numFileEntries-entryOffset,(size_t)chunkSize); 

    //If no entries to read, finished
    if(!nEntriesToRead)
    {
        nEntriesLeft=0;
        return 0;
    }

    char *buffer;   
    try{
        buffer = new char[recordSize];
        outputData.resize(nEntriesToRead);
    }
    catch(std::bad_alloc)
    {
        return RECORDREAD_ERR_NOMEM;
    }

    //Start loading chunks
    f.seekg(entryOffset*recordSize);

    for(size_t ui=0; ui<nEntriesToRead; ui++)
    {
        //Read one record
        if(!f.read(buffer,recordSize))
        {
            delete[] buffer;
            return RECORDREAD_BAD_FILEREAD;
        }
        //call reader function though func pointer
        if(!(*recordReader)(buffer,(char *)(&outputData[ui])))
        {
            delete[] buffer;
            return RECORDREAD_BAD_RECORD;
        }
    }

    nEntriesLeft=numFileEntries-(chunkOffset*chunkSize+nEntriesToRead);

    delete[] buffer;    
    return 0;
}




void EPOS_ENTRY::getIonHit(IonHit &h) const
{
    h[0]=x;
    h[1]=y;
    h[2]=z;
    h[3]=massToCharge;
}

IonHit EPOS_ENTRY::getIonHit() const
{
    IonHit h;
    getIonHit(h);
    return h;
}

const char *getPosFileErrString(unsigned int errCode)
{
    const char *posErrStrings[] = { "",
                    "Memory allocation failure on POS load",
                    "Error opening file for load",
                    "Pos file size appears to have non-integer number of entries",
                    "Pos file empty",
                    "Error reading from pos file (after open)",
                    "Error - Found NaN in pos file",
                    ""
    };

    static_assert( ARRAYSIZE(posErrStrings)-1 == POS_FILE_ENUM_END, "POS Error string size and enum mismatched");

    return posErrStrings[errCode];
}

unsigned int loadPosFile(vector<IonHit> &posIons,const char *posFile)
{


    //buffersize must be a power of two and at least sizeof(IONHIT)
    const unsigned int BUFFERSIZE=8192;
    char *buffer=new char[BUFFERSIZE];
    
    if(!buffer)
        return POS_ALLOC_FAIL;

    //open pos file
    std::ifstream CFile(posFile,std::ios::binary);

    if(!CFile)
    {
        delete[] buffer;
        return POS_OPEN_FAIL;
    }
    
    CFile.seekg(0,std::ios::end);
    std::ifstream::pos_type fileSize=CFile.tellg();
    CFile.seekg(0,std::ios::beg);
    
    if(!fileSize)
    {
        delete[] buffer;
        return POS_SIZE_EMPTY_ERR;
    }
    
    //calculate the number of points stored in the POS file
    IonHit hit;
    unsigned int pointCount=0;
    //regular case
    std::ifstream::pos_type curBufferSize=BUFFERSIZE;
    
    if(fileSize % sizeof(IONHIT))
    {
        delete[] buffer;
        return POS_SIZE_MODULUS_ERR;    
    }
    
    while(fileSize < curBufferSize)
        curBufferSize = curBufferSize >> 1;
    
    try
    {
        posIons.resize(fileSize/(sizeof(IONHIT)));
    }
    catch(std::bad_alloc)
    {
        delete[] buffer;
        return POS_ALLOC_FAIL;
    }


    do
    {
        while(CFile.tellg() <= fileSize-curBufferSize)
        {
            ASSERT(curBufferSize >= sizeof(IONHIT));
            CFile.read(buffer,curBufferSize);
            if(!CFile.good())
            {
                delete[] buffer;
                return POS_READ_FAIL;
            }

            IONHIT *hitStruct;
            hitStruct = (IONHIT *)buffer; 


            unsigned int ui;
            for(ui=0; ui<curBufferSize; ui+=(sizeof(IONHIT)))
            {
                hit.setPos(Point3D(hitStruct->pos));
                hit.setMassToCharge(hitStruct->massToCharge);
                //Data bytes stored in pos files are big
                //endian. flip as required
                #ifdef __LITTLE_ENDIAN__
                    hit.switchEndian(); 
                #endif
                
                if(hit.hasNaN())
                {
                    delete[] buffer;
                    return POS_NAN_LOAD_ERROR;  
                }

                posIons[pointCount] = hit;
            
                pointCount++;
                hitStruct++;
            }   
                
        }

        curBufferSize = curBufferSize >> 1 ;
    }while(curBufferSize >= sizeof(IONHIT));

    ASSERT((unsigned int)CFile.tellg() == fileSize);
    ASSERT(pointCount*sizeof(IONHIT) == fileSize);
    ASSERT(pointCount == posIons.size());
    delete[] buffer;

    return 0;
}

unsigned int loadPosFile(vector<IonHit> &posIons, const char *posFile, unsigned int limitCount)
{


    //open pos file
    std::ifstream CFile(posFile,std::ios::binary);

    if(!CFile)
        return POS_OPEN_FAIL;
    
    CFile.seekg(0,std::ios::end);
    auto fileSize=CFile.tellg();

    if(!fileSize)
        return POS_SIZE_EMPTY_ERR;
    
    CFile.seekg(0,std::ios::beg);
    
    const unsigned int IONHIT_SIZE=16;  
    if(fileSize % IONHIT_SIZE)
        return POS_SIZE_MODULUS_ERR;    

    auto maxIons =fileSize/IONHIT_SIZE;
    limitCount=std::min(limitCount,(unsigned int)maxIons);

    //If we are going to load the whole file, don't use a sampling method to do it.
    if(limitCount == maxIons)
    {
        //Close the file
        CFile.close();
        //Try opening it using the normal functions
        return loadPosFile(posIons,posFile);
    }

    //Use a sampling method to load the pos file
    std::vector<size_t> ionsToLoad;
    
    //Init GSL random generator
    gsl_rng *rng = AtomProbe::randGen.getRng();
    try
    {
        //Allocate space
        posIons.resize(limitCount);

        randomIndices(ionsToLoad,limitCount,maxIons,rng);
        
    }
    catch(std::bad_alloc)
    {
        return POS_ALLOC_FAIL;
    }


    //sort ions to be in increasing size
    //NOTE: I tried to use a functor here to get progress
    // It was not stable with parallel sort 
    std::sort(ionsToLoad.begin(),ionsToLoad.end());


    //TODO: probably not very nice to the disk drive. would be better to
    //scan ahead for contiguous data regions, and load that where possible.
    //Or switch between different algorithms based upon ionsToLoad.size()/  
    std::ios::pos_type  nextIonPos;
    
    char buffer[IONHIT_SIZE];
    for(size_t ui=0;ui<ionsToLoad.size(); ui++)
    {
        nextIonPos =  ionsToLoad[ui]*IONHIT_SIZE;
        
        if(CFile.tellg() !=nextIonPos )
            CFile.seekg(nextIonPos);
    

        CFile.read(buffer,IONHIT_SIZE);
        
        if(!CFile.good())
            return POS_READ_FAIL;

        IONHIT *hitStruct;
        hitStruct = (IONHIT*)buffer;
        posIons[ui].setPos(Point3D(hitStruct->pos));
        posIons[ui].setMassToCharge(hitStruct->massToCharge);
        
        //Data bytes stored in pos files are big
        //endian. flip as required
        #ifdef __LITTLE_ENDIAN__
            posIons[ui].switchEndian(); 
        #endif
        
        if(posIons[ui].hasNaN())
            return POS_NAN_LOAD_ERROR;  
    
        
    }

    return 0;
}

unsigned int savePosFile(const vector<IonHit> &ionVec, const char *filename, bool append)
{
    auto flags = std::ios::binary;
    if(append)
        flags|=std::ios::app;
        
    std::ofstream CFile(filename,flags);
    float floatBuffer[4];

    if(!CFile)
        return 1;

    if(append)
        CFile.seekp(std::ios::end);

    
    for(unsigned int ui=0; ui<ionVec.size(); ui++)
    {

        for(unsigned int uj=0;uj<4;uj++)
        {
            floatBuffer[uj]=ionVec[ui][uj];
#ifdef __LITTLE_ENDIAN__
            floatSwapBytes(floatBuffer+uj);
#endif
        }
        CFile.write((char *)floatBuffer,4*sizeof(float));


        if(!CFile)
            return 1;
    }
    return 0;
}

unsigned int savePosFile(const vector<Point3D> &ptVec, float mass, const char *filename, bool append)
{
    auto flags = std::ios::binary;
    if(append)
        flags|=std::ios::app;

    std::ofstream CFile(filename,flags);
    float floatBuffer[4];

    if(!CFile)
        return 1;

    if(append)
        CFile.seekp(std::ios::end);
    
#ifdef __LITTLE_ENDIAN__
    floatSwapBytes(&mass);
#endif
    
    for(unsigned int ui=0; ui<ptVec.size(); ui++)
    {

        for(unsigned int uj=0;uj<3;uj++)
        {
            floatBuffer[uj]=ptVec[ui][uj];
#ifdef __LITTLE_ENDIAN__
            floatSwapBytes(floatBuffer+uj);
#endif
        }
        CFile.write((char *)floatBuffer,3*sizeof(float));
        CFile.write((char *)&mass,1*sizeof(float));

        if(!CFile)
            return 1;
    }
    return 0;
}

//Load a TAPSIM Binfile
unsigned int loadTapsimBinFile(vector<IonHit> &posIons, const char *posfile)
{
    ifstream f(posfile,std::ios::binary); 

    if(!f)
        return TAPSIM_OPEN_FAIL;

    f.seekg(0,std::ios::end);
    size_t fileSize;
    fileSize = f.tellg();
    f.seekg(0,std::ios::beg);

    //TAPSIM's binary geometery input format is not totally clearly documented
    // but an example is provided. So best efforts are us.

    std::string str;
    getline(f,str);

    vector<string> s;
    splitStrsRef(str.c_str(),' ',s);

    const unsigned int HEADER_ENTRIES=4;
    if(s.size() != HEADER_ENTRIES)
        return TAPSIM_FILE_FORMAT_FAIL;
    
    bool binaryMode;
    

    if(s[0] == "BINARY")
    {
        //Payload is in binary format
        binaryMode=true;
    }
    else if(s[0] == "ASCII")
    {
        //payload is in ascii format
        binaryMode=false;
    }
    else
        return TAPSIM_FILE_FORMAT_FAIL;

    //maximum size of a TAPSIM record, in bytes
    unsigned int MAX_RECORD_SIZE=18;
    char *buffer=new char[MAX_RECORD_SIZE];
    if(binaryMode)
    {
        size_t numEntries;
        if(stream_cast(numEntries,s[1]))
        {
            delete[] buffer;
            return TAPSIM_FILE_FORMAT_FAIL;
        }
        
        //ID specifies atom type
        bool idsPresent;
        //There is an inconsistency between the example files given and the documentation
        // The documentation says the ids are optional. The examples have the option disabled, 
        // but provide the IDs anyway - this changes the payload size.
        if(s[2] == "1" || s[2] == "0" ) 
            idsPresent=true;
        else
        {
            delete[] buffer;
            return TAPSIM_FILE_FORMAT_FAIL;
        }
        

        size_t numbersPresent;
        if(s[3] == "1" )
            numbersPresent=true;
        else if (s[3] == "0") 
            numbersPresent=false;
        else
        {
            delete[] buffer;
            return TAPSIM_FILE_FORMAT_FAIL;
        }


        //Contrary to docuemntation, ids appear to be "short", not "unsigned int"   
        size_t recordSize=( (numbersPresent ? 1:0 )*sizeof(unsigned int) + 
            3*sizeof(float) + (idsPresent ?1:0)*sizeof(short));

        //Check that payload size matches expected size
        if(fileSize -f.tellg() !=recordSize*numEntries)
        {
            delete[] buffer;
            return TAPSIM_FILE_FORMAT_FAIL;
        }

        IonHit h;
        posIons.resize(numEntries);
        for(size_t ui=0;ui<numEntries; ui++)
        {
            f.read(buffer,recordSize);

            //Transfer position data
#ifdef __LITTLE_ENDIAN__
            h.setPos(Point3D((float*)(buffer)));
#elif __BIG_ENDIAN__
            static_assert(false,"Big endian"); //TODO: IMPLEMENT ME
#endif
            //assign the ID as the mass to charge
            h.setMassToCharge(*( (short*)(buffer+12) ));
        
            posIons[ui]=h;
        }

#pragma omp parallel for
        for(size_t ui=0;ui<numEntries;ui++)
        {
            //rescale back to nm, rather than in m
            posIons[ui].setPos(posIons[ui].getPos()*1e9);
        }
    }
    else
    {
        //TODO: IMPLEMENT ME
        delete[] buffer;
        return TAPSIM_FILE_FORMAT_FAIL;
    }

    delete[] buffer;
    return 0;
}


//write a TAPSIM binary file
//FIXME: Improved error reporting
unsigned int saveTapsimBin( const vector<IonHit> &posIons,std::ostream &f)
{
    std::string data,tmp;
    data="BINARY ";
    stream_cast(tmp,posIons.size());
    data+=tmp + " 0 0\n"; 
    
    f.write(data.c_str(),data.size());  
    

    ASSERT(sizeof(short) == 2);
    for(size_t ui=0;ui<posIons.size();ui++)
    {
        float floats[3];
        short id;

        posIons[ui].getPos().copyValueArr(floats);

        //shrink from "nm" to "m" as natural unit
        for(size_t uj=0;uj<3;uj++)
        {
#ifdef __BIG_ENDIAN__
            floatSwapBytes(floats[uj]);
#endif
            floats[uj]*=1e-9;
        }

        id=posIons[ui].getMassToCharge();
        f.write((const char *)floats,sizeof(float)*3);
        f.write((const char *)&id,sizeof(short));

    }

    return 0;
}

unsigned int saveTapsimBin( const vector<IonHit> &posIons,const char *filename)
{
    ofstream f(filename,std::ios::binary);
    if(!f)
        return 1;

    return saveTapsimBin(posIons,f);
}
bool readEposRecord(const char *src, const char *dest)
{
    static_assert(sizeof(float) == 4,"size of float not 4");
    
    EPOS_ENTRY *entry=(EPOS_ENTRY*)dest;
    float *srcPtr = (float*)src;

    //Perform byte swapping and check nan/inf   
    for(size_t ui=0;ui<9;ui++)
    {
        floatSwapBytes(srcPtr+ui);
        //Disallow NaNs/inf
        if( std::isnan(srcPtr[ui]) || std::isinf(srcPtr[ui]))
            return false;

        
    }

    //Map the buffer locations to the locations within the struct
    // remember that you cannot simply do a memcpy - this won't be portable 
    // due to struct alignment
    entry->x= srcPtr[0];
    entry->y= srcPtr[1];
    entry->z= srcPtr[2];
    entry->massToCharge= srcPtr[3];
    entry->timeOfFlight= srcPtr[4];
    entry->voltDC= srcPtr[5];
    entry->voltPulse= srcPtr[6];
    entry->xDetector= srcPtr[7];
    entry->yDetector= srcPtr[8];

    
    //FIXME: Byte swapping!
    
    int4SwapBytes((int32_t*)(src+(4*9)));
    int4SwapBytes((int32_t*)(src+(4*10)));
    entry->deltaPulse = *((int32_t*)(src+(4*9)));
    entry->hitMultiplicity = *((int32_t*)(src+(4*10)));
    
    
    
    return true;
}


size_t loadEposFile(std::vector<EPOS_ENTRY> &outData,const char *filename)
{
    static_assert(sizeof(EPOS_ENTRY) >= EPOS_RECORD_SIZE, "EPOS_RECORD_SIZE");

    size_t errcode;
    errcode=fixedRecordReader(filename,readEposRecord, EPOS_RECORD_SIZE,outData);
    return errcode;
}

size_t chunkLoadEposFile(std::vector<EPOS_ENTRY> &outData, const char *filename, 
    unsigned int nChunksToRead, unsigned int chunkOffset, unsigned int &nChunksLeft)
{
    static_assert(sizeof(EPOS_ENTRY) >= EPOS_RECORD_SIZE, "EPOS_RECORD_SIZE");

    size_t errcode;
    errcode=fixedRecordChunkReader(filename,readEposRecord, EPOS_RECORD_SIZE,
        outData,nChunksToRead,chunkOffset,nChunksLeft);
    return errcode;
}

const char *OPS_ENUM_ERRSTRINGS[] =
{
    "",
        "\"C\" line (exp. parameters) not in expected format",
    "\"I\"  line (Detector parameters) not in expected format",
    "\"-\" line (Time-of-flights) not in expected format",
    "\"V\" line (voltage data) not in expected format",
    "Error interpreting \"V\" line (voltage data) data",
    "Missing beta value from \"V\" line - needs to be in voltage, or in system header (\"C\" line",
    "Incorrect number of \"V\" line (voltage data) data entries",
    "Unknown linetype in data file",
    "\"P\" line (detector channels) not in expected format",
    "Unable to interpret channels line",
    "Incorrect number of events in \"S\" (hit position) line",
    "Unable to interpret event data on \"S\" (hit position) line",
    "Unable to parse \"S\" (hit position) line event count (eg 8 in S8 ...)",
    "\"S\" (hit position) line data not preceded by TOF data, as it should have been",
    "Duplicate system data/experiment setup (\"C\") entry found -- there can only be one",
    "Duplicate detector (\"I\") data entry found -- there can only be one",
    "Trailing \"-\" line found -- should have be followed by a \"S\" line, but wasn't.",
    "Duplicate\"-\" line found -- should have be followed by a \"S\" line, but wasn't.",
    "Unable to open file",
    "unable to read file, after opening",
    "Abort requested"
};

//Function to read POSAP "OPS" files, which contain data
//on voltage, positioned hits and time of flights
//as well as some instrument parameters
unsigned int readPosapOps(const char *file, 
        THREEDAP_EXPERIMENT &data, 
            unsigned int &badLine, unsigned int &progress, 
            std::atomic<bool> &wantAbort, unsigned int nDelayLines,bool strictMode)
{

    static_assert(ARRAYSIZE(OPS_ENUM_ERRSTRINGS) == OPSREADER_ENUM_END, "OPS error strings and error code enum sizes mismatched");

    const char COMMENT_MARKER='*';

    badLine=0;

    size_t totalFilesize;
    if(!getFilesize(file,totalFilesize))
        return OPSREADER_OPEN_ERR;

    std::ifstream f(file);

    if(!f)
        return OPSREADER_OPEN_ERR;


    enum
    {
        OPS_LINE_DASH,
        OPS_LINE_OTHER_OR_NONE
    };

    bool haveSystemParams=false;
    bool haveDetectorSize=false;

    unsigned int lastLine=OPS_LINE_OTHER_OR_NONE;   

    vector<string> strVec;
    while(!f.eof())
    {
        std::string strLine;
        getline(f,strLine);
        badLine++;


        if(f.eof())
            break;

        if(!f.good())
            return OPSREADER_READ_ERR;

        //If user wants to abort, please do so
        if(wantAbort)
            return OPSREADER_ABORT_ERR;

        //Remove any ws delimiters
        strLine=stripWhite(strLine);

        //Is a blank or comment? Don't process line
        if(!strLine.size() || strLine[0] == COMMENT_MARKER)
            continue;

        //Check to see what type of line this is
        strVec.clear();

        splitStrsRef(strLine.c_str(),"\t ",strVec);
        stripZeroEntries(strVec);
        switch(strLine[0])
        {
            //"C" line contains the system parameters
            // flight path, pulse coupling coefficients (alpha, beta -- see Sebastien et al) and time zero offset.
            // "New method for the calibration of three-dimensional atom-probe mass"
            // Review of Scientific Instruments, 2001
            case 'C':
            {
                if(haveSystemParams)
                    return OPSREADER_FORMAT_DUPLICATE_SYSDATA;

                haveSystemParams = true;

                //Should have something eg
                // C 123 1.01 0.8 -123 
                if(strVec.size() != 5)
                    return OPSREADER_FORMAT_CLINE_ERR;

                if(stream_cast( data.params.flightPath,strVec[1]))
                    return OPSREADER_FORMAT_CLINE_ERR;

                if(stream_cast( data.params.alpha,strVec[2]))
                    return OPSREADER_FORMAT_CLINE_ERR;
                
                if(stream_cast( data.params.beta,strVec[3]))
                    return OPSREADER_FORMAT_CLINE_ERR;

                if(stream_cast( data.params.tZero,strVec[4]))
                    return OPSREADER_FORMAT_CLINE_ERR;

                lastLine=OPS_LINE_OTHER_OR_NONE;
                break;
            }
            //"I" Line contains detector radius
            case 'I':
            {
                if(haveDetectorSize)
                    return OPSREADER_FORMAT_DUPLICATE_DETECTORSIZE;

                haveDetectorSize=true;

                //Should have I or IR
                //I - No reflectron
                //IR - reflectron
                //and a floating pt number
                if(strVec.size() != 2)
                    return OPSREADER_FORMAT_DETECTORLINE_ERR;

                if(stream_cast(data.params.detectorRadius,strVec[1]))
                    return OPSREADER_FORMAT_DETECTORLINE_ERR;
                lastLine=OPS_LINE_OTHER_OR_NONE;

                break;
            }
            //"S" line contains the hit TOFs, and hit frequency
            case 'S':
            {
                //OK, so this is the actual data -- each "S"" line should be
                //preceded by a "-" line.
                if(lastLine !=OPS_LINE_DASH)
                {
                    if(strictMode)
                        return OPSREADER_FORMAT_SLINE_PREFIX_ERR;
                    else
                    {
                        continue; //OK, this is kinda a confusing situation, lets just get to the next line, which we hope is an S.
                    }
                }



                //If we have less positions than TOFs, then something is wrong 
                //-- we have to discard the previous "-" line
                // could have a wierd line that makes no sense (for eg, I found this):
                //
                //-151 13
                //S0
                // I assume that this means that 
                if((strVec.size()-1) /(nDelayLines+1)< data.eventData.back().size())
                {
                    ASSERT(data.eventData.size())
                    //We have to discard the previous hit sequence
                    data.eventData.pop_back();
                    lastLine=OPS_LINE_OTHER_OR_NONE;
                    break;
                }


                //OK, so the event numbers are coded like
                //"S8", where S indicates it is an event
                //line, and the next numbers indicate the
                //number of observed events
                strVec[0]= strVec[0].substr(1);

                unsigned int numEvents;
                if(stream_cast(numEvents,strVec[0]))
                {
                    if(strictMode)
                        return OPSREADER_FORMAT_SLINE_FORMAT_ERR;
                    else
                    {
                        data.eventData.pop_back();
                        lastLine=OPS_LINE_OTHER_OR_NONE;
                        break;
                    }
                }


                if(numEvents != (strVec.size()-1)/(nDelayLines+1))
                    return OPSREADER_FORMAT_SLINE_EVENTCOUNT_ERR;

                
                vector<SINGLE_HIT> &curHits=data.eventData.back();
            
                //The TOF must be preset by a "-" line, as checked above.
                //The contents of this "S" line must have at least enough elements
                //to cover the TOFs. The other hits are assigned a zero TOF
                if(curHits.size()> numEvents)
                {
                    if(strictMode)
                        return OPSREADER_FORMAT_SLINE_EVENTCOUNT_ERR;
                    else
                    {
                        data.eventData.pop_back();
                        lastLine=OPS_LINE_OTHER_OR_NONE;
                        break;
                    }
                }


                SINGLE_HIT h;
                h.tof=0;
                curHits.resize(numEvents,h);
                vector<float> tofs;

                //Retrieve the TOF values, so we can redirect them
                //as the TOF values above's correspondence map
                //is stored at the end of the hit sequence, one 
                //entry per hit
                //zero means untimed, otherwise TOF
                tofs.resize(curHits.size());
                for(unsigned int ui=0;ui<curHits.size();ui++)
                    tofs[ui]=curHits[ui].tof;


                //Read up until the TOF timing map
                bool everythingOK=true;
                unsigned int timingIndex=curHits.size()*2+1;
                for(unsigned int ui=1;ui<timingIndex;ui+=2)
                {
                    //Push the current hit into a vector 
                    unsigned int offset;
                    offset=ui/2;

                    if(stream_cast(curHits[offset].x,strVec[ui]))
                    {
                        if(strictMode)
                            return OPSREADER_FORMAT_SLINE_EVENTDATA_ERR;
                        else
                        {
                            everythingOK=false;
                            break;
                        }

                    }
                    
                    if(stream_cast(curHits[offset].y,strVec[ui+1]))
                    {
                        if(strictMode)
                            return OPSREADER_FORMAT_SLINE_EVENTDATA_ERR;
                        else
                        {
                            everythingOK=false;
                            break;
                        }
                    }
                }

                if(!everythingOK)
                {
                    data.eventData.pop_back();
                    lastLine=OPS_LINE_OTHER_OR_NONE;
                    break;
                }   

                vector<unsigned int> timeMap;
                //Read the TOF timing map
                for(unsigned int ui=timingIndex; ui<strVec.size();ui++)
                {
                    unsigned int mapEntry;
                    if(stream_cast(mapEntry,strVec[ui]))
                    {
                        if(strictMode)
                            return OPSREADER_FORMAT_SLINE_EVENTDATA_ERR;
                        else
                        {
                            everythingOK=false;
                            break;

                        }
                    }

                    timeMap.push_back(mapEntry);
                }

                if(!everythingOK)
                {
                    data.eventData.pop_back();
                    lastLine=OPS_LINE_OTHER_OR_NONE;
                    break;
                }   
                ASSERT(timeMap.size() == curHits.size());

                //use the TOF timing map to assign TOF values
                //to observed TOF entries from the "-" line
                for(unsigned int ui=0;ui<curHits.size();ui++)
                {
                    if(timeMap[ui])
                        curHits[ui].tof=tofs[timeMap[ui]-1];
                    else
                        curHits[ui].tof=0;
                }
                
                lastLine=OPS_LINE_OTHER_OR_NONE;
                break;
            }
            //"-" line contains the detector hit positions
            case '-':
            {
                if(lastLine==OPS_LINE_DASH)
                    return OPSREADER_FORMAT_DOUBLEDASH;

                if(strVec.size() < 2)
                    return OPSREADER_FORMAT_LINE_DASH_ERR;


                size_t pulseDelta;
                if(stream_cast(pulseDelta,strVec[0].substr(1)))
                    return OPSREADER_FORMAT_LINE_DASH_ERR;

                //First bit of string is the number of pulses
                //preceding this event
                if(data.eventPulseNumber.empty())
                    data.eventPulseNumber.push_back(pulseDelta+1);
                else
                {
                    data.eventPulseNumber.push_back(
                        data.eventPulseNumber.back()+pulseDelta+1);
                }

                vector<float> tofs;
                for(unsigned int ui=1; ui< strVec.size(); ui++)
                {
                    float timeOfFlight;
                    if(stream_cast(timeOfFlight,strVec[ui]))
                        return OPSREADER_FORMAT_LINE_DASH_ERR;

                    tofs.push_back(timeOfFlight);
                }

                vector<SINGLE_HIT> s;

                s.resize(tofs.size());
                for(unsigned int ui=0;ui<tofs.size();ui++)
                    s[ui].tof =tofs[ui]; 
                data.eventData.push_back(s);

                lastLine=OPS_LINE_DASH;
                break;
            }
            //\"V\" voltage data
            case 'V':
            {

                //V line may have
                if(! (strVec.size() == 4 || strVec.size() == 3) )
                    return OPSREADER_FORMAT_LINE_VOLTAGE_DATACOUNT_ERR;

                VOLTAGE_DATA vDat;
                if(stream_cast(vDat.voltage,strVec[1]))
                    return OPSREADER_FORMAT_LINE_VOLTAGE_DATA_ERR;

                float pulseV;
                if(stream_cast(pulseV,strVec[2]))
                    return OPSREADER_FORMAT_LINE_VOLTAGE_DATA_ERR;
                vDat.pulseVolt=pulseV;

                if(strVec.size() ==4)
                {
                    if(stream_cast(vDat.beta,strVec[3]))
                        return OPSREADER_FORMAT_LINE_VOLTAGE_DATA_ERR;
                }
                else // 3 items
                {
                    //If we have not seen the system parameter's beta, then we cannot compute a good effective voltage
                    if(!haveSystemParams)
                        return OPSREADER_FORMAT_LINE_VOLTAGE_NOBETA;
                    vDat.beta = data.params.beta;
                }

                vDat.nextHitGroupOffset=data.eventData.size();
                data.voltageData.push_back(vDat);
                lastLine=OPS_LINE_OTHER_OR_NONE;
                break;
            }
            case 'P':
            {
                if(strVec.size() != 2)
                    return OPSREADER_FORMAT_CHANNELS_ERR;
                
                if(stream_cast(data.params.detectorChannels,strVec[1]))
                    return OPSREADER_CHANNELS_DATA_ERR;
                break;
            }
            case 'T':
            {
                //This is a timing,  temperature and pressure line.
                // ignore.
                //FIXME: We should read these and report them
                break;
            }
            default:
                return OPSREADER_FORMAT_LINETYPE_ERR;
        }

        progress=(float)f.tellg()/(float)totalFilesize*100.0f;
    }

    if(lastLine==OPS_LINE_DASH)
    {
        //Damn, we have a problem. We saw a "-" line,
        //but failed to find the next "S" line, which should have been there.
        //Something is wrong.
        if(strictMode)
            return OPSREADER_FORMAT_TRAILING_DASH_ERR;
        
        data.eventData.pop_back();
        
    }

    return 0;
}


unsigned int loadATOFile(const char *fileName, vector<ATO_ENTRY> &ions, unsigned int forceEndian)
{
    //open pos file
    std::ifstream CFile(fileName,std::ios::binary);

    if(!CFile)
        return ATO_OPEN_FAIL;

    //Get the filesize
    CFile.seekg(0,std::ios::end);
    size_t fileSize=CFile.tellg();


    //There are differences in the format, unfortunately.
    // Gault et al, Atom Probe Microscopy says 
    // - there are 14 entries of 4 bytes, 
    //   totalling "44" bytes - which cannot be correct. They however,
    //   say that the XYZ is added later.
    // - Header is 2 32 binary
    // - File is serialised as little-endian
    // - Various incompatible versions exist. Unclear how to distinguish. 

    // Larson et al say that 
    // - there are 14 fields
    // - Header byte 0x05 (0-indexed) is version number, and only version 3 is outlined
    // - Pulsenumber can be subject to FP aliasing (bad storage, occurs for values > ~16.7M ), 
    //  - Aliasing errors must be handled, if reading this field
    // - File is in big-endian


    //In summary, we assume there are 14 entries, 4 bytes each, after an 8 byte header.
    // we assume that the endian-ness must be auto-detected somehow, as no sources
    // agree on file endian-ness. If we cannot detect it, we assume little endian

    //Known versions of file : 2, 5

    //Header (8 bytes), record 14 entries, 4 bytes each

    const size_t ATO_HEADER_SIZE=8;
    const size_t ATO_RECORD_SIZE = 14*4;
    const size_t ATO_MIN_FILESIZE = 8 + ATO_RECORD_SIZE;

    if(fileSize < ATO_MIN_FILESIZE)
        return ATO_EMPTY_FAIL;
    
    
    //calculate the number of points stored in the POS file
    IonHit hit;
    size_t pointCount=0;
    if((fileSize - ATO_HEADER_SIZE)  % (ATO_RECORD_SIZE))
        return ATO_SIZE_ERR;    


    //Check that the version number, stored at offxet 0x05 (1-indexed), is 3.
    CFile.seekg(4);
    unsigned int versionByte;
    CFile.read((char*)&versionByte,sizeof(unsigned int));

    //Assume that we can have a new version that doesn't affect the readout
    // assume that earlier versions are compatible. This means, for a random byte
    // in a random length (modulo) file, 
    // we have a 1-5/255 chance of rejection from this test, and a 1/56 chance of
    // rejection from filesize, giving a ~0.02% chance of incorrect acceptance.
    // Known versions : 2, 3, 5
    if(!versionByte || versionByte > 5)
        return ATO_VERSIONCHECK_ERR;


    pointCount = (fileSize-ATO_HEADER_SIZE)/ATO_RECORD_SIZE;
    
    try
    {
        ions.resize(pointCount);
    }
    catch(std::bad_alloc)
    {
        return ATO_MEM_ERR;
    }


    //Heuristic to detect endianness.
    // - Randomly sample 100 pts from file, and check to see if, when interpreted either way,
    // there are any NaN
    //   


    bool endianFlip;
    
    if(forceEndian)
    {
        ASSERT(forceEndian < 3);
#ifdef __LITTLE_ENDIAN__
        endianFlip=(forceEndian == 2);
#elif defined(__BIG_ENDIAN__)
        endianFlip=(forceEndian == 1);
#endif
    }
    else
    {
        //Auto-detect endianness from file content
        size_t numToCheck=std::min(pointCount,(size_t)100);
    
        //Indicies of points to check
        vector<unsigned int> randomNumbers;
        gsl_rng *rng = randGen.getRng();
        randomIndices(randomNumbers,pointCount,
                    numToCheck,rng);

        //Make the traverse in ascending order
        std::sort(randomNumbers.begin(),randomNumbers.end());

        //One for no endian-flip, one for flip
        bool badFloat[2]={ false,false };
        //Track the presence of unreasonably large numbers
        bool veryLargeNumber[2] = { false,false };
    
        //Skip through several records, looking for bad float data,
        auto buffer = new float[ATO_RECORD_SIZE/4];
        for(size_t ui=0;ui<numToCheck;ui++)
        {
            size_t offset;
            offset=randomNumbers[ui];
            
            CFile.seekg(ATO_HEADER_SIZE + ATO_RECORD_SIZE*offset);
            CFile.read((char*)buffer,ATO_RECORD_SIZE);

            const unsigned int BYTES_TO_CHECK[] = { 0,1,2,3,5,6,8,9,10 };
            const size_t CHECKBYTES = 9;

            //Check each field for inf/nan presence
            for(size_t uj=0;uj<CHECKBYTES;uj++)
            {
                if(std::isnan(buffer[BYTES_TO_CHECK[uj]]) ||
                    std::isinf(buffer[BYTES_TO_CHECK[uj]]))
                    badFloat[0]=true;

                //Flip the endian-ness
                floatSwapBytes(buffer+BYTES_TO_CHECK[uj]);

                if(std::isnan(buffer[BYTES_TO_CHECK[uj]]) ||
                    std::isinf(buffer[BYTES_TO_CHECK[uj]]))
                    badFloat[1]=true;
                
                //Swap it back
                floatSwapBytes(buffer+BYTES_TO_CHECK[uj]);


            }
        
            
            //Check for some very likely values.
            //Check for large negative masses
            if( buffer[3] < -1000.0f)
                veryLargeNumber[0] = true;

            // unlikely to exceed 1000 kV
            if( fabs(buffer[6]) > 1000.0f || fabs(buffer[10]) > 1000.0f)
                veryLargeNumber[0] = true;

            //Swap and try again
            floatSwapBytes(buffer+3);
            floatSwapBytes(buffer+6);
            floatSwapBytes(buffer+10);
            if( buffer[3] < -1000.0f) 
                veryLargeNumber[1] = true;

            if( fabs(buffer[6]) > 1000.0f || fabs(buffer[10]) > 1000.0f)
                veryLargeNumber[1] = true;
        }

        delete[] buffer;


        //Now summarise the results

        //If we have a disagreement about bad-float-ness,
        // or stupid-number ness.  then choose the good one.
        // Otherwise abandon detection
        if(badFloat[0] != badFloat[1])
        {
            endianFlip=(badFloat[0]);
        }
        else if(veryLargeNumber[0] != veryLargeNumber[1])
        {
            endianFlip=veryLargeNumber[0];
        }
        else
        {
            //Assume little endian
#ifdef __LITTLE_ENDIAN__
            endianFlip= false;
#else
            endianFlip=true;
#endif
        }
    }

    //File records consist of 14 fields, some of which may not be initialised.
    // each being 4-byte IEEE little-endian float
    // It is unknown how to detect initialised state.
    // Field    Data    
    // 0-3      x,y,z,m/c in Angstrom (x,yz) or Da (m/c)
    // 4        clusterID, if set
    //          - Ignore this field, as this information is redundant
    // 5        Approximate Pulse #, due to Float. Pt. limitation 
    // 6        Standing Voltage (kV)
    // 7        TOF (us) (maybe corrected? maybe not?)
    // 8-9      Detector position (cm)
    // 10       Pulse voltage (kV)
    // 11       "Virtual voltage" for reconstruction.
    //          - Ignore this field, as this information is redundant
    // 12,13    Fourier intensity
    //          - Ignore these fields, as this information is redundant
    //Attempt to detect
    CFile.seekg(8);

    auto buffer = new float[ATO_RECORD_SIZE/4];
    size_t curPos=0;    

    if(endianFlip)
    {
        //Read and swap
        while((size_t)CFile.tellg() < fileSize)
        {
            CFile.read((char*)buffer,ATO_RECORD_SIZE);

            for(size_t ui=0;ui<ATO_RECORD_SIZE;ui++)
                floatSwapBytes(buffer+ui);

            float *fBuf=(float*)buffer;
            ions[curPos].x = fBuf[0];
            ions[curPos].y = fBuf[1];
            ions[curPos].z = fBuf[2];
            ions[curPos].mass= fBuf[3];
            ions[curPos].approxPulse= fBuf[5];
            ions[curPos].voltage = fBuf[6];
            ions[curPos].tof= fBuf[7];
            ions[curPos].detectorX = fBuf[8];
            ions[curPos].detectorY = fBuf[9];
            ions[curPos].pulseVoltage= fBuf[10];
            curPos++;

        }
    }
    else
    {
        //read without swapping
        while((size_t)CFile.tellg() < fileSize)
        {   
            CFile.read((char*)buffer,ATO_RECORD_SIZE);

            float *fBuf=(float*)buffer;
            ions[curPos].x = fBuf[0];
            ions[curPos].y = fBuf[1];
            ions[curPos].z = fBuf[2];
            ions[curPos].mass= fBuf[3];
            ions[curPos].approxPulse= fBuf[5];
            ions[curPos].voltage = fBuf[6];
            ions[curPos].tof= fBuf[7];
            ions[curPos].detectorX = fBuf[8];
            ions[curPos].detectorY = fBuf[9];
            ions[curPos].pulseVoltage= fBuf[10];
            curPos++;
            
        }
    }

    delete[] buffer;


    return 0;
}

bool strhas(const char *cpTest, const char *cpPossible)
{
    while(*cpTest)
    {
        const char *search;
        search=cpPossible;
        while(*search)
        {
            if(*search == *cpTest)
                return true;
            search++;
        }
        cpTest++;
    }

    return false;
}

//A routine for loading numeric data from a text file
unsigned int loadTextData(const char *cpFilename, vector<vector<float> > &dataVec,vector<string> &headerVec, 
        const char *delim, bool allowNan, bool allowConvFails, unsigned int headerCount)
{
    const unsigned int BUFFER_SIZE=8192;
    char inBuffer[BUFFER_SIZE];
    
    unsigned int num_fields=0;


#if !defined(WIN32) && !defined(_WIN64)
    if(isNotDirectory(cpFilename) == false)
        return TEXT_ERR_FILE_OPEN;
#endif

    dataVec.clear();
    //Open a file in text mode
    std::ifstream CFile(cpFilename);

    if(!CFile)
        return TEXT_ERR_FILE_OPEN;

    //Drop the headers, if any
    string str;
    vector<string> strVec;  
    bool atHeader=true;

    vector<string> prevStrs;
    while(CFile.good() && !CFile.eof() && atHeader)
    {
        //Grab a line from the file
        if(!CFile.getline(inBuffer,BUFFER_SIZE))
            break;

        if(!CFile.good())
            return TEXT_ERR_FILE_FORMAT;

        prevStrs=strVec;
        //Split the strings around the deliminator c
        splitStrsRef(inBuffer,delim,strVec);        
        stripZeroEntries(strVec);
        
    
        //Skip blank lines or lines that are only spaces
        if(strVec.empty())
            continue;

        num_fields = strVec.size();
        dataVec.resize(num_fields);     
        //Check to see if we are in the header
        if(atHeader)
        {
            //Count down headers if we still 
            if(headerCount  !=-1)
            {
                headerCount--;
                if(!headerCount)
                {
                    atHeader=false;
                    prevStrs=strVec;
                    continue;
                }
            }

            //If we have the right number of fields 
            //we might be not in the header anymore
            if(num_fields >= 1 && strVec[0].size())
            {
                float f;
                //Assume we are not reading the header
                atHeader=false;

                vector<float> values;
                //Confirm by checking all values
                for(unsigned int ui=0; ui<num_fields; ui++)
                {
    
                    //stream_cast will fail in the case "1 2" if there
                    //is a tab delimiter specified. Check for only 02,3 
                    if(strVec[ui].find_first_not_of("0123456789.Ee+-Na") 
                                        != string::npos)
                    {
                        atHeader=true;
                        break;
                    }
                    //If any cast fails
                    //we are in the header
                    if(stream_cast(f,strVec[ui]))
                    {
                        if(!allowNan || strVec[ui] != "NaN") 
                        {
                            atHeader=true;
                            break;
                        }
                    }

                    values.push_back(f);
                }

                if(!atHeader)
                    break;
            }
        }

    }

    //Drop the blank bits from the field
    stripZeroEntries(prevStrs);
    //switch this out as being the header
    if(prevStrs.size() == num_fields)
        std::swap(headerVec,prevStrs);

    if(atHeader)
    {
        //re-wind back to the beginning of the file
        //as we didn't find a header/data split
        CFile.clear(); // clear EOF bit
        CFile.seekg(0,std::ios::beg);
        CFile.getline(inBuffer,BUFFER_SIZE);
        
        
        splitStrsRef(inBuffer,delim,strVec);    
        stripZeroEntries(strVec);
        num_fields=strVec.size();

    }

    float f;
    while(CFile.good() && !CFile.eof())
    {
        if(strhas(inBuffer,"0123456789"))
        {
            //Split the strings around the tab char
            splitStrsRef(inBuffer,delim,strVec);    
            stripZeroEntries(strVec);
            
            //Check the number of fields    
            //=========
            if(strVec.size() != num_fields)
                return TEXT_ERR_FILE_NUM_FIELDS;    
    

            for(unsigned int ui=0; ui<num_fields; ui++)
            {   
                if(strVec[ui] == "NaN" && allowNan)
                    f = NAN;
                else
                {
                    std::stringstream ss;
                    ss.clear();
                    ss.str(strVec[ui]);
                    ss >> f;
                    if(ss.fail())
                    {
                        if(!allowConvFails)
                            return TEXT_ERR_FILE_FORMAT;
                        f=NAN;
                    }
                }
                dataVec[ui].push_back(f);
            
            }
            //=========
            
        }
        //Grab a line from the file
        if(!CFile.getline(inBuffer,BUFFER_SIZE))
            break;
    }

    return 0;
}



}