source: trunk/MagicSoft/Mars/mimage/MHillasExt.cc@ 2040

/* ======================================================================== *\
!
! *
! * 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>
!   Author(s): Wolfgang Wittek 06/2002 <mailto:wittek@mppmu.mpg.de>
!
!   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
//
// WARNING: Before you can use fAsym, fM3Long and fM3Trans you must
//          multiply by the sign of MHillasSrc::fCosDeltaAlpha
//
////////////////////////////////////////////////////////////////////////////
/*
 // fLeakage1 ratio : (photons in most outer ring of pixels) over fSize
 // fLeakage2 ratio : (photons in the 2 outer rings of pixels) over fSize
 //
 // fAsymna   d/(d na) of ( sum(x*q^na)/sum(q^na), sum(y*q^na)/sum(q^na) )
 //           projected onto the major axis
 // fAsym0    (F-B)/(F+B) along the major axis
 */
#include "MHillasExt.h"

#include <fstream.h>
#include <TArrayF.h>

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

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

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

#include "MHillas.h"

ClassImp(MHillasExt);

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

    Reset();
}

// -------------------------------------------------------------------------
//
void MHillasExt::Reset()
{
    fAsym    =  0;
    fM3Long  =  0;
    fM3Trans =  0;
}

// -------------------------------------------------------------------------
//
void MHillasExt::Print(Option_t *) const
{
    *fLog << "Extended Image Parameters (" << GetName() << ")" << 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
//
Int_t MHillasExt::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hil)
{
    //
    //   calculate the additional image parameters
    // --------------------------------------------
    //
    //  loop to get third moments along ellipse axes and two max pixels
    //
    //  For the moments double precision is used to make sure, that
    //  the complex matrix multiplication and sum is evaluated correctly.
    //
    Double_t m3x = 0;
    Double_t m3y = 0;

    const UInt_t npixevt = evt.GetNumPixels();

    Int_t maxpixid = 0;
    Float_t maxpix = 0;

    for (UInt_t i=0; i<npixevt; i++)
    {
        const MCerPhotPix &pix = evt[i];
        if (!pix.IsPixelUsed())
            continue;

        const Int_t pixid = pix.GetPixId();

        const MGeomPix &gpix = geom[pixid];
        const Double_t dx = gpix.GetX() - hil.GetMeanX();      // [mm]
        const Double_t dy = gpix.GetY() - hil.GetMeanY();      // [mm]

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

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

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

        //
        // Now we are working on absolute values of nphot, which
        // must take pixel size into account
        //
        nphot *= geom.GetPixRatio(pixid);

        if (nphot>maxpix)
        {
            maxpix   = nphot;                                  // [1]
            maxpixid = pixid;
            continue;                                          // [1]
        }

        /*
         //
         //  power na for calculating fAsymna;
         //  the value 1.5 was suggested by Thomas Schweizer
         //
         Double_t na = 1.5;

         //
         // get sums for calculating fAsymna
         // the outer pixels are 4 times as big (in area)
         // as the inner pixels !
         //
         const Double_t dummy = pow(nphot, na)/r;

         sna +=     dummy;
         xna += dzx*dummy;

         sna1 += sna/nphot;
         xna1 += xna/nphot;

         //
         // forward-backward asymmetry
         //
         fb += dzx<0 ? -nphot: nphot;
         */
    }

    const MGeomPix &maxp = geom[maxpixid];

    fAsym = (hil.GetMeanX()-maxp.GetX())*hil.GetCosDelta() +
            (hil.GetMeanY()-maxp.GetY())*hil.GetSinDelta();    // [mm]

    /*
     fAsym0    =       fb / GetSize();
     fAsymna   = na * (sna*xna1 - sna1*xna) / (sna*sna);
     */

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

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

    SetReadyToSave();

    return 0;
}

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

    fAsym    = arr.At(0); // [mm] fDist minus dist: center of ellipse, highest pixel
    fM3Long  = arr.At(1); // [mm] 3rd moment (e-weighted) along major axis
    fM3Trans = arr.At(2); // [mm] 3rd moment (e-weighted) along minor axis
}