Changeset 9424 for trunk

Timestamp:

04/16/09 11:47:17 (16 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mbase

Files:

• MSpline3.cc (modified) (1 diff)
• MSpline3.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mbase/MSpline3.cc

@@ -132 +132 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Return the integral in the splines bin i up to x.
+//
+// The TSpline3 in the Interval [fX[i], fX[i+1]] is defined as:
+//
+//      dx = x-fX[i]
+//      y = fY + dx*fB + dx*dx*fC + dx*dx*dx*fD
+//
+// This yields the integral:
+//
+//   int(y) = dx*fY + 1/2*dx*dx*fB + 1/3*dx*dx*dx*fC + 1/4*dx*dx*dx*dx*fD
+//          = dx*(fY + dx*(1/2*fB + dx*(1/3*fC + dx*(1/4*fD))))
+//
+// Which gives for the integral range [fX[i], fX[i]+w]:
+//   int(fX[i]+w)-int(fX[i]) = w*(fY + w*(1/2*fB + w*(1/3*fC + w*(1/4*fD))))
+//
+// and for the integral range [fX[i]+w, fX[i+1]]:
+//   int(fX[i+1])-int(fX[i]+w) = `
+//     W*(fY + W*(1/2*fB + W*(1/3*fC + W*(1/4*fD)))) -
+//     w*(fY + w*(1/2*fB + w*(1/3*fC + w*(1/4*fD))))
+// with
+//     W := fX[i+1]-fX[i]
+//
+Double_t MSpline3::Integral(Int_t i, Double_t x) const
+{
+    Double_t x0, y, b, c, d;
+    const_cast<MSpline3*>(this)->GetCoeff(i, x0, y, b, c, d);
+
+    const Double_t w = x-x0;
+
+    return w*(y + w*(b/2 + w*(c/3 + w*d/4)));
+}
+
+// --------------------------------------------------------------------------
+//
+// Return the integral of the spline's bin i.
+//
+Double_t MSpline3::Integral(Int_t i) const
+{
+    Double_t x, y;
+
+    GetKnot(i+1, x, y);
+
+    return Integral(i, x);
+}
+
+// --------------------------------------------------------------------------
+//
+// Return the integral from a to b
+//
+Double_t MSpline3::Integral(Double_t a, Double_t b) const
+{
+    const Int_t n = FindX(a);
+    const Int_t m = FindX(b);
+
+    Double_t sum = -Integral(n, a);
+
+    for (int i=n; i<=m-1; i++)
+        sum += Integral(i);
+
+    sum += Integral(m, b);
+
+    return sum;
+}
+
+// --------------------------------------------------------------------------
+//
+// Return the integral between Xmin and Xmax
+//
+Double_t MSpline3::Integral() const
+{
+    Double_t sum = 0;
+
+    for (int i=0; i<GetNp()-1; i++)
+        sum += Integral(i);
+
+    return sum;
+}
+
+// --------------------------------------------------------------------------
+//
+// Return the integral between Xmin and Xmax of int( f(x)*sin(x) )
+//
+// The x-axis is assumed to be in degrees
+//
+Double_t MSpline3::IntegralSolidAngle() const
+{
+    const Int_t n = GetNp();
+
+    MArrayD x(n);
+    MArrayD y(n);
+
+    for (int i=0; i<n; i++)
+    {
+        GetKnot(i, x[i], y[i]);
+
+        x[i] *= TMath::DegToRad();
+        y[i] *= TMath::Sin(x[i]);
+    }
+
+    return TMath::TwoPi()*MSpline3(x.GetArray(), y.GetArray(), n).Integral();
+}
+
+
 // FIXME: As soon as TSpline3 allows access to fPoly we can implement
 //        a much faster evaluation of the spline, especially in

trunk/MagicSoft/Mars/mbase/MSpline3.h

@@ -14 +14 @@
     TGraph  *ConvertGraph(const TGraph &s, Float_t freq) const;
     MArrayD &ConvertFunc(const TF1 &f, Float_t freq) const;
+
+    Double_t Integral(Int_t i, Double_t x) const;
+    Double_t Integral(Int_t i) const;
 
 public:
@@ -49 +52 @@
     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
-   Double_t GetXmax() const { return fXmax; }     // Maximum value of abscissa
+    Double_t GetXmin() const { return fXmin; }     // Minimum value of abscissa
+    Double_t GetXmax() const { return fXmax; }     // Maximum value of abscissa
 
-   Int_t GetNp() const { return fNp; }
+    Int_t GetNp() const { return fNp; }
 
-   TH1 *GetHistogram() const { return (TH1*)fHistogram; }
+    TH1 *GetHistogram() const { return (TH1*)fHistogram; }
 
-   ClassDef(MSpline3, 1) // An extension of the TSpline3
+    Double_t Integral(Double_t a, Double_t b) const;
+    Double_t Integral() const;
+
+    Double_t IntegralSolidAngle() const;
+
+    ClassDef(MSpline3, 1) // An extension of the TSpline3
 };