/* ======================================================================== *\
!
! *
! * 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): Markus Gaug  08/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MCalibrationTestCalc
//
//   PreProcess(): Initialize pointers to MHCalibrationTestCam
//               
//   ReInit():     MCalibrationCam::InitSize(NumPixels) is called from MGeomApply (which allocates
//                 memory in a TClonesArray of type MCalibrationChargePix)
//                 Initializes pointer to MBadPixelsCam
//
//   Process():    Nothing to be done, histograms getting filled by MHCalibrationTestCam
//
//   PostProcess(): Print out interpolation results to file
//
//  Input Containers:
//   MHCalibrationTestCam
//   MBadPixelsCam
//   MGeomCam
//
//  Output Containers:
//   none
//  
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrationTestCalc.h"

#include <TSystem.h>
#include <TH1.h>
#include <TF1.h>

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

#include "MParList.h"

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

#include "MHCalibrationTestCam.h"
#include "MHCalibrationTestPix.h"

#include "MBadPixelsCam.h"
#include "MBadPixelsPix.h"

ClassImp(MCalibrationTestCalc);

using namespace std;

const Float_t MCalibrationTestCalc::fgPhotErrLimit = 4.5;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets the pointer to fTestCam and fGeom to NULL
// Sets outputpath to "."
// Sets outputfile to "TestCalibStat.txt"
// Sets fPhotErrLimit to fgPhotErrLimit  
//
// Calls:
// - Clear()
//
MCalibrationTestCalc::MCalibrationTestCalc(const char *name, const char *title)
    : fMaxNumBadPixelsCluster(-1), 
      fBadPixels(NULL), fTestCam(NULL), fGeom(NULL)
{
        
  fName  = name  ? name  : "MCalibrationTestCalc";
  fTitle = title ? title : "Task to output the results of MHCalibrationTestCam ";
  
  SetPhotErrLimit();

  SetOutputPath();
  SetOutputFile();

}


// -----------------------------------------------------------------------------------
//
// The following containers are searched and created if they were not found:
//
//  - MBadPixelsCam
//
Int_t MCalibrationTestCalc::PreProcess(MParList *pList)
{
  
  // 
  // Containers that are created in case that they are not there.
  //
  fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
  if (!fBadPixels)
    {
      *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl;
      return kFALSE;
    }
  
  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Search for the following input containers and abort if not existing:
//  - MGeomCam
//  - MHCalibrationTestCam
// 
//
Bool_t MCalibrationTestCalc::ReInit(MParList *pList )
{

  fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
  if (!fGeom)
    {
      *fLog << err << "No MGeomCam found... aborting." << endl;
      return kFALSE;
    }
  
  fTestCam = (MHCalibrationTestCam*)pList->FindObject("MHCalibrationTestCam");
  if (!fTestCam)
    {
      *fLog << err << "Cannot find MHCalibrationTestCam... aborting" << endl;
      *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationTestCam before..." << endl;
      return kFALSE;
    }


  return kTRUE;
}

// ----------------------------------------------------------------------------------
//  
// Nothing to be done in Process, but have a look at MHCalibrationTestCam, instead
// 
Int_t MCalibrationTestCalc::Process()
{
  return kTRUE;
}

// -----------------------------------------------------------------------
//
// Return if number of executions is null.
//
// Print out some statistics
//
Int_t MCalibrationTestCalc::PostProcess()
{

  //  if (GetNumExecutions()==0)
  //    return kFALSE;

  //
  // Re-direct the output to an ascii-file from now on:
  //
  MLog asciilog;
  asciilog.SetOutputFile(GetOutputFile(),kTRUE);
  SetLogStream(&asciilog);
  //
  // Finalize calibration statistics
  //
  FinalizeCalibratedPhotons();
  FinalizeNotInterpolated();
  CalcMaxNumBadPixelsCluster();


  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Pixel Interpolation status:" << endl;

  if (fGeom->InheritsFrom("MGeomCamMagic"))
    {
      *fLog << " " << setw(7) << "Not interpolateable Pixels: " 
            << Form("%s%3i%s%3i","Inner: ",fNumUninterpolateable[0],
                    " Outer: ",fNumUninterpolateable[1]) << endl;
      *fLog << " " << setw(7) << "Biggest not-interpolateable cluster: " 
            << fMaxNumBadPixelsCluster << endl;
    }
  
  *fLog << endl;  
  SetLogStream(&gLog);

  return kTRUE;
}


// ------------------------------------------------------------------------
//
//
// First loop: Calculate a mean and mean RMS of calibrated photons per area index 
//             Include only MHCalibrationTestPix's which are not empty (not interpolated)
//              
// Second loop: Get weighted mean number of calibrated photons and its RMS
//              excluding those deviating by more than fPhotErrLimit 
//              sigmas from the mean (obtained in first loop). Set 
//              MBadPixelsPix::kDeviatingNumPhots if excluded.
// 
void MCalibrationTestCalc::FinalizeCalibratedPhotons() const
{

  const UInt_t npixels  = fGeom->GetNumPixels();
  const UInt_t nareas   = fGeom->GetNumAreas();
  const UInt_t nsectors = fGeom->GetNumSectors();

  Double_t lowlim      [nareas];
  Double_t upplim      [nareas];
  Double_t areaphotons [nareas], sectorphotons [nsectors];
  Double_t areavars    [nareas], sectorvars    [nsectors];
  Double_t fittedmean  [nareas], fittedsigma   [nareas];
  Int_t   numareavalid[nareas], numsectorvalid[nsectors];

  memset(lowlim        ,0, nareas   * sizeof(Double_t));
  memset(upplim        ,0, nareas   * sizeof(Double_t));
  memset(fittedmean    ,0, nareas   * sizeof(Double_t));
  memset(fittedsigma   ,0, nareas   * sizeof(Double_t));
  memset(areaphotons   ,0, nareas   * sizeof(Double_t));
  memset(areavars      ,0, nareas   * sizeof(Double_t));
  memset(numareavalid  ,0, nareas   * sizeof(Int_t  ));
  memset(sectorphotons ,0, nsectors * sizeof(Double_t));
  memset(sectorvars    ,0, nsectors * sizeof(Double_t));
  memset(numsectorvalid,0, nsectors * sizeof(Int_t  ));
  
  //
  // First loop: Get mean number of calibrated photons and the RMS
  //             The loop is only to recognize later pixels with very deviating numbers
  //
  MHCamera camphotons(*fGeom,"Camphotons","Photons in Camera");

  for (UInt_t i=0; i<npixels; i++)
    {
      
      MHCalibrationTestPix &pix = (MHCalibrationTestPix&)(*fTestCam)[i];
      //
      // We assume that the pixels have been interpolated so far. 
      // The MBadPixelsCam does not give any more information
      //
      if (pix.IsEmpty())
          continue;

      const Double_t nphot = pix.GetMean();
      const Int_t    aidx  = (*fGeom)[i].GetAidx();

      camphotons.Fill(i,nphot);
      camphotons.SetUsed(i);

      areaphotons [aidx] += nphot;
      areavars    [aidx] += nphot*nphot;
      numareavalid[aidx] ++;
    } 

  for (UInt_t aidx=0; aidx<nareas; aidx++)
    {
      if (numareavalid[aidx] == 0)
        {
          *fLog << warn << GetDescriptor() << ": No pixels with valid number of calibrated photons found "
                << "in area index: " << aidx << endl;
          continue;
        }

      if (numareavalid[aidx] == 1)
        areavars[aidx]  = 0.;
      else if (numareavalid[aidx] == 0)
        {
          areaphotons[aidx] = -1.;
          areavars[aidx]    = -1.;
        }
      else
        {
          areavars[aidx]    = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx]) 
                            / (numareavalid[aidx]-1.);
          areaphotons[aidx] = areaphotons[aidx] / numareavalid[aidx];
        }
      
      if (areavars[aidx] < 0.)
        {
          *fLog << warn << GetDescriptor() << ": No pixels with valid variance of calibrated photons found "
                << "in area index: " << aidx << endl;
          continue;
        }
      
      lowlim [aidx] = areaphotons[aidx] - fPhotErrLimit*TMath::Sqrt(areavars[aidx]);
      upplim [aidx] = areaphotons[aidx] + fPhotErrLimit*TMath::Sqrt(areavars[aidx]);

      TArrayI area(1);
      area[0] = aidx;

      TH1D *hist = camphotons.ProjectionS(TArrayI(),area,"_py",100);
      hist->Fit("gaus","Q");
      const Double_t mean  = hist->GetFunction("gaus")->GetParameter(1);
      const Double_t sigma = hist->GetFunction("gaus")->GetParameter(2);
      const Int_t    ndf   = hist->GetFunction("gaus")->GetNDF();

      if (ndf < 2)
        {
          *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of calibrated photons " 
                << "in the camera with area index: " << aidx << endl;
          *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
          *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
          delete hist;
          continue;
        }
      
      const Double_t prob = hist->GetFunction("gaus")->GetProb();

      if (prob < 0.001)
        {
          *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of calibrated photons " 
                << "in the camera with area index: " << aidx << endl;
          *fLog << warn << GetDescriptor() << ": Fit probability " << prob 
                << " is smaller than 0.001 " << endl;
          *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
          delete hist;
          continue;
        }
      
      fittedmean [aidx] = mean;
      fittedsigma[aidx] = sigma;
      
      lowlim  [aidx] = mean  - fPhotErrLimit*sigma;
      upplim  [aidx] = mean  + fPhotErrLimit*sigma;

      *fLog << inf << GetDescriptor() 
            << ": Fitted number of calib. equiv. Cher. photons in area index " << aidx 
            << ": "  << Form("%7.2f +- %6.2f",fittedmean[aidx],fittedsigma[aidx]) << endl;

      delete hist;
    }

  *fLog << endl;

  memset(numareavalid,0,nareas*sizeof(Int_t));
  memset(areaphotons ,0,nareas*sizeof(Double_t));
  memset(areavars    ,0,nareas*sizeof(Double_t));

  //
  // Second loop: Get mean number of calibrated photons and its RMS excluding 
  //              pixels deviating by more than fPhotErrLimit sigma. 
  // 
  for (UInt_t i=0; i<npixels; i++)
    {
      
      MHCalibrationTestPix &pix = (MHCalibrationTestPix&)(*fTestCam)[i];

      const Int_t    aidx   = (*fGeom)[i].GetAidx();
      const Int_t    sector = (*fGeom)[i].GetSector();
      const Float_t  area   = (*fGeom)[i].GetA();
      const Double_t nphot  = pix.GetMean();

      if ( nphot < lowlim[aidx] || nphot > upplim[aidx] )
        {
          *fLog << warn << GetDescriptor() << ": Number of photons: " 
                << Form("%8.2f out of %3.1f sigma limit: ",nphot,fPhotErrLimit)
                << Form("[%8.2f,%8.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
          MBadPixelsPix &bad = (*fBadPixels)[i];
          bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhots );
          bad.SetUnsuitable  ( MBadPixelsPix::kUnsuitableRun     );
          continue;
        }
      
      areavars     [aidx] += nphot*nphot/area/area;
      areaphotons  [aidx] += nphot/area;
      numareavalid [aidx] ++;

      sectorvars    [sector] += nphot*nphot/area/area;
      sectorphotons [sector] += nphot/area;
      numsectorvalid[sector] ++;
    } 

  *fLog << endl;

  for (UInt_t aidx=0; aidx<nareas; aidx++)
    {

      if (numareavalid[aidx] == 1)
        areavars[aidx] = 0.;
      else if (numareavalid[aidx] == 0)
        {
          areaphotons[aidx] = -1.;
          areavars[aidx]    = -1.;
        }
      else
        {
          areavars[aidx] = (areavars[aidx] - areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx]) 
                         / (numareavalid[aidx]-1.);
          areaphotons[aidx] /=  numareavalid[aidx];
        }
      
      if (areavars[aidx] < 0. || areaphotons[aidx] <= 0.)
        {
          *fLog << warn << GetDescriptor() 
                << ": Mean number of photons per area in area index " 
                << aidx << " could not be calculated! Mean: " << areaphotons[aidx] 
                << " Variance: " << areavars[aidx] << endl;
          continue;
        }
      
      *fLog << inf << GetDescriptor() << ": Mean number of equiv. Cher. photons "
            << "per area in area idx " << aidx << ": "  
            << Form("%5.3f+-%5.4f  [ph/mm^2]",areaphotons[aidx],TMath::Sqrt(areavars[aidx])) << endl;
    }

  *fLog << endl;

  for (UInt_t sector=0; sector<nsectors; sector++)
    {
      
      if (numsectorvalid[sector] == 1)
        sectorvars[sector] = 0.;
      else if (numsectorvalid[sector] == 0)
        {
          sectorphotons[sector]  = -1.;
          sectorvars[sector]     = -1.;
        }
      else
        {
          sectorvars[sector] = (sectorvars[sector] 
                               - sectorphotons[sector]*sectorphotons[sector]/numsectorvalid[sector]
                               ) 
                             / (numsectorvalid[sector]-1.);
          sectorphotons[sector] /=  numsectorvalid[sector];
        }
      
      if (sectorvars[sector] < 0. || sectorphotons[sector] <= 0.)
        {
          *fLog << warn << GetDescriptor() 
                << ": Mean number of calibrated photons per area in sector " 
                << sector << " could not be calculated! Mean: " << sectorphotons[sector] 
                << " Variance: " << sectorvars[sector] << endl;
          continue;
        }
      
  
      *fLog << inf << GetDescriptor() << ": Mean number of equiv. Cher. photons "
            << "per area in sector " << sector << ":   "  
            << Form("%5.3f+-%5.4f  [ph/mm^2]",sectorphotons[sector],TMath::Sqrt(sectorvars[sector])) << endl;
    }

  return;
}


// -----------------------------------------------------------------------------------------------
//
// Print out statistics about not interpolated pixels
// 
void MCalibrationTestCalc::FinalizeNotInterpolated() 
{
  
  const TArrayI &arr = fTestCam->GetNotInterpolateablePixels();
  const Int_t areas  = fGeom->GetNumAreas();

  TArrayI *newarr[areas];

  for (Int_t aidx=0; aidx<areas; aidx++)
    newarr[aidx] = new TArrayI(0);

  fNumUninterpolateable.Set(areas);

  for (Int_t i=0; i<arr.GetSize(); i++)
    {
      const Int_t id   = arr[i];
      const Int_t aidx = (*fGeom)[id].GetAidx();
      const Int_t size = newarr[aidx]->GetSize();
      newarr[aidx]->Set(size+1);
      newarr[aidx]->AddAt(id,size);
      fNumUninterpolateable[aidx]++;
    }

  Int_t num = 0;

  for (Int_t aidx = 0; aidx<areas; aidx++)
    {
      *fLog << endl;
      *fLog << " " << setw(7)
            << "Not interpolated pixels by in area idx " << aidx << ": ";
      for (Int_t i=0; i<newarr[aidx]->GetSize(); i++)
        {
          *fLog << newarr[aidx]->At(i) << " ";
          num++;
        }
      *fLog << endl;
    }


  
  *fLog << " " << setw(7) << num << " total not interpolateable pixels " << endl;
  
}

void MCalibrationTestCalc::CalcMaxNumBadPixelsCluster()
{

  const TArrayI &arr = fTestCam->GetNotInterpolateablePixels();
  const Int_t   size = arr.GetSize();
  
  if (size == 0)
    {
      fMaxNumBadPixelsCluster = 0; 
      return;
    }
  
  if (size == 1)
    {
      fMaxNumBadPixelsCluster = 1; 
      return;
    }

  TArrayI knownpixels(0);
  Int_t clustersize    = 1;
  Int_t oldclustersize = 0;
  //
  // Loop over the not-interpolateable pixels:
  //
  for (Int_t i=0; i<size; i++)
    {

      const Int_t id   = arr[i];
      const Int_t knownsize = knownpixels.GetSize();
      knownpixels.Set(knownsize+1);
      knownpixels[knownsize] = id;
      LoopNeighbours(arr, knownpixels, clustersize, id);
      if (clustersize > oldclustersize)
        oldclustersize = clustersize;
      clustersize = 1;
    }

  fMaxNumBadPixelsCluster = oldclustersize; 

}


void MCalibrationTestCalc::LoopNeighbours( const TArrayI &arr, TArrayI &knownpixels, Int_t &clustersize, const Int_t idx )
{
  
  const MGeomPix &pix = (*fGeom)[idx];
  const Byte_t neighbours = pix.GetNumNeighbors();

  // 
  // Loop over the next neighbours: 
  // - Check if they are already in the list of known pixels
  // - If not, call loopneighbours for the rest
  // - grow clustersize for those
  //
  for (Int_t i=0;i<neighbours;i++)
    {
      const Int_t newid = pix.GetNeighbor(i);
      Bool_t known = kFALSE;

      for (Int_t j=knownpixels.GetSize()-1;j>=0;j--)
        if (newid == knownpixels.At(j))
          {
            known = kTRUE;
            break;
          }
      if (known)
        continue;

      for (Int_t k=0;k<arr.GetSize();k++)
        if (newid == arr.At(k))
          {
            // This is an unknown, new pixel in the cluster!!
            clustersize++;
            const Int_t knownsize = knownpixels.GetSize();
            knownpixels.Set(knownsize+1);
            knownpixels[knownsize] = newid;
            LoopNeighbours(arr, knownpixels, clustersize, newid);
          }
    }
}





// --------------------------------------------------------------------------
//
// Set the path for output file
// 
void MCalibrationTestCalc::SetOutputPath(TString path)
{
  fOutputPath = path;
  if (fOutputPath.EndsWith("/"))
    fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
}

void MCalibrationTestCalc::SetOutputFile(TString file)
{ 
  fOutputFile        = file; 
}

// --------------------------------------------------------------------------
//
// Get the output file
// 
const char* MCalibrationTestCalc::GetOutputFile()
{
  return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
}

const Int_t MCalibrationTestCalc::GetNumUninterpolateable(const Int_t aidx) const
{
  if (aidx < 0)
    return -1;

  return aidx > fNumUninterpolateable.GetSize() ? -1 : fNumUninterpolateable[aidx]; 
}