source: trunk/MagicSoft/Mars/mimage/MHillas.cc@ 7044

/* ======================================================================== *\
!
! *
! * 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-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// 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
//
// Version 3:
// ----------
// fNumCorePixels  moved to MNewImagePar
// fNumUsedPixels  moved to MNewImagePar
// fCosDelta       added
// fSinDelte       added
//
/////////////////////////////////////////////////////////////////////////////
#include "MHillas.h"

#include <TArrayF.h>

#include <TLine.h>
#include <TEllipse.h>
#include <TVector2.h>

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

#include "MSignalPix.h"
#include "MSignalCam.h"

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MHillas);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MHillas::MHillas(const char *name, const char *title)
{
    fName  = name  ? name  : "MHillas";
    fTitle = title ? title : "Storage container for image parameters of one event";

    Reset();
}

// --------------------------------------------------------------------------
//
// 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  =  0;
    fMeanY  =  0;
}

// --------------------------------------------------------------------------
//
//  return the mean position as TVector2(meanx, meany)
//
TVector2 MHillas::GetMean() const
{
    return TVector2(fMeanX, fMeanY);
}

// --------------------------------------------------------------------------
//
//  Analytical estimation of length of border line:
//    U = pi*(a+b - sqrt(a*b))
//
//  GetBorderLine() [mm]
//
Double_t MHillas::GetBorderLine() const
{
    const Double_t a = fWidth+fLength;
    const Double_t s = fWidth*fLength;

    return TMath::Pi()*(1.5*a - TMath::Sqrt(s));
}

// --------------------------------------------------------------------------
//
//  Analytical estimation of length of border line:
//    A = pi*a*b
//
//  GetArea() [mm^2]
//
Double_t MHillas::GetArea() const
{
    return TMath::Pi()*fWidth*fLength;
}

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

// --------------------------------------------------------------------------
//
// Print the hillas Parameters to *fLog, depending on the geometry in
// units of deg
//
void MHillas::Print(const MGeomCam &geom) const
{
    Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg;

    if (atg<0)
        atg += 180;

    *fLog << all;
    *fLog << GetDescriptor() << endl;
    *fLog << " - Length         [deg] = " << fLength*geom.GetConvMm2Deg() << endl;
    *fLog << " - Width          [deg] = " << fWidth *geom.GetConvMm2Deg() << endl;
    *fLog << " - Delta          [deg] = " << fDelta*kRad2Deg << endl;
    *fLog << " - Size           [1]   = " << fSize   << " #CherPhot"   << endl;
    *fLog << " - Meanx          [deg] = " << fMeanX *geom.GetConvMm2Deg() << endl;
    *fLog << " - Meany          [deg] = " << fMeanY *geom.GetConvMm2Deg() << endl;
    *fLog << " - atg(y/x)       [deg] = " << atg     << endl;
}

// --------------------------------------------------------------------------
//
// Paint the ellipse corresponding to the parameters
//
void MHillas::Paint(Option_t *opt)
{
    static const Float_t sOffsetW=40.0;
    static const Float_t sOffsetL=266.0;

    if (fLength<0 || fWidth<0)
        return;

    const Float_t l = fLength+sOffsetL;
    const Float_t w = fWidth +sOffsetW;

    TLine line;
    line.SetLineWidth(1);
    line.SetLineColor(kRed);

    // Length line
    line.PaintLine(-l*fCosDelta+fMeanX, -l*fSinDelta+fMeanY,
                   +l*fCosDelta+fMeanX, +l*fSinDelta+fMeanY);


    // Width line
    line.PaintLine(+w*fSinDelta+fMeanX, -w*fCosDelta+fMeanY,
                   -w*fSinDelta+fMeanX, +w*fCosDelta+fMeanY);

    TEllipse e(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta*kRad2Deg+180);
    e.SetLineWidth(2);
    e.SetLineColor(kGreen);
    e.Paint();
}

// --------------------------------------------------------------------------
//
// 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.
// Returns:
//   0  no error
//   1  number of pixels in event == 0 (special MC events)
//   2  size==0
//   3  number of used pixel < 3
//   4  CorrXY == 0
//
Int_t MHillas::Calc(const MGeomCam &geom, const MSignalCam &evt, Int_t island)
{
    const UInt_t numpix  = evt.GetNumPixels();

    //
    // sanity check 1 (special MC events)
    //
    if (numpix==0)
        return 1;

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

    UInt_t numused = 0;
/*
    MSignalPix *pix = 0;

    TIter Next(evt);
    while ((pix=(MSignalPix*)Next()))*/
    for (UInt_t i=0; i<numpix; i++)
    {
        MSignalPix &pix = evt[i];
        if (!pix.IsPixelUsed())
            continue;

        if (island>=0 && pix.GetIdxIsland()!=island)
            continue;

        const MGeomPix &gpix = geom[i/*pix->GetPixId()*/];

        const Float_t nphot = pix.GetNumPhotons();

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

        numused++;
    }

    //
    // sanity check 2
    //
    if (fSize==0)
        return 2;

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

    //
    // sanity check 3
    //
    if (numused<3)
        return 3;

    //
    // 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]

    //Next.Reset();
    //while ((pix=(MSignalPix*)Next()))
    //{
    for (UInt_t i=0; i<numpix; i++)
    {
        const MSignalPix &pix = evt[i];
        if (!pix.IsPixelUsed())
            continue;

        if (island>=0 && pix.GetIdxIsland()!=island)
            continue;

        const MGeomPix &gpix = geom[i/*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 happen not very often, because fSize>0)
    // we cannot calculate Length and Width.
    // In reallity it is almost impossible to have a distribution of
    // cerenkov photons in the used pixels which is exactly symmetric
    // along one of the axis.
    // It seems to be less than 0.1% of all events.
    //
    if (corrxy==0)
        return 4;

    //
    // 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
    //
    // Rem: I tested replacing sqrt() by hypot() but they exactly
    //      consume the same amount of time
    //
    const Double_t d0    = corryy - corrxx;
    const Double_t d1    = corrxy*2;
    const Double_t d2    = d0 + TMath::Sqrt(d0*d0 + d1*d1);
    const Double_t tand  = d2 / d1;
    const Double_t tand2 = tand*tand;

    fDelta = TMath::ATan(tand);

    const Double_t s2 = tand2+1;
    const Double_t s  = TMath::Sqrt(s2);

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

    const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/fSize;
    const Double_t axis2 = (tand2*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 : TMath::Sqrt(axis1);  // [mm]
    fWidth  = axis2<0 ? 0 : TMath::Sqrt(axis2);  // [mm]

    SetReadyToSave();

    return 0;
}

// --------------------------------------------------------------------------
//
// This function is ment for special usage, please never try to set
// values via this function
//
void MHillas::Set(const TArrayF &arr)
{
    if (arr.GetSize() != 6)
        return;

    fLength = arr.At(0);  // [mm]        major axis of ellipse
    fWidth  = arr.At(1);  // [mm]        minor axis of ellipse
    fDelta  = arr.At(2);  // [rad]       angle of major axis with x-axis
    fSize   = arr.At(3);  // [#CerPhot]  sum of content of all pixels (number of Cherenkov photons)
    fMeanX  = arr.At(4);  // [mm]        x-coordinate of center of ellipse
    fMeanY  = arr.At(5);  // [mm]        y-coordinate of center of ellipse
}