source: trunk/MagicSoft/Mars/mimage/MNewImagePar.cc@ 7384

/* ======================================================================== *\
!
! *
! * 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): Wolfgang Wittek 03/2003 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Thomas Bretz            <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MNewImagePar
//
// Storage Container for new image parameters
//
//    Float_t fLeakage1;             // (photons in most outer ring of pixels) over fSize
//    Float_t fLeakage2;             // (photons in the 2 outer rings of pixels) over fSize
//    Float_t fInnerLeakage1;        // (photons in most outer rings of inner pixels) over fInnerSize
//    Float_t fInnerLeakage2;        // (photons in the 2 outer rings of inner pixels) over fInnerSize
//    Float_t fInnerSize;            //
//
//    Float_t fConc;                 // [ratio] concentration ratio: sum of the two highest pixels / fSize
//    Float_t fConc1;                // [ratio] concentration ratio: sum of the highest pixel / fSize
//
//    Float_t fUsedArea;             // Area of pixels which survived the image cleaning
//    Float_t fCoreArea;             // Area of core pixels
//    Short_t fNumUsedPixels;        // Number of pixels which survived the image cleaning
//    Short_t fNumCorePixels;        // number of core pixels
//    Short_t fNumHGSaturatedPixels; // number of pixels with saturating hi-gains
//    Short_t fNumSaturatedPixels;   // number of pixels with saturating lo-gains
//
// Version 2:
// ----------
//  - added fNumSaturatedPixels
// 
// Version 3:
// ----------
//  - added fNumHGSaturatedPixels
//  - added fInnerLeakage1
//  - added fInnerLeakage2
//  - added fInnerSize
//  - added fUsedArea
//  - added fCoreArea
//
// Version 4:
// ----------
//  - moved cleaning/island independant parameters to MImagePar:
//    + removed fNumHGSaturatedPixels
//    + removed fNumSaturatedPixels
//
//
/////////////////////////////////////////////////////////////////////////////
#include "MNewImagePar.h"

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

#include "MHillas.h"

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

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

ClassImp(MNewImagePar);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MNewImagePar::MNewImagePar(const char *name, const char *title)
{
    fName  = name  ? name  : "MNewImagePar";
    fTitle = title ? title : "New image parameters";

    Reset();
}

// --------------------------------------------------------------------------
//
void MNewImagePar::Reset()
{
    fLeakage1 = -1;
    fLeakage2 = -1;

    fInnerLeakage1 = -1;
    fInnerLeakage2 = -1;
    fInnerSize = -1;

    fConc  = -1;
    fConc1 = -1;

    fNumUsedPixels = -1;
    fNumCorePixels = -1;

    fUsedArea = -1;
    fCoreArea = -1;
}

// --------------------------------------------------------------------------
//
//  Calculation of new image parameters
//
void MNewImagePar::Calc(const MGeomCam &geom, const MSignalCam &evt,
                        const MHillas &hillas, Int_t island)
{
    fNumUsedPixels = 0;
    fNumCorePixels = 0;

    fUsedArea = 0;
    fCoreArea = 0;

    fInnerSize = 0;

    Double_t edgepix1 = 0;
    Double_t edgepix2 = 0;

    Double_t edgepixin1 = 0;
    Double_t edgepixin2 = 0;

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

    const Bool_t ismagiclike =
        geom.GetNumPixels() == 577 &&
        geom.GetNumAreas()  == 2   &&
        geom.GetPixRatio(396) > geom.GetPixRatio(397);

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

        // Check for requested islands
        if (island>=0 && pix.GetIdxIsland()!=island)
            continue;

        // Get geometry of pixel
        //const Int_t pixid = pix->GetPixId();
        const MGeomPix &gpix = geom[i/*pixid*/];

        // count used and core pixels
        if (pix.IsPixelCore())
        {
            fNumCorePixels++;
            fCoreArea += gpix.GetA();
        }

        // count used pixels
        fNumUsedPixels++;
        fUsedArea += gpix.GetA();

        Double_t nphot = pix.GetNumPhotons();

        //
        // count photons in outer rings of camera
        //
        if (gpix.IsInOutermostRing())
           edgepix1 += nphot;
        if (gpix.IsInOuterRing())
           edgepix2 += nphot;

        const Double_t ratio = geom.GetPixRatio(i);
        if (TMath::Nint(ratio)==1) // Means this is a small (= inner pixel)
        {
            fInnerSize += nphot;

            // Do calculation of "inner leakage" only for MAGIC-like geometry,
            // i.e., 577 pixels, pixels of 2 different areas, inner part
            // from pixel 0 to 396:
            if (ismagiclike)
            {
                if(i > 270) // last two "rings" of inner pixels
                {
                    edgepixin2 += nphot;
                    if(i > 330)  // last "ring" of inner pixels
                        edgepixin1 += nphot;
                }
            }
        }

        //
        // Now convert nphot from absolute number of photons or phe to signal
        // density (divide by pixel area), to find the pixel with highest signal
        // density:
        //
        nphot *= ratio;

        // Look for signal density in two highest pixels:
        if (nphot>maxpix1)
        {
            maxpix2  = maxpix1;
            maxpix1  = nphot;                            // [1]
            continue;                                    // [1]
        }

        if (nphot>maxpix2)
            maxpix2 = nphot;                             // [1]
    }
    
    fInnerLeakage1 = edgepixin1 / fInnerSize;
    fInnerLeakage2 = edgepixin2 / fInnerSize;

    fLeakage1 = edgepix1 / hillas.GetSize();
    fLeakage2 = edgepix2 / hillas.GetSize();

    // FIXME?: in case the pixel with highest signal density is an outer pixel,
    // the value of fConc (ratio of signal in two highest pixels to SIZE) should
    // rather be 2*fConc1, under the simplest assumption that the light density 
    // inside the outer (large) pixel is uniform.
    fConc  = (maxpix1+maxpix2)/hillas.GetSize();         // [ratio]
    fConc1 = maxpix1/hillas.GetSize();                   // [ratio]

    SetReadyToSave();
} 

