source: trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc@ 9213

/* ======================================================================== *\
!
! *
! * 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 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
               Qi Zhe,      06/2007 <mailto:qizhe@astro.uni-wuerzburg.de>

!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MCorsikaRunHeader
//
// Root storage container for the RUN HEADER information
//
//
//  RAW DATA FORMAT VERSION
//  =======================
//
//  Format Version 8:
//  -----------------
//   + fNumBytesPerSample;
//   + fFreqSampling;
//   + fNumSignificantBits;
//   * changes in MRawCrateHeader
//
//  Format Version 7:
//  -----------------
//   - unused
//
//  Format Version 6:
//  -----------------
//   + added CameraVersion
//   + added TelescopeNumber
//   + added ObservationMode
//   + added dummies for TelescopeRa/Dec
//
//  Format Version 5:
//  -----------------
//   - now the sub millisecond information of the time is valid and decoded
//     which enhances the precision from 51.2us to 200ns
//
//  Format Version 4:
//  -----------------
//   - added support for pixels with negative IDs
//
//  Format Version 3:
//  -----------------
//   - ???
//
//  Format Version 2:
//  -----------------
//   - removed mjd from data
//   - added start time
//   - added stop  time
//
//
//  MCorsikaRunHeader CLASS VERSION
//  ===========================
//
//  Format Version 6:
//  -----------------
//   - added fNumBytesPerSample;
//   - added fFreqSampling;
//   - added fNumSignificantBits;
//
//  Class Version 5:
//  -----------------
//   - for compatibility with newer camera versions
//
//  Class Version 4:
//  -----------------
//   - added fCameraVersion
//   - added fTelescopeNumber
//   - changed length of fProjectName to 101
//   - changed length of fSourceName  to 81
//
//  Class Version 3:
//  ----------------
//   - enhanced SourceName and ProjectName by one character, because
//     without telling us the guranteed trailing \0-character has
//     skipped
//
//  Class Version 2:
//  ----------------
//   - removed fMJD, fYear, fMonth, fDay
//   - added fRunStart
//   - added fRunStop
// 
//  Class Version 1:
//  ----------------
//   - first implementation
//
////////////////////////////////////////////////////////////////////////////

#include "MCorsikaRunHeader.h"

#include <fstream>
#include <iomanip>

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

ClassImp(MCorsikaRunHeader);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. Creates array which stores the pixel assignment.
//
//
MCorsikaRunHeader::MCorsikaRunHeader(const char *name, const char *title)
    : fNumObsLevel(0),
    fZdMin(0), fZdMax(-1), fAzMin(0), fAzMax(0),
fViewConeInnerAngle(0), fViewConeOuterAngle(-1)
{
    fName  = name  ? name  : "MCorsikaRunHeader";
    fTitle = title ? title : "Raw Run Header Information";
}

// --------------------------------------------------------------------------
//
// Read in one run header from the binary file
//
Bool_t MCorsikaRunHeader::ReadEvt(istream& fin)
{
    char runh[4];
    fin.read(runh, 4);
    if (memcmp(runh, "RUNH", 4))
    {
        *fLog << err << "ERROR - Wrong identifier: RUNH expected." << endl;
        return kFALSE;
    }

    Float_t f[68];
    fin.read((char*)f, 68);

    fRunNumber = TMath::Nint(f[0]);
    fNumEvents = 0;

    fRunStart.SetCorsikaTime(f[1]);

    fProgramVersion = f[2];          //FIXME: INT???
    fNumObsLevel    = TMath::Nint(f[3]);

    memset(fObsLevel, 0, 10*4);
    memcpy(fObsLevel, f+4, fNumObsLevel*4);

    fSlopeSpectrum  = f[14];
    fEnergyMin      = f[15];
    fEnergyMax      = f[16];

    fin.seekg(1020, ios::cur);     // skip the remaining data of this block

    // -------------------- Read first event header -------------------

    // f[76] Cherenkov flag:
    //        bit(1) : CERENKOV option compiled in
    //        bit(2) : IACT option compiled in
    //        bit(3) : CEFFIC option compiled in
    //        bit(4) : ATMEXT option compiled in
    //        bit(5) : ATMEXT option used with refaction enabled
    //        bit(6) : VOLUMEDET option compiled in
    //        bit(7) : CURVED option compiled in
    //        bit(9) : SLATN option compiled in
    //        11-21  : table number for externam athmosphere (but<1024)
    // f[78] Curved athmosphere? (0=flat, 1=curved)

    // f[80] lower edge of theta in °
    // f[81] upper edge of theta in °
    // f[82] lower edge of phi   in °
    // f[83] upper edge of phi   in °
    // f[84] cherenkov bunch size

    // f[93] flag for additinal muon information of particle output file

    // f[95] Cherenkov bandwidth lower end in nm
    // f[96] Cherenkov bandwidth upper end in nm

    // f[97] Numbr i of uses of each cherenkov event
    // f[145] Muon multiple scattering flag
    // f[152] !! inner angle of view cone (°)
    // f[153] !! outer angle of view cone (°)
    // f[156] altitude of horizontal shower axis

    char evth[4];
    fin.read(evth, 4);
    if (memcmp(evth, "EVTH", 4))
    {
        *fLog << err << "ERROR - Wrong identifier: EVTH expected." << endl;
        return kFALSE;
    }

    Float_t g[273];
    fin.read((char*)&g, 273*4);
    if (fin.eof())
        return kFALSE;

    fin.seekg(-274*4, ios::cur);

    const Int_t n = TMath::Nint(g[96]);
    if (n!=1)
    {
        *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl;
        return kFALSE;
    }

    //fImpactMax = g[86];

    fZdMin = g[79];
    fZdMax = g[80];
    fAzMin = 180-g[81];
    fAzMax = 180-g[82];

    fWavelengthMin = g[95];
    fWavelengthMax = g[96];

    fViewConeInnerAngle = g[151];
    fViewConeOuterAngle = g[152];

    return kTRUE;;
}

Bool_t MCorsikaRunHeader::ReadEvtEnd(istream& fin)
{
    char runh[4];
    fin.read(runh, 4);
    if (memcmp(runh, "RUNE", 4))
    {
        *fLog << err << "ERROR - Wrong identifier: RUNE expected." << endl;
        return kFALSE;
    }

    Float_t f[2];
    fin.read((char*)f, 2*4);

    const UInt_t runnum = TMath::Nint(f[0]);
    if (runnum!=fRunNumber)
    {
        *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
        *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
        return kFALSE;
    }

    fNumEvents = TMath::Nint(f[1]);

    fin.seekg(270*4, ios::cur);     // skip the remaining data of this block

    return kTRUE;
}

Bool_t MCorsikaRunHeader::SeekEvtEnd(istream &fin)
{
    // 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 kFALSE;
}

// --------------------------------------------------------------------------
//
// print run header information on *fLog. The option 'header' supresses
// the pixel index translation table.
//
void MCorsikaRunHeader::Print(Option_t *t) const
{
    *fLog << all << endl;
    *fLog << "Run Number:  " << fRunNumber << "  (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl;
    if (fNumEvents>0)
        *fLog << "Num Events:  " << fNumEvents << endl;
    *fLog << "Obs Level:  ";
    for (Byte_t i=0; i<fNumObsLevel; i++)
        *fLog << " " << fObsLevel[i]/100. << "m";
    *fLog << endl;
    *fLog << "Spectrum:    Slope=" << fSlopeSpectrum << "  (" << fEnergyMin << "GeV-" << fEnergyMax << "GeV)" <<  endl;

    if (fViewConeOuterAngle>0)
        *fLog << "ViewCone:    " << fViewConeInnerAngle << "deg-" << fViewConeOuterAngle << "deg" << endl;

    if (fZdMax>=0)
    {
        *fLog << "Zd/Az:       " << fZdMin << "deg";
        if (fZdMin==fZdMax)
            *fLog << " (fixed)";
        else
            *fLog << "-" << fZdMax << "deg";
        *fLog << " / " << fAzMin << "deg";
        if (fAzMin==fAzMax)
            *fLog << " (fixed)";
        else
            *fLog << "-" << fAzMax << "deg";
        *fLog << endl;
    }
}