Changeset 3148 for trunk/MagicSoft/Mars/mtools

Timestamp:

02/14/04 11:48:43 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtools

Files:

• MCubicCoeff.cc (modified) (7 diffs)
• MCubicSpline.cc (modified) (11 diffs)
• MCubicSpline.h (modified) (1 diff)

trunk/MagicSoft/Mars/mtools/MCubicCoeff.cc

@@ -24 +24 @@
 
 //////////////////////////////////////////////////////////////////////////////
-//                                                                          // 
-//  Cubic Spline Interpolation                                              //
-//                                                                          //
+//
+//  Cubic Spline Interpolation
+//
 //////////////////////////////////////////////////////////////////////////////
-
 #include "MCubicCoeff.h"
 
-#include "TMath.h"
+#include <TMath.h>
 
 #include "MLog.h"
@@ -45 +44 @@
 // where x is the independent variable, 2 and 3 are exponents
 //
-
 MCubicCoeff::MCubicCoeff(Double_t x, Double_t xNext, Double_t y, Double_t yNext, 
                          Double_t a, Double_t b, Double_t c) : 
@@ -51 +49 @@
 {
     fH = fXNext - fX;
-    if(!EvalMinMax())
-    {
-        gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl;
-        fMin = 0;
-        fMax = 0;
-    }
+    if (EvalMinMax())
+        return;
+
+    gLog << warn << "Failed to eval interval Minimum and Maximum, returning zeros" << endl;
+    fMin = 0;
+    fMax = 0;
 }
 
@@ -66 +64 @@
 Double_t MCubicCoeff::Eval(Double_t x)
 {
-    Double_t dx = x - fX;
-    return (fY+dx*(fC+dx*(fB+dx*fA)));
+    const Double_t dx = x - fX;
+    return fY + dx*(fC + dx*(fB + dx*fA));
 }
 
@@ -79 +77 @@
 Bool_t MCubicCoeff::EvalMinMax()
 {
-    fMin = fMax = fY;
-    fAbMin = fAbMax = fX;
+    fMin = fY;
+    fMax = fY;
+
+    fAbMin = fX;
+    fAbMax = fX;
+
     if (fYNext < fMin)
     {
-        fMin = fYNext;
+        fMin   = fYNext;
         fAbMin = fXNext;
     }
     if (fYNext > fMax)
     {
-        fMax = fYNext;
+        fMax   = fYNext;
         fAbMax = fXNext;
     }
-    Double_t tempMinMax;
-    Double_t delta = 4.0*fB*fB - 12.0*fA*fC;
-        if (delta >= 0.0 && fA != 0.0)
-        {
-            Double_t sqrtDelta = sqrt(delta);
-            Double_t xPlus  = (-2.0*fB + sqrtDelta)/(6.0*fA);
-            Double_t xMinus = (-2.0*fB - sqrtDelta)/(6.0*fA);
-            if (xPlus >= 0.0 && xPlus <= fH)
-            {
-                tempMinMax = this->Eval(fX+xPlus);
-                if (tempMinMax < fMin)
-                {
-                    fMin = tempMinMax;
-                    fAbMin = fX + xPlus;
-                }
-                if (tempMinMax > fMax)
-                {
-                    fMax = tempMinMax;
-                    fAbMax = fX + xPlus;
-                }
-            }
-            if (xMinus >= 0.0 && xMinus <= fH)
-            {
-                tempMinMax = this->Eval(fX+xMinus);
-                if (tempMinMax < fMin)
-                {
-                    fMin = tempMinMax;
-                    fAbMin = fX + xMinus;
-                }
-                if (tempMinMax > fMax)
-                {
-                    fMax = tempMinMax;
-                    fAbMax = fX + xMinus;
-                }
-            }
-        }
-/* if fA is zero then we have only one possible solution */
-        else if (fA == 0.0 && fB != 0.0)
-        {
-            Double_t xSolo = - (fC/(2.0*fB));
-            if (xSolo >= 0.0 && xSolo <= fH)
-            {
-                tempMinMax = this->Eval(fX+xSolo);
-                if (tempMinMax < fMin)
-                {
-                    fMin = tempMinMax;
-                    fAbMin = fX + xSolo;
-                }
-                if (tempMinMax > fMax)
-                {
-                    fMax = tempMinMax;
-                    fAbMax = fX + xSolo;
-                }
-            }
-        }
+
+    const Double_t delta = fB*fB*4 - fA*fC*12;
+    if (delta >= 0 && fA != 0)
+    {
+        const Double_t sqrtDelta = TMath::Sqrt(delta);
+
+        const Double_t xPlus  = (-fB*2 + sqrtDelta)/(fA*6);
+        const Double_t xMinus = (-fB*2 - sqrtDelta)/(fA*6);
+
+        if (xPlus >= 0 && xPlus <= fH)
+        {
+            const Double_t tempMinMax = Eval(fX+xPlus);
+            if (tempMinMax < fMin)
+            {
+                fMin = tempMinMax;
+                fAbMin = fX + xPlus;
+            }
+            if (tempMinMax > fMax)
+            {
+                fMax = tempMinMax;
+                fAbMax = fX + xPlus;
+            }
+        }
+        if (xMinus >= 0 && xMinus <= fH)
+        {
+            const Double_t tempMinMax = Eval(fX+xMinus);
+            if (tempMinMax < fMin)
+            {
+                fMin = tempMinMax;
+                fAbMin = fX + xMinus;
+            }
+            if (tempMinMax > fMax)
+            {
+                fMax = tempMinMax;
+                fAbMax = fX + xMinus;
+            }
+        }
+        return kTRUE;
+    }
+
+    /* if fA is zero then we have only one possible solution */
+    if (fA == 0 && fB != 0)
+    {
+        const Double_t xSolo = -fC/(fB*2);
+
+        if (xSolo < 0 || xSolo > fH)
+            return kTRUE;
+
+        const Double_t tempMinMax = Eval(fX+xSolo);
+        if (tempMinMax < fMin)
+        {
+            fMin = tempMinMax;
+            fAbMin = fX + xSolo;
+        }
+        if (tempMinMax > fMax)
+        {
+            fMax = tempMinMax;
+            fAbMax = fX + xSolo;
+        }
+        return kTRUE;
+    }
+
     return kTRUE;
 }
