Changeset 8545

Timestamp:

06/11/07 16:49:41 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mextralgo/MExtralgoSpline.cc (modified) (4 diffs)
• mextralgo/MExtralgoSpline.h (modified) (6 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -19 +19 @@
                                                  -*-*- END OF LINE -*-*-
 
- 2007/06/06 Markus Meyer
-
-   * Cosy/tpoint/:
-     - added tpoint files from Jan. 2007 to June 2007 
+ 2007/06/11 Thomas Bretz
+
+   * sponde.cc:
+     - added check for validity of resource file
+
+   * mbase/MMath.cc:
+     - small speed improvement to calclation of three solutions
+       for the third order pol.
+     - for a second order pol. set x1 and x2 if it has only one 
+       solution
+
+   * mbase/MMath.h:
+     - speed improvement using ::cbrt instead of pow(x, 1/3)
+
+   * mcalib/MCalibrationChargeCalc.cc:
+     - improved output
+
+   * mextralgo/MExtralgoSpline.cc:
+     - speed improvement by using a look up table for often used
+       and identical coefficients
+     - use MMath::SolvePol2 to get the null-points of the first
+       derivative (EvalDerivEq0)
+     - removed a lot of old an obsolete comments
 
 

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

@@ -18 +18 @@
 !   Author(s): Markus Gaug 09/2004 <mailto:markus@ifae.es>
 !
-!   Copyright: MAGIC Software Development, 2002-2006
+!   Copyright: MAGIC Software Development, 2002-2007
 !
 !
@@ -53 +53 @@
 
 #include "../mbase/MMath.h"
+#include "../mbase/MArrayF.h"
 
 using namespace std;
@@ -72 +73 @@
         return;
 
+    // Look up table for coefficients
+    static MArrayF lut;
+
+    // If the lut is not et large enough resize and reclaculate
+    if (fNum>(Int_t)lut.GetSize())
+    {
+        lut.Set(fNum);
+
+        lut[0] = 0.;
+        for (Int_t i=1; i<fNum-1; i++)
+            lut[i] = -1.0/(lut[i-1] + 4);
+    }
+
+    // Calculate the coefficients used to get reproduce the first and
+    // second derivative.
     fDer1[0] = 0.;
-    fDer2[0] = 0.;
-
     for (Int_t i=1; i<fNum-1; i++)
     {
-        const Float_t pp = fDer2[i-1] + 4.;
-
-        fDer2[i] = -1.0/pp;
-
         const Float_t d1 = fVal[i+1] - 2*fVal[i] + fVal[i-1];
-        fDer1[i] = (6.0*d1-fDer1[i-1])/pp;
+        fDer1[i] = (fDer1[i-1]-d1)*lut[i];
     }
 
     fDer2[fNum-1] = 0.;
-
     for (Int_t k=fNum-2; k>=0; k--)
-        fDer2[k] = fDer2[k]*fDer2[k+1] + fDer1[k];
-
-    for (Int_t k=fNum-2; k>=0; k--)
-        fDer2[k] /= 6.;
+        fDer2[k] = lut[k]*fDer2[k+1] + fDer1[k];
+}
+
+// --------------------------------------------------------------------------
+//
+// Return the two results x1 and x2 of f'(x)=0 for the third order
+// polynomial (spline) in the interval i. Return the number of results.
+// (0 if the fist derivative does not have a null-point)
+//
+Int_t MExtralgoSpline::EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const
+{
+    const Double_t difder = fDer2[i+1]-fDer2[i];
+    const Double_t difval = fVal[i+1] -fVal[i];
+
+    return MMath::SolvePol2(3*difder, 6*fDer2[i], difval-2*fDer2[i]-fDer2[i+1], x1, x2);
 }
 
@@ -114 +134 @@
     const Double_t c = fVal[i+1]-fVal[i] -2*fDer2[i]-fDer2[i+1];
     const Double_t d = fVal[i] - y;
+
+    // If the first derivative is nowhere==0 and it is increasing
+    // in one point, and the value we search is outside of the
+    // y-interval... it cannot be there
+    // if (c>0 && (d>0 || fVal[i+1]<y) && b*b<3*c*a)
+    //     return -2;
 
     Double_t x1, x2, x3;

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

@@ -30 +30 @@
 
     Float_t fHeightTm;
-
-//    Float_t fResolution;
 
     // Result
@@ -117 +115 @@
     */
 
+    Int_t EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const;
+/*
     inline void EvalDerivEq0(const Int_t i, Float_t &rc1, Float_t &rc2) const
     {
-        /* --- ORIGINAL CODE ---
+        // --- ORIGINAL CODE ---
         Double_t sumder = fDer2[i]+fDer2[i+1];
         Double_t difder = fDer2[i]-fDer2[i+1];
@@ -128 +128 @@
         Double_t denom = -3*(fDer2[i+1]-fDer2[i]);
 
-        rc1 = -(3*fDer2[i] + sqt3)/denom;
-        rc2 = -(3*fDer2[i] - sqt3)/denom;
-         */
-
+        rc1 = (3*fDer2[i] + sqt3)/denom;
+        rc2 = (3*fDer2[i] - sqt3)/denom;
+
+        // --- NEW CODE ---
         Double_t sumder = fDer2[i]+fDer2[i+1];
         Double_t difder = fDer2[i]-fDer2[i+1];
@@ -139 +139 @@
         Double_t sqt3  = sqt1+sqt2<0 ? 0 : sqrt((sqt1 + sqt2)/3);
 
-        rc1 = -(fDer2[i] + sqt3)/difder;
-        rc2 = -(fDer2[i] - sqt3)/difder;
-    }
+        rc1 = (fDer2[i] + sqt3)/difder;
+        rc2 = (fDer2[i] - sqt3)/difder;
+    }*/
 
     // Calculate the "Stammfunktion" of the Eval-function
@@ -355 +355 @@
         // Find analytical maximum in the bin i in the interval [min,max[
 
-        Float_t x1, x2;
-        EvalDerivEq0(i, x1, x2);
+        Double_t x1, x2;
+        if (!EvalDerivEq0(i, x1, x2))
+            return kFALSE;
+
         // const Float_t x1 = EvalDerivEq0S1(i);
         // const Float_t x2 = EvalDerivEq0S2(i);
@@ -479 +481 @@
     void SetExtractionType(ExtractionType_t typ)     { fExtractionType = typ; }
     void SetHeightTm(Float_t h)                      { fHeightTm = h; }
-        //    void SetResolution(Float_t res)                  { fResolution=res; }
 
     Float_t GetTime() const      { return fTime; }