source: trunk/Mars/mcorsika/MCorsikaFormat.cc@ 9937

#include <fstream>
#include <errno.h>


#include "MLogManip.h"

#include "MCorsikaFormat.h"


using namespace std;

MCorsikaFormat * CorsikaFormatFactory(MLog * log, const char * fileName)
{
    ifstream * fileIn = new ifstream(fileName);

    const Bool_t noexist = !(*fileIn);
    if (noexist)
    {
        *log << err << "Cannot open file " << fileName << ": ";
        *log << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl;
        delete fileIn;
        fileIn = NULL;
        return NULL;
    }


    if (!fileIn->is_open())
        {
        *log << err << "Failed to open file " << fileName << endl;
        delete fileIn;
        fileIn = NULL;
        return NULL;
        }
    
    char buffer[4];
    fileIn->read(buffer, 4);
    fileIn->seekg(-4, ios::cur);

    int * syncmarker = reinterpret_cast<int*>(buffer);
    if (memcmp(buffer, "RUNH", 4) == 0)
        return new MCorsikaFormatRaw(log, fileIn);
         
    else if (*syncmarker == -736130505)
        return new MCorsikaFormatEventIO(log, fileIn);

    *log << err << "File " << fileName << 
            " is neither a corsica file, nor an eventio file" << endl;
    delete fileIn;
    fileIn = NULL;

    return NULL;
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

Bool_t MCorsikaFormat::ReadData(Int_t numValues, Float_t * buffer, 
                                Int_t minSeekValues)
{
    fPrevPos = fIn->tellg();

    fIn->read((char*)buffer, numValues * sizeof(Float_t));

    if (numValues < minSeekValues)
        // skip the remaining data of this block
        fIn->seekg( (minSeekValues - numValues) * 4, ios::cur);

    return !fIn->fail();

}

void MCorsikaFormat::UnreadLastData()
{
    fIn->seekg(fPrevPos, ios::beg);
}

void MCorsikaFormat::StorePos()
{
    fPos = fIn->tellg();
//*fLog << all << "storePos: " << fPos << endl;
}

void MCorsikaFormat::ResetPos()
{
    fIn->seekg(fPos, ios::beg);
}


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////


MCorsikaFormat::~MCorsikaFormat() 
{
    delete fIn;
}


// --------------------------------------------------------------------------
//
// Reads the next 4 bytes which should be the id. (for example "RUNH", RUNE".
// "EVTH")
// Returns kTRUE if the functions finds the id
//         kFALSE if the functions does not find the "id" as the next 4 
//                bytes in the file.
// After this call the file position pointer points just after the 4 bytes
// of the id.
//
Bool_t MCorsikaFormatRaw::SeekNextBlock(const char * id, unsigned short type)
{
    char blockHeader[5];
    fIn->read(blockHeader, 4);

    if (memcmp(blockHeader, id, 4))
    {
        blockHeader[4] = 0;
        if (strcmp(id, "EVTH") != 0 || strcmp(blockHeader, "RUNE") != 0 )
            {
            // at the end of a file we are looking for the next Event header,
            // but find the end of a run. This is expected, therefor no error
            // message.
            *fLog << err << "ERROR - Wrong identifier: " << id << " expected." 
                  << " But read " << blockHeader << " from file." << endl;
            }
        return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
void MCorsikaFormatRaw::UnreadLastHeader() 
{
    fIn->seekg(-4, ios::cur);
}

Bool_t MCorsikaFormatRaw::SeekEvtEnd()
{

    // Search subblockwise backward (Block: 5733*4 = 21*273*4)
    for (int i=1; i<22; i++)
    {
        fIn->seekg(-i*273*4, ios::end);

        char runh[4];
        fIn->read(runh, 4);

        if (!memcmp(runh, "RUNE", 4))
        {
            fIn->seekg(-4, ios::cur);
            return kTRUE;
        }
    }
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//                                                                           
// Returns one event (7 Float values) after the other until the EventEnd     
// block is found.                                                           
// If a read error occurred, the readError is set to kTRUE and the function  
// returns kFALSE;                                                           
Bool_t MCorsikaFormatRaw::GetNextEvent(Float_t ** buffer, Bool_t & readError) 
{
    static Float_t data[273];
    static Float_t * next = data + 273;

    if (next == data + 273)
        {
        // read next block of events
        if (!ReadData(273, data))
            {
            readError = kTRUE;
            return kFALSE;
            }

        if (!memcmp(data, "EVTE", 4))
            {
            // we found the end of list of events
            UnreadLastData();
            return kFALSE;
            }

        next = data;
        }

    *buffer = next;
    next += 7;

    return kTRUE;

}


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

// --------------------------------------------------------------------------
//
// Jumps from one top level object to the next until if finds the object with 
// correct type and correct id. The id is identifier of the raw corsika block,
// for example "RUNH", RUNE", "EVTH"
//
// Returns kTRUE if the functions finds the type / id
//         kFALSE if the functions does not find the type / id.
//
// After this call the file position pointer points just after the 4 bytes
// of the id.
//
Bool_t MCorsikaFormatEventIO::SeekNextBlock(const char * id, unsigned short type)
{
    int blockHeader[7];

    while (1)
    {
        // we read - synchronisation marker
        //         - type / version field
        //         - identification field
        //         - length
        //         - unknown field
        //         - id (first 4 bytes of data field)
        fIn->read((char*)blockHeader, 6 * sizeof(int));

        if (fIn->eof())
            {
            *fLog << err << "ERROR - Missing identifier: " << id  <<
                   " type: " << type << endl;
            return kFALSE;
            }

        unsigned short fileType = blockHeader[1] & 0xFFFF;
        if (type == fileType                       )
            break;

         if (type == 1202 && fileType == 1210)
            // we are looking for a event header, but found a runEnd. 
            // This will happen at the end of a run. We stop here looking for
            // the next event header
            return kFALSE;

        // a unknown block, we jump to the next one
//*fLog << "unknown: " <<  id << " type: " << fileType << " sub-blocks: " <<  (blockHeader[3]>>29);
//*fLog <<  " length: " << (blockHeader[3] & 0x3fffffff) <<   "  pos: " << fIn->tellg() << endl;
        int length = blockHeader[3] & 0x3fffffff;
        
        fIn->seekg(length - 2 * (Int_t)(sizeof(int)), ios::cur);    
    }


    return kTRUE;
}

// --------------------------------------------------------------------------
//
void MCorsikaFormatEventIO::UnreadLastHeader() 
{
    fIn->seekg( -6 * (Int_t)(sizeof(int)), ios::cur);
}

// --------------------------------------------------------------------------
//
Bool_t MCorsikaFormatEventIO::SeekEvtEnd()
{
    if (fRunePos != streampos(0))
        {
        fIn->seekg(fRunePos, ios::beg);
        return kTRUE;
        }
    else
        {
        // it is the first time we are looking for the RUNE block

        // is the RUNE block at the very end of the file?
        std::streampos currentPos = fIn->tellg();

        fIn->seekg(-32, ios::end);
        unsigned int blockHeader[4];
        fIn->read((char*)blockHeader, 4 * sizeof(int));
        if ( blockHeader[0]               == 0xd41f8a37 &&
            (blockHeader[1] & 0xffff)     == 1210       &&
            (blockHeader[3] & 0x3fffffff) == 16)
            {
            // this seams to be a RUNE (1210)  block
            fIn->seekg( -4 * (Int_t)(sizeof(int)), ios::cur);
            fRunePos = fIn->tellg();
            return kTRUE;
            }

        // we do not find a RUNE block at the end of the file
        // we have to search in the file
        fIn->seekg(currentPos, ios::beg);
        if (SeekNextBlock("RUNE", 1210))
            {
            UnreadLastHeader();
            fRunePos = fIn->tellg();
            return kTRUE;
            }
        }
    return kFALSE;
}

// --------------------------------------------------------------------------
//                                                                           
// Returns one event (7 Float values) after the other until the EventEnd     
// block is found.                                                           
// If a read error occurred, the readError is set to kTRUE and the function  
// returns kFALSE;                                                           
Bool_t MCorsikaFormatEventIO::GetNextEvent(Float_t ** buffer, 
                                           Bool_t & readError) 
{

    static Float_t * data = NULL;
    static Float_t * next;
    static Int_t   topLevelLength = 0;
    static Int_t   eventLength = 0;


    while (eventLength == 0)
        {
        delete [] data;
        data = NULL;

        if (topLevelLength == 0)
            {
            if (!NextTopLevelBlock(topLevelLength, readError))
                return kFALSE;
            }

        if (!NextEventBlock(eventLength, readError))
            return kFALSE;
        topLevelLength -= eventLength + 3 * sizeof(int);

        // read next block of events
        data = new Float_t [eventLength / sizeof(Float_t)];
        if (!ReadData(eventLength / sizeof(Float_t), data, 0))
            {
            delete [] data;
            data = NULL;
            readError = kTRUE;
            return kFALSE;
            }
        next = data + 3;        
        eventLength -= 3 * sizeof(Float_t);
        }

    eventLength -= 8 * sizeof(Float_t);
    *buffer = next;
    next += 8;


    return kTRUE;
}
// --------------------------------------------------------------------------
//                                                                           
// Looks for the next Block with type 1204 and return kTRUE.                 
// The function also stops moving forward in the file, if it finds a         
// EventEnd block (1209). In this case kFALSE is returned                    
Bool_t MCorsikaFormatEventIO::NextTopLevelBlock(Int_t & length, 
                                                Bool_t & readError)
{
    Int_t blockHeader[4];

    while (1)
    {
        // we read - synchronisation marker
        //         - type / version field
        //         - identification field
        //         - length
        fIn->read((char*)blockHeader, 4 * sizeof(Int_t));

        if (fIn->eof())
            {
            *fLog << err << "ERROR - Missing identifier: 1204 or 1209" << endl;
            readError = kTRUE;
            return kFALSE;
            }

        length = blockHeader[3] & 0x3fffffff;
        unsigned short fileType = blockHeader[1] & 0xFFFF;
        if (fileType == 1204)
            return kTRUE;

        if (fileType == 1209)
            {
            // we found an eventEnd block, reset file pointer
            fIn->seekg( -4 * (Int_t)(sizeof(Int_t)), ios::cur);
            length = 0;
            return kFALSE;
            }


        // a unknown block, we jump to the next one
        fIn->seekg(length, ios::cur);    
    }

    return kTRUE;
}

Bool_t MCorsikaFormatEventIO::NextEventBlock(Int_t & length, 
                                             Bool_t & readError)
{
    Int_t blockHeader[3];

    // we read - synchronisation marker
    //         - type / version field
    //         - identification field
    //         - length
    fIn->read((char*)blockHeader, 3 * sizeof(Int_t));

    if (fIn->eof())
        {
        *fLog << err << "ERROR - Missing identifier: 1205" << endl;
        readError = kTRUE;
        return kFALSE;
        }
    
    unsigned short fileType = blockHeader[0] & 0xFFFF;
    if (fileType != 1205)
        {
        *fLog << err << "ERROR - Unexpected type: " << fileType << "expected 1205" << endl;
        readError = kTRUE;
        return kFALSE;
        }

    unsigned short version = (blockHeader[0] >> 20) & 0x0FFF;
    if (version != 0)
        {
        *fLog << err << "ERROR - Unexpected version: " << version << "expected: 0" << endl;
        readError = kTRUE;
        return kFALSE;
        }

    length = blockHeader[2] & 0x3fffffff;

    return kTRUE;

}