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

/* ======================================================================== *\
!
! *
! * 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;


// --------------------------------------------------------------------------
//
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;
    }

    uint32_t magic = 0;
    fileIn->read((char*)&magic, 4);

    uint32_t blocklength = 0;

    // This seems to be a new corsika binary identifier
    if (magic==kBlockLengthRaw || magic==kBlockLengthThin)
    {
        blocklength = magic;
        fileIn->read((char*)&magic, 4);
    }

    // Jump back that the "RUNH" is can be read again
    fileIn->seekg(-4, ios::cur);

    if (magic==kRUNH)
    {
        gLog << inf2 << "Corsika RAW" << (blocklength==kBlockLengthThin?"(+THIN)":"") << " format detected." << endl;
        return new MCorsikaFormatRaw(fileIn, blocklength);
    }

    if (magic==kSyncMarker)
    {
        gLog << inf2 << "Corsika EventIO format detected." << endl;
        return new MCorsikaFormatEventIO(fileIn);
    }

    gLog << err << "File " << fileName << Form(" (%08x) ", magic);
    gLog << "is neither a CORSIKA raw nor EventIO file" << endl;

    delete fileIn;

    return NULL;
}

Bool_t MCorsikaFormat::Read(void *ptr, Int_t i) const
{
   fIn->read((char*)ptr, i);
   return !fIn->fail();

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

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


// --------------------------------------------------------------------------
//
// After a call to this function, the file pointer is located after the   
// header of the block. As the event block has no header it is located    
// at the beginning of this block, which is as the beginning of the data  
Bool_t MCorsikaFormatRaw::NextBlock(Int_t   readState,
                                    Int_t & blockType, 
                                    Int_t & blockVersion,
                                    Int_t & blockIdentifier, 
                                    Int_t & blockLength) const
{
    if (fBlockLength)
    {
        // In the new corsika format each block
        // starts and end with the block length
        const size_t position = fIn->tellg()%(fBlockLength+8);

        // Whenever we are prior to the end of such a block,
        // we read and check the end and start tag of this
        // and the following block
        if (position==fBlockLength+4)
        {
            uint32_t h[2];
            fIn->read((char*)h, 8);
            if (fIn->eof())
                return kFALSE;

            if (h[0]!=fBlockLength || h[1]!=fBlockLength)
            {
                gLog << err << "ERROR - Block length missing at the end or beginning of a block." << endl;
                return kERROR;
            }
        }
    }

    uint32_t blockHeader = 0;
    fIn->read((char*)&blockHeader, 4);
    if (fIn->eof())
        return kFALSE;

    blockVersion    = 0;
    blockIdentifier = 0;
    blockLength     = fBlockLength/21 - 4;  // in bytes

    switch(blockHeader)
      {
      case kRUNH:
         blockType = 1200;
         break;

      case kRUNE:
         blockType = 1210;
         break;

      case kEVTH:
         if (readState != 10)
            blockType = 1202;     // event header (where readstate := 2)
         else
            {
            blockType = 1105;              // raw data
            fIn->seekg(-4, ios::cur);
            blockLength += 4;
            }
         break;

      case kEVTE:
         blockType = 1209;
         break;

      default:          // the events, they don't have a specific header
         blockType = 1105;                  // raw data
         fIn->seekg(-4, ios::cur);
         blockLength += 4;
      }            
    return kTRUE;
}
// --------------------------------------------------------------------------
//
Bool_t MCorsikaFormatRaw::SeekEvtEnd()
{
    const uint32_t step   =  fBlockLength/21;
    const uint32_t offset = (fBlockLength?4:0);

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

        uint32_t magic = 0;
        fIn->read((char*)&magic, 4);

        if (magic==kRUNE)
        {
//            fIn->seekg(-4, ios::cur);
            return kTRUE;
        }
    }

    return kFALSE;
}
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

Bool_t MCorsikaFormatEventIO::NextBlock(Int_t   readState,
                                        Int_t & blockType, 
                                        Int_t & blockVersion,
                                        Int_t & blockIdentifier, 
                                        Int_t & blockLength) const
{
// we read - synchronisation markerif not subBlock
//         - type / version field
//         - identification field
//         - length
//         - unknown field
//         - id (first 4 bytes of data field)

   if (readState == 4)
//   if (subBlock)
      {
      // this is a sub-block
      int blockHeader[3];
      fIn->read((char*)blockHeader, 3 * sizeof(int));

      if (fIn->eof())
         return kFALSE;


      blockType       = blockHeader[0] & 0xFFFF;
      blockVersion    = (blockHeader[0] & 0xFFF00000) >> 20;
      blockIdentifier = blockHeader[1];
      blockLength     = blockHeader[2] & 0x3FFFFFFF;
      }
   else
      {
       int blockHeader[4];
       fIn->read((char*)blockHeader, 4 * sizeof(int));

       if (fIn->eof())
           return kFALSE;


       blockType       = blockHeader[1] & 0xFFFF;
       blockVersion    = (blockHeader[1] & 0xFFF00000) >> 20;
       blockIdentifier = blockHeader[2];
       blockLength     = blockHeader[3] & 0x3FFFFFFF;

       if (blockType == 1200  || blockType == 1210 ||
           blockType == 1202  || blockType == 1209    )
           {
           // read the "old" corsika header plus additional 4 unknown byte
           char tmp[8];
           fIn->read(tmp, 8);
           blockLength -= 8;
           }
    
      }
    return kTRUE;
}

// --------------------------------------------------------------------------
//
Bool_t MCorsikaFormatEventIO::SeekEvtEnd()
{

    // the RUNE block it at the very end of the file.
    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( 8, ios::cur);
        return kTRUE;
    }

    return kFALSE;
}