source: trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc@ 4235

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtractTimeFastSpline
//
//   Fast arrival Time extractor using a cubic spline algorithm of Numerical Recipes. 
//   It returns the position of the half maximum between absolute maximum 
//   and pedestal of the spline that interpolates the FADC slices.
//
//   The precision of the half-maximum searches can be chosen by:
//   SetPrecision().
//
//   The precision of the maximum-finder is fixed to 0.025 FADC units.
// 
//////////////////////////////////////////////////////////////////////////////
#include "MExtractTimeFastSpline.h"

#include "MPedestalPix.h"

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


ClassImp(MExtractTimeFastSpline);

using namespace std;

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

  fName  = name  ? name  : "MExtractTimeFastSpline";
  fTitle = title ? title : "Calculate photons arrival time using a fast spline";

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

}

MExtractTimeFastSpline::~MExtractTimeFastSpline()
{
  
  if (fHiGainFirstDeriv)
    delete fHiGainFirstDeriv;
  if (fLoGainFirstDeriv)
    delete fLoGainFirstDeriv;
  if (fHiGainSecondDeriv)
    delete fHiGainSecondDeriv;
  if (fLoGainSecondDeriv)
    delete fLoGainSecondDeriv;
  
}


// --------------------------------------------------------------------------
//
// SetRange: 
//
// Calls:
// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
// - Deletes x, if not NULL
// - Creates x according to the range
//
void MExtractTimeFastSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
{

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

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

  if (range < 2)
    {
      *fLog << warn << GetDescriptor()
            << Form("%s%2i%s%2i%s",": Hi-Gain Extraction range [",(int)fHiGainFirst,","
                    ,fHiGainLast,"] too small, ") << endl;
      *fLog << warn << GetDescriptor()
            << " will move higher limit to obtain 4 slices " << endl;
      SetRange(fHiGainFirst, fHiGainLast+4-range,fLoGainFirst,fLoGainLast);
      range = fHiGainLast - fHiGainFirst + 1;
    }
  

  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;

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

    }
}


