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

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

#include <TEllipse.h>

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

#include "MLog.h"

ClassImp(MHillas)

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

    // FIXME: Initialization of values missing
}

MHillas::~MHillas()
{
    Clear();
}

void MHillas::Print(Option_t *)
{
    *fLog << "Hillas Parameter:" << endl;
    *fLog << " - Alpha  = " << fAlpha  << endl;
    *fLog << " - Width  = " << fWidth  << endl;
    *fLog << " - Length = " << fLength << endl;
    *fLog << " - Dist   = " << fDist   << endl;
}

void MHillas::Draw(Option_t *)
{
    //
    // 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() )
    //

    Clear();

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

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

void MHillas::Clear(Option_t *)
{
    if (!fEllipse)
        return;

    delete fEllipse;

    fEllipse = NULL;
}

void MHillas::Calc(MGeomCam &geom, MCerPhotEvt &evt)
{
    const UInt_t nevt = evt.GetNumPixels();

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

    fSize = 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]
    }

    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]

    //
    // resolve four-fold ambiguity of solution
    //
    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
    //
    if (axis1 < 0) axis1=0;
    if (axis2 < 0) axis2=0;

    //
    // check the rotation of the axis (maybe turn by 90ø)
    //
    const int   rotation = axis1<axis2;

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

    fAlpha  = rotation ?
        fabs(atan((-xmean*c - ymean*s)/(c*ymean - s*xmean))) :
        fabs(atan(( ymean*c - xmean*s)/(c*xmean + s*ymean))) ;

        // [deg]

    fAlpha *= kRad2Deg;                                  // [deg]

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

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