/* ======================================================================== *\
!
! *
! * 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): Robert Wagner,   10/2002 <mailto:magicsoft@rwagner.de>
!   Author(s): Wolfgang Wittek, 01/2003 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Thomas Bretz,    04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MSigmabar                                                               //
//                                                                         //
// This is the storage container to hold information about                 //
// "average" pedestal sigmas                                               //
//                                                                         //
// In calculating averages all sigmas are normalized to the area of pixel 0//
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MSigmabar.h"

#include <TMath.h>

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

#include "MParList.h"

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

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

#include "MPedPhotCam.h"
#include "MPedPhotPix.h"

ClassImp(MSigmabar);

using namespace std;

// --------------------------------------------------------------------------
//
MSigmabar::MSigmabar(const char *name, const char *title)
{
    fName  = name  ? name  : "MSigmabar";
    fTitle = title ? title : "Storage container for Sigmabar";
}

// --------------------------------------------------------------------------
//
MSigmabar::~MSigmabar()
{
  // do nothing special.
}

void MSigmabar::Reset()
{
    fSigmabar      = -1;
    fInnerPixels   = -1;
    fOuterPixels   = -1;
    fSigmabarInner = -1;
    fSigmabarOuter = -1;

    memset(fSigmabarSector, 0, sizeof(fSigmabarSector));
    memset(fSigmabarInnerSector, 0, sizeof(fSigmabarInnerSector));
    memset(fSigmabarOuterSector, 0, sizeof(fSigmabarOuterSector));
}

// --------------------------------------------------------------------------
//
// Actual calculation of sigmabar. This is done for each of the six sectors
// separately due to their possibly different HV behavior. Also inner and
// outer pixels are treated separately.
//
// Preliminary! Works for CT1 test, for real MAGIC crosschecks still have
// to be done. Also implementation details will be updated, like 
// determination of sector to which a respective pixel belongs
//
Float_t MSigmabar::Calc(const MGeomCam &geom, const MPedPhotCam &ped, 
                        const MCerPhotEvt &evt)
{
    Int_t innerPixels[6];
    Int_t outerPixels[6];
    Float_t innerSquaredSum[6];
    Float_t outerSquaredSum[6];

    memset(innerPixels, 0, sizeof(innerPixels));
    memset(outerPixels, 0, sizeof(outerPixels));
    memset(innerSquaredSum, 0, sizeof(innerSquaredSum));
    memset(outerSquaredSum, 0, sizeof(outerSquaredSum));

    //
    // sum up sigma^2 for each sector, separately for inner and outer region;
    // all pixels are renormalized to the area of pixel 0
    //
    // consider all pixels with Cherenkov photon information
    // and require "Used"
    //

    const UInt_t npix = evt.GetNumPixels();

    //*fLog << "MSigmabar : npix = " << npix << endl;

    for (UInt_t i=0; i<npix; i++)
    {
        MCerPhotPix &cerpix = evt.operator[](i);
        if (!cerpix.IsPixelUsed())
            continue;

	/*
         if ( cerpix.GetNumPhotons() == 0 )
         {
         *fLog << "MSigmabar::Calc(); no.of photons is 0 for used pixel"
         << endl;
         continue;
         }
	*/

        const Int_t idx = cerpix.GetPixId();
        if (idx == 0)
	{
          //*fLog << "MSigmabar : id = 0;  pixel '0' is used, ignore it" 
	  //      << endl;
          continue;
        }

        // ratio is the area of pixel 0
        //          divided by the area of the current pixel
        const Double_t ratio = geom.GetPixRatio(idx);

        const MGeomPix &gpix = geom[idx];

        Int_t sector = (Int_t)(atan2(gpix.GetY(),gpix.GetX())*6 / (TMath::Pi()*2));
        if (sector<0)
            sector+=6;

        // count only those pixels which have a sigma != 0.0
        const Float_t sigma = ped[idx].GetRms();

        if ( sigma <= 0 )
            continue;

        //if (area < 1.5)
        if (ratio > 0.5)
        {
            innerPixels[sector]++;
            innerSquaredSum[sector]+= sigma*sigma * ratio;
        }
        else
        {
            outerPixels[sector]++;
            outerSquaredSum[sector]+= sigma*sigma * ratio;
        }
    }

    fInnerPixels   = 0;
    fOuterPixels   = 0;
    fSigmabarInner = 0;
    fSigmabarOuter = 0;
    for (UInt_t i=0; i<6; i++)
    {
        fSigmabarInner += innerSquaredSum[i];
        fInnerPixels   += innerPixels[i];
        fSigmabarOuter += outerSquaredSum[i];
        fOuterPixels   += outerPixels[i];
    }

    //
    // this is the sqrt of the average sigma^2;
    //
    fSigmabar=fInnerPixels+fOuterPixels<=0?0:sqrt((fSigmabarInner+fSigmabarOuter)/(fInnerPixels+fOuterPixels));

    if (fInnerPixels > 0) fSigmabarInner /= fInnerPixels;
    if (fOuterPixels > 0) fSigmabarOuter /= fOuterPixels;

    //
    // this is the sqrt of the average sigma^2
    // for the inner and outer pixels respectively
    //
    fSigmabarInner = sqrt( fSigmabarInner );
    fSigmabarOuter = sqrt( fSigmabarOuter );

    for (UInt_t i=0; i<6; i++)
    {
        const Double_t ip  = innerPixels[i];
        const Double_t op  = outerPixels[i];
        const Double_t iss = innerSquaredSum[i];
        const Double_t oss = outerSquaredSum[i];

        const Double_t sum = ip + op;

        fSigmabarSector[i]      = sum<=0 ? 0 : sqrt((iss+oss)/sum);
        fSigmabarInnerSector[i] = ip <=0 ? 0 : sqrt(iss/ip);
        fSigmabarOuterSector[i] = op <=0 ? 0 : sqrt(oss/op);
    }

  return fSigmabar;
}

// --------------------------------------------------------------------------
//
void MSigmabar::Print(Option_t *) const
{
  *fLog << all << endl;
  *fLog << "Total number of inner pixels is " << fInnerPixels << endl;
  *fLog << "Total number of outer pixels is " << fOuterPixels << endl;
  *fLog << endl;

  *fLog << "Sigmabar     Overall : " << fSigmabar << "   ";
  *fLog << " Sectors: ";
  for (Int_t i=0;i<6;i++)
      *fLog << fSigmabarSector[i] << ", ";
  *fLog << endl;


  *fLog << "Sigmabar     Inner   : " << fSigmabarInner << "   ";
  *fLog << " Sectors: ";
  for (Int_t i=0;i<6;i++)
      *fLog << fSigmabarInnerSector[i] << ", ";
  *fLog << endl;


  *fLog << "Sigmabar     Outer   : " << fSigmabarOuter << "   ";
  *fLog << " Sectors: ";
  for (Int_t i=0;i<6;i++)
      *fLog << fSigmabarOuterSector[i] << ", ";
  *fLog << endl;

}