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

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

/////////////////////////////////////////////////////////////////////////////
//
// MCorsikaEvtHeader 
//    
/////////////////////////////////////////////////////////////////////////////
#include "MCorsikaEvtHeader.h"

#include <iomanip>
#include <fstream>

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

#include "MMcEvt.hxx"
#include "MCorsikaRunHeader.h"

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

// --------------------------------------------------------------------------
//
//  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, MCorsikaRunHeader &header)
{
    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];

    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[97];
    //fY = f[117];

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

    if (fEvtNumber==1)
    {
        header.fZdMin = f[79];
        header.fZdMax = f[80];
        header.fAzMin = 180-f[81];
        header.fAzMax = 180-f[82];

        header.fViewConeInnerAngle = f[151];
        header.fViewConeOuterAngle = f[152];

        header.Print();
    }

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

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

    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();
}