source: trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.cc@ 9341

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

/////////////////////////////////////////////////////////////////////////////
//
// MCorsikaEvtHeader 
//
// Class Version 2:
// ----------------
//  - UInt_t fParticleID
//
//
/////////////////////////////////////////////////////////////////////////////
#include "MCorsikaEvtHeader.h"

#include <iomanip>
#include <fstream>

#include <TMath.h>

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

//#include "MMcEvt.hxx"

ClassImp(MCorsikaEvtHeader);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. Create the array to store the data.
//
MCorsikaEvtHeader::MCorsikaEvtHeader(const char *name, const char *title)
{
    fName  = name  ? name  : "MCorsikaEvtHeader";
    fTitle = title ? title : "Raw Event Header Information";
}

// --------------------------------------------------------------------------
//
//  Return Hypot(x, y)
//
Double_t MCorsikaEvtHeader::GetImpact() const
{
    return TMath::Hypot(fX, fY);
}

// --------------------------------------------------------------------------
//
//  This member function prints all Data of one Event to *fLog.
//
void MCorsikaEvtHeader::Print(Option_t *o) const
{
    *fLog << all;
    *fLog << "Event Number:              " << dec << fEvtNumber << endl;
//    *fLog << "Particle ID:               " << MMcEvt::GetParticleName(fParticleID) << endl;
    *fLog << "Energy:                    " << fTotalEnergy << "GeV" << endl;
    *fLog << "Starting Altitude:         " << fStartAltitude << "g/cm²" << endl;
    *fLog << "Number of 1st Target:      " << fFirstTargetNum << endl,
    *fLog << "Height of 1st Interaction: " << fFirstInteractionHeight/100. << "m" << endl;
    *fLog << "Momentum X/Y/Z (GeV/c):    " << fMomentumX << "/" << fMomentumY << "/" << fMomentumZ << endl;
    *fLog << "Zenith/Azimuth Angle:      " << fZd*TMath::RadToDeg() << "°/" << fAz*TMath::RadToDeg() << "°" << endl;
    *fLog << "Impact X/Y:                " << fX/100. << "m/" << fY/100. << "m  (r=" << TMath::Hypot(fX, fY)/100. << "m)" << endl;
    *fLog << "Weighted Num Photons:      " << fWeightedNumPhotons << endl;
    *fLog << endl;
}

// --------------------------------------------------------------------------
//
// read the EVTH information from the input stream
// return FALSE if there is no  header anymore, else TRUE
//
Int_t MCorsikaEvtHeader::ReadEvt(std::istream &fin)
{
    char evth[4];
    fin.read(evth, 4);
    if (memcmp(evth, "EVTH", 4))
    {
        fin.seekg(-4, ios::cur);
        return kFALSE;
    }

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

    fEvtNumber  = TMath::Nint(f[0]);
//    fParticleID = TMath::Nint(f[1]);

    fTotalEnergy            = f[2];
    fStartAltitude          = f[3];
    fFirstTargetNum         = f[4];
    fFirstInteractionHeight = f[5];

    // Pointing opposite particle direction
    // (x=north, y=west, z=upwards)
    //fMomentumX = f[6];
    //fMomentumY = f[7];
    //fMomentumZ = f[8];
    // Pointing along particle direction
    // (x=east, y=north, z=upwards)
    fMomentumX =  f[7];
    fMomentumY = -f[6];
    fMomentumZ = -f[8];

    // FIXME: Correct for direction of magnetic field!
    fZd        = f[9];
    fAz        = TMath::Pi()-f[10];

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

    fX =  f[117];   //fX = f[97];
    fY = -f[97];    //fY = f[117];

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

    return !fin.eof();
}


// this member function is for reading the event end block

Bool_t MCorsikaEvtHeader::ReadEvtEnd(std::istream &fin)
{
    //fin.seekg(-1088,ios::cur);

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

    const UInt_t evtnum = TMath::Nint(f[0]);
    if (evtnum!=fEvtNumber)
    {
        *fLog << err << "ERROR - Mismatch in stream: Event number in EVTE (";
        *fLog << evtnum << ") doesn't match EVTH (" << fEvtNumber << ")." << endl;
        return kERROR;
    }

    fWeightedNumPhotons = f[1];

    fin.seekg(1080,ios::cur);

    return !fin.eof();
}