// --------------------------------------------------------------------------
//
// Calculates the arrival time for each pixel 
//
void MExtractTimeFastSpline::FindTimeHiGain(Byte_t *first, Float_t &time, Float_t &dtime, 
                                        Byte_t &sat, const MPedestalPix &ped) const
{
  
  const Int_t range = fHiGainLast - fHiGainFirst + 1;
  const Byte_t *end = first + range;
  Byte_t *p = first;
  Byte_t max    = 0;
  Byte_t maxpos = 0;

  //
  // Check for saturation in all other slices
  //
  while (++p<end)
    {
      if (*p > max)
        {
          max    = *p;
          maxpos =  p-first;
        }

      if (*p >= fSaturationLimit)
        {
          sat++;
          break;
        }
    }
  
  if (sat)
    return;

  if (maxpos < 2)
    return;
  
  Float_t pp;

  p = first;
  fHiGainSecondDeriv[0] = 0.;
  fHiGainFirstDeriv[0]  = 0.;

  for (Int_t i=1;i<range-1;i++)
    {
      pp = fHiGainSecondDeriv[i-1] + 4.;
      fHiGainSecondDeriv[i] = -1.0/pp;
      fHiGainFirstDeriv [i] = *(p+2) - 2.* *(p+1) + *(p);
      fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
      p++;
    }

  fHiGainSecondDeriv[range-1] = 0.;

  for (Int_t k=range-2;k>=0;k--)
    fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
  
  //
  // Now find the maximum  
  //
  Float_t step  = 0.2; // start with step size of 1ns and loop again with the smaller one
  Float_t lower = (Float_t)maxpos-1.;
  Float_t upper = (Float_t)maxpos;
  Float_t abmax = 0.;
  Float_t x     = lower;
  Float_t y     = 0.;
  Float_t a     = 1.;
  Float_t b     = 0.;
  Int_t   klo = maxpos-1;
  Int_t   khi = maxpos;
  Float_t klocont = (Float_t)*(first+klo);
  Float_t khicont = (Float_t)*(first+khi);
  time  = lower;

  //
  // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to 
  // interval maxpos+1.
  //
  while (x<upper+0.1)
    {

      x += step;
      a -= step;
      b += step;

      y = a*klocont
        + b*khicont
        + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
             +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;

      if (y > abmax)
        {
          abmax  = y;
          time   = x;
        }

    }

 if (time > upper-0.1)
    {

      upper = (Float_t)maxpos+1.;
      lower = (Float_t)maxpos;
      x     = lower;
      a     = 1.;
      b     = 0.;
      khi   = maxpos+1;
      klo   = maxpos;
      klocont = (Float_t)*(first+klo);
      khicont = (Float_t)*(first+khi);

      while (x<upper+0.1)
        {

          x += step;
          a -= step;
          b += step;
          
          y = a* klocont
            + b* khicont
            + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
                 +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
          
          if (y > abmax)
            {
              abmax  = y;
              time   = x;
            }
        }
  }
 
  const Float_t up = time+step+0.015;

  x     = time;
  a     = upper - x;
  b     = x - lower;

  step  = 0.04; // step size of 83 ps 

  while (x<up)
    {

      x += step;
      a -= step;
      b += step;
      
      y = a* klocont
        + b* khicont
        + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
             +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
      
      if (y > abmax)
        {
          abmax  = y;
          time   = x;
        }
      
    }

  Float_t lo = time-0.225;

  x     = time;
  a     = upper - x;
  b     = x - lower;
  khi--;
  klo--;
  klocont = (Float_t)*(first+klo);
  khicont = (Float_t)*(first+khi);

  while (x>lo)
    {

      x -= step;
      a += step;
      b -= step;
      
      y = a* klocont
        + b* khicont
        + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
             +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
      
      if (y > abmax)
        {
          abmax  = y;
          time   = x;
        }
      
    }

  const Float_t pedes   = ped.GetPedestal();
  const Float_t halfmax = pedes + (abmax - pedes)/2.;

  //
  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
  // First, find the right FADC slice:
  // 
  klo   = maxpos;
  while (klo > maxpos-4)
    {
      if (*(first+klo) < (Byte_t)halfmax)
        break;
      klo--;
    }

  //
  // Loop from the beginning of the slice upwards to reach the halfmax:
  // With means of bisection:
  // 
  x     = (Float_t)klo;
  a     = 1.;
  b     = 0.;
  klocont = (Float_t)*(first+klo);
  khicont = (Float_t)*(first+klo+1);
  time  = x;

  step = 0.5;
  Bool_t back = kFALSE;

  while (step > fResolution)
    {
      
      if (back)
        {
          x -= step;
          a += step;
          b -= step;
        }
      else
        {
          x += step;
          a -= step;
          b += step;
        }
      
      y = a*klocont
        + b*khicont
        + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
             +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;

      if (y >= halfmax)
        back = kTRUE;
      else
        back = kFALSE;

      step /= 2.;
    }

  time  = (Float_t)fHiGainFirst + x;
  dtime = fResolution;
}


