/* ======================================================================== *\
!
! *
! * 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): Hendrik Bartko, 09/2004 <mailto:hbartko@mppmu.mpg.de> 
!   Author(s): Markus Gaug, 05/2004 <mailto:markus@ifae.es>
!   Author(s): Diego Tescaro, 05/2004 <mailto:tescaro@pd.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtractTimeAndChargeDigitalFilter
//
//   Hendrik has promised to write more documentation
//
//
//   The following variables have to be set by the derived class and 
//   do not have defaults:
//   - fNumHiGainSamples
//   - fNumLoGainSamples
//   - fSqrtHiGainSamples
//   - fSqrtLoGainSamples
//
// Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MPedestalCam
//
// Output Containers:
//   MArrivalTimeCam
//   MExtractedSignalCam
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractTimeAndChargeDigitalFilter.h"

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

#include "MPedestalPix.h"

#include "TFile.h" // remove!!
#include "TH1F.h"  // remove!!
#include "TString.h"

ClassImp(MExtractTimeAndChargeDigitalFilter);

using namespace std;

const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst  =  0;
const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast   = 14;
const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst  =  3;
const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast   = 14;
const Int_t  MExtractTimeAndChargeDigitalFilter::fgHiGainWindowSize  = 6;
const Int_t  MExtractTimeAndChargeDigitalFilter::fgLoGainWindowSize  = 6;
const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolution = 10;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Calls: 
// - SetWindowSize();
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
// - SetBinningResolution();
//
MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title) 
{

  fName  = name  ? name  : "MExtractTimeAndChargeDigitalFilter";
  fTitle = title ? title : "Digital Filter";

  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
  SetWindowSize();
  SetBinningResolution();

  for (int i=0; i<60; i++){
    fw_amp[i]  = 1;
    fw_time[i] = 1;
  }
}


// ---------------------------------------------------------------------------------------
//
// Checks:
// - if a window is bigger than the one defined by the ranges, set it to the available range
// 
// Sets:
// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
//  
void MExtractTimeAndChargeDigitalFilter::SetWindowSize(Int_t windowh, Int_t windowl)
{

  fWindowSizeHiGain = windowh;
  fWindowSizeLoGain = windowl;
  
  const Int_t availhirange = (Int_t)(fHiGainLast-fHiGainFirst+1);

  if (fWindowSizeHiGain > availhirange)
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",fWindowSizeHiGain,
                    " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
      fHiGainLast = fHiGainFirst + fWindowSizeHiGain;
      *fLog << warn << GetDescriptor() 
            << ": Will set the upper range to: " << (int)fHiGainLast << endl;
    }
  
  if (fWindowSizeHiGain < 2) 
    {
      fWindowSizeHiGain = 2;
      *fLog << warn << GetDescriptor() << ": High Gain window size set to two samples" << endl;
    }
  
  if (fLoGainLast != 0 && fWindowSizeLoGain != 0)
    {
      const Int_t availlorange = (Int_t)(fLoGainLast-fLoGainFirst+1);
      
      if (fWindowSizeLoGain > availlorange)
        {
          *fLog << warn << GetDescriptor() 
                << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",fWindowSizeLoGain,
                        " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
          fLoGainLast = fLoGainFirst + fWindowSizeLoGain;
          *fLog << warn << GetDescriptor() 
            << ": Will set the upper range to: " << (int)fLoGainLast << endl;
        }
      
      if (fWindowSizeLoGain<2) 
        {
          fWindowSizeLoGain = 2;
          *fLog << warn << GetDescriptor() << ": Low Gain window size set to two samples" << endl;
        }
    }
  
  fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
  fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
  
}

//----------------------------------------------------------------------------
//
// Read a pre-defined weights file into the class. 
// This is mandatory for the extraction
//
void MExtractTimeAndChargeDigitalFilter::SetWeightsFile(TString filename)
{

  TFile *f = new TFile(filename);
  TH1F  *hw_amp  = (TH1F*)f->Get("hw_amp");
  TH1F  *hw_time = (TH1F*)f->Get("hw_time");

  if (!hw_amp)
    {
      *fLog << warn << GetDescriptor()
            << ": Could not find hw_amp in " << filename << endl;
      return;
    }
  
  if (!hw_time)
    {
      *fLog << warn << GetDescriptor()
            << ": Could not find hw_time in " << filename << endl;
      return;
    }
  
  for (int i=0; i<fBinningResolution*fWindowSizeHiGain;i++)
    {
    fw_amp [i] = hw_amp->GetBinContent(i+1);
    fw_time[i] = hw_time->GetBinContent(i+1);
    //    y[i]=hshape->GetBinContent(i+1);
    //    derivative[i]=hderivative->GetBinContent(i+1);
  
    //    cout << "w_amp[" << i << "] = " << fw_amp[i] << endl;
  }
  delete f;
}

void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Float_t &dsum, 
                                                                 Float_t &time, Float_t &dtime,
                                                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
{

  const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
  
  Float_t time_sum  = 0.;
  Float_t fmax      = 0.;
  Float_t ftime_max = 0.;
  
  Float_t pedes = ped.GetPedestal();

  const Float_t ABoffs = ped.GetPedestalABoffset();
	
  Float_t PedMean[2];
  PedMean[0] = pedes + ABoffs;
  PedMean[1] = pedes - ABoffs;

  //
  // Calculate the sum of the first fWindowSize slices
  //
  Byte_t *p     = ptr;
  Byte_t *max_p = ptr;

  while (p < end-fWindowSizeHiGain)
    {
      sum      = 0.;
      time_sum = 0.;

      //
      // Slide with a window of size fWindowSizeHiGain over the sample 
      // and multiply the entries with the corresponding weights
      //
      Byte_t *p1 = p;
      for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
      {
        const Int_t   idx = fBinningResolution*sample+fBinningResolutionHalf;
        const Float_t pex = (Float_t)*p1-PedMean[(Int_t(p1-ptr)+abflag) & 0x1];
	sum      += fw_amp [idx]*pex; 
	time_sum += fw_time[idx]*pex;
	p1++;
      }

      if (sum>fmax)
      {
	fmax      = sum;
	ftime_max = time_sum;
	max_p     = p;
      }
      p++;
    } /* while (p<ptr+fWindowSizeHiGain-fHiLoLast) */
  
