source: trunk/MagicSoft/Mars/mimage/MConcentration.cc@ 3526

/* ======================================================================== *\
!
! *
! * 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): Harald Kornmayer 1/2001
!   Author(s): Rudolf Bock     10/2001 <mailto:Rudolf.Bock@cern.ch>
!   Author(s): Wolfgang Wittek  6/2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MConcentration
//
// Storage Container for image parameters
//
//    basic image parameters
//
// Version 1:
// ----------
// fLength   [mm]       major axis of ellipse
// fWidth    [mm]       minor axis
// fDelta    [rad]      angle of major axis with x-axis
//                      by definition the major axis is pointing into
//                      the hemisphere x>0, thus -pi/2 < delta < pi/2
// fSize     [#CerPhot] total sum of pixels
// fMeanX    [mm]       x of center of ellipse
// fMeanY    [mm]       y of center of ellipse
//
// Version 2:
// ----------
// fNumCorePixels  number of pixels called core
// fNumUsedPixels  number of pixels which survived the cleaning
//
// Version 3:
// ----------
// fNumCorePixels  moved to MNewImagePar
// fNumUsedPixels  moved to MNewImagePar
// fCosDelta       added
// fSinDelte       added
//
/////////////////////////////////////////////////////////////////////////////
#include "MConcentration.h"

//#include <fstream>

#include <TArrayF.h>
#include <TArrayI.h>
//#include <TEllipse.h>

#include "MHillas.h"

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

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

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

ClassImp(MConcentration);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MConcentration::MConcentration(const char *name, const char *title)
{
    fName  = name  ? name  : "MConcentration";
    fTitle = title ? title : "Storage container for concentrations";

    Reset();
}

// --------------------------------------------------------------------------
//
// Destructor. Deletes the TEllipse if one exists.
//
MConcentration::~MConcentration()
{
}

// --------------------------------------------------------------------------
//
// Initializes the values with defaults. For the default values see the
// source code.
//
void MConcentration::Reset()
{
    for (int i=0; i<9; i++)
        fConc[i] = -1;
}

// --------------------------------------------------------------------------
//
// Print the hillas Parameters to *fLog
//
void MConcentration::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "Concentrations (" << GetName() << ")" << endl;
    //*fLog << " - Length      [mm]  = " << fLength << endl;
}

// --------------------------------------------------------------------------
//
// Calculate the image parameters from a Cherenkov photon event
// assuming Cher.photons/pixel and pixel coordinates are given
// In case you don't call Calc from within an eventloop make sure, that
// you call the Reset member function before.
// Returns:
//   0  no error
//   1  number of pixels < 3
//   2  size==0
//   3  number of used pixel < 3
//   4  CorrXY == 0
//
Int_t MConcentration::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hillas)
{
    const UInt_t npixevt = evt.GetNumPixels();

    Float_t maxpix[9] = {0,0,0,0,0,0,0,0,0};             // [#phot]

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

        // skip unused pixels
        if (!pix.IsPixelUsed())
            continue;

        const Int_t pixid = pix.GetPixId();

        const MGeomPix &gpix = geom[pixid];

        Double_t nphot = pix.GetNumPhotons();

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

        // Compute Concentrations 1-8   

        if(maxpix[0]<=nphot) {
          for(int i=0;i<8;i++)
            maxpix[8-i]=maxpix[7-i];
          maxpix[0]=nphot;
        }
        else 
          for(int i=0;i<8;i++){
            if (nphot<maxpix[7-i]){
              for(int j=0;j<i-1;j++){  
                maxpix[7-j]=maxpix[6-j];                    // [1]
              }
              maxpix[8-i]=nphot;
              break;
            }
          }

    }

    // Compute concentrations from the 8 pixels with higher signal
   
    fConc[0]=maxpix[0];
    for(int i=1;i<8;i++)
      {
        fConc[i]=fConc[i-1]+maxpix[i];
      }
    
    for(int i=0;i<8;i++)
      {
        fConc[i]/=hillas.GetSize();                       // [ratio]
      }

    SetReadyToSave();

}