Changeset 9226 for trunk/MagicSoft/Mars/mextralgo

Timestamp:

01/17/09 14:52:43 (16 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mextralgo

Files:

• ExtralgoIncl.h (modified) (1 diff)
• MExtralgoSpline.cc (modified) (8 diffs)
• MExtralgoSpline.h (modified) (4 diffs)

trunk/MagicSoft/Mars/mextralgo/ExtralgoIncl.h

@@ -1 +1 @@
 #ifndef __CINT__
 
+#include "../mbase/MArrayF.h"
+
 #endif // __CINT__

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

@@ -18 +18 @@
 !   Author(s): Markus Gaug 09/2004 <mailto:markus@ifae.es>
 !
-!   Copyright: MAGIC Software Development, 2002-2008
+!   Copyright: MAGIC Software Development, 2002-2009
 !
 !
@@ -120 +120 @@
 // --------------------------------------------------------------------------
 //
-// Returns the highest x value in [min;max[ at which the spline in
-// the bin i is equal to y
-//
-// min and max are defined to be [0;1]
-//
-// The default for min is 0, the default for max is 1
-// The defaule for y is 0
-//
-Double_t MExtralgoSpline::FindY(Int_t i, Bool_t downwards, Double_t y, Double_t min, Double_t max) const
+// Solve the polynomial
+//
+//    y = a*x^3 + b*x^2 + c*x + d'
+//    0 = a*x^3 + b*x^2 + c*x + d' - y
+//
+// to find y in the i-th bin. Return the result as x1, x2, x3 and the return
+// code from MMath::SolvPol3.
+//
+Int_t MExtralgoSpline::SolvePol3(Int_t i, Double_t y, Double_t &x1, Double_t &x2, Double_t &x3) const
 {
     // y = a*x^3 + b*x^2 + c*x + d'
@@ -145 +145 @@
     //     return -2;
 
+    return MMath::SolvePol3(a, b, c, d, x1, x2, x3);
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns the highest x value in [min;max[ at which the spline in
+// the bin i is equal to y
+//
+// min and max must be in the range [0;1]
+//
+// The default for min is 0, the default for max is 1
+// The default for y is 0
+//
+Double_t MExtralgoSpline::FindYdn(Int_t i, Double_t y, Double_t min, Double_t max) const
+{
     Double_t x1, x2, x3;
-    const Int_t rc = MMath::SolvePol3(a, b, c, d, x1, x2, x3);
-
-    if (downwards==kTRUE)
-    {
-        Double_t x = -1;
-
-        if (rc>0 && x1>=min && x1<max && x1>x)
-            x = x1;
-        if (rc>1 && x2>=min && x2<max && x2>x)
-            x = x2;
-        if (rc>2 && x3>=min && x3<max && x3>x)
-            x = x3;
-
-        return x<0 ? -2 : x+i;
-    }
-    else
-    {
-        Double_t x = 2;
-
-        if (rc>0 && x1>min && x1<=max && x1<x)
-            x = x1;
-        if (rc>1 && x2>min && x2<=max && x2<x)
-            x = x2;
-        if (rc>2 && x3>min && x3<=max && x3<x)
-            x = x3;
-
-        return x>1 ? -2 : x+i;
-    }
-
-    return -2;
+    const Int_t rc = SolvePol3(i, y, x1, x2, x3);
+
+    Double_t x = -1;
+
+    if (rc>0 && x1>=min && x1<max && x1>x)
+        x = x1;
+    if (rc>1 && x2>=min && x2<max && x2>x)
+        x = x2;
+    if (rc>2 && x3>=min && x3<max && x3>x)
+        x = x3;
+
+    return x<0 ? -2 : x+i;
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns the lowest x value in [min;max[ at which the spline in
+// the bin i is equal to y
+//
+// min and max must be in the range [0;1]
+//
+// The default for min is 0, the default for max is 1
+// The default for y is 0
+//
+Double_t MExtralgoSpline::FindYup(Int_t i, Double_t y, Double_t min, Double_t max) const
+{
+    Double_t x1, x2, x3;
+    const Int_t rc = SolvePol3(i, y, x1, x2, x3);
+
+    Double_t x = 2;
+
+    if (rc>0 && x1>min && x1<=max && x1<x)
+        x = x1;
+    if (rc>1 && x2>min && x2<=max && x2<x)
+        x = x2;
+    if (rc>2 && x3>min && x3<=max && x3<x)
+        x = x3;
+
+    return x>1 ? -2 : x+i;
 }
 
@@ -181 +205 @@
 //
 // Search analytically downward for the value y of the spline, starting
-// at x, until x==0. If y is not found -2 is returned.
+// at x, until x==0. If y is not found or out of range -2 is returned.
 //
 Double_t MExtralgoSpline::SearchYdn(Float_t x, Float_t y) const
@@ -192 +216 @@
         return -2;
 
-    Double_t rc = FindY(i, kTRUE, y, 0, x-i);
+    Double_t rc = FindYdn(i, y, 0, x-i);
     while (--i>=0 && rc<0)
-        rc = FindY(i, kTRUE, y);
+        rc = FindYdn(i, y);
 
     return rc;
 }
 
+// --------------------------------------------------------------------------
+//
+// Search analytically upwards for the value y of the spline, starting
+// at x, until x==fNum-1. If y is not found or out of range -2 is returned.
+//
 Double_t MExtralgoSpline::SearchYup(Float_t x, Float_t y) const
 {
@@ -208 +237 @@
         return -2;
 
-    Double_t rc = FindY(i, kFALSE, y, x-i, 1.);
+    Double_t rc = FindYup(i, y, x-i, 1.);
     while (++i<fNum-1 && rc<0)
-        rc = FindY(i, kFALSE, y);
+        rc = FindYup(i, y);
 
     return rc;
@@ -222 +251 @@
 Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const
 {
-    // In the future we will calculate the intgeral analytically.
-    // It has been tested that it gives identical results within
-    // acceptable differences.
-
     // We allow extrapolation of 1/2 slice.
     const Float_t min = fRiseTime;        //-0.5+fRiseTime;
@@ -236 +261 @@
 
     return EvalInteg(pos-fRiseTime, pos+fFallTime);
+}
+
+MArrayF MExtralgoSpline::GetIntegral(bool norm) const
+{
+    MArrayF val(fNum);
+
+    //val[0] = 0;
+
+    Double_t integ = 0;
+    for (int i=0; i<fNum-1; i++)
+    {
+        integ += EvalInteg(i);
+
+        val[i+1] = integ;
+    }
+
+    if (norm)
+        for (int i=0; i<fNum-1; i++)
+            val[i+1] /= val[fNum-1];
+
+    return val;
 }
 

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

@@ -6 +6 @@
 #endif
 
+class MArrayF;
 class TComplex;
 
@@ -75 +76 @@
     }
 
-    Double_t FindY(Int_t i, Bool_t downwards, Double_t y=0, Double_t min=0, Double_t max=1) const;
+    Int_t SolvePol3(Int_t i, Double_t y, Double_t &x1, Double_t &x2, Double_t &x3) const;
+    Double_t FindYdn(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
+    Double_t FindYup(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
+    //Double_t FindY(Int_t i, Bool_t downwards, Double_t y=0, Double_t min=0, Double_t max=1) const;
 
     Int_t EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const;
@@ -153 +157 @@
     }
 
-    // Identical to sum EvalInteg(i, 0, 1) for i=0 to i<b but much faster
+    // Identical to sum of EvalInteg(i, 0, 1) for i=a to i=b-1,
+    // but much faster
+    // It is identical to EvalInteg(fVal[a], fVal[b])
     // Be carefull: NO RANGECHECK!
     inline Double_t EvalInteg(Int_t a, Int_t b) const
@@ -329 +335 @@
     Double_t SearchYdn(Float_t y) const { return SearchYdn(fNum, y); }
     Double_t SearchYup(Float_t y) const { return SearchYup(0,    y); }
+
+    MArrayF GetIntegral(bool norm=true) const;
 };