source: trunk/MagicSoft/Mars/mtools/MCubicSpline.cc@ 3029

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