@@ -157 +166 @@
 // There could be three or one real solutions
 //
-
 Short_t MCubicCoeff::FindCardanRoot(Double_t y, Double_t *x)
 {
-
-    Short_t whichRoot = -1;
-
-    Double_t a = fB/fA;
-    Double_t b = fC/fA;
-    Double_t c = (fY - y)/fA;
-
-    Double_t q = (a*a-3.0*b)/9.0;
-    Double_t r = (2.0*a*a*a-9.0*a*b+27.0*c)/54.0;
-
-    Double_t aOver3 = a/3.0;
-    Double_t r2 = r*r;
-    Double_t q3 = q*q*q;
+    const Double_t a = fB/fA;
+    const Double_t b = fC/fA;
+    const Double_t c = (fY - y)/fA;
+
+    const Double_t q = (a*a - b*3)/9;
+    const Double_t r = (a*a*a*2 - a*b*9 + c*27)/54;
+
+    const Double_t aOver3 = a/3;
+    const Double_t r2 = r*r;
+    const Double_t q3 = q*q*q;
     
     if (r2 < q3) //3 real sol
     {
-        Double_t sqrtQ = TMath::Sqrt(q);
-        Double_t min2SqQ = -2.0*sqrtQ;
-        Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ));
-
-        x[0] = min2SqQ * TMath::Cos(theta/3.0) - aOver3;
-        x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3.0) - aOver3;
-        x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3.0) - aOver3;
+        const Double_t sqrtQ = TMath::Sqrt(q);
+        const Double_t min2SqQ = -sqrtQ*2;
+        const Double_t theta = TMath::ACos(r/(sqrtQ*sqrtQ*sqrtQ));
+
+        x[0] = min2SqQ * TMath::Cos(theta/3) - aOver3;
+        x[1] = min2SqQ * TMath::Cos((theta+TMath::TwoPi())/3) - aOver3;
+        x[2] = min2SqQ * TMath::Cos((theta-TMath::TwoPi())/3) - aOver3;
+
         for (Int_t i = 0; i < 3; i++)
-            if (x[i] >= 0.0 && x[i] <= fH)
+            if (x[i] >= 0 && x[i] <= fH)
             {
-                x[i] = x[i] + fX;
-                whichRoot = i;
-                break;
+                x[i] += fX;
+                return i;
             }
