Changeset 8524 for trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

Timestamp:

05/17/07 12:19:59 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h (modified) (7 diffs)

trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

@@ -119 +119 @@
     inline void EvalDerivEq0(const Int_t i, Float_t &rc1, Float_t &rc2) const
     {
+        /* --- ORIGINAL CODE ---
         Double_t sumder = fDer2[i]+fDer2[i+1];
         Double_t difder = fDer2[i]-fDer2[i+1];
@@ -125 +126 @@
         Double_t sqt2  = difder*(fVal[i+1]-fVal[i]);
         Double_t sqt3  = sqrt(3*sqt1 + 3*sqt2);
-        Double_t denom = 3*(fDer2[i+1]-fDer2[i]);
+        Double_t denom = -3*(fDer2[i+1]-fDer2[i]);
 
         rc1 = -(3*fDer2[i] + sqt3)/denom;
         rc2 = -(3*fDer2[i] - sqt3)/denom;
+         */
+
+        Double_t sumder = fDer2[i]+fDer2[i+1];
+        Double_t difder = fDer2[i]-fDer2[i+1];
+
+        Double_t sqt1  = sumder*sumder - fDer2[i]*fDer2[i+1];
+        Double_t sqt2  = difder*(fVal[i+1]-fVal[i]);
+        Double_t sqt3  = sqt1+sqt2<0 ? 0 : sqrt((sqt1 + sqt2)/3);
+
+        rc1 = -(fDer2[i] + sqt3)/difder;
+        rc2 = -(fDer2[i] - sqt3)/difder;
     }
 
@@ -134 +146 @@
     inline Double_t EvalPrimitive(Int_t i, Float_t x) const
     {
-        /* TO BE CHECKED!
+        Align(i, x);
+
         if (x==0)
-            return 0;
+            return -fDer2[i]/4;
 
         if (x==1)
-            return (fVal[i+1]+fVal[i])/2 - fDer2[i+1]/4;
-            */
-        Align(i, x);
+            return (fVal[i+1] + fVal[i])/2 - fDer2[i+1]/4 - fDer2[i]/2;
 
         const Double_t x2  = x*x;
@@ -149 +160 @@
 
         return x2*fVal[i+1]/2 + (x4/2-x2)*fDer2[i+1]/2 + (x-x2/2)*fVal[i] + (x2/2-x-x14/4)*fDer2[i];
+
     }
 
@@ -167 +179 @@
     // Calculate the intgeral of the Eval-function in
     // bin i from a=[0;1[ to b=[0;1[
-    inline Double_t EvalInteg(Int_t i, Float_t a=0, Float_t b=1) const
+    inline Double_t EvalInteg(Int_t i, Float_t a, Float_t b) const
     {
         return EvalPrimitive(i, b)-EvalPrimitive(i, a);
+    }
+
+    // Identical to EvalInteg(i, 0, 1) but much faster
+    // Be carefull: NO RANGECHECK!
+    inline Double_t EvalInteg(Int_t i) const
+    {
+        return (fVal[i+1] + fVal[i])/2 - (fDer2[i+1] + fDer2[i])/4;
+    }
+
+    // Identical to sum EvalInteg(i, 0, 1) for i=0 to i<b but much faster
+    // Be carefull: NO RANGECHECK!
+    inline Double_t EvalInteg(Int_t a, Int_t b) const
+    {
+       /*
+        Double_t sum = 0;
+        for (int i=a; i<b; i++)
+            sum += EvalInteg(i);
+
+        return sum;
+        */
+        Double_t sum=0;
+        for (const Float_t *ptr=fDer2+a+1; ptr<fDer2+b; ptr++)
+            sum -= *ptr;
+
+        sum -= (fDer2[a]+fDer2[b])/2;
+
+        sum /= 2;
+
+        for (const Float_t *ptr=fVal+a+1; ptr<fVal+b; ptr++)
+            sum += *ptr;
+
+        sum += (fVal[a]+fVal[b])/2;
+
+        return sum;
     }
 
@@ -175 +221 @@
     inline Double_t EvalInteg(Float_t x0, Float_t x1) const
     {
+        // RANGE CHECK MISSING!
+
         const Int_t min = TMath::CeilNint(x0);
         const Int_t max = TMath::FloorNint(x1);
@@ -183 +231 @@
 
         // Sum complete intervals
-        Double_t sum = 0;
-        for (int i=min; i<max; i++)
-            sum += EvalInteg(i);
+        Double_t sum = EvalInteg(min, max);
 
         // Sum the incomplete intervals at the beginning and end