source: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc@ 5339

/* ======================================================================== *\
!
! *
! * 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 <errno.h>
#include <fstream>

#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TString.h>
#include <TMatrix.h>

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

#include "MPedestalPix.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::fgWindowSizeHiGain  = 6;
const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain  = 6;
const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain = 4;
const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain = 4;

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

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

    ReadWeightsFile("");
}

// ---------------------------------------------------------------------------------------
//
// 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)
{

  if (windowh != fgWindowSizeHiGain)
    *fLog << warn << GetDescriptor() 
          << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default." << endl;
  if (windowl != fgWindowSizeLoGain)
    *fLog << warn << GetDescriptor() 
          << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default" << endl;

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

  if (fWindowSizeHiGain > availhirange)
  {
      // Please simplify this!
      *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)
        {
      // Please simplify this!
          *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;
        }
    }
  
  // 
  // 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 = (Float_t)fWindowSizeHiGain/2.5;
  fNumLoGainSamples = (Float_t)fWindowSizeLoGain/2.5;
  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
  
}

// --------------------------------------------------------------------------
//
// ReInit
//
// Calls:
// - MExtractor::ReInit(pList);
// - Creates new arrays according to the extraction range
//
Bool_t MExtractTimeAndChargeDigitalFilter::ReInit(MParList *pList)
{
    if (!MExtractTimeAndCharge::ReInit(pList))
        return kFALSE;

    fHiGainSignal.Set(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
    fLoGainSignal.Set(fLoGainLast - fLoGainFirst + 1);

    fTimeShiftHiGain = (Float_t)fHiGainFirst + 0.5 + 1./fBinningResolutionHiGain;
    fTimeShiftLoGain = (Float_t)fLoGainFirst + 0.5 + 1./fBinningResolutionLoGain;

    return kTRUE;
}

Int_t MExtractTimeAndChargeDigitalFilter::PreProcess(MParList *pList)
{
    *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;

    return MExtractTimeAndCharge::PreProcess(pList);
}

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)
{

  Int_t range = fHiGainLast - fHiGainFirst + 1;
  const Byte_t *end = ptr + range;
  Byte_t *p     = ptr;
  Byte_t maxpos = 0;
  Byte_t max    = 0;
  Int_t count   = 0;

  //
  // Check for saturation in all other slices
  //
  while (p<end)
    {

      fHiGainSignal[count++] = (Float_t)*p;

      if (*p > max)
        {
          max    = *p;
          maxpos =  p-ptr;
        }

      if (*p++ >= fSaturationLimit)
        sat++;
    }
  
  if (fHiLoLast != 0)
    {

      end = logain + fHiLoLast;

      while (logain<end)
        {

          fHiGainSignal[count++] = (Float_t)*logain;
          range++;

          if (*logain > max)
            {
              max    = *logain;
              maxpos =  range;
            }
          
          if (*logain++ >= fSaturationLimit)
            sat++;
        }
    }
  
  //
  // allow one saturated slice 
  //
  if (sat > 0)
    return;

  Float_t time_sum  = 0.;
  Float_t fmax      = 0.;
  Float_t ftime_max = 0.;
  Int_t   max_p     = 0;
  
  const 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
  //
  for (Int_t i=0;i<range-fWindowSizeHiGain;i++)
    {
      sum      = 0.;
      time_sum = 0.;
      
      //
      // 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 = fBinningResolutionHiGain*sample+fBinningResolutionHalfHiGain;
        const Float_t pex = fHiGainSignal[sample+i]-pedmean[(sample+i+abflag) & 0x1];
        sum              += fAmpWeightsHiGain [idx]*pex; 
        time_sum         += fTimeWeightsHiGain[idx]*pex;
      }

      if (sum>fmax)
      {
        fmax      = sum;
        ftime_max = time_sum;
        max_p     = i;
      }
    } /*   for (Int_t i=0;i<range-fWindowSizeHiGain;i++) */

  time = 0;
  if (fmax==0)
      return;

  ftime_max        /= fmax;
  Int_t t_iter      = Int_t(ftime_max*fBinningResolutionHiGain);
  Int_t sample_iter = 0;

  while ( t_iter > fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain )
    {
      if (t_iter > fBinningResolutionHalfHiGain-1)
        {
          t_iter -= fBinningResolutionHiGain;
          max_p--; 
          sample_iter--;
        }
      if (t_iter < -fBinningResolutionHalfHiGain)
        {
          t_iter += fBinningResolutionHiGain;
          max_p++; 
          sample_iter++;
        }
    }
  
  sum = 0.;
  //
  // 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 = fBinningResolutionHiGain*sample + fBinningResolutionHalfHiGain + t_iter;
    const Int_t   ids = max_p + sample;
    const Float_t pex = ids < 0 ? 0. : 
      ( ids > range ? 0. : fHiGainSignal[ids]-pedmean[(ids+abflag) & 0x1]);
    sum              += fAmpWeightsHiGain [idx]*pex; 
    time_sum         += fTimeWeightsHiGain[idx]*pex;
  }

  if (sum == 0)
      return;

  time = max_p + fTimeShiftHiGain /* this shifts the time to the start of the rising edge */
      - ((Float_t)t_iter)/fBinningResolutionHiGain - time_sum/sum;
}

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)
{

  Int_t range = fLoGainLast - fLoGainFirst + 1;
  const Byte_t *end = ptr + range;
  Byte_t *p     = ptr;
  Byte_t maxpos = 0;
  Byte_t max    = 0;
  Int_t count   = 0;

  //
  // Check for saturation in all other slices
  //
  while (p<end)
    {

      fLoGainSignal[count++] = (Float_t)*p;

      if (*p > max)
        {
          max    = *p;
          maxpos =  p-ptr;
        }

      if (*p++ >= fSaturationLimit)
        sat++;
    }
  
  Float_t time_sum  = 0.;
  Float_t fmax      = 0.;
  Float_t ftime_max = 0.;
  Int_t   max_p     = 0;
  
  const 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
  //
  for (Int_t i=0;i<range-fWindowSizeLoGain;i++)
    {
      sum      = 0.;
      time_sum = 0.;
      
      //
      // 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 = fBinningResolutionLoGain*sample+fBinningResolutionHalfLoGain;
        const Float_t pex = fLoGainSignal[sample+i]-pedmean[(sample+i+abflag) & 0x1];
        sum              += fAmpWeightsLoGain [idx]*pex; 
        time_sum         += fTimeWeightsLoGain[idx]*pex;
      }

      if (sum>fmax)
      {
        fmax      = sum;
        ftime_max = time_sum;
        max_p     = i;
      }
    } /*   for (Int_t i=0;i<range-fWindowSizeLoGain;i++) */

  time = 0;
  if (fmax==0)
      return;

  ftime_max        /= fmax;
  Int_t t_iter      = Int_t(ftime_max*fBinningResolutionLoGain);
  Int_t sample_iter = 0;

  while ( t_iter > fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain )
    {
      if (t_iter > fBinningResolutionHalfLoGain-1)
        {
          t_iter -= fBinningResolutionLoGain;
          max_p--; 
          sample_iter--;
        }
      if (t_iter < -fBinningResolutionHalfLoGain)
        {
          t_iter += fBinningResolutionLoGain;
          max_p++; 
          sample_iter++;
        }
    }
  
  sum = 0.;

  //
  // 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 = fBinningResolutionLoGain*sample + fBinningResolutionHalfLoGain + t_iter;
    const Int_t   ids = max_p + sample;
    const Float_t pex = ids < 0 ? 0. : 
      ( ids > range ? 0. : fLoGainSignal[ids]-pedmean[(ids+abflag) & 0x1]);
    sum              += fAmpWeightsLoGain [idx]*pex; 
    time_sum         += fTimeWeightsLoGain[idx]*pex;
  }

  if (sum == 0)
      return;

  time = max_p + fTimeShiftLoGain /* this shifts the time to the start of the rising edge */
      - ((Float_t)t_iter)/fBinningResolutionLoGain - time_sum/sum;
}

// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//   MJPedestal.MExtractor.WindowSizeHiGain: 6
//   MJPedestal.MExtractor.WindowSizeLoGain: 6
//   MJPedestal.MExtractor.BinningResolutionHiGain: 10
//   MJPedestal.MExtractor.BinningResolutionLoGain: 10
//   MJPedestal.MExtractor.WeightsFile: filename
//
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);

    Bool_t rc2 = kFALSE;
    Int_t brh = fBinningResolutionHiGain;
    Int_t brl = fBinningResolutionLoGain;

    if (IsEnvDefined(env, prefix, "BinningResolutionHiGain", print))
    {
        brh = GetEnvValue(env, prefix, brh);
        rc2 = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "BinningResolutionLoGain", print))
    {
        brl = GetEnvValue(env, prefix, brl);
        rc2 = kTRUE;
    }

    if (rc2)
    {
        SetBinningResolution(brh, brl);
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "WeightsFile", print))
    {
        if (!ReadWeightsFile(GetEnvValue(env, prefix, "WeightsFile", "")))
            return kERROR;
        rc = kTRUE;
    }

    return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
}

//----------------------------------------------------------------------------
//
// 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.
//
Bool_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename)
{
    fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain);
    fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain);
    fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain);
    fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain);

    if (filename.IsNull())
    {
        fAmpWeightsHiGain.Reset(1);
        fTimeWeightsHiGain.Reset(1);
        fAmpWeightsLoGain.Reset(1);
        fTimeWeightsLoGain.Reset(1);
        return kTRUE;
    }

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

    *fLog << inf << "Reading weights file " << 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 kFALSE;
            }

            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 kFALSE;
            }

            len = fBinningResolutionHiGain*fWindowSizeHiGain;
            fAmpWeightsHiGain .Set(len);
            fTimeWeightsHiGain.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 kFALSE;
            }

            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 kFALSE;
            }

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

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

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

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

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

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

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

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

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

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

    return kTRUE;
}

