source: trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.cc@ 9251

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz  11/2008 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: Software Development, 2000-2008
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MCorsikaRead
//
//  Input Containers:
//   -/-
//
//  Output Containers:
//   MCorsikaRunHeader
//   MCorsikaEvtHeader
//   MPhotonEvent
//
//////////////////////////////////////////////////////////////////////////////
#include "MCorsikaRead.h"

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

#include <TSystem.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"
#include "MStatusDisplay.h"

#include "MCorsikaRunHeader.h"
#include "MCorsikaEvtHeader.h"

#include "MPhotonEvent.h"

ClassImp(MCorsikaRead);

using namespace std;

/*  ----------- please don't delete and don't care about (Thomas) ------------
#define kBUFSZ 64 //1024*1024*64
#include <iomanip.h>
class bifstream : public istream, public streambuf
{
private:
    char fBuffer[kBUFSZ]; //!
    FILE *fd;

    int sync()
    {
        memset(fBuffer, 0, kBUFSZ);
        return 0; 
    }
    int underflow()
    {
        int sz=fread(fBuffer, kBUFSZ, 1, fd);
        //int sz=fread(fBuffer, 1, kBUFSZ, fd);
        setg(fBuffer, fBuffer, fBuffer+kBUFSZ);

        return sz==1 ? *(unsigned char*)fBuffer : EOF;//EOF;
        //return sz==kBUFSZ ? *(unsigned char*)fBuffer : EOF;//EOF;
    }
public:
    bifstream(const char *name) : istream(this)
    {
        fd = fopen(name, "rb");
        setbuf(fBuffer, kBUFSZ);
    }
};
*/

// --------------------------------------------------------------------------
//
// Default constructor. It tries to open the given file.
//
MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title)
    : fRunHeader(0), fEvtHeader(0), fEvent(0), /*fEvtData(0),*/ fForceMode(kFALSE),
    fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0),
    fIn(0), fParList(0)
{
    fName  = name  ? name  : "MRead";
    fTitle = title ? title : "Read task to read DAQ binary files";

    fFileNames = new TList;
    fFileNames->SetOwner();

    if (fname!=NULL)
        AddFile(fname);
}

// --------------------------------------------------------------------------
//
// Destructor. Delete input stream.
//
MCorsikaRead::~MCorsikaRead()
{
    delete fFileNames;
    if (fIn)
        delete fIn;
}

/*
Byte_t MCorsikaRead::IsFileValid(const char *name)
{
    MZlib fin(name);
    if (!fin)
        return 0;

    Byte_t c[4];
    fin.read((char*)c, 4);
    if (!fin)
        return 0;

    if (c[0]!=0xc0)
        return 0;

    if (c[1]==0xc0)
        return 1;

    if (c[1]==0xc1)
        return 2;

    return 0;
}
*/

// --------------------------------------------------------------------------
//
// Add a new file to a list of files to be processed, Returns the number
// of files added. (We can enhance this with a existance chack and
// wildcard support)
//
Int_t MCorsikaRead::AddFile(const char *fname, Int_t entries)
{
    TNamed *name = new TNamed(fname, "");
    fFileNames->AddLast(name);
    return 1;

}

