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

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

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

#include "MHillas.h"

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

#include "MCerPhotEvt.h"
#include "MCerPhotPix.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 MCerPhotEvt &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]

    MCerPhotPix *pix = 0;

    TIter Next(evt);
    while ((pix=(MCerPhotPix*)Next()))
    {
        // 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[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;

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

        // count inner pixels: To dependent on MAGIC Camera --> FIXME

        if (pixid<397){
            fInnerSize += nphot;
            if(pixid>270){
                edgepixin2 += nphot;
                if(pixid>330)
                    edgepixin1 += nphot;
            }
        }

        // Compute Concentration 1-2
        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();

    fConc  = (maxpix1+maxpix2)/hillas.GetSize();         // [ratio]
    fConc1 = maxpix1/hillas.GetSize();                   // [ratio]

    SetReadyToSave();
} 

// --------------------------------------------------------------------------
//
void MNewImagePar::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "New Image Parameters (" << GetName() << ")" << endl;
    *fLog << " - Leakage1       [1]   = " << fLeakage1             << endl;
    *fLog << " - Leakage2       [1]   = " << fLeakage2             << endl;
    *fLog << " - InnerLeakage1  [1]   = " << fInnerLeakage1        << endl;
    *fLog << " - InnerLeakage2  [1]   = " << fInnerLeakage2        << endl;
    *fLog << " - InnerSize      [#]   = " << fInnerSize            << " CerPhot" << endl;
    *fLog << " - Conc           [1]   = " << fConc                 << " (ratio)" << endl;
    *fLog << " - Conc1          [1]   = " << fConc1                << " (ratio)" << endl;
    *fLog << " - Used Pixels    [#]   = " << fNumUsedPixels        << " Pixels" << endl;
    *fLog << " - Core Pixels    [#]   = " << fNumCorePixels        << " Pixels" << 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 << "New Image Parameters (" << GetName() << ")" << endl;
    *fLog << " - Leakage1       [1]   = " << fLeakage1             << endl;
    *fLog << " - Leakage2       [1]   = " << fLeakage2             << endl;
    *fLog << " - InnerLeakage1  [1]   = " << fInnerLeakage1        << endl;
    *fLog << " - InnerLeakage2  [1]   = " << fInnerLeakage2        << endl;
    *fLog << " - InnerSize      [#]   = " << fInnerSize            << " CerPhot" << endl;
    *fLog << " - Conc           [1]   = " << fConc                 << " (ratio)" << endl;
    *fLog << " - Conc1          [1]   = " << fConc1                << " (ratio)" << endl;
    *fLog << " - Used Pixels    [#]   = " << fNumUsedPixels        << " Pixels" << endl;
    *fLog << " - Core Pixels    [#]   = " << fNumCorePixels        << " Pixels" << endl;
    *fLog << " - Used Area    [deg^2] = " << fUsedArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
    *fLog << " - Core Area    [deg^2] = " << fCoreArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
}