-    }
-    else //only 1 real sol
-    {
-        Double_t s;
-        if (r == 0.0)
-            s = 0.0;
-        else if (r < 0.0)
-            s = TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0);
-        else // r > 0.0
-            s = - TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3),1.0/3.0);
-        if (s == 0.0)
-            x[0] = - aOver3;
-        else
-            x[0] = (s + q/s) - aOver3;
-        if (x[0] >= 0.0 && x[0] <= fH)
-        {
-            x[0] = x[0] + fX;
-            whichRoot = 0;
-        }
-    }
-    return whichRoot;
+        return -1;
+    }
+
+    const Double_t s = r==0 ? 0 : -TMath::Sign(TMath::Power(TMath::Abs(r) + TMath::Sqrt(r2 - q3), 1./3), r);
+
+    x[0] = s==0 ? - aOver3 : (s + q/s) - aOver3;
+
+    if (x[0] < 0 || x[0] > fH)
+        return -1;
+
+    x[0] += fX;
+    return 0;
 }
 
@@ -220 +216 @@
 Bool_t MCubicCoeff :: IsIn(Double_t x)
 {
-    if (x >= fX && x <= fXNext)
-        return kTRUE;
-    else
-        return kFALSE;
-}
+    return x >= fX && x <= fXNext;
+}

trunk/MagicSoft/Mars/mtools/MCubicSpline.cc

@@ -24 +24 @@
 
 //////////////////////////////////////////////////////////////////////////////
-//                                                                          // 
-//  Cubic Spline Interpolation                                              //
-//                                                                          //
+//
+//  Cubic Spline Interpolation
+//
 //////////////////////////////////////////////////////////////////////////////
-
 #include "MCubicSpline.h"
 
-#include "MCubicCoeff.h"
-
-#include "TMath.h"
+#include <TMath.h>
 
 #include "MLog.h"
 #include "MLogManip.h"
 
+#include "MCubicCoeff.h"
+
 ClassImp(MCubicSpline);
 
@@ -48 +47 @@
 //
 MCubicSpline::MCubicSpline(Byte_t *y, Byte_t *x, Bool_t areAllEq,
-                           Int_t n, Double_t begSD, Double_t endSD):
-    fN(n)
+                           Int_t n, Double_t begSD, Double_t endSD)
 {
     Init(y,x,areAllEq,n,begSD,endSD);
@@ -62 +60 @@
 {
     Byte_t x[]={0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0A,0x0B,0x0C,0x0D,0x0E};
-    fN = 15;
     Init(y,x,kTRUE,15,0.0,0.0);
 }
@@ -75 +72 @@
 
 {
-    Double_t *temp = new Double_t[fN-1];
-    Double_t *ysd = new Double_t[fN-1];
-    Double_t p,h;
-    MCubicCoeff *tempCoeff;
-    fCoeff = new TObjArray(fN-1,0);
+    Double_t *temp = new Double_t[n-1];
+    Double_t *ysd  = new Double_t[n-1];
+
+    fCoeff = new TObjArray(n-1,0);
+
+    ysd[0]  =begSD;
+    temp[0] =begSD;
+    ysd[n-1]=endSD;
     
+    Double_t h = x[1]-x[0];
+
     if (areAllEq)
     {
-        h = x[1]-x[0];
-        ysd[0]=temp[0]=begSD;
-        ysd[n-1]=endSD;
         for(Int_t i = 1; i < n-1; i++)
         {
-            p = 0.5*ysd[i-1]+2.0;
-            ysd[i] = (-0.5)/p;
-            temp[i] = (y[i+1]-2*y[i]+y[i-1])/h;
-            temp[i] = (6.0*temp[i]/h-0.5*temp[i-1])/p;
-        }
-        for(Int_t i = n-2; i > 0; i--)
-            ysd[i]=ysd[i]*ysd[i+1]+temp[i];
-        for(Int_t i = 0; i < n-1; i++)
-        {
-            tempCoeff = new MCubicCoeff(x[i],x[i+1],y[i],y[i+1],(ysd[i+1]-ysd[i])/(6*h),
-                                        ysd[i]/2.0,(y[i+1]-y[i])/h-(h*(ysd[i+1]+2*ysd[i]))/6);
-            fCoeff->AddAt(tempCoeff,i);
-        }
-        delete [] temp;
-        delete [] ysd;
+            const Double_t p = ysd[i-1]/2+2;
+
+            ysd[i]  = -0.5/p;
+            temp[i] = (y[i+1] - y[i]*2 + y[i-1])/h;
+            temp[i] = (temp[i]*6/h-temp[i-1]/2)/p;
+        }
     }
     else
     {
-        Double_t sig;
-        ysd[0]=temp[0]=begSD;
-        ysd[n-1]=endSD;
         for(Int_t i = 1; i < n-1; i++)
         {
-            sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
-            p = sig*ysd[i-1]+2.0;
-            ysd[i] = (sig-1.0)/p;
+            const Double_t sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
+
+            const Double_t p = sig*ysd[i-1]+2;
+
+            ysd[i]  = (sig-1.0)/p;
             temp[i] = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
-            temp[i] = (6.0*temp[i]/(x[i+1]-x[i-1])-sig*temp[i-1])/p;
+            temp[i] = (temp[i]*6/(x[i+1]-x[i-1])-sig*temp[i-1])/p;
         }
-        for(Int_t i = n-2; i > 0; i--)
-            ysd[i]=ysd[i]*ysd[i+1]+temp[i];
-        for(Int_t i = 0; i < n-1; i++)
-        {
-            h = x[i+1]-x[i];
-            tempCoeff = new MCubicCoeff(x[i],x[i+1],y[i],y[i+1],(ysd[i+1]-ysd[i])/(6*h),
-                                        ysd[i]/2.0,(y[i+1]-y[i])/h-(h*(ysd[i+1]+2*ysd[i]))/6);
-            fCoeff->AddAt(tempCoeff,i);
-        }
-        delete [] temp;
-        delete [] ysd;
-    }
+    }
+
+    for(Int_t i = n-2; i > 0; i--)
+        ysd[i] = ysd[i]*ysd[i+1] + temp[i];
+
+    for(Int_t i = 0; i < n-1; i++)
+    {
+        if (!areAllEq)
+            h = x[i+1]-x[i];
+
+        MCubicCoeff *c = new MCubicCoeff(x[i], x[i+1], y[i], y[i+1], (ysd[i+1]-ysd[i])/(h*6),
+                                         ysd[i]/2, (y[i+1]-y[i])/h-(h*(ysd[i+1]+ysd[i]*2))/6);
+        fCoeff->AddAt(c, i);
+    }
+
+    delete [] temp;
+    delete [] ysd;
 }
 
