Changeset 9312 for trunk/MagicSoft/Mars/mbase

Timestamp:

02/10/09 20:00:10 (16 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mbase

Files:

• MSpline3.cc (modified) (5 diffs)
• MSpline3.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mbase/MSpline3.cc

@@ -38 +38 @@
 #include <TF1.h>
 
+#include "MArrayD.h"
+
 ClassImp(MSpline3);
 
@@ -48 +50 @@
 MSpline3::MSpline3(const TF1 &f, const char *opt, Double_t valbeg, Double_t valend)
     : TSpline3("MSpline3", f.GetXmin(), f.GetXmax(), &f, f.GetNpx(), opt, valbeg, valend)
+{
+}
+
+MSpline3::MSpline3(const TF1 &f, Double_t freq, const char *opt,Double_t valbeg, Double_t valend)
+    : TSpline3("MSpline3", f.GetXmin()*freq, f.GetXmax()*freq, ConvertFunc(f, freq).GetArray(), f.GetNpx(), opt, valbeg, valend)
 {
 }
@@ -105 +112 @@
 // a discrete binning is done similar to the constructor of TSpline
 //
-TGraph *MSpline3::ConvertFunc(const TF1 &f, Float_t freq) const
+MArrayD &MSpline3::ConvertFunc(const TF1 &f, Float_t freq) const
 {
     const UInt_t npx = f.GetNpx();
@@ -111 +118 @@
     // WARNING: This is a stupid workaround because the TSpline3-
     // constructor takes a pointer as input! It is not thread-safe!
-    static TGraph g;
+    static MArrayD g;
     g.Set(npx);
 
@@ -119 +126 @@
     {
         const Double_t x = f.GetXmin() + i*step;
-        g.SetPoint(i, x*freq, f.Eval(x));
+        g[i] = f.Eval(x);
     }
 
-    return &g;
+    return g;
 }
+
+// FIXME: As soon as TSpline3 allows access to fPoly we can implement
+//        a much faster evaluation of the spline, especially in
+//        special conditions like in MAnalogSignal::AddPulse
+//        This will be the case for root > 5.22/00
+
+/*
+Double_t MSpline3::EvalFast(Double_t x) const
+{
+    // Eval this spline at x
+    const Int_t klow=FindFast(x);
+    return fPoly[klow].Eval(x);
+}
+
+Int_t MSpline3::FindFast(Double_t x) const
+{
+    //
+    // If out of boundaries, extrapolate
+    // It may be badly wrong
+
+    // if (x<=fXmin)
+    //     return 0;
+    //
+    // if (x>=fXmax)
+    //     return fNp-1;
+
+    //
+    // Equidistant knots, use histogramming
+    if (fKstep)
+        return TMath::Min(Int_t((x-fXmin)/fDelta),fNp-1);
+
+    //
+    // Non equidistant knots, binary search
+    Int_t klow = 0;
+    Int_t khig = fNp-1;
+
+    Int_t khalf;
+    while (khig-klow>1)
+        if(x>fPoly[khalf=(klow+khig)/2].X())
+            klow=khalf;
+        else
+            khig=khalf;
+
+    // This could be removed, sanity check
+    //if(!(fPoly[klow].X()<=x && x<=fPoly[klow+1].X()))
+    //    Error("Eval",
+    //          "Binary search failed x(%d) = %f < %f < x(%d) = %f\n",
+    //          klow,fPoly[klow].X(),x,fPoly[klow+1].X());
+
+    return klow;
+}
+*/

trunk/MagicSoft/Mars/mbase/MSpline3.h

@@ -6 +6 @@
 #endif
 
+class MArrayD;
+
 class MSpline3 : public TSpline3
 {
 private:
-    TGraph *ConvertSpline(const TSpline &s, Float_t freq) const;
-    TGraph *ConvertGraph(const TGraph &s, Float_t freq) const;
-    TGraph *ConvertFunc(const TF1 &f, Float_t freq) const;
+    TGraph  *ConvertSpline(const TSpline &s, Float_t freq) const;
+    TGraph  *ConvertGraph(const TGraph &s, Float_t freq) const;
+    MArrayD &ConvertFunc(const TF1 &f, Float_t freq) const;
 
 public:
-   MSpline3(const TGraph &g,
-            const char *opt=0, Double_t valbeg=0, Double_t valend=0)
-        : TSpline3("MSpline3", &g, opt, valbeg, valend)
-   {
-   }
+    MSpline3() { }
 
-   MSpline3(const TGraph &g, Double_t freq,
-            const char *opt=0, Double_t valbeg=0, Double_t valend=0)
-        : TSpline3("MSpline3", ConvertGraph(g, freq), opt, valbeg, valend)
-   {
-   }
+    // None equidistant binning (evaluation a bit slower in case of many bins)
+    MSpline3(const TGraph &g,
+             const char *opt=0, Double_t valbeg=0, Double_t valend=0)
+       : TSpline3("MSpline3", &g, opt, valbeg, valend)
+    {
+    }
 
-   MSpline3(const TSpline &s, Double_t freq,
-            const char *opt=0, Double_t valbeg=0, Double_t valend=0)
-        : TSpline3("MSpline3", ConvertSpline(s, freq), opt, valbeg, valend)
-   {
-   }
-        /*
-   MSpline3(Double_t xmin, Double_t xmax,
-            const TF1 *func, Int_t n, const char *opt=0,
-            Double_t valbeg=0, Double_t valend=0)
-        : TSpline3("MSpline3", xmin, xmax, func, n, opt, valbeg, valend)
-   {
-   }*/
-   MSpline3(const TF1 &f, const char *opt=0,Double_t valbeg=0, Double_t valend=0);
+    MSpline3(const TGraph &g, Double_t freq,
+             const char *opt=0, Double_t valbeg=0, Double_t valend=0)
+       : TSpline3("MSpline3", ConvertGraph(g, freq), opt, valbeg, valend)
+    {
+    }
 
-   MSpline3(const TF1 &f, Double_t freq,
-            const char *opt=0,Double_t valbeg=0, Double_t valend=0)
-        : TSpline3("MSpline3", ConvertFunc(f, freq), opt, valbeg, valend)
-   {
-   }
+    MSpline3(const TSpline &s, Double_t freq,
+             const char *opt=0, Double_t valbeg=0, Double_t valend=0)
+       : TSpline3("MSpline3", ConvertSpline(s, freq), opt, valbeg, valend)
+    {
+    }
 
-   MSpline3(const Double_t x[], const Double_t y[], Int_t n, const char *opt=0,
-            Double_t valbeg=0, Double_t valend=0)
-        : TSpline3("MSpline3", const_cast<Double_t*>(x), const_cast<Double_t*>(y), n, opt, valbeg, valend)
-   {
-   }
+    MSpline3(const Double_t x[], const Double_t y[], Int_t n, const char *opt=0,
+             Double_t valbeg=0, Double_t valend=0)
+       : TSpline3("MSpline3", const_cast<Double_t*>(x), const_cast<Double_t*>(y), n, opt, valbeg, valend)
+    {
+    }
+
+    // Equidistant binning (evaluation a bit faster in case of many bins)
+
+    // FIXME: In principle TF1 can be evaluated histogram like which should be faster
+    MSpline3(const TF1 &f, const char *opt=0,Double_t valbeg=0, Double_t valend=0);
+    MSpline3(const TF1 &f, Double_t freq, const char *opt=0,Double_t valbeg=0, Double_t valend=0);
 
    Double_t GetXmin() const { return fXmin; }     // Minimum value of abscissa
@@ -57 +54 @@
    Int_t GetNp() const { return fNp; }
 
+   TH1 *GetHistogram() const { return fHistogram; }
+
    ClassDef(MSpline3, 1) // An extension of the TSpline3
 };