// --------------------------------------------------------------------------
//
void MNewImagePar::Print(Option_t *) const
{
    *fLog << all;
    *fLog << GetDescriptor() << endl;
    *fLog << " - Leakage1         [1] = " << fLeakage1      << endl;
    *fLog << " - Leakage2         [1] = " << fLeakage2      << endl;
    *fLog << " - InnerLeakage1    [1] = " << fInnerLeakage1 << endl;
    *fLog << " - InnerLeakage2    [1] = " << fInnerLeakage2 << endl;
    *fLog << " - InnerSize      [phe] = " << fInnerSize     << endl;
    *fLog << " - Conc             [1] = " << fConc          << endl;
    *fLog << " - Conc1            [1] = " << fConc1         << endl;
    *fLog << " - Num Used Pixels  [#] = " << fNumUsedPixels << endl;
    *fLog << " - Num Core Pixels  [#] = " << fNumCorePixels << endl;
    *fLog << " - Used Area     [mm^2] = " << fUsedArea      << endl;
    *fLog << " - Core Area     [mm^2] = " << fCoreArea      << endl;
}

// -------------------------------------------------------------------------
//
// Print contents of MNewImagePar to *fLog, depending on the geometry in
// units of deg.
//
void MNewImagePar::Print(const MGeomCam &geom) const
{
    *fLog << all;
    *fLog << GetDescriptor() << endl;
    *fLog << " - Leakage1         [1] = " << fLeakage1      << endl;
    *fLog << " - Leakage2         [1] = " << fLeakage2      << endl;
    *fLog << " - InnerLeakage1    [1] = " << fInnerLeakage1 << endl;
    *fLog << " - InnerLeakage2    [1] = " << fInnerLeakage2 << endl;
    *fLog << " - InnerSize      [phe] = " << fInnerSize     << endl;
    *fLog << " - Conc             [1] = " << fConc          << endl;
    *fLog << " - Conc1            [1] = " << fConc1         << endl;
    *fLog << " - Num Used Pixels  [#] = " << fNumUsedPixels << endl;
    *fLog << " - Num Core Pixels  [#] = " << fNumCorePixels << endl;
    *fLog << " - Used Area    [deg^2] = " << fUsedArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
    *fLog << " - Core Area    [deg^2] = " << fCoreArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
}