source: trunk/MagicSoft/Mars/msim/MPhotonData.h@ 9249

#ifndef MARS_MPhotonData
#define MARS_MPhotonData

#ifndef MARS_MMcEvtBasic
#include "MMcEvtBasic.h"
#endif

#ifndef ROOT_TVector2
#include <TVector2.h>
#endif

#ifndef ROOT_TVector3
#include <TVector3.h>
#endif

#ifndef ROOT_TQuaternion
#include <TQuaternion.h>
#endif

// gcc 3.2
//class ifstream;
#include <iosfwd>
/*
class TH1;
class TH2;
class TH3;
*/
class MCorsikaRunHeader;

class MPhotonData : public TObject //MParContainer
{
private:
    Float_t fPosX;                       // [cm] X (north) at observation level
    Float_t fPosY;                       // [cm] Y (west) at observation level

    Float_t fCosU;                       // [cos x] U direction cosine to x-axis
    Float_t fCosV;                       // [cos y] V direction cosine to y-axis

    Float_t fTime;                       // [ns] Time since first interaction or entrance into atmosphere
    // 17M
    UShort_t fWavelength;                // [nm] Wavelength
    // 19M
    UInt_t   fNumPhotons;                // Number of cherenkov photons ins bunch
    Float_t  fProductionHeight;          // [cm] Height of bunch production
    MMcEvtBasic::ParticleId_t fPrimary;  // Type of emitting particle
    // 22M
    // gzip
    // 25M
    // raw
    // 32M

    Int_t   fTag;    //! A tag for external use
    Float_t fWeight; //! A weight for external use

protected:
    virtual Int_t FillCorsika(Float_t f[7]);
    virtual Int_t FillRfl(Float_t f[8]);

public:
    MPhotonData(/*const char *name=NULL, const char *title=NULL*/);
    //MPhotonData(const MPhotonData &ph);

    // Getter 1D
    Float_t  GetPosX()  const { return fPosX; }
    Float_t  GetPosY()  const { return fPosY; }

    Float_t  GetCosU()  const { return fCosU; }
    Float_t  GetCosV()  const { return fCosV; }
    Double_t GetCosW()  const { return TMath::Sqrt(1 - fCosU*fCosU - fCosV*fCosV); } // cos(Theta)
    Double_t GetTheta() const { return TMath::ACos(GetCosW()); }

    Double_t GetTime()  const { return fTime; }

    // Getter 2D
    TVector2 GetPos2()  const { return TVector2(fPosX, fPosY); }
    TVector2 GetDir2()  const { return TVector2(fCosU, fCosV);}
//    TVector2 GetDir2() const { return TVector2(fCosU, fCosV)/(1-TMath::Hypot(fCosU, fCosV)); }

    // Getter 3D
    TVector3 GetPos3()  const { return TVector3(fPosX, fPosY, 0); }
    TVector3 GetDir3()  const { return TVector3(fCosU, fCosV, -GetCosW()); }

    // Getter 4D                                      // cm   // ns
    TQuaternion GetPosQ() const { return TQuaternion(GetPos3(), fTime); }
    // FIXME: v in air!                                       // m/s  cm? / ns
    TQuaternion GetDirQ() const { return TQuaternion(GetDir3(), 1./(TMath::C()*100/1e9)); }

    // Getter Others
    UShort_t GetWavelength() const { return fWavelength; }
    MMcEvtBasic::ParticleId_t GetPrimary() const { return fPrimary; }

    //virtual Float_t GetWeight() const { return 1; }
    virtual Float_t GetWeight() const { return fWeight; }
    void SetWeight(Float_t w=1) { fWeight=w; }

    // Setter
    void SetPosition(Float_t x, Float_t y)  { fPosX=x; fPosY=y; }
    void SetDirection(Float_t u, Float_t v) { fCosU=u; fCosV=v; }

    void SetPosition(const TVector2 &p)     { fPosX=p.X(); fPosY=p.Y(); }
    void SetDirection(const TVector2 &d)    { fCosU=d.X(); fCosV=d.Y(); }

    void SetPosition(const TQuaternion &p)  { fPosX=p.fVectorPart.X(); fPosY=p.fVectorPart.Y(); fTime=p.fRealPart; }
    void SetDirection(const TQuaternion &d) { fCosU=d.fVectorPart.X(); fCosV=d.fVectorPart.Y(); }

    void SetPrimary(MMcEvtBasic::ParticleId_t p) { fPrimary=p; }
    void SetWavelength(UShort_t wl) { fWavelength=wl; }

    void AddTime(Double_t dt) { fTime += dt; }
    void SetTime(Double_t t)  { fTime  = t; }

    void Copy(TObject &obj) const;

    void SetTag(Int_t tag) { fTag=tag; }
    Int_t GetTag() const { return fTag; }

    //void Clear(Option_t * = NULL);
    void Print(Option_t * = NULL) const;
    //void Draw (Option_t * = NULL);
    Bool_t IsSortable() const { return kTRUE; }
    Int_t  Compare(const TObject *obj) const
    {
        const MPhotonData &d = *static_cast<const MPhotonData*>(obj);
        if (fTime<d.fTime)
            return -1;
        if (fTime>d.fTime)
            return 1;
        return 0;
    }
/*
    void FillRad(TH1 &hist, Float_t scale=1) const;
    void FillRad(TH2 &hist, Double_t x, Float_t scale=1) const;
    void Fill(TH2 &hist, Float_t scale=1) const;
    void Fill(TH3 &hist, Double_t z, Float_t scale=1) const;
 */
    Int_t ReadCorsikaEvt(istream &fin);
    Int_t ReadRflEvt(istream &fin);

    ClassDef(MPhotonData, 1) //Container to store a cherenkov photon bunch from a CORSUKA file
};
/*
class MPhotonDataWeighted : public MPhotonData
{
private:
    Float_t fWeight; // A weight for external use

public:
    Float_t GetWeight() const { return fWeight; }

    ClassDef(MPhotonData, 1)
};
*/
#endif