@@ -134 +128 @@
 {
     fCoeff->Delete();
+    delete fCoeff;
 }
 
@@ -142 +137 @@
 Double_t MCubicSpline :: Eval(Double_t x)
 {
-    for (Int_t i = 0; i < fN-1; i++)
-        if (((MCubicCoeff*)fCoeff->At(i))->IsIn(x))
-            return ((MCubicCoeff*)fCoeff->At(i))->Eval(x);
+    const Int_t n = fCoeff->GetSize()-1;
+    for (Int_t i = 0; i < n; i++)
+    {
+        MCubicCoeff *c = (MCubicCoeff*)fCoeff->UncheckedAt(i);
+        if (c->IsIn(x))
+            return c->Eval(x);
+    }
+
     gLog << warn << "Cannot evaluate Spline at " << x << "; returning 0";
-    return 0.0;
+
+    return 0;
 }
 
@@ -155 +156 @@
 Double_t MCubicSpline :: EvalMax()
 {
-    Double_t temp_max;
-    Double_t max = ((MCubicCoeff*)fCoeff->At(0))->GetMax();
-    for (Int_t i = 1; i < fN-1; i++)
-    {
-        temp_max = ((MCubicCoeff*)fCoeff->At(i))->GetMax();
-        if (temp_max > max)
-            max = temp_max;
-    }
+    Double_t max = -FLT_MAX;
+
+    TIter Next(fCoeff);
+    MCubicCoeff *c;
+    while ((c=(MCubicCoeff*)Next()))
+        max = TMath::Max(max, c->GetMax());
+
     return max;
 }
@@ -172 +172 @@
 Double_t MCubicSpline :: EvalMin()
 {
-    Double_t temp_min;
-    Double_t min = ((MCubicCoeff*)fCoeff->At(0))->GetMin();
-    for (Int_t i = 1; i < fN-1; i++)
-    {
-        temp_min = ((MCubicCoeff*)fCoeff->At(i))->GetMin();
-        if (temp_min < min)
-            min = temp_min;
-    }
+    Double_t min = FLT_MAX;
+
+    TIter Next(fCoeff);
+    MCubicCoeff *c;
+    while ((c=(MCubicCoeff*)Next()))
+        min = TMath::Min(min, c->GetMax());
+
     return min;
 }
