source: trunk/MagicSoft/Mars/manalysis/MHillas.cc@ 842

/* ======================================================================== *\
!
! *
! * 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): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
!   Author(s): Thomas Bretz  12/2000 (tbretz@uni-sw.gwdg.de)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MHillas                                                                 //
//                                                                         //
// Storage Container for the Hillas parameter                              //
//                                                                         //
// FIXME: Here everybody should find an explanation of the parameters      //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MHillas.h"

#include <math.h>
#include <fstream.h>

#include <TEllipse.h>

#include "MCerPhotEvt.h"
#include "MCerPhotPix.h"
#include "MGeomCam.h"

#include "MLog.h"

ClassImp(MHillas)

// --------------------------------------------------------------------------
//
// Default constructor.
//
MHillas::MHillas(const char *name, const char *title) : fEllipse(NULL)
{
    *fName  = name  ? name  : "MHillas";
    *fTitle = title ? title : "Storage container for Hillas parameter of one event";

    Reset();
    // FIXME: Initialization of values missing
}

// --------------------------------------------------------------------------
//
// Destructor. Deletes the TEllipse if one exists.
//
MHillas::~MHillas()
{
    Clear();
}

void MHillas::Reset()
{
    fAlpha  = 0;
    fTheta  = 0;
    fWidth  = 0;
    fLength = 0;
    fSize   = 0;
    fDist   = 0;

    Clear();
}

// --------------------------------------------------------------------------
//
// Print the hillas Parameters to *fLog
//
void MHillas::Print(Option_t *)
{
    *fLog << "Hillas Parameter:" << endl;
    *fLog << " - Alpha  = " << fabs(fAlpha)  << endl;
    *fLog << " - Width  = " << fWidth  << endl;
    *fLog << " - Length = " << fLength << endl;
    *fLog << " - Size   = " << fSize   << endl;
    *fLog << " - Dist   = " << fDist   << endl;
}

// --------------------------------------------------------------------------
//
// call the Paint function of the Ellipse if a TEllipse exists
//
void MHillas::Paint(Option_t *)
{
    if (!fEllipse)
        return;

    fEllipse->Paint();
}

// --------------------------------------------------------------------------
//
// Instead of adding MHillas itself to the Pad
// (s. AppendPad in TObject) we create an ellipse,
// which is added to the Pad by it's Draw function
// You can remove it by deleting the Ellipse Object
// (s. Clear() )
//
void MHillas::Draw(Option_t *)
{
    Clear();

    fEllipse = new TEllipse(cos(fTheta)*fDist, sin(fTheta)*fDist,
                            fLength, fWidth,
                            0, 360, fTheta*kRad2Deg+fAlpha-180);

    fEllipse->SetLineWidth(2);
    fEllipse->Draw();

    /*
     This is from TH1
     TString opt = option;
   opt.ToLower();
   if (gPad && !opt.Contains("same")) {
      //the following statement is necessary in case one attempts to draw
      //a temporary histogram already in the current pad
      if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
      gPad->Clear();
   }
   AppendPad(opt.Data());
   */
}

// --------------------------------------------------------------------------
//
// If a TEllipse object exists it is deleted
//
void MHillas::Clear(Option_t *)
{
    if (!fEllipse)
        return;

    delete fEllipse;

    fEllipse = NULL;
}

// --------------------------------------------------------------------------
//
// Calculate the Hillas parameters from a cerenkov photon event
// (The calcualtion is some kind of two dimentional statistics)
//
Bool_t MHillas::Calc(MGeomCam &geom, MCerPhotEvt &evt)
{
    const UInt_t nevt = evt.GetNumPixels();

    //
    // sanity check
    //
    if (nevt <= 2)
        return kFALSE;

    //
    // calculate mean valu of pixels
    //
    float xmean =0;
    float ymean =0;

    fSize = 0;

    //
    // FIXME! npix should be retrieved from MCerPhotEvt
    //
    UShort_t npix=0;
    for (UInt_t i=0; i<nevt; i++)
    {
        const MCerPhotPix &pix = evt[i];

        if (!pix.IsPixelUsed())
            continue;

        const MGeomPix &gpix = geom[pix.GetPixId()];

        const float nphot = pix.GetNumPhotons();

        fSize += nphot;
        xmean += nphot * gpix.GetX(); // [mm]
        ymean += nphot * gpix.GetY(); // [mm]

        npix++;
    }

    //
    // sanity check
    //
    if (fSize==0 || npix<=2)
        return kFALSE;

    xmean /= fSize; // [mm]
    ymean /= fSize; // [mm]

    //
    // calculate sdev
    //
    float sigmaxx=0;
    float sigmaxy=0;
    float sigmayy=0;

    for (UInt_t i=0; i<nevt; i++)
    {
        const MCerPhotPix &pix = evt[i];

        if (!pix.IsPixelUsed())
            continue;

        const MGeomPix &gpix = geom[pix.GetPixId()];

        const float dx = gpix.GetX() - xmean;
        const float dy = gpix.GetY() - ymean;

        const float nphot = pix.GetNumPhotons();

        sigmaxx += nphot * dx*dx; // [mm^2]
        sigmaxy += nphot * dx*dy; // [mm^2]
        sigmayy += nphot * dy*dy; // [mm^2]
    }

    //
    // check for orientation
    //
    const float theta = atan(sigmaxy/(sigmaxx-sigmayy)*2)/2;

    float c = cos(theta); // [1]
    float s = sin(theta); // [1]

    //
    // calculate the length of the two axis
    //
    float axis1 =  2.0*c*s*sigmaxy + c*c*sigmaxx + s*s*sigmayy; // [mm^2]
    float axis2 = -2.0*c*s*sigmaxy + s*s*sigmaxx + c*c*sigmayy; // [mm^2]

    axis1 /= fSize; // [mm^2]
    axis2 /= fSize; // [mm^2]

    //
    // check for numerical negatives
    // (very small number can get negative by chance)
    //
    if (axis1 < 0) axis1=0;
    if (axis2 < 0) axis2=0;

    //
    // calculate the main Hillas parameters
    //
    // fLength, fWidth describes the two axis of the ellipse
    // fAlpha is the angle between the length-axis and the center
    //    of the camera
    // fDist is the distance between the center of the camera and the
    //    denter of the ellipse
    //
    const int rotation = axis1<axis2;

    fLength = rotation ? sqrt(axis2) : sqrt(axis1);  // [mm]
    fWidth  = rotation ? sqrt(axis1) : sqrt(axis2);  // [mm]

    const float a = c*xmean + s*ymean;
    const float b = c*ymean - s*xmean;

    fAlpha  = rotation ? atan(a/b) : atan(-b/a);     // [rad]
    fAlpha *= kRad2Deg;                              // [deg]

    fDist   = sqrt(xmean*xmean + ymean*ymean);       // [mm]

    fTheta  = atan(ymean/xmean);                     // [rad]
    if (xmean<0) fTheta += kPI;                      // [rad]

    SetHasChanged();

    return kTRUE;
}

void MHillas::AsciiRead(ifstream &fin)
{
    fin >> fAlpha;
    fin >> fTheta;
    fin >> fWidth;
    fin >> fLength;
    fin >> fSize;
    fin >> fDist;
}

void MHillas::AsciiWrite(ofstream &fout) const
{
    fout << fAlpha << " ";
    fout << fTheta << " ";
    fout << fWidth << " ";
    fout << fLength << " ";
    fout << fSize << " ";
    fout << fDist << endl;
}