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

/* ======================================================================== *\
!
! *
! * 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>
!   Author(s): Wolfgang Wittek  6/2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHillas
//
// Storage Container for image parameters
//
//    basic image parameters
//
// Version 1:
// ----------
// fLength   [mm]       major axis of ellipse
// fWidth    [mm]       minor axis
// fDelta    [rad]      angle of major axis with x-axis
//                      by definition the major axis is pointing into
//                      the hemisphere x>0, thus -pi/2 < delta < pi/2
// fSize     [#CerPhot] total sum of pixels
// fMeanX    [mm]       x of center of ellipse
// fMeanY    [mm]       y of center of ellipse
//
// Version 2:
// ----------
// fNumCorePixels  number of pixels called core
// fNumUsedPixels  number of pixels which survived the cleaning
//
/////////////////////////////////////////////////////////////////////////////
#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();
}

// --------------------------------------------------------------------------
//
// Initializes the values with defaults. For the default values see the
// source code.
//
void MHillas::Reset()
{
    fLength = -1;
    fWidth  = -1;
    fDelta  =  0;

    fSize   = -1;
    fMeanX  = -1;
    fMeanY  = -1;

    fNumUsedPixels = -1;
    fNumCorePixels = -1;

    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 << " - Delta    [deg] = " << fDelta*kRad2Deg << endl;
    *fLog << " - Size     [1]   = " << fSize   << " #CherPhot"   << endl;
    *fLog << " - Meanx    [mm]  = " << fMeanX  << endl;
    *fLog << " - Meany    [mm]  = " << fMeanY  << endl;
    *fLog << " - atg(y/x) [deg] = " << atg     << 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
// In case you don't call Calc from within an eventloop make sure, that
// you call the Reset member function before.
//
Bool_t MHillas::Calc(const MGeomCam &geom, const MCerPhotEvt &evt)
{
    const UInt_t npixevt = evt.GetNumPixels();

    //
    // sanity check
    //
    if (npixevt <= 2)
    {
        *fLog << warn << "MHillas::Calc: Event has less than two pixels... skipped." << endl;
        return kFALSE;
    }

    //
    // calculate mean value of pixel coordinates and fSize
    // -----------------------------------------------------
    //
    // Because this are only simple sums of roughly 600 values
    // with an accuracy less than three digits there should not
    // be any need to use double precision (Double_t) for the
    // calculation of fMeanX, fMeanY and fSize.
    //
    fMeanX = 0;
    fMeanY = 0;
    fSize  = 0;

    fNumUsedPixels = 0;
    fNumCorePixels = 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_t nphot = pix.GetNumPhotons();

        if (nphot==0)
            *fLog << warn << GetDescriptor() << ": Pixel #" << pix.GetPixId() << " has no photons." << endl;

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

        if (pix.IsPixelCore())
            fNumCorePixels++;

        fNumUsedPixels++;
    }

    //
    // sanity checks
    //
    if (fSize==0)
    {
        *fLog << warn << GetDescriptor() << ": Event has zero cerenkov photons... skipped." << endl;
        return kFALSE;
    }

    if (fNumUsedPixels<3)
    {
        *fLog << warn << GetDescriptor() << ": Event has less than 3 used pixels... skipped." << endl;
        return kFALSE;
    }

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

    //
    // calculate 2nd moments
    // ---------------------
    //
    // To make sure that the more complicated sum is evaluated as
    // accurate as possible (because it is needed for more
    // complicated calculations (see below) we calculate it using
    // double prcision.
    //
    Double_t corrxx=0;                               // [m^2]
    Double_t corrxy=0;                               // [m^2]
    Double_t 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_t dx = gpix.GetX() - fMeanX;     // [mm]
        const Float_t dy = gpix.GetY() - fMeanY;     // [mm]

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

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

    //
    // If corrxy=0 (which should never happen, because fSize>0) we
    // cannot calculate Length and Width. The calculation failed
    // and returnes kFALSE
    //
    if (corrxy==0)
    {
        *fLog << warn << GetDescriptor() << ": Event has CorrXY==0... skipped." << endl;
        return kFALSE;
    }

    //
    // calculate the basic Hillas parameters: orientation and size of axes
    // -------------------------------------------------------------------
    //
    // fDelta is the angle between the major axis of the ellipse and
    //  the x axis. It is independent of the position of the ellipse
    //  in the camera it has values between -pi/2 and pi/2 degrees
    //
    const Double_t d0   = corryy - corrxx;
    const Double_t d1   = corrxy*2;
    const Double_t d2   = d0 + sqrt(d0*d0 + d1*d1);
    const Double_t tand = d2 / d1;

    fDelta = atan(tand);

    const Double_t s2 = tand*tand+1;
    const Double_t s  = sqrt(s2);

    fCosDelta =  1.0/s;   // need these in derived classes
    fSinDelta = tand/s;   // like MHillasExt

    Double_t axis1 = (tand*tand*corryy + d2 + corrxx)/s2/fSize;
    Double_t axis2 = (tand*tand*corrxx - d2 + corryy)/s2/fSize;

    //
    // fLength^2 is the second moment along the major axis of the ellipse
    // fWidth^2  is the second moment along the minor axis of the ellipse
    //
    // From the algorithm we get: fWidth <= fLength is always true
    //
    // very small numbers can get negative by rounding
    //
    fLength = axis1<0 ? 0 : sqrt(axis1);  // [mm]
    fWidth  = axis2<0 ? 0 : 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;
}
*/