/* ======================================================================== *\
!
! *
! * 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()
//
// Class Version 2:
//  + Byte_t fCheckFlags; // Bit-field to hold the possible check flags
//
// ClassVersionb 3:
//  -  TString fOutputPath;                           //  Path to the output file
//  -  TString fOutputFile;                           //  Name of the output file
//
//
//  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
//
// 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";

  SetCheckFitResults       ( kFALSE );
  SetCheckDeviatingBehavior( kFALSE );
  SetCheckHistOverflow     ( kFALSE );
  SetCheckOscillations     ( kFALSE );

  SetRelTimeResolutionLimit();

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

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

  //
  // 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 resol. less than ", fRelTimeResolutionLimit, " FADC sl. from Mean:"));
  PrintUncalibrated(MBadPixelsPix::kRelTimeOscillating,   
                    "Pixels with changing Rel. Times over time:");
  PrintUncalibrated(MBadPixelsPix::kRelTimeNotFitted,     
                    "Pixels with unsuccesful Gauss fit to the times:");

  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 << "Area   " << setw(4) << aidx << ": Less than 20 pixels with valid time resolution found." << 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 << "Pixel  " << setw(4) << i << ": Deviating time resolution: "
                << Form("%4.2f",res) << " out of range ["
                << Form("%4.2f,%4.2f",lowlim[aidx],upplim[aidx]) << "]" << 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 (IsCheckDeviatingBehavior())
          if (bad.IsUncalibrated(MBadPixelsPix::kDeviatingTimeResolution))
              bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);

      if (IsCheckFitResults())
          if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeNotFitted))
              bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);

      if (IsCheckOscillations())
          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;

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

    const Int_t nareas = fGeom->GetNumAreas();

    TArrayI unsuit(nareas);
    TArrayI unrel(nareas);

    for (int aidx=0; aidx<nareas; aidx++)
    {
        unsuit[aidx] += badcam->GetNumUnsuitable(MBadPixelsPix::kUnsuitableRun, fGeom, aidx);
        unrel[aidx]  += badcam->GetNumUnsuitable(MBadPixelsPix::kUnreliableRun, fGeom, aidx);
        relcam->SetNumUnsuitable(unsuit[aidx], aidx);
        relcam->SetNumUnreliable(unrel[aidx],  aidx);
    }

    if (fGeom->InheritsFrom("MGeomCamMagic"))
    {
        PrintUncalibrated("Uncalibrated Pixels:", unsuit[0], unsuit[1]);
        PrintUncalibrated("Unreliable Pixels:",   unrel[0],  unrel[1]);
    }
}

// -----------------------------------------------------------------------------------------------
//
// Print out statistics about BadPixels of type UncalibratedType_t 
// 
void MCalibrationRelTimeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const 
{
    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++)
    {
        if ((*badcam)[i].IsUncalibrated(typ))
        {
            if (fGeom->GetPixRatio(i) == 1.)
                countinner++;
            else
                countouter++;
        }
    }

    PrintUncalibrated(text, countinner, countouter);
}

void MCalibrationRelTimeCalc::PrintUncalibrated(const char *text, Int_t in, Int_t out) const
{
    *fLog << " " << setfill(' ') << setw(56) << setiosflags(ios::left) << text;
    *fLog << " Inner: " << Form("%3i", in);
    *fLog << " Outer: " << Form("%3i", out);
    *fLog << resetiosflags(ios::left) << endl;

}

// --------------------------------------------------------------------------
//
// MCalibrationRelTimeCam.CheckFitResults: Yes
// MCalibrationRelTimeCam.CheckDeviatingBehavior: Yes
// MCalibrationRelTimeCam.CheckHistOverflow: Yes
// MCalibrationRelTimeCam.CheckOscillations: Yes
// 
Int_t MCalibrationRelTimeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;

    if (IsEnvDefined(env, prefix, "CheckFitResults", print))
    {
        SetCheckFitResults(GetEnvValue(env, prefix, "CheckFitResults", IsCheckFitResults()));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "CheckDeviatingBehavior", print))
    {
        SetCheckDeviatingBehavior(GetEnvValue(env, prefix, "CheckDeviatingBehavior", IsCheckDeviatingBehavior()));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "CheckHistOverflow", print))
    {
        SetCheckHistOverflow(GetEnvValue(env, prefix, "CheckHistOverflow", IsCheckHistOverflow()));
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "CheckOscillations", print))
    {
        SetCheckOscillations(GetEnvValue(env, prefix, "CheckOscillations", IsCheckOscillations()));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "RelTimeResolutionLimit", print))
    {
        SetRelTimeResolutionLimit(GetEnvValue(env, prefix, "RelTimeResolutionLimit", fRelTimeResolutionLimit));
        rc = kTRUE;
    }

    return rc;
}