#if 0
  //
  // If part of the "low-Gain" slices are used, 
  // repeat steps one and two for the logain part until fHiLoLast
  //
  Byte_t *l = logain;
  while (l<logain+fHiLoLast)
    {
      sum      = 0;
      time_sum = 0;

      //
      // Slide with a window of size fWindowSizeHiGain over the sample 
      // and multiply the entries with the corresponding weights
      //
      Byte_t *p1 = l++;
      for (Int_t sample=0; sample < fHiLoLast; sample++)
      {
        const Int_t   idx = fBinningResolution*sample+fBinningResolutionHalf;
        const Float_t pex = (Float_t)*p1-PedMean[(Int_t(p1-logain)+abflag) & 0x1];
	sum      += fw_amp [idx]*pex; 
	time_sum += fw_time[idx]*pex;
	p1++;
      }

      if (sum>fmax)
      {
	fmax      = sum;
	ftime_max = time_sum;
	max_p     = p;
      }

    }
#endif  
  
    if (fmax!=0)
    {
      ftime_max        /= fmax;
      Int_t t_iter      = Int_t(ftime_max*fBinningResolution);
      Int_t sample_iter = 0;
      p                 = max_p;

      while((t_iter)>fBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf)
        {
          if (t_iter>fBinningResolutionHalf-1)
            {
              t_iter -= fBinningResolution;
              p--;
              sample_iter--;
            }
          if (t_iter < -fBinningResolutionHalf)
            {
              t_iter += fBinningResolution;
              p++; 
              sample_iter++;
            }
        }
      
      sum        = 0.;
      Byte_t *p2 = p;
     
      //
      // Slide with a window of size fWindowSizeHiGain over the sample 
      // and multiply the entries with the corresponding weights
      //
      for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
      {
        const Int_t   idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter;
        const Float_t pex = (Float_t)*p2-PedMean[(Int_t(p2-ptr)+abflag) & 0x1];
	sum      += fw_amp [idx]*pex; 
        //        cout << idx << " " << fw_amp[idx] << "  " << pex << "  " << t_iter << "  " << sum << "  " << ((Int_t(p2-ptr)+abflag) & 0x1) << "  " << PedMean[(Int_t(p2-ptr)+abflag) & 0x1] << endl;
	time_sum += fw_time[idx]*pex;
	p2++;
      }

      if (sum != 0)
	time = (Float_t)(max_p - ptr) + 1. + sample_iter 
             - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ; 
      else 
        time = 0.;
    } /* if (max!=0) */ 
    else 
      time=0.;
    
  return;
}

