Changeset 9212 for trunk/MagicSoft/Mars/mextralgo

Timestamp:

01/06/09 13:07:00 (16 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mextralgo

Files:

• MExtralgoSpline.cc (modified) (6 diffs)
• MExtralgoSpline.h (modified) (7 diffs)

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

@@ -18 +18 @@
 !   Author(s): Markus Gaug 09/2004 <mailto:markus@ifae.es>
 !
-!   Copyright: MAGIC Software Development, 2002-2007
+!   Copyright: MAGIC Software Development, 2002-2008
 !
 !
@@ -47 +47 @@
 //   (Note that the variables fDer are not real first derivative coefficients.)
 //
+//   Note, this spline is not optimized to be evaluated many many times, but
+//   it is optimized to be initialized very fast with new values again and
+//   again.
+//
 //////////////////////////////////////////////////////////////////////////////
 #include "MExtralgoSpline.h"
@@ -179 +183 @@
 // at x, until x==0. If y is not found -2 is returned.
 //
-Double_t MExtralgoSpline::SearchY(Float_t x, Float_t y) const
+Double_t MExtralgoSpline::SearchYdn(Float_t x, Float_t y) const
 {
     if (x>=fNum-1)
@@ -185 +189 @@
 
     Int_t i = TMath::FloorNint(x);
+    if (i<0)
+        return -2;
+
     Double_t rc = FindY(i, kTRUE, y, 0, x-i);
     while (--i>=0 && rc<0)
@@ -198 +205 @@
 
     Int_t i = TMath::FloorNint(x);
+    if (i>fNum-2)
+        return -2;
+
     Double_t rc = FindY(i, kFALSE, y, x-i, 1.);
-    while (i++<fNum-1 && rc<0)
+    while (++i<fNum-1 && rc<0)
         rc = FindY(i, kFALSE, y);
 
@@ -287 +297 @@
     // Search downwards for fHeight/2
     // By doing also a search upwards we could extract the pulse width
-    fTime      = SearchY(maxpos, h);
+    fTime      = SearchYdn(maxpos, h);
     fTimeDev   = 0;
     if (width)

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

@@ -58 +58 @@
 
     // Evaluate first derivative of spline in the interval i with x=[0;1[
-    inline Double_t EvalDeriv1(const Float_t x, const Int_t i) const
+    inline Double_t EvalDeriv1(const Int_t i, const Float_t x) const
     {
         // x := [0; 1[
@@ -64 +64 @@
         const Double_t difder = fDer2[i+1]-fDer2[i];
 
-        return 3*difder*x*x + 6*fDer2[i]*x - 2*fDer2[i] - fDer2[i+1] + difval;
+        //return 3*difder*x*x + 6*fDer2[i]*x - 2*fDer2[i] - fDer2[i+1] + difval;
+        return 3*difder*x*x + (6*x - 2)*fDer2[i] - fDer2[i+1] + difval;
     }
 
     // Evaluate second derivative of spline in the interval i with x=[0;1[
-    inline Double_t EvalDeriv2(const Float_t x, const Int_t i) const
+    inline Double_t EvalDeriv2(const Int_t i, const Float_t x) const
     {
         // x := [0; 1[
@@ -75 +76 @@
 
     Double_t FindY(Int_t i, Bool_t downwards, Double_t y=0, Double_t min=0, Double_t max=1) const;
-    Double_t SearchY(Float_t maxpos, Float_t y) const;
-    Double_t SearchYup(Float_t maxpos, Float_t y) const;
 
     Int_t EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const;
@@ -141 +140 @@
 
     // Calculate the intgeral of the Eval-function in
-    // bin i from a=[0;1[ to b=[0;1[
+    // bin i from 0 <= a < b < 1
     inline Double_t EvalInteg(Int_t i, Float_t a, Float_t b) const
     {
@@ -268 +267 @@
             return kFALSE;
 
-        const Bool_t ismax1 = x1>=min && x1<max && EvalDeriv2(x1, i)<0;
-        const Bool_t ismax2 = x2>=min && x2<max && EvalDeriv2(x2, i)<0;
+        const Bool_t ismax1 = x1>=min && x1<max && EvalDeriv2(i, x1)<0;
+        const Bool_t ismax2 = x2>=min && x2<max && EvalDeriv2(i, x2)<0;
 
         if (!ismax1 && !ismax2)
@@ -323 +322 @@
 
     Float_t EvalAt(const Float_t x) const;
+    Float_t Deriv1(const Float_t x) const;
+
+    Double_t SearchYdn(Float_t maxpos, Float_t y) const;
+    Double_t SearchYup(Float_t maxpos, Float_t y) const;
+
+    Double_t SearchYdn(Float_t y) const { return SearchYdn(fNum, y); }
+    Double_t SearchYup(Float_t y) const { return SearchYup(0,    y); }
 };
 
@@ -335 +341 @@
 }
 
+inline Float_t MExtralgoSpline::Deriv1(const Float_t x) const
+{
+    Int_t   i = TMath::FloorNint(x);
+    Float_t f = x-i;
+
+    Align(i, f);
+
+    return EvalDeriv1(i, f);
+}
+
 #endif