/* ======================================================================== *\
!
! *
! * 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 "MGeomCam.h"
#include "MGeomPix.h"

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

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


ClassImp(MCalibrationRelTimeCalc);

using namespace std;

const Float_t MCalibrationRelTimeCalc::fgRelTimeRelErrLimit      = 15.;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets all pointers to NULL
// 
// Initializes:
// - fRelTimeRelErrLimit to fgRelTimeRelErrimit
//
// Calls:
// - Clear()
//
MCalibrationRelTimeCalc::MCalibrationRelTimeCalc(const char *name, const char *title)
    : fBadPixels(NULL), fCam(NULL), fGeom(NULL)
{
        
  fName  = name  ? name  : "MCalibrationRelTimeCalc";
  fTitle = title ? title : "Task to finalize the relative time calibration";
  
  SetRelTimeRelErrLimit();
 
  Clear();
  
}

// --------------------------------------------------------------------------
//
// Sets:
// - all flags to kFALSE
//
void MCalibrationRelTimeCalc::Clear(const Option_t *o)
{
    SkipHiLoGainCalibration( kFALSE );    
}


// -----------------------------------------------------------------------------------
//
// The following containers are searched and created if they were not found:
//
//  - MBadPixelsCam
//
Int_t MCalibrationRelTimeCalc::PreProcess(MParList *pList)
{
  

  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
//  - MCalibrationRelTimeCam
// 
// 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;
    }
  
  fCam = (MCalibrationRelTimeCam*)pList->FindObject("MCalibrationRelTimeCam");
  if (!fCam)
    {
      *fLog << err << "Cannot find MCalibrationRelTimeCam... aborting" << endl;
      *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationRelTimeCam before..." << endl;
      return kFALSE;
    }

  
  UInt_t npixels     = fGeom->GetNumPixels();
  
  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*fCam)  [i];
      MBadPixelsPix         &bad = (*fBadPixels)[i];
      
      pix.SetPixId(i);
      
      if (bad.IsBad())
        {
          pix.SetExcluded();
          continue;
        }
      
    }

  return kTRUE;
}

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

// -----------------------------------------------------------------------
//
// Return if number of executions is null.
//
// 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::PostProcess()
{

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

  //
  // First loop over pixels, call FinalizePedestals and FinalizeRelTimes
  //
  FinalizeRelTimes();

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

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

  fCam      ->SetReadyToSave();
  fBadPixels->SetReadyToSave();

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

  PrintUncalibrated(MBadPixelsPix::kDeviatingTimeResolution,    
                    Form("%s%2.1f%s","Time resolution less than ",fRelTimeRelErrLimit," sigma 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 fRelTimeRelErrLimit mean 
//              sigmas from the mean (obtained in first loop). Set 
//              MBadPixelsPix::kDeviatingTimeResolution if excluded.
// 
void MCalibrationRelTimeCalc::FinalizeRelTimes()
{

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

  Float_t lowlim      [nareas];
  Float_t upplim      [nareas];
  Float_t areasum     [nareas];
  Float_t areasum2    [nareas];
  Int_t   numareavalid[nareas];

  memset(lowlim        ,0, nareas   * sizeof(Float_t));
  memset(upplim        ,0, nareas   * sizeof(Float_t));
  memset(areasum       ,0, nareas   * sizeof(Float_t));
  memset(areasum2      ,0, nareas   * sizeof(Float_t));
  memset(numareavalid  ,0, nareas   * sizeof(Int_t  ));

  //
  // 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&)(*fCam)[i];
      MBadPixelsPix          &bad = (*fBadPixels)[i];
      
      if (pix.IsExcluded())
        continue;

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

      if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
        continue;
      
      const Float_t res   = pix.GetTimePrecision();
      const Int_t   aidx  = (*fGeom)[i].GetAidx();

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


  for (UInt_t aidx=0; aidx<nareas; aidx++)
    {
      if (numareavalid[aidx] == 0)
        {
          *fLog << warn << GetDescriptor() << ": No pixels with valid time resolution found "
                << "in area index: " << aidx << endl;
          continue;
        }
 
      *fLog << err << areasum2[aidx] << "  " << areasum[aidx] << "  " << numareavalid[aidx] << endl;

      // 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]  = areasum [aidx] - fRelTimeRelErrLimit*areasum2[aidx];
      upplim  [aidx]  = areasum [aidx] + fRelTimeRelErrLimit*areasum2[aidx];
    }
  


  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*fCam)[i];
      MBadPixelsPix          &bad = (*fBadPixels)[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 << GetDescriptor() << ": Deviating time resolution: " 
            << Form("%4.2f",res) << " out of accepted limits: [" 
                << Form("%4.2f%s%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()
{
  
  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    {
      
      MBadPixelsPix    &bad    = (*fBadPixels)[i];
      MCalibrationPix  &pix    = (*fCam)[i];

      if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingTimeResolution))
        bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun   );
 
      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 fCam
//
void MCalibrationRelTimeCalc::FinalizeUnsuitablePixels()
{
  
  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Rel. Times Calibration status:" << endl;
  *fLog << dec << setfill(' ');

  const Int_t nareas = fGeom->GetNumAreas();

  Int_t counts[nareas];
  memset(counts,0,nareas*sizeof(Int_t));

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

  for (Int_t aidx=0; aidx<nareas; aidx++)
    fCam->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;

  memset(counts,0,nareas*sizeof(Int_t));

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

  for (Int_t aidx=0; aidx<nareas; aidx++)
    fCam->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 
{
  
  UInt_t countinner = 0;
  UInt_t countouter = 0;
  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    {
      MBadPixelsPix &bad = (*fBadPixels)[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;
}