Skip to content
Snippets Groups Projects
OdfIO.cpp 7.87 KiB
Newer Older
Brown's avatar
Brown committed
#include "OdfIO.h"

namespace sammy{
    
    OdfIO::OdfIO() {}

    double OdfIO::readFloat(std::istream &fs, bool unix){
    
        unsigned char tmp1,tmp2,tmp3,tmp4;
        int i=0;
        float value;
       
        if(!unix){ // wacky VMS
            fs.read( (char*)&tmp1, 1 ); // 1 byte each
            fs.read( (char*)&tmp2, 1 );
            fs.read( (char*)&tmp3, 1 );
            fs.read( (char*)&tmp4, 1 );
            if((tmp2 & 0x7f) == 0){
                if(tmp1+tmp3+tmp4 != 0) {
                    tmp1=tmp2=tmp3=0;
                }   
                tmp2=0;
            }
            else tmp2--;
            
            i=(tmp2 << 24) + (tmp1 << 16) + (tmp4 << 8) + tmp3;
        }
        else{
            fs.read( (char*)&tmp1, 1 );
            fs.read( (char*)&tmp2, 1 );
            fs.read( (char*)&tmp3, 1 );
            fs.read( (char*)&tmp4, 1 );
            i=(tmp4 << 24) + (tmp3 << 16) + (tmp2 << 8 ) +tmp1;
        }
        // grab int pointer, cast to float, dereference
        value = *reinterpret_cast<float *>( &i );
        return (double)value;
    }

    int OdfIO::readInteger(std::istream &fs){
        size_t intsize = sizeof(int);
        int val;
        fs.read( (char*)&val, intsize );
        return val;
    }

    std::vector<double> OdfIO::getSection(int start, int readSize, int section, 
                                          std::istream &in){
        in.seekg(0,std::ios::beg);
        std::vector<double> dataSec;
        int block,skipped,channel,max;
        size_t intsize = sizeof(int);

        dataSec.resize(readSize);

        block = start/125;
        channel = start - block * 125;
        block = head["ifb"] + block * head["numSection"];
        block--; 
        skipped=( block + (section-1) )*512;     // skip 512 bytes per block
        skipped += 4;                            // account for control word
        skipped += channel*4;                    // and add channels not to read
        in.ignore(skipped);                      // jump skipped-many bytes, alt: in.seekg(skipped,std::ios::cur);

        max = std::min(readSize,125-channel);

        int read = 0;
        for(int i=0;i<max;i++){
            if(      head["mode"] == 3 ) dataSec[read] = readFloat(in,head["unix"]);
            else if( head["mode"] == 1 ) dataSec[read] = readInteger(in);
            read++;
        }
        in.ignore(intsize*2); // skip over two ints in the block which contain no data

        while( read<readSize ){
            skipped=(head["numSection"]-1)*512;  // skip over other sections
            in.ignore(skipped);
            in.ignore(intsize);          // discard control word

            for( int i=0; i<125 && read<readSize; i++ ) {
                if(     head["mode"] == 3) dataSec[read] = readFloat(in,head["unix"]);
                else if(head["mode"] == 1) dataSec[read] = readInteger(in);
                read++;
            }
            in.ignore(intsize*2); // skip over two ints in the block which contain no data
        }

        return dataSec;
    }

    std::vector<std::vector<double>> OdfIO::getAllSections(std::istream &in){
        // --- go to start of file ---
        in.seekg(0, std::ios::beg);

        std::vector<std::vector<double>> data;
        data.resize(head["numSection"]);

        // --- grab each section ---
        for( int section=1; section<=head["numSection"]; ++section ){
            data[section-1].resize(head["numChan"]);
            data[section-1] = getSection(0,head["numChan"],section,in);
        }

        return data;
    }

