Index: trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc
===================================================================
--- trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc	(revision 3971)
+++ trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc	(revision 3971)
@@ -0,0 +1,654 @@
+/* ======================================================================== *\
+!
+! *
+! * 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.025; // 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  = 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  = x;
+  dtime = fResolution;
+}
+
+
Index: trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.h
===================================================================
--- trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.h	(revision 3971)
+++ trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.h	(revision 3971)
@@ -0,0 +1,42 @@
+#ifndef MARS_MExtractTimeFastSpline
+#define MARS_MExtractTimeFastSpline
+
+#ifndef MARS_MExtractTime
+#include "MExtractTime.h"
+#endif
+
+class MPedestalPix;
+class MExtractTimeFastSpline : public MExtractTime
+{
+
+  static const Byte_t  fgHiGainFirst;    // Default for fHiGainFirst  (now set to: 2)
+  static const Byte_t  fgHiGainLast;     // Default for fHiGainLast   (now set to: 14)
+  static const Byte_t  fgLoGainFirst;    // Default for fLOGainFirst  (now set to: 3)
+  static const Byte_t  fgLoGainLast;     // Default for fLoGainLast   (now set to: 14)
+  static const Float_t fgResolution;     // Default for fResolution   (now set to: 0.003)
+
+  Float_t *fHiGainFirstDeriv;
+  Float_t *fLoGainFirstDeriv;  
+  Float_t *fHiGainSecondDeriv;
+  Float_t *fLoGainSecondDeriv;  
+
+  Float_t fResolution;                        // The time resolution in FADC units
+  
+  void FindTimeHiGain(Byte_t *first, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const;
+  void FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const;
+  
+public:
+
+  MExtractTimeFastSpline(const char *name=NULL, const char *title=NULL);
+  ~MExtractTimeFastSpline();
+    
+  void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0);      
+  void SetResolution(Float_t f=fgResolution)     { fResolution = f;  }
+  
+  ClassDef(MExtractTimeFastSpline, 0)   // Task to Extract the Arrival Times using a Fast Spline
+};
+
+#endif
+
+
+
Index: trunk/MagicSoft/Mars/msignal/Makefile
===================================================================
--- trunk/MagicSoft/Mars/msignal/Makefile	(revision 3962)
+++ trunk/MagicSoft/Mars/msignal/Makefile	(revision 3971)
@@ -43,4 +43,5 @@
            MExtractTime.cc \
            MExtractTimeSpline.cc \
+           MExtractTimeFastSpline.cc \
            MExtractTimeHighestIntegral.cc \
            MArrivalTime.cc \
Index: trunk/MagicSoft/Mars/msignal/SignalLinkDef.h
===================================================================
--- trunk/MagicSoft/Mars/msignal/SignalLinkDef.h	(revision 3962)
+++ trunk/MagicSoft/Mars/msignal/SignalLinkDef.h	(revision 3971)
@@ -23,4 +23,5 @@
 #pragma link C++ class MExtractTime+;
 #pragma link C++ class MExtractTimeSpline+;
+#pragma link C++ class MExtractTimeFastSpline+;
 #pragma link C++ class MExtractTimeHighestIntegral+;
 #pragma link C++ class MArrivalTimeCam;