//----------------------------------------------------------------------------
//
// Create the weights file
// Beware that the shape-histogram has to contain the pulse starting at bin 1
//
Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi,
                                                            TH1F *shapelo, TH2F *autocorrlo )
{

  const Int_t nbinshi = shapehi->GetNbinsX();
  Float_t binwidth    = shapehi->GetBinWidth(1);

  TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"),
                                Form("%s%s",shapehi->GetTitle()," derivative"),
                                nbinshi,
                                shapehi->GetBinLowEdge(1),
                                shapehi->GetBinLowEdge(nbinshi)+binwidth);

  //
  // Calculate the derivative of shapehi
  // 
  for (Int_t i = 1; i<nbinshi+1;i++)
    {
      derivativehi->SetBinContent(i,
                                ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth));
      derivativehi->SetBinError(i,
                              (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1)
                                    +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth));
    }
  
  //
  // normalize the shapehi, such that the integral for fWindowSize slices is one!
  //
  Float_t sum    = 0;
  Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain);
  lasttemp       = lasttemp > nbinshi ? nbinshi : lasttemp;
  
  for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) {
    sum += shapehi->GetBinContent(i);
  }
  sum /= fBinningResolutionHiGain;

  shapehi->Scale(1./sum);
  derivativehi->Scale(1./sum);

  //
  // read in the noise auto-correlation function:
  //
  TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain);

  for (Int_t i=0; i<fWindowSizeHiGain; i++){
    for (Int_t j=0; j<fWindowSizeHiGain; j++){
      Bhi[i][j]=autocorrhi->GetBinContent(i+1,j+1); //+fSignalStartBinHiGain +fSignalStartBinHiGain
    }
  }  
  Bhi.Invert();

  const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain;
  fAmpWeightsHiGain.Set(nsizehi);
  fTimeWeightsHiGain.Set(nsizehi);  
  
  //
  // Loop over relative time in one BinningResolution interval
  //
  Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);

  for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++)
    {
  
      TMatrix g(fWindowSizeHiGain,1);
      TMatrix gT(1,fWindowSizeHiGain);
      TMatrix d(fWindowSizeHiGain,1);
      TMatrix dT(1,fWindowSizeHiGain);
    
      for (Int_t count=0; count < fWindowSizeHiGain; count++){
      
        g[count][0]=shapehi->GetBinContent(start
                                         +fBinningResolutionHiGain*count+i); 
        gT[0][count]=shapehi->GetBinContent(start
                                          +fBinningResolutionHiGain*count+i);
        d[count][0]=derivativehi->GetBinContent(start
                                              +fBinningResolutionHiGain*count+i);
        dT[0][count]=derivativehi->GetBinContent(start
                                               +fBinningResolutionHiGain*count+i);
      }
    
      TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g));
      Float_t   denom = m_denom[0][0];  // ROOT thinks, m_denom is still a matrix
      
      TMatrix m_first = dT*(Bhi*d);       // ROOT thinks, m_first is still a matrix
      Float_t   first = m_first[0][0]/denom;
      
      TMatrix m_last  = gT*(Bhi*d);       // ROOT thinks, m_last  is still a matrix
      Float_t   last  = m_last[0][0]/denom;
      
      TMatrix m1 = gT*Bhi;
      m1 *= first;
      
      TMatrix m2 = dT*Bhi; 
      m2 *=last;
      
      TMatrix w_amp = m1 - m2;
      
      TMatrix m_first1 = gT*(Bhi*g);
      Float_t   first1 = m_first1[0][0]/denom;
      
      TMatrix m_last1  = gT*(Bhi*d);
      Float_t   last1  = m_last1 [0][0]/denom;
      
      TMatrix m11 = dT*Bhi; 
      m11 *=first1;
      
      TMatrix m21 = gT*Bhi;
      m21 *=last1;
      
      TMatrix w_time= m11 - m21; 
      
      for (Int_t count=0; count < fWindowSizeHiGain; count++)
        {
          const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1;
          fAmpWeightsHiGain [idx] = w_amp [0][count];
          fTimeWeightsHiGain[idx] = w_time[0][count];
        }
      
    } // end loop over i

  //
  // Low Gain histograms
  //
  TH1F *derivativelo = NULL;
  if (shapelo)
    {
      const Int_t nbinslo  = shapelo->GetNbinsX();
      binwidth = shapelo->GetBinWidth(1);
      
      derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"),
                              Form("%s%s",shapelo->GetTitle()," derivative"),
                              nbinslo,
                              shapelo->GetBinLowEdge(1),
                              shapelo->GetBinLowEdge(nbinslo)+binwidth);
      
      //
      // Calculate the derivative of shapelo
      // 
      for (Int_t i = 1; i<nbinslo+1;i++)
        {
          derivativelo->SetBinContent(i,
                                      ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth));
          derivativelo->SetBinError(i,
                                    (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1)
                                          +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth));
        }
      
      //
      // normalize the shapelo, such that the integral for fWindowSize slices is one!
      //
      sum      = 0;
      lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain);
      lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp;
      
      for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++) 
        sum += shapelo->GetBinContent(i);

      sum /= fBinningResolutionLoGain;
      
      shapelo->Scale(1./sum);
      derivativelo->Scale(1./sum);
      
      //
      // read in the noise auto-correlation function:
      //
      TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain);
      
      for (Int_t i=0; i<fWindowSizeLoGain; i++){
        for (Int_t j=0; j<fWindowSizeLoGain; j++){
          Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain);
        }
      }  
      Blo.Invert();

      const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain;
      fAmpWeightsLoGain.Set(nsizelo);
      fTimeWeightsLoGain.Set(nsizelo);  
  
      //
      // Loop over relative time in one BinningResolution interval
      //
      Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain;
      
      for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++)
        {
          
          TMatrix g(fWindowSizeLoGain,1);
          TMatrix gT(1,fWindowSizeLoGain);
          TMatrix d(fWindowSizeLoGain,1);
          TMatrix dT(1,fWindowSizeLoGain);
          
          for (Int_t count=0; count < fWindowSizeLoGain; count++){
            
            g[count][0] = shapelo->GetBinContent(start
                                             +fBinningResolutionLoGain*count+i); 
            gT[0][count]= shapelo->GetBinContent(start
                                              +fBinningResolutionLoGain*count+i);
            d[count][0] = derivativelo->GetBinContent(start
                                                  +fBinningResolutionLoGain*count+i);
            dT[0][count]= derivativelo->GetBinContent(start
                                                   +fBinningResolutionLoGain*count+i);
          }
          
          TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g));
          Float_t   denom = m_denom[0][0];  // ROOT thinks, m_denom is still a matrix
          
          TMatrix m_first = dT*(Blo*d);       // ROOT thinks, m_first is still a matrix
          Float_t   first = m_first[0][0]/denom;
          
          TMatrix m_last  = gT*(Blo*d);       // ROOT thinks, m_last  is still a matrix
          Float_t   last  = m_last[0][0]/denom;
          
          TMatrix m1 = gT*Blo;
          m1 *= first;
          
          TMatrix m2 = dT*Blo; 
          m2 *=last;
          
          TMatrix w_amp = m1 - m2;
          
          TMatrix m_first1 = gT*(Blo*g);
          Float_t   first1 = m_first1[0][0]/denom;
          
          TMatrix m_last1  = gT*(Blo*d);
          Float_t   last1  = m_last1 [0][0]/denom;
          
          TMatrix m11 = dT*Blo; 
          m11 *=first1;
          
          TMatrix m21 = gT*Blo;
          m21 *=last1;
          
          TMatrix w_time= m11 - m21; 
          
          for (Int_t count=0; count < fWindowSizeLoGain; count++)
            {
              const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1;
              fAmpWeightsLoGain [idx] = w_amp [0][count];
              fTimeWeightsLoGain[idx] = w_time[0][count];
            }
      
        } // end loop over i
    }
  
  ofstream fn(filename.Data());

  fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl;
  fn << "# (Amplitude)  (Time) " << endl;

  for (Int_t i=0; i<nsizehi; i++)
    fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl;

  fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl;
  fn << "# (Amplitude)  (Time) " << endl;

  for (Int_t i=0; i<nsizehi; i++)
    fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl;

  delete derivativehi;
  if (derivativelo)
    delete derivativelo;

  return kTRUE;
}

void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
{
    *fLog << all;
    *fLog << GetDescriptor() << ":" << endl;

    *fLog << " Time Shift HiGain:  " << fTimeShiftHiGain << endl;
    *fLog << " Time Shift LoGain:  " << fTimeShiftLoGain << endl;
    *fLog << " Window Size HiGain: " << fWindowSizeHiGain << endl;
    *fLog << " Window Size LoGain: " << fWindowSizeLoGain << endl;
    *fLog << " Binning Res HiGain: " << fBinningResolutionHiGain << endl;
    *fLog << " Binning Res LoGain: " << fBinningResolutionHiGain << endl;
  
    MExtractTimeAndCharge::Print(o);
}