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

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Reiner Rohlfs 2010
!   Author(s): Thomas Bretz  2010 <mailto:thomas.bretz@epfl.ch>
!
!   Copyright: Software Development, 2000-2010
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MCorsikaFormat
//
//////////////////////////////////////////////////////////////////////////////
#include "MCorsikaFormat.h"

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

#include "MLogManip.h"

using namespace std;

const unsigned int MCorsikaFormat::kSyncMarker = 0xd41f8a37;

// --------------------------------------------------------------------------
//
MCorsikaFormat *MCorsikaFormat::CorsikaFormatFactory(const char * fileName)
{
    ifstream * fileIn = new ifstream(fileName);

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

    char buffer[5]="\0\0\0\0";
    fileIn->read(buffer, 4);
    fileIn->seekg(-4, ios::cur);

    if (strcmp(buffer, "RUNH") == 0)
        return new MCorsikaFormatRaw(fileIn);

    if (*reinterpret_cast<unsigned int*>(buffer) == kSyncMarker)
        return new MCorsikaFormatEventIO(fileIn);

    gLog << err << "File " << fileName <<
            " is neither a CORSIKA raw nor EventIO file" << endl;
    delete fileIn;

    return NULL;
}

void MCorsikaFormat::Read(void *ptr, Int_t i) const
{
    fIn->read((char*)ptr, i);
}

// --------------------------------------------------------------------------
//
Bool_t MCorsikaFormat::Eof() const
{
    return fIn->eof();
}

// --------------------------------------------------------------------------
//
streampos MCorsikaFormat::GetCurrPos() const
{
    return fIn->tellg();
}

// --------------------------------------------------------------------------
//
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() const
{
    fIn->seekg(fPrevPos, ios::beg);
}

// --------------------------------------------------------------------------
//
void MCorsikaFormat::StorePos()
{
    fPos = fIn->tellg();
}

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

