Index: fact/tools/pyscripts/pyfact/fits.h
===================================================================
--- fact/tools/pyscripts/pyfact/fits.h	(revision 13319)
+++ fact/tools/pyscripts/pyfact/fits.h	(revision 13319)
@@ -0,0 +1,718 @@
+#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"
+#define ___err___ err
+#define ___all___ all
+#endif
+
+#define HAVE_ZLIB
+#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;
+            string unit;
+        };
+
+        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());
+                const string unit = Get<string>("TUNIT"+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, unit };
+
+                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(bool display_all=false) const
+        {
+            for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++)
+            {
+                if (!display_all &&
+                    (it->first.substr(0, 6)=="TTYPE" ||
+                     it->first.substr(0, 6)=="TFORM" ||
+                     it->first.substr(0, 6)=="TUNIT" ||
+                     it->first=="TFIELDS"  ||
+                     it->first=="XTENSION" ||
+                     it->first=="NAXIS"    ||
+                     it->first=="BITPIX"   ||
+                     it->first=="PCOUNT"   ||
+                     it->first=="GCOUNT")
+                   )
+                    continue;
+
+                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 << "| ";
+                gLog << it->second.num << 'x';
+                switch (it->second.type)
+                {
+                case 'L': gLog << "bool(8)";    break;
+                case 'B': gLog << "byte(8)";    break;
+                case 'I': gLog << "short(16)";  break;
+                case 'J': gLog << "int(32)";    break;
+                case 'K': gLog << "int(64)";    break;
+                case 'E': gLog << "float(32)";  break;
+                case 'D': gLog << "double(64)"; break;
+                }
+                gLog << ": " << it->first.second << " [" << it->second.unit << "]" << 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>();
+        }
+
+            // Values of keys are always signed
+        template<typename T>
+            T Get(const string &key, const string &deflt) const
+        {
+            const map<string,Entry>::const_iterator it = keys.find(key);
+            return it==keys.end() ? deflt :it->second.Get<T>();
+        }
+
+        size_t GetN(const string &key) const
+        {
+            const Columns::const_iterator it = cols.find(key);
+            if (it==cols.end()){
+                ostringstream str;
+                str << "Key '" << key << "' not found." << endl;
+                str << "Possible keys are:" << endl;
+                for ( Columns::const_iterator it=cols.begin() ; it != cols.end(); ++it){
+                    str << it->first << endl;
+                }
+#ifdef __EXCEPTIONS
+                throw runtime_error(str.str());
+#else
+                gLog << ___err___ << "ERROR - " << str.str() << endl;
+                return 0;
+#endif
+            }
+            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);
+    }
+
+    template<typename T>
+        T Get(const string &key, const string &deflt) const
+    {
+        return fTable.Get<T>(key, deflt);
+    }
+
+    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