source: trunk/MagicSoft/Mars/manalysis/MHillasExt.cc@ 1211

/* ======================================================================== *\
!
! *
! * 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): Rudolf Bock     10/2001 <mailto:Rudolf.Bock@cern.ch>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHillasExt
//
// Storage Container for extended image parameters
//
//    extended image parameters
// fConc     ratio of sum of two highest pixels over fSize
// fConc1    ratio of highest pixel over fSize
// fAsym     distance from highest pixel to center, projected onto major axis
// fM3Long   third moment along major axis
// fM3Trans  third moment along minor axis
//
////////////////////////////////////////////////////////////////////////////
#include "MHillasExt.h"

#include <fstream.h>

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

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

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

ClassImp(MHillasExt);

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

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

// -------------------------------------------------------------------------
//
void MHillasExt::Reset()
{
    MHillas::Reset();

    fConc    = 0;
    fConc1   = 0; 
    fAsym    = 0;
    fM3Long  = 0;
    fM3Trans = 0;
}

// -------------------------------------------------------------------------
//
void MHillasExt::Print(Option_t *) const
{
    MHillas::Print();

    *fLog << "Extended Image Parameters (" << GetName() << ")" << endl;
    *fLog << " - Conc      = "        << fConc    << " (ratio)" << endl;
    *fLog << " - Conc1     = "        << fConc1   << " (ratio)" << endl;
    *fLog << " - Asymmetry = "        << fAsym    << " mm" << endl;
    *fLog << " - 3rd Moment Long  = " << fM3Long  << " mm" << endl;
    *fLog << " - 3rd Moment Trans = " << fM3Trans << " mm" << endl;
}

// -------------------------------------------------------------------------
//
//  calculation of additional parameters based on the camera geometry
// and the cerenkov photon event
//
Bool_t MHillasExt::Calc(const MGeomCam &geom, const MCerPhotEvt &evt)
{
    if (!MHillas::Calc(geom, evt))
        return kFALSE;

    //
    //   calculate the additional image parameters
    // --------------------------------------------
    //
    //   loop to get third moments along ellipse axes and two max pixels
    //
    float m3x = 0;
    float m3y = 0;

    float maxpix1 = 0;                                               // [#phot]
    float maxpix2 = 0;                                               // [#phot]

    float maxpixx = 0;                                               // [mm]
    float maxpixy = 0;                                               // [mm]

    const UInt_t npixevt = evt.GetNumPixels();

    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() - GetMeanX();                   // [mm]
        const float dy = gpix.GetY() - GetMeanY();                   // [mm]

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

        const float dzx =  GetCosDelta()*dx + GetSinDelta()*dy;      // [mm]
        const float dzy = -GetSinDelta()*dx + GetCosDelta()*dy;      // [mm]

        m3x += nphot * dzx*dzx*dzx;                                  // [mm^3]
        m3y += nphot * dzy*dzy*dzy;                                  // [mm^3]

        if (nphot>maxpix1)
        {
            maxpix2 = maxpix1;
            maxpix1 = nphot;                                         // [1]
            maxpixx = dx + GetMeanX();                               // [mm]
            maxpixy = dy + GetMeanY();                               // [mm]
            continue;                                                // [1]
        }

        if (nphot>maxpix2)
            maxpix2 = nphot;                                         // [1]
    }

    fAsym  = (GetMeanX()-maxpixx)*GetCosDelta() + (GetMeanY()-maxpixy)*GetSinDelta(); // [mm]
    fConc  = (maxpix1+maxpix2)/GetSize();                            // [ratio]
    fConc1 = maxpix1/GetSize();                                      // [ratio]

    //
    // Third moments along axes get normalized
    //
    m3x /= GetSize();
    m3y /= GetSize();

    fM3Long  = m3x<0 ? -pow(fabs(m3x), 1./3) : pow(fabs(m3x), 1./3); // [mm]
    fM3Trans = m3y<0 ? -pow(fabs(m3y), 1./3) : pow(fabs(m3y), 1./3); // [mm]

    SetReadyToSave();

    return kTRUE;
}

// -------------------------------------------------------------------------
//
void MHillasExt::AsciiRead(ifstream &fin)
{
    MHillas::AsciiRead(fin);

    fin >> fConc;
    fin >> fConc1;
    fin >> fAsym;
    fin >> fM3Long;
    fin >> fM3Trans;
}

// -------------------------------------------------------------------------
//
void MHillasExt::AsciiWrite(ofstream &fout) const
{
    MHillas::AsciiWrite(fout);

    fout << " ";
    fout << fConc   << " ";
    fout << fConc1  << " ";
    fout << fAsym   << " ";
    fout << fM3Long << " ";
    fout << fM3Trans;
}