/* ======================================================================== *\
!
! *
! * 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>
!              Wolfgang Wittek 01/2003 <mailto:wittek@mppmu.mpg.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 "MPedestalCam.h"
#include "MPedestalPix.h"
#include "MGeomPix.h"
#include "MCerPhotEvt.h"
#include "MCerPhotPix.h"

ClassImp(MSigmabar);

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

    fCalcPixNum=kTRUE;
}

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

// --------------------------------------------------------------------------
//
// 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 MPedestalCam &ped, 
                        const MCerPhotEvt &evt)
{
  fSigmabar      = 0.0;
  fSigmabarInner = 0.0;
  fSigmabarOuter = 0.0;


  Float_t innerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  Float_t outerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  Int_t innerPixels[6]       = {0,0,0,0,0,0};
  Int_t outerPixels[6]       = {0,0,0,0,0,0};

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

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

      Int_t j = cerpix.GetPixId();
      Double_t area = geom.GetPixRatio(j);

      const MGeomPix    &gpix = geom[j];
      const MPedestalPix &pix =  ped[j]; 

      //angle = 6.0*atan2(gpix.GetX(),gpix.GetY()) / (2.0*TMath::Pi()) - 0.5;
      //if (angle<0.0) angle+=6.0;

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

      Int_t currentSector=(Int_t)angle;
       
      // count only those pixels which have a sigma != 0.0 
      Float_t sigma = pix.GetMeanRms();

      if ( sigma != 0.0 )
      {  
        if (area < 1.5)
        {
          innerPixels[currentSector]++;
          innerSquaredSum[currentSector]+= sigma*sigma / area;
        }
        else
        {
          outerPixels[currentSector]++;
          outerSquaredSum[currentSector]+= sigma*sigma / area;
        }
      }
  }
 
  fSigmabarInner=0; fInnerPixels=0;
  fSigmabarOuter=0; fOuterPixels=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;
  //
  if ( (fInnerPixels+fOuterPixels) > 0)
    fSigmabar=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++) {
    fSigmabarSector[i]         = 0.0;
    fSigmabarInnerSector[i]    = 0.0;
    fSigmabarOuterSector[i]    = 0.0;

    if ( (innerPixels[i]+outerPixels[i]) > 0)
      fSigmabarSector[i] = sqrt(  (innerSquaredSum[i] + outerSquaredSum[i])
                                / (innerPixels[i]     + outerPixels[i]    ) );
    if ( innerPixels[i] > 0)
      fSigmabarInnerSector[i] = innerSquaredSum[i] / innerPixels[i];
    if ( outerPixels[i] > 0)
      fSigmabarOuterSector[i] = outerSquaredSum[i] / outerPixels[i];


    fSigmabarInnerSector[i] = sqrt( fSigmabarInnerSector[i] );
    fSigmabarOuterSector[i] = sqrt( fSigmabarOuterSector[i] );
  }
    
  return fSigmabar;
}

// --------------------------------------------------------------------------
//
void MSigmabar::Print(Option_t *) const
{
  *fLog << 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;

}
// --------------------------------------------------------------------------