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

/* ======================================================================== *\
!
! *
! * 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>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MNewImagePar
//
// Storage Container for new image parameters
//
// fLeakage1 ratio: (photons in most outer ring of pixels) over fSize
// fLeakage2 ratio: (photons in the 2 outer rings of pixels) over fSize
// fNumSaturatedPixels: number of pixels in which at least one slice
//                      of the low gain FADC was saturated.
//
// Version 2:
// ----------
//  - added fNumSaturatedPixels
// 
// Version 3:
// ----------
//  - added fNumHGSaturatedPixels
//  - added fInnerLeakage1
//  - added fInnerLeakage2
//  - added fInnerSize
//  - added fUsedArea
//  - added fCoreArea
// 
// 
/////////////////////////////////////////////////////////////////////////////
#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";
}

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

    fNumSaturatedPixels = -1;
    fNumHGSaturatedPixels = -1;
}

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

    fUsedArea = 0;
    fCoreArea = 0;

    fNumSaturatedPixels = 0;
    fNumHGSaturatedPixels = 0;

    fInnerSize = 0;

    const UInt_t npixevt = evt.GetNumPixels();

    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]

    for (UInt_t i=0; i<npixevt; i++)
    {
        const MCerPhotPix &pix = evt[i];

        // count saturated pixels
        if (pix.IsPixelHGSaturated())
            fNumHGSaturatedPixels++;
        if (pix.IsPixelSaturated())
            fNumSaturatedPixels++;

        // skip unused pixels
        if (!pix.IsPixelUsed())
            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 Concetration 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 << " - Conc           [1]   = " << fConc          << " (ratio)" << endl;
    *fLog << " - Conc1          [1]   = " << fConc1         << " (ratio)" << endl;
    *fLog << " - Used Pixels    [#]   = " << fNumUsedPixels << " Pixels" << endl;
    *fLog << " - Core Pixels    [#]   = " << fNumCorePixels << " Pixels" << endl;
    *fLog << " - Sat. Pixels (HG) [#]  = " << fNumHGSaturatedPixels << " Pixels" << endl;
    *fLog << " - Sat. Pixels (LG) [#]  = " << fNumSaturatedPixels << " Pixels" << endl;
}