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

//////////////////////////////////////////////////////////////////////////////
//
//   MCalibrationRelTimeCalc
//
//   Task to finalize the relative time calibration obtained 
//   from the fit results to the summed FADC slice distributions delivered by 
//   MCalibrationRelTimeCam, calculated and filled by MHCalibrationRelTimeCam, 
//
//   PreProcess(): Initialize pointers to MCalibrationRelTimeCam, 
//               
//   ReInit():     Initializes pointer to MBadPixelsCam
//
//   Process():    Nothing to be done, histograms getting filled by MHCalibrationChargeCam
//
//   PostProcess(): - FinalizeRelTimes()
//                  - FinalizeBadPixels()
//
//  Input Containers:
//   MCalibrationRelTimeCam
//   MBadPixelsCam
//   MGeomCam
//
//  Output Containers:
//   MCalibrationRelTimeCam
//   MBadPixelsCam
//
//  
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrationRelTimeCalc.h"

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

#include "MParList.h"

#include "MStatusDisplay.h"

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

#include "MCalibrationIntensityRelTimeCam.h"
#include "MCalibrationRelTimeCam.h"
#include "MCalibrationRelTimePix.h"

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


ClassImp(MCalibrationRelTimeCalc);

using namespace std;

const Float_t MCalibrationRelTimeCalc::fgRelTimeResolutionLimit = 1.0;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets all pointers to NULL
// 
// Initializes:
// - fRelTimeResolutionLimit to fgRelTimeResolutionimit
// - fOutputPath        to "."
// - fOutputFile        to "TimeCalibStat.txt"
//
// Calls:
// - Clear()
//
MCalibrationRelTimeCalc::MCalibrationRelTimeCalc(const char *name, const char *title)
    : fGeom(NULL), fFlags(0)
{
        
  fName  = name  ? name  : "MCalibrationRelTimeCalc";
  fTitle = title ? title : "Task to finalize the relative time calibration";
  
  SetRelTimeResolutionLimit();
  SetOutputPath();
  SetOutputFile("");
 
  Clear();
  
}

// --------------------------------------------------------------------------
//
// Sets:
// - all flags to kFALSE
// - all pointers to NULL
//
void MCalibrationRelTimeCalc::Clear(const Option_t *o)
{

  fIntensBad  = NULL;
  fBadPixels  = NULL;
  fCam        = NULL;
  fIntensCam  = NULL;

}

// --------------------------------------------------------------------------
//
// Search for the following input containers and abort if not existing:
//  - MGeomCam
//  - MCalibrationIntensityRelTimeCam or MCalibrationRelTimeCam 
//  - MBadPixelsIntensityCam or MBadPixelsCam
// 
// It defines the PixId of every pixel in:
//
// - MCalibrationRelTimeCam 
//
// It sets all pixels in excluded which have the flag fBadBixelsPix::IsBad() set in:
// 
// - MCalibrationRelTimePix
//
Bool_t MCalibrationRelTimeCalc::ReInit(MParList *pList )
{

  fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
  if (!fGeom)
    {
      *fLog << err << "No MGeomCam found... aborting." << endl;
      return kFALSE;
    }
  
  fIntensBad = (MBadPixelsIntensityCam*)pList->FindObject(AddSerialNumber("MBadPixelsIntensityCam"));
  if (fIntensBad)
    *fLog << inf << "Found MBadPixelsIntensityCam ... " << endl;
  else
    {
      fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
      if (!fBadPixels)
        {
          *fLog << err << "Cannot find MBadPixelsCam ... abort." << endl;
          return kFALSE;
        }
    }

  fIntensCam = (MCalibrationIntensityRelTimeCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityRelTimeCam"));
  if (fIntensCam)
    *fLog << inf << "Found MCalibrationIntensityRelTimeCam ... " << endl;
  else
    {
      fCam = (MCalibrationRelTimeCam*)pList->FindObject(AddSerialNumber("MCalibrationRelTimeCam"));
      if (!fCam)
        {
          *fLog << err << "Cannot find MCalibrationRelTimeCam ... abort." << endl;
          *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationRelTimeCam before..." << endl;
          return kFALSE;
        }
    }
  
  UInt_t npixels     = fGeom->GetNumPixels();
  
  MCalibrationRelTimeCam *relcam = fIntensCam 
    ? (MCalibrationRelTimeCam*)fIntensCam->GetCam()  : fCam;
  MBadPixelsCam          *badcam    = fIntensBad 
    ? (MBadPixelsCam*)        fIntensBad->GetCam()  : fBadPixels;

  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
      MBadPixelsPix          &bad =                          (*badcam)[i];
      
      if (bad.IsBad())
        {
          pix.SetExcluded();
          continue;
        }

      if (IsDebug())
        pix.SetDebug();
    }

  return kTRUE;
}

