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

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz    12/2000 <mailto:tbretz@uni-sw.gwdg.de>
!   Author(s): Harald Kornmayer 1/2001
!   Author(s): Rudolf Bock     10/2001 <mailto:Rudolf.Bock@cern.ch>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHillas
//
// Storage Container for image parameters
//
//    basic image parameters
// fLength   major axis of ellipse
// fWidth    minor axis
// fDelta    angle of major axis wrt x-axis
// fSize     total sum of pixels
// fMeanX    x of center of ellipse
// fMeanY    y of center of ellipse
//
/////////////////////////////////////////////////////////////////////////////
#include "MHillas.h"

#include <fstream.h>

#include <TEllipse.h>

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

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

#include "MLog.h"
#include "MLogManip.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 image parameters of one event";

    Reset();
    // FIXME: (intelligent) initialization of values missing

    fEllipse = new TEllipse;
}

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

// --------------------------------------------------------------------------
//
void MHillas::Reset()
{
    fLength = 0;
    fWidth  = 0;
    fDelta  = 0;
    fSize   = 0;
    fMeanX  = 0;
    fMeanY  = 0;

    Clear();
}

// --------------------------------------------------------------------------
//
// Print the hillas Parameters to *fLog
//
void MHillas::Print(Option_t *) const
{
    Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg;

    if (atg<0)
        atg += 180;

    *fLog << all;
    *fLog << "Basic Image Parameters (" << GetName() << ")" << endl;
    *fLog << " - Length   [mm]  = " << fLength << endl;
    *fLog << " - Width    [mm]  = " << fWidth  << endl;
    *fLog << " - Meanx    [mm]  = " << fMeanX  << endl;
    *fLog << " - Meany    [mm]  = " << fMeanY  << endl;
    *fLog << " - Delta    [deg] = " << fDelta*kRad2Deg << endl;
    *fLog << " - atg(y/x) [deg] = " << atg     << endl;
    *fLog << " - Size     [1]   = " << fSize   << " #CherPhot"   << endl;
}

/*
// -----------------------------------------------------------
//
// call the Paint function of the Ellipse if a TEllipse exists
//
void MHillas::Paint(Option_t *)
{
     fEllipse->SetLineWidth(2);
     fEllipse->PaintEllipse(fMeanX, fMeanY, fLength, fWidth,
                            0, 360, fDelta*kRad2Deg+180);
}
*/

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

    Clear();

    fEllipse = new TEllipse(fMeanX, fMeanY, fLength, fWidth,
                            0, 360, fDelta*kRad2Deg+180);

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

    /*
     fEllipse->SetPhimin();
     fEllipse->SetPhimax();
     fEllipse->SetR1(fLength);
     fEllipse->SetR2(fWidth);
     fEllipse->SetTheta(fDelta*kRad2Deg+180);
     fEllipse->SetX1(fMeanX);
     fEllipse->SetY1(fMeanY);

     fEllipse->SetLineWidth(2);
     fEllipse->PaintEllipse(fMeanX, fMeanY, fLength, fWidth,
                            0, 360, fDelta*kRad2Deg+180);

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

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

    delete fEllipse;

    fEllipse = NULL;
}


// --------------------------------------------------------------------------
//
// Calculate the image parameters from a Cherenkov photon event
// assuming Cher.photons/pixel and pixel coordinates are given
//
Bool_t MHillas::Calc(const MGeomCam &geom, const MCerPhotEvt &evt)
{
    // 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)
    //

    const UInt_t npixevt = evt.GetNumPixels();

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

    //
    // calculate mean value of pixel coordinates and fSize
    // -----------------------------------------------------
    //
    fMeanX = 0;
    fMeanY = 0;
    fSize  = 0;

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

        if (!pix.IsPixelUsed())
            continue;

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

        const float nphot = pix.GetNumPhotons();

        fSize  += nphot;                             // [counter]
        fMeanX += nphot * gpix.GetX();               // [mm]
        fMeanY += nphot * gpix.GetY();               // [mm]

        npix++;
    }

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

    fMeanX /= fSize;                                 // [mm]
    fMeanY /= fSize;                                 // [mm]

    //
    // calculate 2nd moments
    // -------------------
    //
    float corrxx=0;                                  // [m^2]
    float corrxy=0;                                  // [m^2]
    float corryy=0;                                  // [m^2]

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

        if (!pix.IsPixelUsed())
            continue;

        const MGeomPix &gpix = geom[pix.GetPixId()];
        const float dx = gpix.GetX() - fMeanX;       // [mm]
        const float dy = gpix.GetY() - fMeanY;       // [mm]

        const float nphot = pix.GetNumPhotons();     // [#phot]

        corrxx += nphot * dx*dx;                     // [mm^2]
        corrxy += nphot * dx*dy;                     // [mm^2]
        corryy += nphot * dy*dy;                     // [mm^2]
    }

    //
    // calculate the basic Hillas parameters: orientation and size of axes
    // -------------------------------------------------------------------
    //
    const float d = (corryy - corrxx)/(corrxy*2);

    fDelta = atan(d + sqrt(d*d + 1));

    fCosDelta = cos(fDelta);   // need these in derived classes
    fSinDelta = sin(fDelta);   // like MHillasExt

    float axis1 = ( fCosDelta*fSinDelta*corrxy*2 + fCosDelta*fCosDelta*corrxx + fSinDelta*fSinDelta*corryy) / fSize; // [mm^2]
    float axis2 = (-fCosDelta*fSinDelta*corrxy*2 + fSinDelta*fSinDelta*corrxx + fCosDelta*fCosDelta*corryy) / fSize; // [mm^2]
 
    // very small numbers can get negative by rounding
    if (axis1 < 0) axis1=0;
    if (axis2 < 0) axis2=0; 

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

    SetReadyToSave();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
void MHillas::AsciiRead(ifstream &fin)
{
    fin >> fLength;
    fin >> fWidth;
    fin >> fDelta;
    fin >> fSize;
    fin >> fMeanX;
    fin >> fMeanY;
}

// --------------------------------------------------------------------------
/*
void MHillas::AsciiWrite(ofstream &fout) const
{
    fout << fLength << " ";
    fout << fWidth  << " ";
    fout << fDelta  << " ";
    fout << fSize   << " ";
    fout << fMeanX  << " ";
    fout << fMeanY;
}
*/