Index: trunk/Mars/mcore/fits.h
===================================================================
--- trunk/Mars/mcore/fits.h	(revision 11806)
+++ trunk/Mars/mcore/fits.h	(revision 11806)
@@ -0,0 +1,657 @@
+#ifndef MARS_fits
+#define MARS_fits
+
+#include <stdint.h>
+
+#include <map>
+#include <string>
+#include <sstream>
+#include <algorithm>
+
+#ifdef __EXCEPTIONS
+#include <stdexcept>
+#endif
+
+#ifdef __CINT__
+#define off_t size_t
+#endif
+
+#ifndef __MARS__
+#include <vector>
+#include <iomanip>
+#include <iostream>
+#define gLog cerr
+#define err ""
+#define all ""
+#else
+#include "MLog.h"
+#include "MLogManip.h"
+#endif
+
+#ifdef HAVE_ZLIB
+#include "izstream.h"
+#else
+#include <fstream>
+#define izstream ifstream
+#warning Support for zipped FITS files disabled.
+#endif
+
+namespace std
+{
+
+class fits : public izstream
+{
+public:
+    struct Entry
+    {
+        char   type;
+        string value;
+        string comment;
+
+        template<typename T>
+            T Get() const
+        {
+            T t;
+
+            istringstream str(value);
+            str >> t;
+
+            return t;
+        }
+    };
+
+    struct Table
+    {
+        off_t offset;
+
+        string name;
+        size_t bytes_per_row;
+        size_t num_rows;
+        size_t num_cols;
+
+        struct Column
+        {
+            size_t offset;
+            size_t num;
+            size_t size;
+            char   type;
+        };
+
+        typedef map<string, Entry>  Keys;
+        typedef map<string, Column> Columns;
+
+        Columns cols;
+        Keys    keys;
+
+        string Trim(const string &str, char c=' ') const
+        {
+            // Trim Both leading and trailing spaces
+            const size_t pstart = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces
+            const size_t pend   = str.find_last_not_of(c);  // Find the first character position from reverse af
+
+            // if all spaces or empty return an empty string
+            if (string::npos==pstart || string::npos==pend)
+                return string();
+
+            return str.substr(pstart, pend-pstart+1);
+        }
+
+        bool Check(const string &key, char type, const string &value="") const
+        {
+            const Keys::const_iterator it = keys.find(key);
+            if (it==keys.end())
+            {
+                ostringstream str;
+                str << "Key '" << key << "' not found.";
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << err << "ERROR - " << str.str() << endl;
+                return false;
+#endif
+            }
+
+            if (it->second.type!=type)
+            {
+                ostringstream str;
+                str << "Wrong type for key '" << key << "': expected " << type << ", found " << it->second.type << ".";
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << err << "ERROR - " << str.str() << endl;
+                return false;
+#endif
+            }
+
+            if (!value.empty() && it->second.value!=value)
+            {
+                ostringstream str;
+                str << "Wrong value for key '" << key << "': expected " << value << ", found " << it->second.value << ".";
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << err << "ERROR - " << str.str() << endl;
+                return false;
+#endif
+            }
+
+            return true;
+        }
+
+        Keys ParseBlock(const vector<string> &vec) const
+        {
+            map<string,Entry> rc;
+
+            for (unsigned int i=0; i<vec.size(); i++)
+            {
+                const string key = Trim(vec[i].substr(0,8));
+                // Keywords without a value, like COMMENT / HISTORY
+                if (vec[i].substr(8,2)!="= ")
+                    continue;
+
+                char type = 0;
+
+                string com;
+                string val = Trim(vec[i].substr(10));
+                if (val[0]=='\'')
+                {
+                    // First skip all '' in the string
+                    size_t p = 1;
+                    while (1)
+                    {
+                        const size_t pp = val.find_first_of('\'', p);
+                        if (pp==string::npos)
+                            break;
+
+                        p = val[pp+1]=='\'' ? pp+2 : pp+1;
+                    }
+
+                    // Now find the comment
+                    const size_t ppp = val.find_first_of('/', p);
+
+                    // Set value, comment and type
+                    com  = ppp==string::npos ? "" : Trim(val.substr(ppp+1));
+                    val  = Trim(val.substr(1, p-2));
+                    type = 'T';
+                }
+                else
+                {
+                    const size_t p = val.find_first_of('/');
+
+                    com = Trim(val.substr(p+2));
+                    val = Trim(val.substr(0, p));
+
+                    if (val.empty() || val.find_first_of('T')!=string::npos || val.find_first_of('F')!=string::npos)
+                        type = 'B';
+                    else
+                        type = val.find_last_of('.')==string::npos ? 'I' : 'F';
+                }
+
+                const Entry e = { type, val, com };
+                rc[key] = e;
+            }
+
+            return rc;
+        }
+
+        Table() : offset(0) { }
+        Table(const vector<string> &vec, off_t off) :
+            offset(off), keys(ParseBlock(vec))
+        {
+            if (!Check("XTENSION", 'T', "BINTABLE") ||
+                !Check("NAXIS",    'I', "2")        ||
+                !Check("BITPIX",   'I', "8")        ||
+                !Check("PCOUNT",   'I', "0")        ||
+                !Check("GCOUNT",   'I', "1")        ||
+                !Check("EXTNAME",  'T')             ||
+                !Check("NAXIS1",   'I')             ||
+                !Check("NAXIS2",   'I')             ||
+                !Check("TFIELDS",  'I'))
+                return;
+
+            bytes_per_row = Get<size_t>("NAXIS1");
+            num_rows      = Get<size_t>("NAXIS2");
+            num_cols      = Get<size_t>("TFIELDS");
+
+            size_t bytes = 0;
+            for (size_t i=1; i<=num_cols; i++)
+            {
+                ostringstream num;
+                num << i;
+
+                if (!Check("TTYPE"+num.str(), 'T') ||
+                    !Check("TFORM"+num.str(), 'T'))
+                    return;
+
+                const string id  = Get<string>("TTYPE"+num.str());
+                const string fmt = Get<string>("TFORM"+num.str());
+
+                istringstream sin(fmt);
+                int n = 0;
+                sin >> n;
+                if (!sin)
+                    n = 1;
+
+                const char type = fmt[fmt.length()-1];
+
+                size_t size = 0;
+                switch (type)
+                {
+                    // We could use negative values to mark floats
+                    // otheriwse we could just cast them to int64_t?
+                case 'L': size = 1; break; // logical
+                // case 'X': size = n; break; // bits (n=number of bytes needed to contain all bits)
+                case 'B': size = 1; break; // byte
+                case 'I': size = 2; break; // short
+                case 'J': size = 4; break; // int
+                case 'K': size = 8; break; // long long
+                case 'E': size = 4; break; // float
+                case 'D': size = 8; break; // double
+                // case 'C': size =  8; break; // complex float
+                // case 'M': size = 16; break; // complex double
+                // case 'P': size =  8; break; // array descriptor (32bit)
+                // case 'Q': size = 16; break; // array descriptor (64bit)
+                default:
+                    {
+                        ostringstream str;
+                        str << "FITS format TFORM='" << fmt << "' not yet supported.";
+#ifdef __EXCEPTIONS
+                        throw runtime_error(str.str());
+#else
+                        gLog << err << "ERROR - " << str.str() << endl;
+                        return;
+#endif
+                    }
+                }
+
+                const Table::Column col = { bytes, n, size, type };
+
+                cols[id] = col;
+                bytes  += n*size;
+            }
+
+            if (bytes!=bytes_per_row)
+            {
+#ifdef __EXCEPTIONS
+                throw runtime_error("Column size mismatch");
+#else
+                gLog << err << "ERROR - Column size mismatch" << endl;
+                return;
+#endif
+            }
+
+            name = Get<string>("EXTNAME");
+        }
+
+        void PrintKeys() const
+        {
+            for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++)
+                gLog << all << setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << endl;
+        }
+
+        void PrintColumns() const
+        {
+            typedef map<pair<size_t, string>, Column> Sorted;
+
+            Sorted sorted;
+
+            for (Columns::const_iterator it=cols.begin(); it!=cols.end(); it++)
+                sorted[make_pair(it->second.offset, it->first)] = it->second;
+
+            for (Sorted::const_iterator it=sorted.begin(); it!=sorted.end(); it++)
+                gLog << all << setw(6) << it->second.offset << "| " << it->first.second << " [" << it->second.type << '=' << it->second.num << 'x' << it->second.size << "byte]" << endl;
+        }
+
+        operator bool() const { return !name.empty(); }
+
+        bool HasKey(const string &key) const
+        {
+            return keys.find(key)!=keys.end();
+        }
+
+            // Values of keys are always signed
+        template<typename T>
+            T Get(const string &key) const
+        {
+            const map<string,Entry>::const_iterator it = keys.find(key);
+            if (it==keys.end())
+            {
+                ostringstream str;
+                str << "Key '" << key << "' not found." << endl;
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << err << "ERROR - " << str.str() << endl;
+                return 0;
+#endif
+            }
+            return it->second.Get<T>();
+        }
+
+        size_t GetN(const string &key) const
+        {
+            const Columns::const_iterator it = cols.find(key);
+            return it==cols.end() ? 0 : it->second.num;
+        }
+    };
+
+private:
+    Table fTable;
+
+    typedef pair<void*, Table::Column> Address;
+    typedef vector<Address> Addresses;
+    //map<void*, Table::Column> fAddresses;
+    Addresses fAddresses;
+
+    vector<char> fBufferRow;
+    vector<char> fBufferDat;
+
+    size_t fRow;
+
+    vector<string> ReadBlock(vector<string> &vec)
+    {
+        bool endtag = false;
+        for (int i=0; i<36; i++)
+        {
+            char c[81];
+            c[80] = 0;
+            read(c, 80);
+            if (!good())
+                break;
+
+            if (c[0]==0)
+                return vector<string>();
+
+            string str(c);
+
+//            if (!str.empty())
+//                cout << setw(2) << i << "|" << str << "|" << (endtag?'-':'+') << endl;
+
+            if (str=="END                                                                             ")
+            {
+                endtag = true;
+
+                // Make sure that no empty vector is returned
+                if (vec.size()%36==0)
+                    vec.push_back(string("END     = '' / "));
+            }
+
+            if (endtag)
+                continue;
+
+            vec.push_back(str);
+        }
+
+        return vec;
+    }
+
+    string Compile(const string &key, int16_t i=-1)
+    {
+        if (i<0)
+            return key;
+
+        ostringstream str;
+        str << key << i;
+        return str.str();
+    }
+
+public:
+    fits(const string &fname) : izstream(fname.c_str())
+    {
+        char simple[10];
+        read(simple, 10);
+        if (!good())
+            return;
+
+        if (memcmp(simple, "SIMPLE  = ", 10))
+        {
+            clear(rdstate()|ios::badbit);
+#ifdef __EXCEPTIONS
+            throw runtime_error("File is not a FITS file.");
+#else
+            gLog << err << "ERROR - File is not a FITS file." << endl;
+            return;
+#endif
+        }
+
+        seekg(0);
+
+        while (good())
+        {
+            vector<string> block;
+            while (1)
+            {
+                // FIXME: Set limit on memory consumption
+                ReadBlock(block);
+                if (!good())
+                {
+                    clear(rdstate()|ios::badbit);
+#ifdef __EXCEPTIONS
+                    throw runtime_error("FITS file corrupted.");
+#else
+                    gLog << err << "ERROR - FITS file corrupted." << endl;
+                    return;
+#endif
+                }
+
+                if (block.size()%36)
+                    break;
+            }
+
+            if (block.size()==0)
+                break;
+
+            if (block[0].substr(0, 9)=="SIMPLE  =")
+                continue;
+
+            if (block[0].substr(0, 9)=="XTENSION=")
+            {
+                // FIXME: Check for table name
+
+                fTable = Table(block, tellg());
+                fRow   = (size_t)-1;
+
+                //fTable.PrintKeys();
+
+                if (!fTable)
+                {
+                    clear(rdstate()|ios::badbit);
+                    return;
+                }
+
+                fBufferRow.resize(fTable.bytes_per_row);
+                fBufferDat.resize(fTable.bytes_per_row);
+
+                /*
+                 // Next table should start at:
+                 const size_t size = fTable.bytes_per_row*fTable.num_rows;
+                 const size_t blks = size/(36*80);
+                 const size_t rest = size%(36*80);
+
+                seekg((blks+(rest>0?1:0))*(36*80), ios::cur);
+                if (!good())
+                     gLog << err << "File seems to be incomplete (less data than expected from header)." << endl;
+
+                fRow   = fTable.num_rows;
+                 */
+
+                break;
+            }
+        }
+    }
+
+    void ReadRow(size_t row)
+    {
+        // if (row!=fRow+1) // Fast seeking is ensured by izstream
+        seekg(fTable.offset+row*fTable.bytes_per_row);
+
+        fRow = row;
+
+        read(fBufferRow.data(), fBufferRow.size());
+        //fin.clear(fin.rdstate()&~ios::eofbit);
+    }
+
+    template<size_t N>
+        void revcpy(char *dest, const char *src, int num)
+    {
+        const char *pend = src + num*N;
+        for (const char *ptr = src; ptr<pend; ptr+=N, dest+=N)
+            reverse_copy(ptr, ptr+N, dest);
+    }
+
+    bool GetRow(size_t row)
+    {
+        ReadRow(row);
+        if (!good())
+            return good();
+
+        for (Addresses::const_iterator it=fAddresses.begin(); it!=fAddresses.end(); it++)
+        {
+            const Table::Column &c = it->second;
+
+            const char *src = fBufferRow.data() + c.offset;
+            char *dest = reinterpret_cast<char*>(it->first);
+
+            // Let the compiler do some optimization by
+            // knowing the we only have 1, 2, 4 and 8
+            switch (c.size)
+            {
+            case 1: memcpy   (dest, src, c.num*c.size); break;
+            case 2: revcpy<2>(dest, src, c.num);        break;
+            case 4: revcpy<4>(dest, src, c.num);        break;
+            case 8: revcpy<8>(dest, src, c.num);        break;
+            }
+        }
+
+        return good();
+    }
+
+    bool GetNextRow()
+    {
+        return GetRow(fRow+1);
+    }
+
+    bool SkipNextRow()
+    {
+        seekg(fTable.offset+(++fRow)*fTable.bytes_per_row);
+        return good();
+    }
+
+    static bool Compare(const Address &p1, const Address &p2)
+    {
+        return p1.first>p2.first;
+    }
+
+    template<typename T>
+    bool SetPtrAddress(const string &name, T *ptr, size_t cnt)
+    {
+        if (fTable.cols.count(name)==0)
+        {
+            ostringstream str;
+            str <<"SetPtrAddress('" << name << "') - Column not found." << endl;
+#ifdef __EXCEPTIONS
+            throw runtime_error(str.str());
+#else
+            gLog << err << "ERROR - " << str.str() << endl;
+            return false;
+#endif
+        }
+
+        if (sizeof(T)!=fTable.cols[name].size)
+        {
+            ostringstream str;
+            str << "SetPtrAddress('" << name << "') - Element size mismatch: expected "
+                << fTable.cols[name].size << " from header, got " << sizeof(T) << endl;
+#ifdef __EXCEPTIONS
+            throw runtime_error(str.str());
+#else
+            gLog << err << "ERROR - " << str.str() << endl;
+            return false;
+#endif
+        }
+
+        if (cnt!=fTable.cols[name].num)
+        {
+            ostringstream str;
+            str << "SetPtrAddress('" << name << "') - Element count mismatch: expected "
+                << fTable.cols[name].num << " from header, got " << cnt << endl;
+#ifdef __EXCEPTIONS
+            throw runtime_error(str.str());
+#else
+            gLog << err << "ERROR - " << str.str() << endl;
+            return false;
+#endif
+        }
+
+        // if (fAddresses.count(ptr)>0)
+        //     gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
+
+        //fAddresses[ptr] = fTable.cols[name];
+        fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
+        sort(fAddresses.begin(), fAddresses.end(), Compare);
+        return true;
+    }
+
+    template<class T>
+    bool SetRefAddress(const string &name, T &ptr)
+    {
+        return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T));
+    }
+
+    template<typename T>
+    bool SetVecAddress(const string &name, vector<T> &vec)
+    {
+        return SetPtrAddress(name, vec.data(), vec.size());
+    }
+
+    template<typename T>
+        T Get(const string &key) const
+    {
+        return fTable.Get<T>(key);
+    }
+
+    bool SetPtrAddress(const string &name, void *ptr)
+    {
+        if (fTable.cols.count(name)==0)
+        {
+            ostringstream str;
+            str <<"SetPtrAddress('" << name << "') - Column not found." << endl;
+#ifdef __EXCEPTIONS
+            throw runtime_error(str.str());
+#else
+            gLog << err << "ERROR - " << str.str() << endl;
+            return false;
+#endif
+        }
+
+        // if (fAddresses.count(ptr)>0)
+        //     gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl;
+
+        //fAddresses[ptr] = fTable.cols[name];
+        fAddresses.push_back(make_pair(ptr, fTable.cols[name]));
+        sort(fAddresses.begin(), fAddresses.end(), Compare);
+        return true;
+    }
+
+    bool     HasKey(const string &key) const { return fTable.HasKey(key); }
+    int64_t  GetInt(const string &key) const { return fTable.Get<int64_t>(key); }
+    uint64_t GetUInt(const string &key) const { return fTable.Get<uint64_t>(key); }
+    double   GetFloat(const string &key) const { return fTable.Get<double>(key); }
+    string   GetStr(const string &key) const { return fTable.Get<string>(key); }
+
+    size_t GetN(const string &key) const
+    {
+        return fTable.GetN(key);
+    }
+
+    size_t GetNumRows() const { return fTable.num_rows; }
+    size_t GetRow() const { return fRow; }
+
+    operator bool() const { return fTable && fTable.offset!=0; }
+
+    void PrintKeys() const { fTable.PrintKeys(); }
+    void PrintColumns() const { fTable.PrintColumns(); }
+};
+
+};
+#endif
Index: trunk/Mars/mcore/izstream.h
===================================================================
--- trunk/Mars/mcore/izstream.h	(revision 11806)
+++ trunk/Mars/mcore/izstream.h	(revision 11806)
@@ -0,0 +1,176 @@
+#ifndef MARS_izstream
+#define MARS_izstream
+
+#include <istream>
+#include <streambuf>
+
+#ifdef __MARS__
+#ifndef ROOT_TObject
+#include <TObject.h> // Needed for ClassDef
+#endif
+#endif
+
+#ifdef __CINT__
+typedef void *gzFile;
+#else
+#include <zlib.h>
+#endif
+
+class izstream : public std::streambuf, public std::istream
+{
+private:
+    static const int fgBufferSize = 3276804; // maximum size of a fact event + 4
+
+    gzFile fFile;                 // file handle for compressed file
+    char   fBuffer[fgBufferSize]; // data buffer
+
+    int underflow()
+    {
+        if (gptr() && gptr()<egptr())
+            return * reinterpret_cast<unsigned char *>(gptr());
+
+        if (!is_open())
+            return EOF;
+
+        // gptr()-eback(): if more than four bytes are already flushed
+        const int iputback = gptr()-eback()>4 ? 4 : gptr()-eback();
+
+        // Copy the last four bytes flushed into the putback area
+        memcpy(fBuffer+(4-iputback), gptr()-iputback, iputback);
+
+        // Fill the buffer starting at the current file position and reset buffer
+        // pointers by calling setg
+        const int num = gzread(fFile, fBuffer+4, fgBufferSize-4);
+        if (num <= 0) // ERROR or EOF
+            return EOF;
+
+        // reset buffer pointers
+        setg(fBuffer+(4-iputback), fBuffer+4, fBuffer+4+num);
+
+        // return next character
+        return *reinterpret_cast<unsigned char *>(gptr());
+    }
+
+
+public:
+    izstream() : std::istream(this), fFile(0)
+    {
+        setg(fBuffer+4, fBuffer+4, fBuffer+4);
+    }
+    izstream(const char *name) : std::istream(this), fFile(0)
+    {
+        setg(fBuffer+4, fBuffer+4, fBuffer+4);
+        open(name);
+    }
+    ~izstream() { izstream::close(); }
+
+    int is_open() { return fFile!=0; }
+
+    // --------------------------------------------------------------------------
+    //
+    // Open a file by name. Test if it is open like for an ifstream
+    // It doesn't matter whether the file is gzip compressed or not.
+    //
+    void open(const char* name)
+    {
+        if (is_open())
+        {
+            clear(rdstate()|std::ios::failbit);
+            return;
+        }
+
+        fFile = gzopen(name, "rb");
+        if (fFile == 0)
+        {
+            clear(rdstate()|std::ios::failbit);
+            return;
+        }
+    }
+    // --------------------------------------------------------------------------
+    //
+    // Close an open file.
+    //
+    void close()
+    {
+        if (!is_open())
+            return;
+
+        if (gzclose(fFile) != Z_OK)
+            clear(rdstate()|std::ios::failbit);
+
+        fFile = 0;
+    }
+
+    std::streambuf::pos_type seekoff(std::streambuf::off_type offset, std::ios_base::seekdir dir,
+                                     std::ios_base::openmode = std::ios_base::in)
+    {
+        // Using a switch instead results in:
+        //  In member function `virtual std::streampos izstream::seekoff(long int, std::_Ios_Seekdir, std::_Ios_Openmode)':
+        //  warning: enumeration value `_M_ios_seekdir_end' not handled in switch
+        //  warning: case value `0' not in enumerated type `_Ios_Seekdir'
+        //  warning: case value `1' not in enumerated type `_Ios_Seekdir'
+        //  warning: case value `2' not in enumerated type `_Ios_Seekdir'
+
+        if (dir==std::ios::end)
+        {
+            clear(rdstate()|std::ios::failbit);
+            return EOF;
+        }
+
+        // We only do relative seeking to avoid unnecessary decompression
+        // of the whole file
+        if (dir==std::ios::beg)
+            offset -= tellg();
+
+        // Calculate future position in streambuffer
+        const char *ptr = gptr()+offset;
+
+        // This is the number of bytes still available in the buffer
+        const size_t sbuf = egptr()-gptr();
+
+        // Check if the new position will still be in the buffer
+        // In this case the target data was already decompressed.
+        if (ptr>=eback() && ptr<egptr())
+        {
+            // Absolute position in z-stream
+            const z_off_t zpos = gztell(fFile)-sbuf; //gzseek(fFile, 0, SEEK_CUR);
+
+            gbump(offset);
+
+            return zpos+offset;
+        }
+
+        const streampos pos = gzseek(fFile, offset-sbuf, SEEK_CUR);
+
+        // Buffer is empty - force refilling
+        setg(fBuffer+4, fBuffer+4, fBuffer+4);
+
+        return pos<0 ? streampos(EOF) : pos;
+
+        /*
+         // SEEK_END not supported by zlib
+         if (dir==ios::end)
+         {
+             // Position in z-stream
+             const z_off_t zpos = gzseek(fFile, offset, SEEK_END);
+             if (zpos<0)
+                 return EOF;
+
+             return fill_buffer()==EOF ? EOF : zpos;
+         }
+         */
+        return EOF;
+    }
+
+    std::streambuf::pos_type seekpos(std::streambuf::pos_type pos,
+                                     std::ios_base::openmode = std::ios_base::in)
+    {
+        return seekoff(pos, std::ios::beg);
+    }
+
+#ifdef __MARS__
+    ClassDef(izstream, 0) // A C++ wrapper to istream zlib files
+#endif
+};
+
+#endif
