source: trunk/Mars/mbase/MSpline3.cc @ 19273

Last change on this file since 19273 was 19128, checked in by tbretz, 16 months ago
Added a Scale function.
File size: 7.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18!   Author(s): Thomas Bretz  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20!   Copyright: Software Development, 2000-2009
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27//  MSpline3
28//
29// This is a extension of TSpline3. In addition to TSpline3 it allows access
30// to Xmin, Xman and Np. The construction is a bit simplified because no
31// title hase to be given (it can be given later by SetTitle anyway)
32// and is provides constructors which allow to scale the x-values by
33// pre-defined multiplier (e.g. frequency) to create the spline.
34//
35//////////////////////////////////////////////////////////////////////////////
36#include "MSpline3.h"
37
38#include <TF1.h>
39#include <TMath.h>
40
41#include "MArrayD.h"
42
43ClassImp(MSpline3);
44
45using namespace std;
46
47// --------------------------------------------------------------------------
48//
49//  Constructor.
50//
51MSpline3::MSpline3(const TF1 &f, const char *opt, Double_t valbeg, Double_t valend)
52    : TSpline3("MSpline3", f.GetXmin(), f.GetXmax(), &f, f.GetNpx(), opt, valbeg, valend)
53{
54}
55
56MSpline3::MSpline3(const TF1 &f, Double_t freq, const char *opt,Double_t valbeg, Double_t valend)
57    : TSpline3("MSpline3", f.GetXmin()*freq, f.GetXmax()*freq, ConvertFunc(f, freq).GetArray(), f.GetNpx(), opt, valbeg, valend)
58{
59}
60
61// --------------------------------------------------------------------------
62//
63//  This is a helper to convert the x-values by multiplying with freq
64// before initializing the spline
65//
66TGraph *MSpline3::ConvertSpline(const TSpline &s, Float_t freq) const
67{
68    const UInt_t npx = s.GetNpx();
69
70    // WARNING: This is a stupid workaround because the TSpline3-
71    // constructor takes a pointer as input! It is not thread-safe!
72    static TGraph g;
73    g.Set(npx);
74
75    for (UInt_t i=0; i<npx; i++)
76    {
77        Double_t x, y;
78        s.GetKnot(i, x, y);
79        g.SetPoint(i, x*freq, y);
80    }
81
82    return &g;
83}
84
85// --------------------------------------------------------------------------
86//
87//  This is a helper to convert the x-values by multiplying with freq
88// before initializing the spline
89//
90TGraph *MSpline3::ConvertGraph(const TGraph &s, Float_t freq) const
91{
92    const UInt_t npx = s.GetN();
93
94    // WARNING: This is a stupid workaround because the TSpline3-
95    // constructor takes a pointer as input! It is not thread-safe!
96    static TGraph g;
97    g.Set(npx);
98
99    for (UInt_t i=0; i<npx; i++)
100    {
101        Double_t x, y;
102        s.GetPoint(i, x, y);
103        g.SetPoint(i, x*freq, y);
104    }
105
106    return &g;
107}
108
109// --------------------------------------------------------------------------
110//
111//  This is a helper to convert the x-values by multiplying with freq
112// before initializing the spline. The conversion from the function to
113// a discrete binning is done similar to the constructor of TSpline
114//
115MArrayD &MSpline3::ConvertFunc(const TF1 &f, Float_t freq) const
116{
117    const UInt_t npx = f.GetNpx();
118
119    // WARNING: This is a stupid workaround because the TSpline3-
120    // constructor takes a pointer as input! It is not thread-safe!
121    static MArrayD g;
122    g.Set(npx);
123
124    const Double_t step = (f.GetXmax()-f.GetXmin())/(npx-1);
125
126    for (UInt_t i=0; i<npx; ++i)
127    {
128        const Double_t x = f.GetXmin() + i*step;
129        g[i] = f.Eval(x);
130    }
131
132    return g;
133}
134
135// --------------------------------------------------------------------------
136//
137// Return the integral in the splines bin i up to x.
138//
139// The TSpline3 in the Interval [fX[i], fX[i+1]] is defined as:
140//
141//      dx = x-fX[i]
142//      y = fY + dx*fB + dx*dx*fC + dx*dx*dx*fD
143//
144// This yields the integral:
145//
146//   int(y) = dx*fY + 1/2*dx*dx*fB + 1/3*dx*dx*dx*fC + 1/4*dx*dx*dx*dx*fD
147//          = dx*(fY + dx*(1/2*fB + dx*(1/3*fC + dx*(1/4*fD))))
148//
149// Which gives for the integral range [fX[i], fX[i]+w]:
150//   int(fX[i]+w)-int(fX[i]) = w*(fY + w*(1/2*fB + w*(1/3*fC + w*(1/4*fD))))
151//
152// and for the integral range [fX[i]+w, fX[i+1]]:
153//   int(fX[i+1])-int(fX[i]+w) = `
154//     W*(fY + W*(1/2*fB + W*(1/3*fC + W*(1/4*fD)))) -
155//     w*(fY + w*(1/2*fB + w*(1/3*fC + w*(1/4*fD))))
156// with
157//     W := fX[i+1]-fX[i]
158//
159Double_t MSpline3::IntegralBin(Int_t i, Double_t x) const
160{
161    Double_t x0, y, b, c, d;
162    const_cast<MSpline3*>(this)->GetCoeff(i, x0, y, b, c, d);
163
164    const Double_t w = x-x0;
165
166    return w*(y + w*(b/2 + w*(c/3 + w*d/4)));
167}
168
169// --------------------------------------------------------------------------
170//
171// Return the integral of the spline's bin i.
172//
173Double_t MSpline3::IntegralBin(Int_t i) const
174{
175    Double_t x, y;
176
177    GetKnot(i+1, x, y);
178
179    return IntegralBin(i, x);
180}
181
182// --------------------------------------------------------------------------
183//
184// Return the integral from a to b
185//
186Double_t MSpline3::Integral(Double_t a, Double_t b) const
187{
188    const Int_t n = FindX(a);
189    const Int_t m = FindX(b);
190
191    Double_t sum = -IntegralBin(n, a);
192
193    for (int i=n; i<=m-1; i++)
194        sum += IntegralBin(i);
195
196    sum += IntegralBin(m, b);
197
198    return sum;
199}
200
201// --------------------------------------------------------------------------
202//
203// Return the integral between Xmin and Xmax
204//
205Double_t MSpline3::Integral() const
206{
207    Double_t sum = 0;
208
209    for (int i=0; i<GetNp()-1; i++)
210        sum += IntegralBin(i);
211
212    return sum;
213}
214
215// --------------------------------------------------------------------------
216//
217// Return the integral between Xmin and Xmax of int( f(x)*sin(x) )
218//
219// The x-axis is assumed to be in degrees
220//
221Double_t MSpline3::IntegralSolidAngle() const
222{
223    const Int_t n = GetNp();
224
225    MArrayD x(n);
226    MArrayD y(n);
227
228    for (int i=0; i<n; i++)
229    {
230        GetKnot(i, x[i], y[i]);
231
232        x[i] *= TMath::DegToRad();
233        y[i] *= TMath::Sin(x[i]);
234    }
235
236    return TMath::TwoPi()*MSpline3(x.GetArray(), y.GetArray(), n).Integral();
237}
238
239
240// FIXME: As soon as TSpline3 allows access to fPoly we can implement
241//        a much faster evaluation of the spline, especially in
242//        special conditions like in MAnalogSignal::AddPulse
243//        This will be the case for root > 5.22/00
244
245/*
246Double_t MSpline3::EvalFast(Double_t x) const
247{
248    // Eval this spline at x
249    const Int_t klow=FindFast(x);
250    return fPoly[klow].Eval(x);
251}
252
253Int_t MSpline3::FindFast(Double_t x) const
254{
255    //
256    // If out of boundaries, extrapolate
257    // It may be badly wrong
258
259    // if (x<=fXmin)
260    //     return 0;
261    //
262    // if (x>=fXmax)
263    //     return fNp-1;
264
265    //
266    // Equidistant knots, use histogramming
267    if (fKstep)
268        return TMath::Min(Int_t((x-fXmin)/fDelta),fNp-1);
269
270    //
271    // Non equidistant knots, binary search
272    Int_t klow = 0;
273    Int_t khig = fNp-1;
274
275    Int_t khalf;
276    while (khig-klow>1)
277        if(x>fPoly[khalf=(klow+khig)/2].X())
278            klow=khalf;
279        else
280            khig=khalf;
281
282    // This could be removed, sanity check
283    //if(!(fPoly[klow].X()<=x && x<=fPoly[klow+1].X()))
284    //    Error("Eval",
285    //          "Binary search failed x(%d) = %f < %f < x(%d) = %f\n",
286    //          klow,fPoly[klow].X(),x,fPoly[klow+1].X());
287
288    return klow;
289}
290*/
291
292void MSpline3::Scale(double scale)
293{
294    //return fY + dx*fB + dx*dx*fC + dx*dx*dx*fD;
295
296    for (int i=0; i<fNp; i++)
297    {
298        fPoly[i].B() *= scale;
299        fPoly[i].C() *= scale;
300        fPoly[i].D() *= scale;
301        fPoly[i].Y() *= scale;
302    }
303}
Note: See TracBrowser for help on using the repository browser.