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

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

#include <math.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) :
    fAlpha(0), fTheta(0), fWidth(0), fLength(0), fSize(0), fDist(0), 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  = " << fabs(fAlpha)  << endl;
    *fLog << " - Width  = " << fWidth  << endl;
    *fLog << " - Length = " << fLength << endl;
    *fLog << " - Size   = " << fSize   << endl;
    *fLog << " - Dist   = " << fDist   << endl;
}

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

    fEllipse->Paint();
}

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

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

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

    delete fEllipse;

    fEllipse = NULL;
}

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

    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]

    return kTRUE;
}