/* ======================================================================== *\
!
! *
! * 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
//
// The reading of automatic weights files (color, type) can be switched
// off using EnableAutomaticWeights(kFALSE).
//
// An empty name or "-" as the weights file name is a synonym for
// setting all weights to 1
//
// The digital filter can even do a little extrapolation around the
// extraction window. For real Pulses having their maximum around
// slice 0 the extracted time is similar to a gaussian around 0 with
// sigma 0.4.
//
//
// Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MPedestalCam
//   [MCalibrationPattern]
//
// Output Containers:
//   MArrivalTimeCam
//   MExtractedSignalCam
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractTimeAndChargeDigitalFilter.h"

#include <errno.h>
#include <fstream>

#include <TRandom.h>

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

#include "MParList.h"

#include "MRawRunHeader.h"
#include "MCalibrationPattern.h"
#include "MExtractedSignalCam.h"
#include "MExtralgoDigitalFilter.h"

ClassImp(MExtractTimeAndChargeDigitalFilter);

using namespace std;

const Byte_t  MExtractTimeAndChargeDigitalFilter::fgHiGainFirst             =  0;
const Byte_t  MExtractTimeAndChargeDigitalFilter::fgHiGainLast              = 16;
const Int_t   MExtractTimeAndChargeDigitalFilter::fgLoGainFirst             =  1;
const Byte_t  MExtractTimeAndChargeDigitalFilter::fgLoGainLast              = 14;
const Int_t   MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
const Int_t   MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain            =  0.95;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Calls: 
// - SetWindowSize();
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
// - SetBinningResolution();
//
// Sets all weights to 1.
//
MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title) 
    : fBinningResolutionHiGain(fgBinningResolutionHiGain),
    fBinningResolutionLoGain(fgBinningResolutionLoGain),
    fAutomaticWeights(kTRUE)
{
    fName  = name  ? name  : "MExtractTimeAndChargeDigitalFilter";
    fTitle = title ? title : "Digital Filter";

    SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
    SetWindowSize(3, 5);
    SetOffsetLoGain(fgOffsetLoGain); 
}

// ---------------------------------------------------------------------------------------
//
// 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
//
// This function might be used to turn the digital filter into a
// sliding window extractor by setting the filename to NULL
//
void MExtractTimeAndChargeDigitalFilter::SetWindowSize(Int_t windowh, Int_t windowl)
{
    if (windowh > fHiGainLast-fHiGainFirst+1)
    {
        *fLog << err << "ERROR - The new hi-gain window size exceeds the extraction range." << endl;
        return;
    }

    if (windowl > (int)fLoGainLast-fLoGainFirst+1)
    {
        *fLog << err << "ERROR - The new hi-gain window size exceeds the extraction range." << endl;
        return;
    }

    fWindowSizeHiGain = windowh;
    fWindowSizeLoGain = windowl;

  /*
  const Int_t availhirange = (Int_t)();

  if (fWindowSizeHiGain > availhirange)
  {
      *fLog << warn << GetDescriptor() << ": Hi Gain window size: " << Form("%2i",fWindowSizeHiGain);
      *fLog << " is bigger than available range: [" << Form("%2i", (int)fHiGainFirst);
      *fLog << "," << Form("%21", (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() << ": Lo Gain window size: " << Form("%2i",fWindowSizeLoGain);
            *fLog << " is bigger than available range: [" << Form("%2i", (int)fLoGainFirst);
            *fLog << "," << Form("%21", (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;
        }
    }*/
  // 
  // We need here the effective number of samples which is about 2.5 in the case of a window
  // size of 6. The exact numbers have to be found still.
  //
  fNumHiGainSamples  = fWindowSizeHiGain;
  fNumLoGainSamples  = fWindowSizeLoGain;
  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
}