// -----------------------------------------------------------------------
//
// Return if number of executions is null.
//
Int_t MCalibrationRelTimeCalc::PostProcess()
{

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

  return Finalize();
}

// -----------------------------------------------------------------------
//
// First loop over pixels, average areas and sectors, call: 
//  - FinalizeRelTimes()
// for every entry. Count number of valid pixels in loop and return kFALSE
// if there are none (the "Michele check").
//
// Call FinalizeBadPixels()
//
// Call MParContainer::SetReadyToSave() for fCam
//
// Print out some statistics
//
Int_t MCalibrationRelTimeCalc::Finalize()
{
  
  //
  // First loop over pixels, call FinalizePedestals and FinalizeRelTimes
  //
  FinalizeRelTimes();

  //
  // Finalize Bad Pixels
  // 
  FinalizeBadPixels();

  //
  // Re-direct the output to an ascii-file from now on:
  //
  MLog *asciilog = fOutputFile.IsNull() ? 0 : new MLog;
  if (asciilog)
  {
      asciilog->SetOutputFile(GetOutputFile(),kTRUE);
      SetLogStream(asciilog);
  }

  //
  // Finalize calibration statistics
  //
  FinalizeUnsuitablePixels();

  if (fIntensCam)
    fIntensCam->SetReadyToSave();
  else
    fCam      ->SetReadyToSave();

  if (fIntensBad)
    fIntensBad->SetReadyToSave();
  else
    fBadPixels->SetReadyToSave();

  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Errors statistics:" << endl;  

  PrintUncalibrated(MBadPixelsPix::kDeviatingTimeResolution,    
                    Form("%s%2.1f%s","Time resolution less than ",fRelTimeResolutionLimit," FADC slices from Mean:   "));
  PrintUncalibrated(MBadPixelsPix::kRelTimeOscillating,   
                    "Pixels with changing Rel. Times   over time:      ");
  PrintUncalibrated(MBadPixelsPix::kRelTimeNotFitted,     
                    "Pixels with unsuccesful Gauss fit to the times:   ");

  if (asciilog)
  {
      SetLogStream(&gLog);
      delete asciilog;
  }

  return kTRUE;
}


// ----------------------------------------------------------------------------------------------------
//
//
// First loop: Calculate a mean and mean RMS of time resolution per area index 
//             Include only pixels which are not MBadPixelsPix::kUnsuitableRun or 
//             MBadPixelsPix::kUnreliableRun (see FinalizeBadPixels())
//              
// Second loop: Exclude those deviating by more than fRelTimeResolutionLimit FADC slices
//              from the mean (obtained in first loop). Set 
//              MBadPixelsPix::kDeviatingTimeResolution if excluded.
// 
void MCalibrationRelTimeCalc::FinalizeRelTimes()
{

  MCalibrationRelTimeCam *relcam = fIntensCam 
    ? (MCalibrationRelTimeCam*)fIntensCam->GetCam()  : fCam;
  MBadPixelsCam         *badcam    = fIntensBad 
    ? (MBadPixelsCam*)        fIntensBad->GetCam()  : fBadPixels;

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

  TArrayF lowlim       (nareas);
  TArrayF upplim       (nareas);
  TArrayF areasum      (nareas);
  //  Float_t areasum2     [nareas];
  TArrayI numareavalid (nareas);
  TArrayI useunreliable(nareas);

  //
  // Apero loop: Count number of unreliable pixels:
  //     
  for (UInt_t i=0; i<npixels; i++)
    {
      MBadPixelsPix &bad = (*badcam)[i];      
      const Int_t  aidx  = (*fGeom)[i].GetAidx();

      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        continue;

      if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
        continue;

      numareavalid[aidx] ++;
    }
  
  for (UInt_t aidx=0; aidx<nareas; aidx++)
    if (numareavalid[aidx] < 100)
      useunreliable[aidx] = 1;

  numareavalid.Reset();
  //
  // First loop: Get mean time resolution the RMS
  //             The loop is only to recognize later pixels with very deviating numbers
  //
  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
      MBadPixelsPix          &bad =                          (*badcam)[i];
      
      if (pix.IsExcluded())
        continue;

      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        continue;

      const Int_t   aidx  = (*fGeom)[i].GetAidx();

      if (!useunreliable[aidx])
        if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
          continue;

      const Float_t res   = pix.GetTimePrecision();

      areasum     [aidx] += res;
      //      areasum2    [aidx] += res*res;
      numareavalid[aidx] ++;
    }


  for (UInt_t aidx=0; aidx<nareas; aidx++)
    {
      if (numareavalid[aidx] < 20)
        {
          *fLog << warn << GetDescriptor() << ": Less than 20 pixels with valid time resolution found "
                << "in area index: " << aidx << endl;
          continue;
        }
 
      // Calculate the rms out of sum2:
      /*
      areasum2[aidx]  = (areasum2[aidx] - areasum[aidx]*areasum[aidx]/numareavalid[aidx]);
      areasum2[aidx] /= (numareavalid[aidx]-1.);
      */
      areasum [aidx] /= numareavalid[aidx];
      lowlim  [aidx]  = 0.;
      upplim  [aidx]  = areasum [aidx] + fRelTimeResolutionLimit;
      
    }
  *fLog << endl;  


  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
      MBadPixelsPix          &bad =                          (*badcam)[i];
      
      if (pix.IsExcluded())
        continue;
      
      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        continue;
      
      const Float_t res    = pix.GetTimePrecision();
      const Int_t   aidx   = (*fGeom)[i].GetAidx();
      
      if ( res < lowlim[aidx] || res > upplim[aidx] )
        {
          *fLog << warn << "Deviating time resolution: "
            << Form("%4.2f",res) << " out of accepted limits ["
                << Form("%4.2f,%4.2f",lowlim[aidx],upplim[aidx]) << "] in pixel " << i << endl;
          bad.SetUncalibrated( MBadPixelsPix::kDeviatingTimeResolution);
          pix.SetExcluded();
        }
    }
}


