source: trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc@ 4658

/* ======================================================================== *\
!
! *
! * 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 analyzing 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       05/2004 <mailto:markus@ifae.es> 
!
!   Copyright: MAGIC Software Development, 2002-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtractTimeAndChargeSpline
//
//   Fast Spline extractor using a cubic spline algorithm of Numerical Recipes. 
//   It returns the integral below the interpolating spline. 
// 
//   Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) 
//         to modify the ranges. Ranges have to be an even number. In case of odd 
//         ranges, the last slice will be reduced by one.
//
//  Defaults are: 
// 
//   fHiGainFirst =  fgHiGainFirst =  3 
//   fHiGainLast  =  fgHiGainLast  =  14
//   fLoGainFirst =  fgLoGainFirst =  3 
//   fLoGainLast  =  fgLoGainLast  =  14
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractFixedWindowSpline.h"

#include "MExtractedSignalCam.h"

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

ClassImp(MExtractFixedWindowSpline);

using namespace std;

const Byte_t  MExtractFixedWindowSpline::fgHiGainFirst  = 2;
const Byte_t  MExtractFixedWindowSpline::fgHiGainLast   = 14;
const Byte_t  MExtractFixedWindowSpline::fgLoGainFirst  = 3;
const Byte_t  MExtractFixedWindowSpline::fgLoGainLast   = 14;
// --------------------------------------------------------------------------
//
// Default constructor.
//
// Calls: 
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
// 
MExtractFixedWindowSpline::MExtractFixedWindowSpline(const char *name, const char *title) 
    : fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL),
      fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL)
{

  fName  = name  ? name  : "MExtractFixedWindowSpline";
  fTitle = title ? title : "Signal Extractor for a fixed FADC window using a fast spline";

  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
}

MExtractFixedWindowSpline::~MExtractFixedWindowSpline()
{
  
  if (fHiGainFirstDeriv)
    delete [] fHiGainFirstDeriv;
  if (fLoGainFirstDeriv)
    delete [] fLoGainFirstDeriv;
  if (fHiGainSecondDeriv)
    delete [] fHiGainSecondDeriv;
  if (fLoGainSecondDeriv)
    delete [] fLoGainSecondDeriv;
  
}

// --------------------------------------------------------------------------
//
// SetRange: 
//
// Checks: 
// - if the window defined by (fHiGainLast-fHiGainFirst-1) are odd, subtract one
// - if the window defined by (fLoGainLast-fLoGainFirst-1) are odd, subtract one
// - if the Hi Gain window is smaller than 2, set fHiGainLast to fHiGainFirst+1
// - if the Lo Gain window is smaller than 2, set fLoGainLast to fLoGainFirst+1
// 
// Calls:
// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
// 
// Sets:
// - fNumHiGainSamples to: (Float_t)(fHiGainLast-fHiGainFirst+1)
// - fNumLoGainSamples to: (Float_t)(fLoGainLast-fLoGainFirst+1)
// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)  
//  
void MExtractFixedWindowSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
{

  const Byte_t windowhi = hilast-hifirst+1;
  const Byte_t windowlo = lolast-lofirst+1;
  
  const Byte_t whieven = windowhi & ~1;
  const Byte_t wloeven = windowlo & ~1;

  if (whieven != windowhi)
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s%2i",": Hi Gain window size has to be even, set last slice from "
                    ,(int)hilast," to ",(int)(hilast-1)) << endl;
      hilast -= 1;
    }
  
  if (wloeven != windowlo)
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s%2i",": Lo Gain window size has to be even, set last slice from "
                    ,(int)lolast," to ",(int)(lolast-1)) << endl;
      lolast -= 1;
    }
  
  if (whieven<2) 
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s%2i",": Hi Gain window is smaller than 2 FADC sampes, set last slice from" 
                    ,(int)hilast," to ",(int)(hifirst+1)) << endl;
      hilast = hifirst+1;
    }
  
  if (wloeven<2) 
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s%2i",": Lo Gain window is smaller than 2 FADC sampes, set last slice from" 
                    ,(int)lolast," to ",(int)(lofirst+1)) << endl;
      lolast = lofirst+1;        
    }


  MExtractor::SetRange(hifirst,hilast,lofirst,lolast);

  //
  // Very important: Because the spline interpolates between the slice, 
  //                 the number of samples for the pedestal subtraction 
  //                 is now 1 less than with e.g. MExtractFixedWindow
  //
  fNumHiGainSamples = (Float_t)(fHiGainLast-fHiGainFirst);
  fNumLoGainSamples = (Float_t)(fLoGainLast-fLoGainFirst);  

  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);  
  
}

// --------------------------------------------------------------------------
//
// ReInit
//
// Calls:
// - MExtractor::ReInit(pList);
// - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
//                                fLoGainFirst, fLoGainLast, fNumLoGainSamples);
// 
// Deletes all arrays, if not NULL
// Creates new arrays according to the extraction range
//
Bool_t MExtractFixedWindowSpline::ReInit(MParList *pList)
{

  if (!MExtractor::ReInit(pList))
    return kFALSE;

  fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
                              fLoGainFirst, fLoGainLast, fNumLoGainSamples);

  if (fHiGainFirstDeriv)
    delete [] fHiGainFirstDeriv;
  if (fLoGainFirstDeriv)
    delete [] fLoGainFirstDeriv;
  if (fHiGainSecondDeriv)
    delete [] fHiGainSecondDeriv;
  if (fLoGainSecondDeriv)
    delete [] fLoGainSecondDeriv;
  
  Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;

  fHiGainFirstDeriv  = new Float_t[range];
  memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
  fHiGainSecondDeriv = new Float_t[range];
  memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));

  range = fLoGainLast - fLoGainFirst + 1;

  fLoGainFirstDeriv  = new Float_t[range];
  memset(fLoGainFirstDeriv,0,range*sizeof(Float_t));
  fLoGainSecondDeriv = new Float_t[range];
  memset(fLoGainSecondDeriv,0,range*sizeof(Float_t));

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// FindSignalHiGain:
//
// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
// - Sum up contents of *ptr
// - If *ptr is greater than fSaturationLimit, raise sat by 1
// 
// - If fHiLoLast is not 0, loop also from logain to (logain+fHiLoLast)
// - Sum up contents of logain
// - If *logain is greater than fSaturationLimit, raise sat by 1
//
void MExtractFixedWindowSpline::FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat) const
{
  
  const Byte_t *end = ptr + fHiGainLast - fHiGainFirst;
  Int_t range = fHiGainLast - fHiGainFirst + fHiLoLast + 1;
  
  Float_t pp;
  Int_t   i = 0;

  Int_t summ = 0;
  sum = (Float_t)*ptr++/2.;
  //
  // Check for saturation in all other slices
  //
  while (ptr<end)
    {

      summ += *ptr;
      i++;

      pp = fHiGainSecondDeriv[i-1] + 4.;
      fHiGainSecondDeriv[i] = -1.0/pp;
      fHiGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
      fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;

      if (*ptr++ >= fSaturationLimit)
        sat++;
    }
  
  if (*ptr++ >= fSaturationLimit)
    sat++;

  if (fHiLoLast == 0)
    {
      sum += (Float_t)*ptr/2.;
      fHiGainSecondDeriv[++i] = 0.;      

    }
  else
    {
      summ += *ptr;
      i++;

      pp = fHiGainSecondDeriv[i-1] + 4.;
      fHiGainSecondDeriv[i] = -1.0/pp;
      fHiGainFirstDeriv [i] = *(logain) - 2.* *(ptr) + *(ptr-1);
      fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;

      end    = logain + fHiLoLast - 1;

      while (logain<end)
        {

          summ += *logain;
          i++;
          
          pp = fHiGainSecondDeriv[i-1] + 4.;
          fHiGainSecondDeriv[i] = -1.0/pp;
          fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(ptr-1);
          fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;

          if (*logain++ >= fSaturationLimit)
            sat++;
          
        }
      sum += (Float_t)*logain/2;
      fHiGainSecondDeriv[++i] = 0.;
    }

  for (Int_t k=range-2;k>0;k--)
    {
      fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
      sum += 0.25*fHiGainSecondDeriv[k];
    }
  
  sum += (Float_t)summ;
}

// --------------------------------------------------------------------------
//
// FindSignalLoGain:
//
// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
// - Sum up contents of *ptr
// - If *ptr is greater than fSaturationLimit, raise sat by 1
// 
void MExtractFixedWindowSpline::FindSignalLoGain(Byte_t *ptr, Float_t &sum, Byte_t &sat) const
{
  
  const Byte_t *end = ptr + fLoGainLast - fLoGainFirst;
  Int_t range = fLoGainLast - fLoGainFirst + 1;
  
  Float_t pp;
  Int_t   i = 0;

  fLoGainSecondDeriv[0] = 0.;
  fLoGainFirstDeriv[0]  = 0.;

  Int_t summ = 0;
  sum = (Float_t)*ptr++/2.;
  //
  // Check for saturation in all other slices
  //
  while (ptr<end)
    {

      summ += *ptr;
      i++;

      pp = fLoGainSecondDeriv[i-1] + 4.;
      fLoGainSecondDeriv[i] = -1.0/pp;
      fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
      fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;

      if (*ptr++ >= fSaturationLimit)
        sat++;
    }
  
  if (*ptr++ >= fSaturationLimit)
    sat++;

  sum += (Float_t)*ptr/2.;
  fLoGainSecondDeriv[++i] = 0.;      
  
  for (Int_t k=range-2;k>0;k--)
    {
      fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
      sum += 0.25*fLoGainSecondDeriv[k];
    }
  
  sum += (Float_t)summ;
}