/* ======================================================================== *\
!
! *
! * 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)
    : 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();

  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Not interpolated pixels statistics:" << endl;  

  PrintNotInterpolated();
  
  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;
        }

      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 << " Fitted number of photons in area index " << aidx 
            << ": "  << Form("%6.2f +- %6.2f",fittedmean[aidx],fittedsigma[aidx]) << endl;

      delete hist;
    }

  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];
      MBadPixelsPix        &bad = (*fBadPixels)[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 << "Deviating number of calibrated photons: " 
                << Form("%4.2f",nphot) << " out of accepted limits: [" 
                << Form("%4.2f%s%4.2f",lowlim[aidx],",",upplim[aidx]) << "] in pixel " << i << endl;
          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] ++;
    } 

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

      areavars   [aidx] -=  areaphotons[aidx]*areaphotons[aidx]/numareavalid[aidx];
      areaphotons[aidx] /=  numareavalid[aidx];
      areavars   [aidx] /=  numareavalid[aidx]-1.;
      
      if (areavars[aidx] <= 0. || areaphotons[aidx] <= 0.)
        {
          *fLog << warn << " Mean number of photons per area in area index " 
                << aidx << " cannot be calculated! Mean: " << areaphotons[aidx] 
                << " Variance: " << areavars[aidx] << endl;
          continue;
        }
      
      *fLog << " Mean number of photons per area in area index " << aidx 
            << ": "  << Form("%5.4f +- %5.4f",areaphotons[aidx],TMath::Sqrt(areavars[aidx])) << endl;
    }

  for (UInt_t sector=0; sector<nsectors; sector++)
    {

      sectorvars   [sector] -=  sectorphotons[sector]*sectorphotons[sector]/numsectorvalid[sector];
      sectorphotons[sector] /=  numsectorvalid[sector];
      sectorvars   [sector] /=  numsectorvalid[sector]-1.;

      if (sectorvars[sector] <= 0. || sectorphotons[sector] <= 0.)
        {
          *fLog << warn << " Mean number of calibrated photons per area from sector " 
                << sector << " cannot be calculated! " << endl;
          continue;
        }
      
      *fLog << " Mean number of photons per area in sector " << sector
            << ": "  << Form("%5.4f +- %5.4f",sectorphotons[sector],TMath::Sqrt(sectorvars[sector])) << endl;
    }

  return;
}


// -----------------------------------------------------------------------------------------------
//
// Print out statistics about not interpolated pixels
// 
void MCalibrationTestCalc::PrintNotInterpolated() const 
{
  
  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);

  TArrayI numtot;
  numtot.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);
      numtot[aidx]++;
    }

  Int_t num = 0;

  for (Int_t aidx = 0; aidx<areas; aidx++)
    {
      *fLog << " " << setw(7) 
            << numtot[aidx] << " not interpolateable pixels in area index " << aidx << endl;
      *fLog << " " << setw(7)
            << "Pixel software idx: ";
      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;
  
}

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