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

namespace sammy{
    
    odfIO::odfIO() {}

    double odfIO::readFloat(std::ifstream &fs, bool unix){
    
        unsigned char tmp1,tmp2,tmp3,tmp4;
        int i=0;
        float value;
        float * dummy_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 );
            // std::cout << "1,2,3,4 = " << (int)tmp1<<","<< (int)tmp2<<","<< (int)tmp3<<","<< (int)tmp4<<","<<std::endl;
            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;   
        } 
        dummy_value=(float*) &i;
        value = *dummy_value;
        return value;
    }

    std::vector<double> odfIO::getSection(int start, int readSize, int section, 
                                          header &hd, std::ifstream &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 = hd.ifb + block * hd.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(      hd.mode == 3 ) dataSec[read] = readFloat(in,hd.unix);
            else if( hd.mode == 1 ) dataSec[read] = readFloat(in,hd.unix); // needs to read integer not float
            read++;
        }
        in.ignore(intsize*2); // skip over two ints in the block which contain no data

        while( read<readSize ){
            skipped=(hd.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(     hd.mode == 3) dataSec[read] = readFloat(in,hd.unix);
                else if(hd.mode == 1) dataSec[read] = readFloat(in,hd.unix); // needs to read integer not float
                read++;
            }
            in.ignore(intsize*2); // skip over two ints in the block which contain no data
        }

        return dataSec;
    }

    void odfIO::readHeader(std::ifstream &in, header &head){
        size_t intsize = sizeof(int);
        int bytes_per_block = 512;
        int bytes_per_word = 4;
        // comments based on OPRODF manual: Jack Craven, report no. ORNL/CSD/TM-45, 1978
        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
        

        /* the following doesn't work properly yet */
        // ------------------------------------
        // --- Read crunch table
        // ------------------------------------
        // energy index table. each word contains the largest energy for each n blocks in the data set. N=(NDWRDS/125)+1
        head.ndtbl.resize(100);
        for( int i=0;i<100;++i ) in.read( (char*)&head.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);
        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(header &head){
        std::cout << "mode       = " << head.mode       << std::endl       
                  << "source     = " << head.source     << std::endl     
                  << "irun       = " << head.irun       << std::endl       
                  << "ncblks     = " << head.ncblks     << std::endl     
                  << "ncwrds     = " << head.ncwrds     << std::endl     
                  << "nsblks     = " << head.nsblks     << std::endl     
                  << "nswrds     = " << head.nswrds     << std::endl     
                  << "ncstrt     = " << head.ncstrt     << std::endl     
                  << "ncntrs     = " << head.ncntrs     << std::endl     
                  << "nxstrt     = " << head.nxstrt     << std::endl     
                  << "nxwrds     = " << head.nxwrds     << std::endl 
                  << "npblks     = " << head.npblks     << std::endl 
                  << "npwrds     = " << head.npwrds     << std::endl 
                  << "ndtype     = " << head.ndtype     << std::endl 
                  << "numSection = " << head.numSection << std::endl 
                  << "ifb        = " << head.ifb        << std::endl 
                  << "numChan    = " << head.numChan    << std::endl 
                  << "zan        = " << head.zan        << std::endl 
                  << "awr        = " << head.awr        << std::endl 
                  << "mat        = " << head.mat        << std::endl 
                  << "mf         = " << head.mf         << std::endl 
                  << "mt         = " << head.mt         << std::endl 
                  << "varswt     = " << head.varswt     << std::endl  
                  << "dctswt     = " << head.dctswt     << std::endl  
                  << "strt       = " << head.strt       << std::endl  
                  << "ndwend     = " << head.ndwend     << std::endl; 

        // for( int i=0;i<ndtbl.size();++i ){
        //     std::cout << ndtbl[i] << " ";
        //     if( (i+1)%10==0 ) std::cout << std::endl;
        // }
    }

}