// -----------------------------------------------------------------------------------
//
// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
// - MBadPixelsPix::kRelTimeIsPedestal
// - MBadPixelsPix::kRelTimeErrNotValid 
// - MBadPixelsPix::kRelTimeRelErrNotValid 
//
// - Call MCalibrationPix::SetExcluded() for the bad pixels
//
void MCalibrationRelTimeCalc::FinalizeBadPixels()
{
  
  MCalibrationRelTimeCam *relcam = fIntensCam 
    ? (MCalibrationRelTimeCam*)fIntensCam->GetCam()  : fCam;
  MBadPixelsCam         *badcam    = fIntensBad 
    ? (MBadPixelsCam*)        fIntensBad->GetCam()  : fBadPixels;

  for (Int_t i=0; i<badcam->GetSize(); i++)
    {
      
      MBadPixelsPix          &bad =                          (*badcam)[i];
      MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];

      if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingTimeResolution))
        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun   );
 
      if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeNotFitted))
        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun   );
 
      if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeOscillating))
        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun   );
 
      if (bad.IsUnsuitable(   MBadPixelsPix::kUnsuitableRun    ))
        pix.SetExcluded();

    }
}


// -----------------------------------------------------------------------------------------------
//
// - Print out statistics about BadPixels of type UnsuitableType_t 
// - store numbers of bad pixels of each type in fIntensCam or fCam
//
void MCalibrationRelTimeCalc::FinalizeUnsuitablePixels()
{
  
  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Rel. Times Calibration status:" << endl;
  *fLog << dec << setfill(' ');

  MCalibrationRelTimeCam *relcam = fIntensCam 
    ? (MCalibrationRelTimeCam*)fIntensCam->GetCam()  : fCam;
  MBadPixelsCam         *badcam    = fIntensBad 
    ? (MBadPixelsCam*)        fIntensBad->GetCam()  : fBadPixels;

  const Int_t nareas = fGeom->GetNumAreas();

  TArrayI counts(nareas);

  for (Int_t i=0; i<badcam->GetSize(); i++)
    {
      MBadPixelsPix &bad = (*badcam)[i];
      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        {
          const Int_t aidx = (*fGeom)[i].GetAidx();
          counts[aidx]++;
        }
    }

  for (Int_t aidx=0; aidx<nareas; aidx++)
    relcam->SetNumUnsuitable(counts[aidx], aidx);

  if (fGeom->InheritsFrom("MGeomCamMagic"))
    *fLog << " " << setw(7) << "Uncalibrated Pixels:            " 
          << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;

  counts.Reset();

  for (Int_t i=0; i<badcam->GetSize(); i++)
    {
      MBadPixelsPix &bad = (*badcam)[i];
      if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
        {
          const Int_t aidx = (*fGeom)[i].GetAidx();
          counts[aidx]++;
        }
    }

  for (Int_t aidx=0; aidx<nareas; aidx++)
    relcam->SetNumUnreliable(counts[aidx], aidx);

  *fLog << " " << setw(7) << "Unreliable Pixels:              "
        << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;

}

// -----------------------------------------------------------------------------------------------
//
// Print out statistics about BadPixels of type UncalibratedType_t 
// 
void MCalibrationRelTimeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const 
{
  
  MBadPixelsCam         *badcam    = fIntensBad 
    ? (MBadPixelsCam*)        fIntensBad->GetCam()  : fBadPixels;

  UInt_t countinner = 0;
  UInt_t countouter = 0;
  for (Int_t i=0; i<badcam->GetSize(); i++)
    {
      MBadPixelsPix &bad = (*badcam)[i];
      if (bad.IsUncalibrated(typ))
        {
          if (fGeom->GetPixRatio(i) == 1.)
            countinner++;
          else
            countouter++;
        }
    }

  *fLog << " " << setw(7) << text  
        << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl;
}

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

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