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

/* ======================================================================== *\
!
! *
! * 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    //
//        - using boooleans for fIsPixelUsed, fIsPixelCore, ... is rather  //
//          slow because you have to loop over all pixels in any loop.     //
//          There could be a huge speed improvement using Hash Tables      //
//          (linked lists, see THashTable and THashList, too)              //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MHillas.h"

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

#include <TEllipse.h>

#include "MGeomPix.h"
#include "MGeomCam.h"

#include "MCerPhotPix.h"
#include "MCerPhotEvt.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: (intelligent) 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 *) const
{
    *fLog << "Hillas Parameter:" << endl;
    *fLog << " - Alpha  = " << fabs(fAlpha) << "°"    << endl;
    *fLog << " - Width  = " << fWidth  << " mm"       << endl;
    *fLog << " - Length = " << fLength << " mm"       << endl;
    *fLog << " - Size   = " << fSize   << " #CerPhot" << endl;
    *fLog << " - Dist   = " << fDist   << " mm"       << 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 *opt)
{
    Clear();

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

    fEllipse->SetLineWidth(2);
    fEllipse->Draw();
    //AppendPad(opt);

    /*
     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)
//
//   FIXME: MHillas::Calc is rather slow at the moment because it loops
//          unnecessarily over all pixels in all its loops (s.MImgCleanStd)
//          The speed could be improved very much by using Hash-Tables
//          (linked lists, see THashTable and THashList, too)
//
Bool_t MHillas::Calc(const MGeomCam &geom, const 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]

    SetReadyToSave();

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