@@ -189 +188 @@
 Double_t MCubicSpline :: EvalAbMax()
 {
-    Double_t temp_max;
-    Double_t abMax = ((MCubicCoeff*)fCoeff->At(0))->GetAbMax();
-    Double_t max = ((MCubicCoeff*)fCoeff->At(0))->GetMax();
-    for (Int_t i = 1; i < fN-1; i++)
-    {
-        temp_max = ((MCubicCoeff*)fCoeff->At(i))->GetMax();
-        if (temp_max > max)
-        {
-            max = temp_max;
-            abMax = ((MCubicCoeff*)fCoeff->At(i))->GetAbMax();
-        }
-    }
-    return abMax;
+    Double_t max = -FLT_MAX;
+
+    TIter Next(fCoeff);
+
+    MCubicCoeff *c;
+    MCubicCoeff *cmax=0;
+
+    while ((c=(MCubicCoeff*)Next()))
+    {
+        const Double_t temp = c->GetMax();
+        if (temp <= max)
+            continue;
+
+        max = temp;
+        cmax = c;
+    }
+
+    return cmax ? cmax->GetAbMax() : -FLT_MAX;
 }
 
@@ -210 +214 @@
 Double_t MCubicSpline :: EvalAbMin()
 {
-    Double_t temp_min;
-    Double_t abMin = ((MCubicCoeff*)fCoeff->At(0))->GetAbMin();
-    Double_t min = ((MCubicCoeff*)fCoeff->At(0))->GetMin();
-    for (Int_t i = 1; i < fN-1; i++)
-    {
-        temp_min = ((MCubicCoeff*)fCoeff->At(i))->GetMin();
-        if (temp_min < min)
-        {
-            min = temp_min;
-            abMin = ((MCubicCoeff*)fCoeff->At(i))->GetAbMin();
-        }
-    }
-    return abMin;
+    Double_t min = FLT_MAX;
+
+    TIter Next(fCoeff);
+
+    MCubicCoeff *c;
+    MCubicCoeff *cmin=0;
+
+    while ((c=(MCubicCoeff*)Next()))
+    {
+        const Double_t temp = c->GetMax();
+        if (temp >= min)
+            continue;
+
+        min = temp;
+        cmin = c;
+    }
+
+    return cmin ? cmin->GetAbMax() : FLT_MAX;
 }
 
@@ -233 +242 @@
 Double_t MCubicSpline :: FindVal(Double_t y, Double_t x0, Char_t direction = 'l')
 {
-    Short_t whichRoot;
-    Double_t tempRoot;
-    Double_t *roots = new Double_t[3];
-
-    for (Int_t i = 0; i < fN-1; i++)
-    {
-        if(((MCubicCoeff*)fCoeff->At(i))->IsIn(x0))
-        {
-            if(direction == 'l')
-            {
-                for (Int_t j = i; j >= 0; j--)
-                {
-                    whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
-                    if (whichRoot >= 0 )
-                    {
-                        tempRoot = roots[whichRoot];
-                        delete [] roots;
-                        return tempRoot;
-                    }
-                }
-            }
-            if(direction == 'r')
-            {
-                for (Int_t j = i; j < fN-1; j++)
-                {
-                    whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
-                    if (whichRoot >= 0)
-                    {
-                        tempRoot = roots[whichRoot];
-                        delete [] roots;
-                        return tempRoot;
-                    }
-                }
-            }
-        }
-    }
+    Double_t roots[3] = { 0, 0, 0 };
+
+    const Int_t n = fCoeff->GetSize()-1;
+
+    for (Int_t i = 0; i < n; i++)
+    {
+        if (!((MCubicCoeff*)fCoeff->At(i))->IsIn(x0))
+            continue;
+
+        switch (direction)
+        {
+        case 'l':
+            for (Int_t j = i; j >= 0; j--)
+            {
+                const Int_t whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
+                if (whichRoot >= 0 )
+                    return roots[whichRoot];
+            }
+            break;
+
+        case 'r':
+            for (Int_t j = i; j < n; j++)
+            {
+                const Int_t whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
+                if (whichRoot >= 0)
+                    return roots[whichRoot];
+            }
+            break;
+        }
+    }
+
     gLog << warn << "Nothing found calling MCubicSpline :: FindVal(), returning 0" << endl;
-    return 0.0;
-}
+
+    return 0;
+}

trunk/MagicSoft/Mars/mtools/MCubicSpline.h

@@ -16 +16 @@
  private:
     TObjArray *fCoeff; //array of the coefficients
-    Int_t fN; //number of points
+
     void Init(Byte_t *y, Byte_t *x, Bool_t areAllEq, Int_t n, Double_t begSD, Double_t endSD);