Bool_t MCorsikaRead::ReadEvtEnd()
{
    if (!fRunHeader->SeekEvtEnd(*fIn))
    {
        *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl;
        if (!fForceMode)
            return kFALSE;
    }

    if (!fRunHeader->ReadEvtEnd(*fIn))
    {
        *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl;
        if (!fForceMode)
            return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// This opens the next file in the list and deletes its name from the list.
//
Int_t MCorsikaRead::OpenNextFile(Bool_t print)
{
    //
    // open the input stream and check if it is really open (file exists?)
    //
    if (fIn)
        delete fIn;
    fIn = NULL;

    //
    // Check for the existance of a next file to read
    //
    TObject *file = fFileNames->At(fNumFile);
    if (!file)
        return kFALSE;

    //
    // open the file which is the first one in the chain
    //
    const char *name = file->GetName();

    const char *expname = gSystem->ExpandPathName(name);
    fIn = new ifstream(expname);

    const Bool_t noexist = !(*fIn);
    if (noexist)
    {
        *fLog << err << "Cannot open file " << expname << ": ";
        *fLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl;
    }
    else
    {
        *fLog << inf << "Open file: '" << name << "'" << endl;

        if (fDisplay)
        {
            TString txt = GetFileName();
            txt += " @ ";
            txt += GetNumExecutions()-1;
            fDisplay->SetStatusLine2(txt);
        }
    }

    delete [] expname;

    if (noexist)
        return kERROR;

    fNumFile++;

    //
    // Read RUN HEADER (see specification) from input stream
    //
    if (!fRunHeader->ReadEvt(*fIn))
        return kERROR;
//    if (!fEvtHeader->ReadRunHeader(*fIn, *fRunHeader))
//        return kERROR;

    const streampos pos = fIn->tellg();
    if (!ReadEvtEnd())
        return kERROR;
    fIn->seekg(pos, ios::beg);


    fNumEvents += fRunHeader->GetNumEvents();
    fRunHeader->SetReadyToSave();

    //
    // Print Run Header
    //  We print it after the first event was read because
    //  we still miss information which is stored in the event header?!?
    if (print)
        fRunHeader->Print();

    if (!fParList)
        return kTRUE;

    //
    // Search for MTaskList
    //
    MTask *tlist = (MTask*)fParList->FindObject("MTaskList");
    if (!tlist)
    {
        *fLog << err << dbginf << "MTaskList not found... abort." << endl;
        return kERROR;
    }

    //
    // A new file has been opened and new headers have been read.
    //  --> ReInit tasklist
    //
    return tlist->ReInit(fParList) ? kTRUE : kERROR;
}

// --------------------------------------------------------------------------
//
// Return file name of current file.
//
TString MCorsikaRead::GetFullFileName() const
{
    const TObject *file = fFileNames->At(fNumFile-1);
    return file ? file->GetName() : "";
}

// --------------------------------------------------------------------------
//
// Restart with the first file
//
Bool_t MCorsikaRead::Rewind()
{
    fNumFile=0;
    fNumEvents=0;
    return OpenNextFile()==kTRUE;
}

Bool_t MCorsikaRead::CalcNumTotalEvents()
{
    fNumTotalEvents = 0;

    Bool_t rc = kTRUE;

    while (1)
    {
        switch (OpenNextFile(kFALSE))
        {
        case kFALSE:
            break;
        case kERROR:
            rc = kFALSE;
            break;
        case kTRUE:
            if (!ReadEvtEnd())
            {
                rc = kFALSE;
                break;
            }

            fNumTotalEvents += fRunHeader->GetNumEvents();
            continue;
        }
        break;
    }

    if (fIn)
        delete fIn;
    fIn = NULL;

    return rc;
}

// --------------------------------------------------------------------------
//
// The PreProcess of this task checks for the following containers in the
// list:
//   MCorsikaRunHeader <output>   if not found it is created
//   MCorsikaEvtHeader <output>   if not found it is created
//   MCorsikaEvtData <output>     if not found it is created
//   MCorsikaCrateArray <output>  if not found it is created
//   MCorsikaEvtTime <output>     if not found it is created (MTime)
//
// If all containers are found or created the run header is read from the
// binary file and printed.  If the Magic-Number (file identification)
// doesn't match we stop the eventloop.
//
// Now the EvtHeader and EvtData containers are initialized.
//
Int_t MCorsikaRead::PreProcess(MParList *pList)
{
    //
    // open the input stream
    // first of all check if opening the file in the constructor was
    // successfull
    //
    fParList = 0;

    if (!OpenStream())
        return kFALSE;

    //
    //  check if all necessary containers exist in the Parameter list.
    //  if not create one and add them to the list
    //
    fRunHeader = (MCorsikaRunHeader*)pList->FindCreateObj("MCorsikaRunHeader");
    if (!fRunHeader)
        return kFALSE;

    fEvtHeader = (MCorsikaEvtHeader*)pList->FindCreateObj("MCorsikaEvtHeader");
    if (!fEvtHeader)
        return kFALSE;

    fEvent = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent");
    if (!fEvent)
        return kFALSE;

    *fLog << inf << "Calculating number of total events..." << flush;
    if (!CalcNumTotalEvents())
        return kFALSE;
    *fLog << inf << " " << fNumTotalEvents << " found." << endl;

    fNumFile=0;
    fNumEvents=0;

    fParList = pList;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Read a single event from the stream
//
Bool_t MCorsikaRead::ReadEvent(istream &fin)
{
    //
    // Read in the next EVENT HEADER (see specification),
    // if there is no next event anymore stop eventloop
    //
    Int_t rc = fEvtHeader->ReadEvt(fin); //read event header block
    if (!rc)
        return kFALSE;

    rc = fEvent->ReadCorsikaEvt(fin);

    /*
    // Check if we are allowed to stop reading in the middle of the data
    if (rc==kFALSE && !fForce)
    {
        *fLog << err;
        *fLog << "ERROR - End of file in the middle of the data stream!" << endl;
        *fLog << "        Set the force-option to force processing." << endl;
        return kERROR;
    }
    */

    return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fin) : rc;
}

// --------------------------------------------------------------------------
//
// The Process reads one event from the binary file:
//  - The event header is read
//  - the run header is read
//  - all crate information is read
//  - the raw data information of one event is read
//
Int_t MCorsikaRead::Process()
{
    while (1)
    {
        if (fIn)
        {
            // Read a single event from file
            const Bool_t rc = ReadEvent(*fIn);

            // kFALSE means: end of file (try next one)
            if (rc!=kFALSE)
                return rc;

            if (!fRunHeader->ReadEvtEnd(*fIn))
                if (!fForceMode)
                    return kERROR;
        }

        //
        // If an event could not be read from file try to open new file
        //
        const Int_t rc = OpenNextFile();
        if (rc!=kTRUE)
            return rc;
    }
    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Close the file. Check whether the number of read events differs from
//  the number the file should containe (MCorsikaRunHeader). Prints a warning
//  if it doesn't match.
//
Int_t MCorsikaRead::PostProcess()
{
    //
    // Sanity check for the number of events
    //
    if (fNumEvents==GetNumExecutions()-1 || GetNumExecutions()==0)
        return kTRUE;

    *fLog << warn << dec;
    *fLog << "Warning - number of read events (" << GetNumExecutions()-1;
    *fLog << ") doesn't match number in run header(s) (";
    *fLog << fNumEvents << ")." << endl;

    return kTRUE;
}