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

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

/////////////////////////////////////////////////////////////////////////////
//
// 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)
    : fX(0), fY(0)
{
    fName  = name  ? name  : "MCorsikaEvtHeader";
    fTitle = title ? title : "Raw Event Header Information";
}

// --------------------------------------------------------------------------
//
// Return Impact (distance of ground incident point from telescope axis)
//
// Distance d between a point q and a line u (via p):
//   d = | ( vec(q) - vec(p) ) x vec(u) |  /  | vec(u) |
//
// If p = 0
//
//   ==>   d = | vec(q) x vec(u) |  /  | vec(u) |
//          w :=  q = (x/y/0)
//          v :=  u = (Alt/Az)
//
// For Alt/Az = Theta/Phi:
//
//  x = r sin(theta) cos(phi)
//  y = r sin(theta) sin(phi)
//  z = r cos(theta)
//
//                   (   -q3*u2    )   ( -cos(theta)*y                                     )
// vec(q) x vec(u) = (    q3*u1    ) = (  cos(theta)*x                                     )
//                   ( q1*u2-q2*u1 )   ( sin(theta) cos(phi) * y - sin(theta) sin(phi) * x )
//
//   ==> d = sqrt( cos(theta)^2 (x^2  + y^2 ) +
//                 sin(theta)^2 ( cos(phi) y + sin(phi) x)^2 )
//
Double_t MCorsikaEvtHeader::GetImpact() const
{
    const Double_t c = TMath::Cos(fZd);
    const Double_t s = TMath::Sin(fZd);

    const Double_t p = TMath::Cos(fAz)*fX - TMath::Sin(fAz)*fY;

    return TMath::Sqrt(c*c*(fX*fX + fY*fY) + s*s* p*p);
}

// --------------------------------------------------------------------------
//
//  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" << UTF8::kSquare << 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() << UTF8::kDeg << "/" << fAz*TMath::RadToDeg() << UTF8::kDeg << 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]; // FIXME: Could be negative

    // 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];  // Does this minus make sense?!

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