// --------------------------------------------------------------------------
//
// Executing MExtractTimeAndCharge::PreProcess and searching for
// MCalibrationPattern
//
Int_t MExtractTimeAndChargeDigitalFilter::PreProcess(MParList *pList)
{
    if (!MExtractTimeAndCharge::PreProcess(pList))
        return kFALSE;

    fCalibPattern = (MCalibrationPattern*)pList->FindObject("MCalibrationPattern");
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// The weights are determined using GetAutimaticWeights().
//
// kFALSE is returned if it returned an error.
// kTRUE  is returned if no new weights were set.
//
// If new weights are set
//  fNumHiGainSamples
//  fNumLoGainSamples
//  fSqrtHiGainSamples
//  fSqrtLoGainSamples
// and
//  fSignals->SetUsedFADCSlices(...)
// is updated accordingly.
//
Bool_t MExtractTimeAndChargeDigitalFilter::GetWeights()
{
    switch (GetAutomaticWeights())
    {
    case kERROR: // An error occured
        return kFALSE;
    case kFALSE: // No new weights set
        return kTRUE;
    }

    //
    // We need here the effective number of samples. In pricipal the number
    // is different depending on the weights used and must be set
    // event by event.
    //
    fNumHiGainSamples  = fAmpWeightsHiGain.GetSum()/fBinningResolutionHiGain;
    fNumLoGainSamples  = fAmpWeightsLoGain.GetSum()/fBinningResolutionLoGain;
    fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
    fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);

    // From MExtractTimeAndCharge::ReInit
    if (fSignals)
        fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples,
                                    fLoGainFirst, fLoGainLast, fNumLoGainSamples);
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// InitArrays
//
// Gets called in the ReInit() and initialized the arrays
//
Bool_t MExtractTimeAndChargeDigitalFilter::InitArrays(Int_t n)
{
    if (!fRunHeader)
        return kFALSE;

    return GetWeights();
}

// --------------------------------------------------------------------------
//
// Check if reading a new weights file is necessary because the calibration
// pattern has changed. (Cannot be done in ReInit, because at this time
// the calibration pattern is not available.
// Then process the event.
//
Int_t MExtractTimeAndChargeDigitalFilter::Process()
{
    // Change Weights if the calibration patter changes
    if (!GetWeights())
        return kERROR;

    // Process event
    return MExtractTimeAndCharge::Process();
}

// --------------------------------------------------------------------------
//
// Apply the digital filter algorithm to the high-gain slices.
//
void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain2(const Float_t *ptr, Int_t num,
                                                                  Float_t &sum, Float_t &dsum,
                                                                  Float_t &time, Float_t &dtime,
                                                                  Byte_t sat, Int_t maxpos) const
{
    // Do some handling if maxpos is last slice!

    MExtralgoDigitalFilter df(fBinningResolutionHiGain, fWindowSizeHiGain,
                              fAmpWeightsHiGain.GetArray(),
                              fTimeWeightsHiGain.GetArray(),
                              fPulseHiGain.GetArray());
    df.SetData(num, ptr);

    if (IsNoiseCalculation())
    {
        sum = df.ExtractNoise();
        return;
    }

    df.Extract(/*maxpos*/);
    df.GetSignal(sum, dsum);
    df.GetTime(time, dtime);
}

void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain2(const Float_t *ptr, Int_t num,
                                                                  Float_t &sum, Float_t &dsum,
                                                                  Float_t &time, Float_t &dtime,
                                                                  Byte_t sat, Int_t maxpos) const
{
    MExtralgoDigitalFilter df(fBinningResolutionLoGain, fWindowSizeLoGain,
                              fAmpWeightsLoGain.GetArray(),
                              fTimeWeightsLoGain.GetArray(),
                              fPulseLoGain.GetArray());

    df.SetData(num, ptr);

    if (IsNoiseCalculation())
    {
        sum = df.ExtractNoise();
        return;
    }

    df.Extract(/*maxpos*/);
    df.GetSignal(sum, dsum);
    df.GetTime(time, dtime);
}


// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//   MJPedestal.MExtractor.WeightsFile: filename
//   MJPedestal.MExtractor.AutomaticWeights: off
//
Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{

  Bool_t rc = kFALSE;

  if (IsEnvDefined(env, prefix, "AutomaticWeights", print))
  {
      EnableAutomaticWeights(GetEnvValue(env, prefix, "AutomaticWeights", fAutomaticWeights));
      rc = kTRUE;
  }

  if (IsEnvDefined(env, prefix, "WeightsFile", print))
    {
        SetNameWeightsFile(GetEnvValue(env, prefix, "WeightsFile", ""));
        rc = kTRUE;
    }
  
  return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
}

//----------------------------------------------------------------------------
//
// If automatic weights are requested, no default weights (name.IsNull())
// are requested, fRunHeader is available and fRunHeader->IsMonteCarloRun()
// is true prepend "MC_" in from of the name.
//
// return poth+name;
//
TString MExtractTimeAndChargeDigitalFilter::CompileWeightFileName(TString path, const TString &name) const
{
    if (fAutomaticWeights && !name.IsNull() && fRunHeader && fRunHeader->IsMonteCarloRun())
        path += "MC_";

    path += name;

    return path;
}

//----------------------------------------------------------------------------
//
// Read a pre-defined weights file into the class. 
// This is mandatory for the extraction
//
// If filenname is empty, then all weights will be set to 1.
//
// Returns:
//  kTRUE:  new weights set
//  kFALSE: no  weights set
//  kERROR: error
//
Int_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename, TString path)
{
    if (filename.IsNull())
    {
        fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain);
        fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain);
        fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain);
        fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain);

        fAmpWeightsHiGain.Reset(1);
        fTimeWeightsHiGain.Reset(1);
        fAmpWeightsLoGain.Reset(1);
        fTimeWeightsLoGain.Reset(1);
        return kTRUE;
    }

    // Add "MC_" in front of the filename if necessary
    filename = CompileWeightFileName(path, filename);

    //filename = MJob::ExpandPath(filename);

    if (fNameWeightsFileSet==filename)
        return kFALSE; // No file read

    ifstream fin(filename.Data());
    if (!fin)
    {
        *fLog << err << GetDescriptor() << ": ERROR - Cannot open file " << filename << ": ";
        *fLog << strerror(errno) << endl;
        return kERROR;
    }

    *fLog << all << GetDescriptor() << ": Reading weights in " << filename << "..." << flush;

    Int_t len = 0;
    Int_t cnt = 0;
    Int_t line = 0;
    Bool_t hi = kFALSE;
    Bool_t lo = kFALSE;

    TString str;

    while (1)
    {
        str.ReadLine(fin);
        if (!fin)
            break;

        line++;

        if (str.Contains("# High Gain Weights:"))
        {
            if (hi)
            {
                *fLog << err << "ERROR - 'High Gain Weights' found twice in line #" << line << "." << endl;
                return kERROR;
            }

            if (2!=sscanf(str.Data(), "# High Gain Weights: %2i %2i", &fWindowSizeHiGain, &fBinningResolutionHiGain))
            {
                *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
                *fLog << str << endl;
                return kERROR;
            }

            len = fBinningResolutionHiGain*fWindowSizeHiGain;
            fAmpWeightsHiGain .Set(len);
            fTimeWeightsHiGain.Set(len);
            fPulseHiGain.Set(len);
            hi = kTRUE;
            continue;
        }

        if (str.Contains("# Low Gain Weights:"))
        {
            if (lo)
            {
                *fLog << err << "ERROR - 'Lo Gain Weights' found twice in line #" << line << "." << endl;
                return kERROR;
            }

            if (2!=sscanf(str.Data(),"# Low Gain Weights: %2i %2i", &fWindowSizeLoGain, &fBinningResolutionLoGain))
            {
                *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
                *fLog << str << endl;
                return kERROR;
            }

            len = fBinningResolutionLoGain*fWindowSizeLoGain;
            fAmpWeightsLoGain .Set(len);
            fTimeWeightsLoGain.Set(len);
            fPulseLoGain.Set(len);
            lo = kTRUE;
            continue;
        }

        // Handle lines with comments
        if (str.Contains("#"))
            continue;

        // Nothing found so far
        if (len == 0)
            continue;

        if (3!=sscanf(str.Data(), "%f %f %f",
                      lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt],
                      lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt],
                      lo ? &fPulseLoGain[cnt] : &fPulseHiGain[cnt]))
        {
            *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
            *fLog << str << endl;
            return kERROR;
        }

        if (++cnt == len)
        {
            len = 0;
            cnt = 0;
        }
    }

    if (cnt != len)
    {
        *fLog << err << "ERROR - Size mismatch in weights file " << filename << endl;
        return kERROR;
    }

    if (!hi)
    {
        *fLog << err << "ERROR - No correct header found in weights file " << filename << endl;
        return kERROR;
    }

    *fLog << "done." << endl;

    *fLog << inf << " File contains " << fWindowSizeHiGain << " hi-gain slices ";
    *fLog << "with a resolution of " << fBinningResolutionHiGain << endl;

    *fLog << inf << " File contains " << fWindowSizeLoGain << " lo-gain slices ";
    *fLog << "with a resolution of " << fBinningResolutionLoGain << endl;

    //CalcBinningResArrays();

    switch (fWindowSizeHiGain)
      {
      case 4:
	SetResolutionPerPheHiGain(0.036);
	break;
      case 6:
	SetResolutionPerPheHiGain(0.021);
	break;
      default:
	*fLog << warn << "Could not set the high-gain extractor resolution per phe for window size " 
	      << fWindowSizeHiGain << endl;
      }

    switch (fWindowSizeLoGain)
      {
      case 4:
	SetResolutionPerPheLoGain(0.005);
	break;
      case 6:
	SetResolutionPerPheLoGain(0.004);
	break;
      default:
	*fLog << warn << "Could not set the low-gain extractor resolution per phe for window size " 
	      << fWindowSizeLoGain << endl;
      }

    fNameWeightsFileSet = filename;

    return kTRUE;
}