void MCorsikaFormat::Rewind() const
{
    fIn->seekg(0, 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.
//
Int_t MCorsikaFormatRaw::SeekNextBlock(const char * id, unsigned short type) const
{
    char blockHeader[5]="\0\0\0\0";
    fIn->read(blockHeader, 4);

    if (strcmp(blockHeader, id)==0)
        return kTRUE;

    // 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.
    if (strcmp(blockHeader, "RUNE")==0)
        return kFALSE;

    gLog << err << "ERROR - Wrong identifier: " << id << " expected.";
    gLog << " But read " << blockHeader << " from file." << endl;

    return kERROR;
}

// --------------------------------------------------------------------------
//
void MCorsikaFormatRaw::UnreadLastHeader() const
{
    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[5]="\0\0\0\0";
        fIn->read(runh, 4);

        if (!strcmp(runh, "RUNE"))
        {
            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;
//
Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, Int_t)
{
    static Float_t data[273];
    static Float_t * next = data + 273;

    if (next == data + 273)
    {
        // read next block of events
        if (!ReadData(273, data))
            return kERROR;

        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.
//
Int_t MCorsikaFormatEventIO::SeekNextBlock(const char * id, unsigned short type) const
{
    cout << "Seek " << type << endl;
    while (1)
    {
        // we read - synchronisation marker
        //         - type / version field
        //         - identification field
        //         - length
        //         - unknown field
        //         - id (first 4 bytes of data field)
        int blockHeader[4];
        fIn->read((char*)blockHeader, 4 * sizeof(int));

        if (fIn->eof())
        {
            gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - Unexpected end-of-file seeking " << id << " (" << type << ")." << endl;
            return kERROR;
        }

        const unsigned short objType = blockHeader[1] & 0xFFFF;
        cout << "objType=" << objType << endl;
        if (type == objType)
        {
            if (!id)
                break;

            char c[9];
            fIn->read(c, 8);

            if (memcmp(c+4, id, 4)==0)
                break;

            c[8] = 0;
            gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - '" << c+4 << "' doesn't match expected one '" << id << "' in object " << type << endl;
            return kFALSE;
        }

        // 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
        if (type == 1202 && objType == 1210)
            return kFALSE;

        // a unknown block, we jump to the next one
        const int length = blockHeader[3] & 0x3fffffff;
        fIn->seekg(length, ios::cur);
    }


    return kTRUE;
}

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

// --------------------------------------------------------------------------
//
Bool_t MCorsikaFormatEventIO::SeekEvtEnd()
{
    if (fRunePos != streampos(0))
    {
        fIn->seekg(fRunePos, ios::beg);
        return kTRUE;
    }

    // it is the first time we are looking for the RUNE block

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

    fIn->seekg(-32, ios::end);

    unsigned int blockHeader[4];
    fIn->read((char*)blockHeader, 4 * sizeof(int));

    if ( blockHeader[0]               == kSyncMarker &&
        (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)!=kTRUE)
        return kFALSE;

    UnreadLastHeader();
    fRunePos = fIn->tellg();

    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;                                                           
//
Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, Int_t telescope)
{
    static Float_t *data = NULL;

    static int lengthPhotonData = 0;
    static int lengthEvent      = 0;

    if (lengthPhotonData>0)
    {
        cout << "Return Bunch l2=" << lengthPhotonData << endl;

        lengthPhotonData -= 8 * sizeof(Float_t);
        *buffer += 8; // Return the pointer to the start of the current photon data
        return kTRUE;
    }

    if (data)
    {
        delete [] data;
        data = NULL;
        cout << "Return end-of-tel LE=" << lengthEvent << endl;
        return kFALSE;
    }

    if (lengthEvent==0)
    {
        cout << "Searching 1204... " << flush;
        // The length of the block is stored and we use it to determine
        // when the next top level block is expected
        const Int_t rc = NextEventObject(lengthEvent);
        if (rc==kERROR)
            return kERROR;
        if (!rc)
        {
            cout << "skip EVTE." << endl;
            return kFALSE;
        }

        cout << "found new array." << endl;
    }

    cout << "Searching 1205..." << flush;

    // Look for the next event photon bunch (object type 1205)
    const Int_t tel = NextPhotonObject(lengthPhotonData);
    if (tel<0)
        return kERROR;

    lengthEvent -= lengthPhotonData+12; // Length of data + Length of header

    lengthPhotonData -= 12;     // Length of data before the photon bunches
    fIn->seekg(12, ios::cur);   // Skip this data

    cout << " found (tel=" << tel << ",LEN=" << lengthEvent << ",len=" << lengthPhotonData << ")... " << flush;

    if (lengthPhotonData==0)
    {
        cout << "empty!" << endl;
        return kFALSE;
    }

    const UInt_t cnt = lengthPhotonData / sizeof(Float_t);
    // Read next object (photon bunch)
    data = new Float_t[cnt];
    if (!ReadData(cnt, data, 0))
    {
        delete [] data;
        data = NULL;
        return kERROR;
    }

    cout << "return!" << endl;

    lengthPhotonData -= 8 * sizeof(Float_t);
    *buffer = data; // Return the pointer to the start of the current photon data

    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
//
Int_t MCorsikaFormatEventIO::NextEventObject(Int_t &length) const
{
    while (1)
    {
        // we read - synchronisation marker
        //         - type / version field
        //         - identification field
        //         - length
        UInt_t blockHeader[4];
        fIn->read((char*)blockHeader, 4 * sizeof(Int_t));

        if (fIn->eof())
        {
            gLog << err << "MCorsikaFormatEventIO::NextEventObject: ERROR - Unexpected end-of-file." << endl;
            return kERROR;
        }

        // For speed reasons we skip the check of the synchronization marker

        // Decode the object type
        const unsigned short objType = blockHeader[1] & 0xFFFF;

        cout << "objType=" << objType << endl;

        // Decode length of block
        length = blockHeader[3] & 0x3fffffff;

        // Check if we found the next array (reuse / impact parameter)
        // blockHeader[2] == array number (reuse)
        if (objType==1204)
            return kTRUE;

        // we found an eventEnd block, reset file pointer
        if (objType==1209)
            break;

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

    fIn->seekg( -4 * (Int_t)(sizeof(Int_t)), ios::cur);
    length = 0;

    return kFALSE;
}

// --------------------------------------------------------------------------
//
Int_t MCorsikaFormatEventIO::NextPhotonObject(Int_t &length) const
{
    UInt_t blockHeader[3];

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

    if (fIn->eof())
    {
        gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected end-of-file.." << endl;
        return -1;
    }

    const unsigned short objType = blockHeader[0] & 0xFFFF;
    if (objType != 1205)
    {
        gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object type " << objType << " (expected 1205)." << endl;
        return -1;
    }

    const unsigned short version = (blockHeader[0] >> 20) & 0x0FFF;
    if (version != 0)
    {
        gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object version " << version << " (expected 0 for object type 1205)." << endl;
        return -1;
    }

    length = blockHeader[2] & 0x3fffffff;

    // blockHeader[1] == 1000 x array number (reuse) + telescope number
    return blockHeader[1] % 1000;
}