void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum, 
                                                                 Float_t &time, Float_t &dtime,
                                                                 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
{

  const Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;
  
  Float_t time_sum  = 0;
  Float_t fmax      = 0;
  Float_t ftime_max = 0;
  
  Float_t pedes = ped.GetPedestal();

  const Float_t ABoffs = ped.GetPedestalABoffset();
	
  Float_t PedMean[2];
  PedMean[0] = pedes + ABoffs;
  PedMean[1] = pedes - ABoffs;

  //
  // Calculate the sum of the first fWindowSize slices
  //
  sat           = 0;
  Byte_t *p     = ptr;
  Byte_t *max_p = ptr;

  while (p < end-fWindowSizeLoGain)
    {
      sum      = 0;
      time_sum = 0;

      //
      // Slide with a window of size fWindowSizeLoGain over the sample 
      // and multiply the entries with the corresponding weights
      //
      Byte_t *p1 = p++;
      for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
      {
        const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf;
        const Int_t adx = (Int_t(p1-ptr)+abflag) & 0x1;
	sum      += fw_amp [idx]*(*p1-PedMean[adx]); 
	time_sum += fw_time[idx]*(*p1-PedMean[adx]);
	p1++;
      }

      if (sum>fmax)
      {
	fmax      = sum;
	ftime_max = time_sum;
	max_p     = p;
      }
    } /* while (p<end-fWindowSizeLoGain) */
  
    if (fmax!=0)
    {
      ftime_max        /= fmax;
      Int_t t_iter      = Int_t(ftime_max*fBinningResolution);
      Int_t sample_iter = 0;
      p                 = max_p;

      while((t_iter)>fBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf)
        {
          if (t_iter>fBinningResolutionHalf-1)
            {
              t_iter -= fBinningResolution;
              p--;
              sample_iter--;
            }
          if (t_iter < -fBinningResolutionHalf)
            {
              t_iter += fBinningResolution;
              p++; 
              sample_iter++;
            }
        }
      
      sum       = 0;
      Byte_t *p2 = p;
     
      //
      // Slide with a window of size fWindowSizeLoGain over the sample 
      // and multiply the entries with the corresponding weights
      //
      for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
      {
        const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf+t_iter;
        const Int_t adx = (Int_t(p2-ptr)+abflag) & 0x1;
	sum      += fw_amp [idx]*(*p2-PedMean[adx]); 
	time_sum += fw_time[idx]*(*p2-PedMean[adx]);
	p2++;
      }

      if (sum != 0)
	time = (Float_t)(max_p - ptr) + 1. + sample_iter 
             - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ; 
      else 
        time = 0.;
    } /* if (max!=0) */ 
    else 
      time=0.;
    
  return;
}

Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
  
  Byte_t hw = fWindowSizeHiGain;
  Byte_t lw = fWindowSizeLoGain;
  Bool_t rc = kFALSE;
  
  if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
    {
      hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
      rc = kTRUE;
    }
  if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
    {
      lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
      rc = kTRUE;
    }
  
  if (rc)
    SetWindowSize(hw, lw);
    
  if (IsEnvDefined(env, prefix, "BinningResolution", print))
    {
      SetBinningResolution(GetEnvValue(env, prefix, "BinningResolution", fBinningResolution));
      rc = kTRUE;
    }
  
  rc = MExtractor::ReadEnv(env, prefix, print) ? kTRUE : rc;
  
  return rc;
}