//----------------------------------------------------------------------------
//
// The default (+ prepending possible "MC_") is read for:
//
//   - RunType: Pedestal (independant of fAutomaticWeights)
//   - fAutomaticWeights disabled
//
//  if fAutomaticWeights enabled:
//   - fNameWeightsFile.IsNull()
//   - !fCalibPattern
//   - fCalibPattern->GetPulserColor()==MCalibrationCam::kNONE
//
// If automatic weights are enabled, the case above didn't take place and
// fNameWeightsFile starts with "calibration_weights_"
//   - the color (blue, UV) is replaced by the appropriate one
//     taken from the calibration pattern
//
// In most cases a debug output is printed. Further output about the color
// determination can be switched on with debug level > 5;
//
// Returns:
//  kFALSE: No new weights set
//  kTRUE:  New weights set
//  kERROR: Error
//
Int_t MExtractTimeAndChargeDigitalFilter::GetAutomaticWeights()
{
    const Ssiz_t pos = fNameWeightsFile.Last('/')+1;
    const Ssiz_t len = fNameWeightsFile.Length();

    // Split file name in path and name
    TString path = fNameWeightsFile(0, pos>=0?pos:len);
    TString name = fNameWeightsFile(pos>=0?pos:0, len);

    // Remove trailing "MC_" for automatic weights
    if (fAutomaticWeights && name.BeginsWith("MC_"))
        name.Remove(0, 3);

    // In case of a pedetsal run no calibration pattern can be available
    // the default weights are always used.
    if (fRunHeader->GetRunType()==MRawRunHeader::kRTPedestal)
    {
        *fLog << dbg << "Pedestal file... using default weights: " << fNameWeightsFile << endl;
        return ReadWeightsFile(name, path);
    }

    // If automatic weights are switched off use default weights
    if (!fAutomaticWeights)
    {
        *fLog << dbg << "Automatic weights switched off... using default weights: " << fNameWeightsFile << endl;
        return ReadWeightsFile(name, path);
    }

    // If automatic weights are switched on but no filename is given raise error
    if (fNameWeightsFile.IsNull())
    {
        *fLog << err << "ERROR - Cannot get automatic weights without default filename." << endl;
        return kERROR;
    }

    // If this is no pedestal run, automatic weights are requested and a
    // filename for the weights file is given pedestal-extraction from
    // cosmics data is assumed.
    if (!fCalibPattern)
    {
        *fLog << dbg << "No decoded calibration pattern available... using default weights: " << fNameWeightsFile << endl;
        return ReadWeightsFile(name, path);
    }

    const Bool_t debug = gLog.GetDebugLevel()>5;

    // If no calibration pattern is available do not change the
    // current weighs or current weights file name.
    if (fCalibPattern->GetPulserColor()==MCalibrationCam::kNONE)
    {
        // If we are extracting data and the calibration pattern is kNONE
        // we assume that it is a data file without interleaved events
        // and calibration pattern information available.
        if ((fRunHeader->GetRunType()!=MRawRunHeader::kRTData && !fRunHeader->IsMonteCarloRun()) || debug)
            *fLog << dbg << "No calibration color set so far... guessing default: " << fNameWeightsFile << endl;

        return ReadWeightsFile(name, path);
    }

    if (debug)
    {
        *fLog << dbg << endl;
        *fLog << underline << GetDescriptor() << endl;
        *fLog << " Trying to get automatic weight for " << fNameWeightsFile << endl;
        *fLog << " Run type: ";
    }

    if (name.BeginsWith("calibration_weights_") && fCalibPattern)
    {
        if (debug)
            *fLog << " Calibration with color " << fCalibPattern->GetPulserColorStr() << ", setting ";
        switch (fCalibPattern->GetPulserColor())
        {
        case MCalibrationCam::kBLUE:  // 2
        case MCalibrationCam::kGREEN: // 1
            if (debug)
                *fLog << "blue/green, ";
            name.ReplaceAll("UV", "blue");
            break;

        case MCalibrationCam::kUV:    // 3
        case MCalibrationCam::kCT1:   // 0
            if (debug)
                *fLog << "UV/CT1, ";
            name.ReplaceAll("blue", "UV");
            break;
        case MCalibrationCam::kNONE:
            break;
        default: // kNone + etc
            *fLog << err << "ERROR - Cannot get automatic weights for " << fCalibPattern->GetPulserColorStr() << endl;
            return kERROR;
        }
    }

    return ReadWeightsFile(name, path);
}