    void OdfIO::readHeader(std::istream &in, int unix){
        size_t intsize = sizeof(int);
        int bytes_per_block = 512;
        int bytes_per_word = 4;
        head["unix"] = unix;

        // comments based on OPRODF manual: Jack Craven, report no. ORNL/CSD/TM-45, 1978
        //  ** keep in mind that changes may have happened between that report and
        //  ** the writing of the files!!! For example it seems that 36 bit ints are no 
Brown's avatar
Brown committed
        //  ** longer used (32 bit now). 
Brown's avatar
Brown committed
        in.read( (char*)&head["mode"], intsize ); // 0,1,3 for 18-bit ints, 36-bit ints, floating-point vals
        in.read( (char*)&head["source"], intsize ); // 0,1,2 for SEL, CSISRS, or ENDF/B data
        in.read( (char*)&head["irun"], intsize ); // run number
        in.read( (char*)&head["ncblks"], intsize ); // starting block no. for the comment section
        in.read( (char*)&head["ncwrds"], intsize ); // number of words in comment section
        in.read( (char*)&head["nsblks"], intsize ); // starting block for scalar section
        in.read( (char*)&head["nswrds"], intsize ); // number of words in scalar section
        in.read( (char*)&head["ncstrt"], intsize ); // starting word number of DAQ scalar and counter section in scalar section
        in.read( (char*)&head["ncntrs"], intsize ); // number of words in DAQ scalar and counter section
        in.read( (char*)&head["nxstrt"], intsize ); // starting word number of DAQ variable section in scalar section
        in.read( (char*)&head["nxwrds"], intsize ); // number of words in DAQ variable section
        in.read( (char*)&head["npblks"], intsize ); // starting block number of parameter section
        in.read( (char*)&head["npwrds"], intsize ); // number of words in parameter section
        in.read( (char*)&head["ndtype"], intsize ); // 0,1 for data section described by table in parameter section, data section is dependent on dataset 1 (i.e. datasets 2-x correspond to energy in dataset 1)
        in.read( (char*)&head["numSection"], intsize ); // number of sections in each dataset
        in.read( (char*)&head["ifb"], intsize ); // 
        in.read( (char*)&head["numChan"], intsize ); // number of channels in each section
        in.read( (char*)&head["zan"], intsize ); // ENDF/B charge/mass id
        in.read( (char*)&head["awr"], intsize ); // ENDF/B ratio of nuclear mass (iso,ele,molec,mix) to neutron
        in.read( (char*)&head["mat"], intsize ); // ENDF/B MAT number
        in.read( (char*)&head["mf"], intsize ); // ENDF/B file number
        in.read( (char*)&head["mt"], intsize ); // ENDF/B reaction type
        in.read( (char*)&head["varswt"], intsize ); // 0,1 for dataset 1 decreasing, and increasing in value (ignored if ndtype is 0)
        in.read( (char*)&head["dctswt"], intsize ); // 0,1 for not dead-time corrected, and dead-time corrected
        in.read( (char*)&head["strt"], intsize ); // starting word number of each dataset in data section (mode 0 only)
        in.read( (char*)&head["ndwend"], intsize ); // ending word number of last word written in block NDBEND

        
        /* TODO: the following doesn't work properly yet, we ignore for now */
        // ------------------------------------
        // --- Read crunch table
        // ------------------------------------
        // energy index table. each word contains the largest energy for each n blocks in the data set. N=(NDWRDS/125)+1
        ndtbl.resize(100);
        for( int i=0;i<100;++i ) in.read( (char*)&ndtbl[i], intsize );
        // ------------------------------------
        // --- comments section ---
        // ------------------------------------
        int bytes_to_comments = (head["ncblks"]) * bytes_per_block;
        std::vector<unsigned char> comments;
        int bytes_in_comments = head["ncwrds"]*bytes_per_word;
        comments.resize( bytes_in_comments );

        in.seekg(bytes_to_comments,std::ios::beg);
        in.ignore(bytes_per_word);

        for( int i=0;i<bytes_in_comments;++i ){
            in.read( (char*)&comments[i], 1 );
            comments[i] = comments[i];
        }
        // ------------------------------------
        // --- scalars section ---
        // ------------------------------------
        int bytes_to_scalars = (head["nsblks"]) * bytes_per_block;
    }

    void OdfIO::printHeader(){
        for( auto member : head ){
            std::cout << member.first << " = " << head[member.first] << std::endl;
        }

    }

}