source: branches/Mars_McMismatchStudy/mtools/MCubicSpline.cc@ 19401

Last change on this file since 19401 was 3390, checked in by raducci, 21 years ago
*** empty log message ***
File size: 6.6 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): Sebastian Raducci 01/2004 <mailto:raducci@fisica.uniud.it>
19!
20! Copyright: MAGIC Software Development, 2001-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// Cubic Spline Interpolation
28//
29//////////////////////////////////////////////////////////////////////////////
30#include "MCubicSpline.h"
31
32#include <TMath.h>
33
34#include "MLog.h"
35#include "MLogManip.h"
36
37#include "MCubicCoeff.h"
38
39ClassImp(MCubicSpline);
40
41using namespace std;
42
43//---------------------------------------------------------------------------
44//
45// Contructor
46//
47//
48MCubicSpline::MCubicSpline(const Byte_t *y, const Byte_t *x, Bool_t areAllEq,
49 Int_t n, Double_t begSD, Double_t endSD)
50{
51 Init(y,x,areAllEq,n,begSD,endSD);
52}
53
54//---------------------------------------------------------------------------
55//
56// Constructor for FADC slice (only the FADC counts are needed)
57//
58//
59MCubicSpline::MCubicSpline(const Byte_t *y)
60{
61 const Byte_t x[]={0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0A,0x0B,0x0C,0x0D,0x0E};
62 Init(y,x,kTRUE,15,0.0,0.0);
63}
64
65//---------------------------------------------------------------------------
66//
67// Constructors common part
68//
69//
70void MCubicSpline::Init(const Byte_t *y, const Byte_t *x, Bool_t areAllEq,
71 Int_t n, Double_t begSD, Double_t endSD)
72
73{
74 Double_t *temp = new Double_t[n];
75 Double_t *ysd = new Double_t[n];
76
77 fCoeff = new TObjArray(n-1,0);
78
79 ysd[0] =begSD;
80 temp[0] =begSD;
81 ysd[n-1]=endSD;
82
83 Double_t h = x[1]-x[0];
84
85 if (areAllEq)
86 {
87 for(Int_t i = 1; i < n-1; i++)
88 {
89 const Double_t p = ysd[i-1]/2+2;
90
91 ysd[i] = -0.5/p;
92 temp[i] = (y[i+1] - y[i]*2 + y[i-1])/h;
93 temp[i] = (temp[i]*6/h-temp[i-1]/2)/p;
94 }
95 }
96 else
97 {
98 for(Int_t i = 1; i < n-1; i++)
99 {
100 const Double_t sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
101
102 const Double_t p = sig*ysd[i-1]+2;
103
104 ysd[i] = (sig-1.0)/p;
105 temp[i] = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
106 temp[i] = (temp[i]*6/(x[i+1]-x[i-1])-sig*temp[i-1])/p;
107 }
108 }
109
110 for(Int_t i = n-2; i > 0; i--)
111 ysd[i] = ysd[i]*ysd[i+1] + temp[i];
112
113 for(Int_t i = 0; i < n-1; i++)
114 {
115 if (!areAllEq)
116 h = x[i+1]-x[i];
117
118 MCubicCoeff *c = new MCubicCoeff(x[i], x[i+1], y[i], y[i+1], (ysd[i+1]-ysd[i])/(h*6),
119 ysd[i]/2, (y[i+1]-y[i])/h-(h*(ysd[i+1]+ysd[i]*2))/6);
120 fCoeff->AddAt(c, i);
121 }
122
123 delete [] temp;
124 delete [] ysd;
125}
126
127MCubicSpline::~MCubicSpline()
128{
129 fCoeff->Delete();
130 delete fCoeff;
131}
132
133//---------------------------------------------------------------------------
134//
135// Evaluate the spline at a given point
136//
137Double_t MCubicSpline :: Eval(Double_t x)
138{
139 const Int_t n = fCoeff->GetSize();
140 for (Int_t i = 0; i < n; i++)
141 {
142 MCubicCoeff *c = (MCubicCoeff*)fCoeff->UncheckedAt(i);
143 if (c->IsIn(x))
144 return c->Eval(x);
145 }
146
147 gLog << warn << "Cannot evaluate Spline at " << x << "; returning 0";
148
149 return 0;
150}
151
152//----------------------------------------------------------------------------
153//
154// Search for max
155//
156Double_t MCubicSpline :: EvalMax()
157{
158 Double_t max = -FLT_MAX;
159
160 TIter Next(fCoeff);
161 MCubicCoeff *c;
162 while ((c=(MCubicCoeff*)Next()))
163 max = TMath::Max(max, c->GetMax());
164
165 return max;
166}
167
168//----------------------------------------------------------------------------
169//
170// Search for min
171//
172Double_t MCubicSpline :: EvalMin()
173{
174 Double_t min = FLT_MAX;
175
176 TIter Next(fCoeff);
177 MCubicCoeff *c;
178 while ((c=(MCubicCoeff*)Next()))
179 min = TMath::Min(min, c->GetMin());
180
181 return min;
182}
183
184//----------------------------------------------------------------------------
185//
186// Search for abscissa of the max
187//
188Double_t MCubicSpline :: EvalAbMax()
189{
190 Double_t max = -FLT_MAX;
191
192 TIter Next(fCoeff);
193
194 MCubicCoeff *c;
195 MCubicCoeff *cmax=0;
196
197 while ((c=(MCubicCoeff*)Next()))
198 {
199 const Double_t temp = c->GetMax();
200 if (temp <= max)
201 continue;
202
203 max = temp;
204 cmax = c;
205 }
206
207 return cmax ? cmax->GetAbMax() : -FLT_MAX;
208}
209
210//----------------------------------------------------------------------------
211//
212// Search for abscissa of the min
213//
214Double_t MCubicSpline :: EvalAbMin()
215{
216 Double_t min = FLT_MAX;
217
218 TIter Next(fCoeff);
219
220 MCubicCoeff *c;
221 MCubicCoeff *cmin=0;
222
223 while ((c=(MCubicCoeff*)Next()))
224 {
225 const Double_t temp = c->GetMin();
226 if (temp >= min)
227 continue;
228
229 min = temp;
230 cmin = c;
231 }
232
233 return cmin ? cmin->GetAbMin() : FLT_MAX;
234}
235
236//----------------------------------------------------------------------------
237//
238// Finds the abscissa where the spline reaches y starting from x0 going in
239// direction direction
240// You have to give as input a starting point and a direction ("l" or "r")
241//
242Double_t MCubicSpline :: FindVal(Double_t y, Double_t x0, Char_t direction = 'l')
243{
244 Double_t roots[3] = { 0, 0, 0 };
245
246 const Int_t n = fCoeff->GetSize()-1;
247
248 for (Int_t i = 0; i < n; i++)
249 {
250 if (!((MCubicCoeff*)fCoeff->At(i))->IsIn(x0))
251 continue;
252
253 switch (direction)
254 {
255 case 'l':
256 for (Int_t j = i; j >= 0; j--)
257 {
258 const Int_t whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
259 if (whichRoot >= 0 )
260 return roots[whichRoot];
261 }
262 break;
263
264 case 'r':
265 for (Int_t j = i; j < n; j++)
266 {
267 const Int_t whichRoot = ((MCubicCoeff*)fCoeff->At(j))->FindCardanRoot(y, roots);
268 if (whichRoot >= 0)
269 return roots[whichRoot];
270 }
271 break;
272 }
273 }
274
275 //gLog << warn << "Nothing found calling MCubicSpline :: FindVal(), returning 0" << endl;
276
277 return 0;
278}
Note: See TracBrowser for help on using the repository browser.