// --------------------------------------------------------------------------
//
// Calculates the arrival time for each pixel 
//
void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, 
                                        Byte_t &sat, const MPedestalPix &ped) const
{
  
  const Int_t range = fLoGainLast - fLoGainFirst + 1;
  const Byte_t *end = first + range;
  Byte_t *p = first;
  Byte_t max    = 0;
  Byte_t maxpos = 0;

  //
  // Check for saturation in all other slices
  //
  while (++p<end)
    {
      if (*p > max)
        {
          max    = *p;
          maxpos =  p-first;
        }

      if (*p >= fSaturationLimit)
        {
          sat++;
          break;
        }
    }
  
  if (sat)
    return;

  if (maxpos < 2)
    return;
  
  Float_t pp;

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

  for (Int_t i=1;i<range-1;i++)
    {
      pp = fLoGainSecondDeriv[i-1] + 4.;
      fLoGainSecondDeriv[i] = -1.0/pp;
      fLoGainFirstDeriv [i] = *(p+2) - 2.* *(p+1) + *(p);
      fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
      p++;
    }

  fLoGainSecondDeriv[range-1] = 0.;

  for (Int_t k=range-2;k>=0;k--)
    fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
  
  //
  // Now find the maximum  
  //
  Float_t step  = 0.2; // start with step size of 1ns and loop again with the smaller one
  Float_t lower = (Float_t)maxpos-1.;
  Float_t upper = (Float_t)maxpos;
  Float_t abmax = 0.;
  Float_t x     = lower;
  Float_t y     = 0.;
  Float_t a     = 1.;
  Float_t b     = 0.;
  Int_t   klo = maxpos-1;
  Int_t   khi = maxpos;
  Float_t klocont = (Float_t)*(first+klo);
  Float_t khicont = (Float_t)*(first+khi);
  time  = lower;

  //
  // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to 
  // interval maxpos+1.
  //
  while (x<upper+0.1)
    {

      x += step;
      a -= step;
      b += step;

      y = a*klocont
        + b*khicont
        + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
             +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;

      if (y > abmax)
        {
          abmax  = y;
          time   = x;
        }

    }

 if (time > upper-0.1)
    {

      upper = (Float_t)maxpos+1.;
      lower = (Float_t)maxpos;
      x     = lower;
      a     = 1.;
      b     = 0.;
      khi   = maxpos+1;
      klo   = maxpos;
      klocont = (Float_t)*(first+klo);
      khicont = (Float_t)*(first+khi);

      while (x<upper+0.1)
        {

          x += step;
          a -= step;
          b += step;
          
          y = a* klocont
            + b* khicont
            + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
                 +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
          
          if (y > abmax)
            {
              abmax  = y;
              time   = x;
            }
        }
  }
 
  const Float_t up = time+step+0.015;

  x     = time;
  a     = upper - x;
  b     = x - lower;

  step  = 0.025; // step size of 165 ps 

  while (x<up)
    {

      x += step;
      a -= step;
      b += step;
      
      y = a* klocont
        + b* khicont
        + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
             +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
      
      if (y > abmax)
        {
          abmax  = y;
          time   = x;
        }
      
    }

  Float_t lo = time-0.225;

  x     = time;
  a     = upper - x;
  b     = x - lower;
  khi--;
  klo--;
  klocont = (Float_t)*(first+klo);
  khicont = (Float_t)*(first+khi);

  while (x>lo)
    {

      x -= step;
      a += step;
      b -= step;
      
      y = a* klocont
        + b* khicont
        + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
             +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
      
      if (y > abmax)
        {
          abmax  = y;
          time   = x;
        }
      
    }

  const Float_t pedes   = ped.GetPedestal();
  const Float_t halfmax = pedes + (abmax - pedes)/2.;

  //
  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
  // First, find the right FADC slice:
  // 
  klo   = maxpos;
  while (klo > maxpos-4)
    {
      if (*(first+klo) < (Byte_t)halfmax)
        break;
      klo--;
    }

  //
  // Loop from the beginning of the slice upwards to reach the halfmax:
  // With means of bisection:
  // 
  x     = (Float_t)klo;
  a     = 1.;
  b     = 0.;
  klocont = (Float_t)*(first+klo);
  khicont = (Float_t)*(first+klo+1);
  time  = x;

  step = 0.5;
  Bool_t back = kFALSE;

  while (step > fResolution)
    {
      
      if (back)
        {
          x -= step;
          a += step;
          b -= step;
        }
      else
        {
          x += step;
          a -= step;
          b += step;
        }
      
      y = a*klocont
        + b*khicont
        + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
             +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;

      if (y >= halfmax)
        back = kTRUE;
      else
        back = kFALSE;

      step /= 2.;
    }

  time  = (Float_t)fLoGainFirst + x;
  dtime = fResolution;
}