//----------------------------------------------------------------------------
//
// Print the setup of the digital filter extraction used. Use
//  the option "weights" if you want to print also all weights.
//
void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
{
    if (IsA()==Class())
        *fLog << GetDescriptor() << ":" << endl;

    MExtractTimeAndCharge::Print(o);
    *fLog << " Window Size HiGain: " << setw(2) << fWindowSizeHiGain        << "  LoGain: " << setw(2) << fWindowSizeLoGain << endl;
    *fLog << " Binning Res HiGain: " << setw(2) << fBinningResolutionHiGain << "  LoGain: " << setw(2) << fBinningResolutionHiGain << endl;
    *fLog << " Weights File desired: " << (fNameWeightsFile.IsNull()?"-":fNameWeightsFile.Data()) << endl;
    if (!fNameWeightsFileSet.IsNull())
        *fLog << " Weights File set:     " << fNameWeightsFileSet << endl;

    TString opt(o);
    if (!opt.Contains("weights"))
        return;

    *fLog << endl;
    *fLog << inf << "Using the following weights: " << endl;
    *fLog << "Hi-Gain:" << endl;
    for (Int_t i=0; i<fBinningResolutionHiGain*fWindowSizeHiGain; i++)
        *fLog << " " << fAmpWeightsHiGain[i] << " \t " << fTimeWeightsHiGain[i] << endl;

    *fLog << "Lo-Gain:" << endl;
    for (Int_t i=0; i<fBinningResolutionLoGain*fWindowSizeLoGain; i++)
        *fLog << " " << fAmpWeightsLoGain[i] << " \t " << fTimeWeightsLoGain[i] << endl;
}