]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/STEERBase/AliCheb3D.h
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliCheb3D.h
CommitLineData
0eea9d4d 1// Author: ruben.shahoyan@cern.ch 09/09/2006
40389866 2
0eea9d4d 3////////////////////////////////////////////////////////////////////////////////
4// //
5// AliCheb3D produces the interpolation of the user 3D->NDimOut arbitrary //
6// function supplied in "void (*fcn)(float* inp,float* out)" format //
7// either in a separate macro file or as a function pointer. //
8// Only coefficients needed to guarantee the requested precision are kept. //
9// //
10// The user-callable methods are: //
11// To create the interpolation use: //
12// AliCheb3D(const char* funName, // name of the file with user function //
13// or //
14// AliCheb3D(void (*ptr)(float*,float*),// pointer on the user function //
15// Int_t DimOut, // dimensionality of the function's output //
16// Float_t *bmin, // lower 3D bounds of interpolation domain //
17// Float_t *bmax, // upper 3D bounds of interpolation domain //
18// Int_t *npoints, // number of points in each of 3 input //
19// // dimension, defining the interpolation grid //
20// Float_t prec=1E-6); // requested max.absolute difference between //
21// // the interpolation and any point on grid //
22// //
23// To test obtained parameterization use the method //
24// TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0); //
25// it will compare the user output of the user function and interpolation //
26// for idim-th output dimension and fill the difference in the supplied //
27// histogram. If no histogram is supplied, it will be created. //
28// //
29// To save the interpolation data: //
30// SaveData(const char* filename, Bool_t append ) //
31// write text file with data. If append is kTRUE and the output file already //
32// exists, data will be added in the end of the file. //
33// Alternatively, SaveData(FILE* stream) will write the data to //
34// already existing stream. //
35// //
36// To read back already stored interpolation use either the constructor //
37// AliCheb3D(const char* inpFile); //
38// or the default constructor AliCheb3D() followed by //
39// AliCheb3D::LoadData(const char* inpFile); //
40// //
41// To compute the interpolation use Eval(float* par,float *res) method, with //
42// par being 3D vector of arguments (inside the validity region) and res is //
43// the array of DimOut elements for the output. //
44// //
45// If only one component (say, idim-th) of the output is needed, use faster //
46// Float_t Eval(Float_t *par,int idim) method. //
47// //
48// void Print(option="") will print the name, the ranges of validity and //
49// the absolute precision of the parameterization. Option "l" will also print //
50// the information about the number of coefficients for each output //
51// dimension. //
52// //
53// NOTE: during the evaluation no check is done for parameter vector being //
54// outside the interpolation region. If there is such a risk, use //
55// Bool_t IsInside(float *par) method. Chebyshev parameterization is not //
56// good for extrapolation! //
57// //
58// For the properties of Chebyshev parameterization see: //
59// H.Wind, CERN EP Internal Report, 81-12/Rev. //
60// //
61////////////////////////////////////////////////////////////////////////////////
62
63
5406439e 64#ifndef ALICHEB3D_H
65#define ALICHEB3D_H
66
0eea9d4d 67#include <TNamed.h>
0eea9d4d 68#include <TObjArray.h>
99adacae 69#include "AliCheb3DCalc.h"
70
0eea9d4d 71class TString;
72class TSystem;
73class TRandom;
5406439e 74class TH1;
75class TMethodCall;
76class TRandom;
77class TROOT;
78class stdio;
79
40389866 80
81
0eea9d4d 82class AliCheb3D: public TNamed
83{
84 public:
40389866 85 AliCheb3D();
86 AliCheb3D(const AliCheb3D& src);
87 AliCheb3D(const char* inpFile);
88 AliCheb3D(FILE* stream);
0eea9d4d 89 //
90#ifdef _INC_CREATION_ALICHEB3D_
f69ec3f4 91 AliCheb3D(const char* funName, Int_t DimOut, const Float_t *bmin, const Float_t *bmax, Int_t *npoints, Float_t prec=1E-6, const Float_t* precD=0);
92 AliCheb3D(void (*ptr)(float*,float*), Int_t DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec=1E-6, const Float_t* precD=0);
93 AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Int_t *npX,Int_t *npY,Int_t *npZ, Float_t prec=1E-6, const Float_t* precD=0);
94 AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Float_t prec=1E-6, Bool_t run=kTRUE, const Float_t* precD=0);
0eea9d4d 95#endif
96 //
97 ~AliCheb3D() {Clear();}
98 //
40389866 99 AliCheb3D& operator=(const AliCheb3D& rhs);
1d18ebe0 100 void Eval(const Float_t *par, Float_t *res);
101 Float_t Eval(const Float_t *par,int idim);
102 void Eval(const Double_t *par, Double_t *res);
103 Double_t Eval(const Double_t *par,int idim);
40389866 104 //
5406439e 105 void EvalDeriv(int dimd, const Float_t *par, Float_t *res);
106 void EvalDeriv2(int dimd1, int dimd2, const Float_t *par,Float_t *res);
107 Float_t EvalDeriv(int dimd, const Float_t *par, int idim);
108 Float_t EvalDeriv2(int dimd1,int dimd2, const Float_t *par, int idim);
109 void EvalDeriv3D(const Float_t *par, Float_t dbdr[3][3]);
110 void EvalDeriv3D2(const Float_t *par, Float_t dbdrdr[3][3][3]);
111 void Print(const Option_t* opt="") const;
1d18ebe0 112 Bool_t IsInside(const Float_t *par) const;
113 Bool_t IsInside(const Double_t *par) const;
ff66b122 114 //
0eea9d4d 115 AliCheb3DCalc* GetChebCalc(int i) const {return (AliCheb3DCalc*)fChebCalc.UncheckedAt(i);}
116 Float_t GetBoundMin(int i) const {return fBMin[i];}
117 Float_t GetBoundMax(int i) const {return fBMax[i];}
40389866 118 Float_t* GetBoundMin() const {return (float*)fBMin;}
119 Float_t* GetBoundMax() const {return (float*)fBMax;}
0eea9d4d 120 Float_t GetPrecision() const {return fPrec;}
121 void ShiftBound(int id,float dif);
122 //
123 void LoadData(const char* inpFile);
124 void LoadData(FILE* stream);
125 //
126#ifdef _INC_CREATION_ALICHEB3D_
7f42fc51 127 void InvertSign();
1d18ebe0 128 int* GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec, Int_t npCheck=30);
f69ec3f4 129 void EstimateNPoints(float prec, int gridBC[3][3],Int_t npd1=30,Int_t npd2=30,Int_t npd3=30);
0eea9d4d 130 void SaveData(const char* outfile,Bool_t append=kFALSE) const;
131 void SaveData(FILE* stream=stdout) const;
132 //
133 void SetUsrFunction(const char* name);
134 void SetUsrFunction(void (*ptr)(float*,float*));
da7cd221 135 void EvalUsrFunction(const Float_t *x, Float_t *res);
0eea9d4d 136 TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0);
5406439e 137 static Int_t CalcChebCoefs(const Float_t *funval,int np, Float_t *outCoefs, Float_t prec=-1);
0eea9d4d 138#endif
139 //
140 protected:
5406439e 141 void Clear(const Option_t* option = "");
f69ec3f4 142 void SetDimOut(const int d, const float* prec=0);
5406439e 143 void PrepareBoundaries(const Float_t *bmin,const Float_t *bmax);
0eea9d4d 144 //
145#ifdef _INC_CREATION_ALICHEB3D_
146 void EvalUsrFunction();
147 void DefineGrid(Int_t* npoints);
148 Int_t ChebFit(); // fit all output dimensions
149 Int_t ChebFit(int dmOut);
2572efdf 150 void SetPrecision(float prec) {fPrec = prec;}
0eea9d4d 151#endif
152 //
1d18ebe0 153 Float_t MapToInternal(Float_t x,Int_t d) const; // map x to [-1:1]
154 Float_t MapToExternal(Float_t x,Int_t d) const {return x/fBScale[d]+fBOffset[d];} // map from [-1:1] to x
155 Double_t MapToInternal(Double_t x,Int_t d) const; // map x to [-1:1]
156 Double_t MapToExternal(Double_t x,Int_t d) const {return x/fBScale[d]+fBOffset[d];} // map from [-1:1] to x
40389866 157 //
0eea9d4d 158 protected:
159 Int_t fDimOut; // dimension of the ouput array
160 Float_t fPrec; // requested precision
161 Float_t fBMin[3]; // min boundaries in each dimension
162 Float_t fBMax[3]; // max boundaries in each dimension
163 Float_t fBScale[3]; // scale for boundary mapping to [-1:1] interval
164 Float_t fBOffset[3]; // offset for boundary mapping to [-1:1] interval
165 TObjArray fChebCalc; // Chebyshev parameterization for each output dimension
166 //
167 Int_t fMaxCoefs; //! max possible number of coefs per parameterization
168 Int_t fNPoints[3]; //! number of used points in each dimension
169 Float_t fArgsTmp[3]; //! temporary vector for coefs caluclation
0eea9d4d 170 Float_t * fResTmp; //! temporary vector for results of user function caluclation
171 Float_t * fGrid; //! temporary buffer for Chebyshef roots grid
172 Int_t fGridOffs[3]; //! start of grid for each dimension
173 TString fUsrFunName; //! name of user macro containing the function of "void (*fcn)(float*,float*)" format
174 TMethodCall* fUsrMacro; //! Pointer to MethodCall for function from user macro
175 //
f69ec3f4 176 static const Float_t fgkMinPrec; // smallest precision
177 //
1d18ebe0 178 ClassDef(AliCheb3D,2) // Chebyshev parametrization for 3D->N function
0eea9d4d 179};
180
0eea9d4d 181//__________________________________________________________________________________________
1d18ebe0 182inline Bool_t AliCheb3D::IsInside(const Float_t *par) const
0eea9d4d 183{
184 // check if the point is inside of the fitted box
1d18ebe0 185 for (int i=3;i--;) if (fBMin[i]>par[i] || par[i]>fBMax[i]) return kFALSE;
0eea9d4d 186 return kTRUE;
187}
188
0eea9d4d 189//__________________________________________________________________________________________
1d18ebe0 190inline Bool_t AliCheb3D::IsInside(const Double_t *par) const
191{
192 // check if the point is inside of the fitted box
193 for (int i=3;i--;) if (fBMin[i]>par[i] || par[i]>fBMax[i]) return kFALSE;
194 return kTRUE;
195}
196
197//__________________________________________________________________________________________
198inline void AliCheb3D::Eval(const Float_t *par, Float_t *res)
0eea9d4d 199{
200 // evaluate Chebyshev parameterization for 3d->DimOut function
201 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
202 for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp);
203 //
204}
1d18ebe0 205//__________________________________________________________________________________________
206inline void AliCheb3D::Eval(const Double_t *par, Double_t *res)
207{
208 // evaluate Chebyshev parameterization for 3d->DimOut function
209 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
210 for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp);
211 //
212}
213
214//__________________________________________________________________________________________
215inline Double_t AliCheb3D::Eval(const Double_t *par, int idim)
216{
217 // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
218 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
219 return GetChebCalc(idim)->Eval(fArgsTmp);
220 //
221}
0eea9d4d 222
223//__________________________________________________________________________________________
1d18ebe0 224inline Float_t AliCheb3D::Eval(const Float_t *par, int idim)
0eea9d4d 225{
226 // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
227 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
228 return GetChebCalc(idim)->Eval(fArgsTmp);
229 //
230}
231
40389866 232//__________________________________________________________________________________________
5406439e 233inline void AliCheb3D::EvalDeriv3D(const Float_t *par, Float_t dbdr[3][3])
40389866 234{
235 // return gradient matrix
236 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
237 for (int ib=3;ib--;) for (int id=3;id--;) dbdr[ib][id] = GetChebCalc(ib)->EvalDeriv(id,fArgsTmp)*fBScale[id];
238}
239
1cf34ee8 240//__________________________________________________________________________________________
5406439e 241inline void AliCheb3D::EvalDeriv3D2(const Float_t *par, Float_t dbdrdr[3][3][3])
1cf34ee8 242{
243 // return gradient matrix
244 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
245 for (int ib=3;ib--;) for (int id=3;id--;)for (int id1=3;id1--;)
246 dbdrdr[ib][id][id1] = GetChebCalc(ib)->EvalDeriv2(id,id1,fArgsTmp)*fBScale[id]*fBScale[id1];
247}
248
40389866 249//__________________________________________________________________________________________
5406439e 250inline void AliCheb3D::EvalDeriv(int dimd, const Float_t *par, Float_t *res)
40389866 251{
252 // evaluate Chebyshev parameterization derivative for 3d->DimOut function
253 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
254 for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];;
255 //
256}
257
1cf34ee8 258//__________________________________________________________________________________________
5406439e 259inline void AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t *par, Float_t *res)
1cf34ee8 260{
261 // evaluate Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function
262 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
263 for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2];
264 //
265}
266
40389866 267//__________________________________________________________________________________________
5406439e 268inline Float_t AliCheb3D::EvalDeriv(int dimd, const Float_t *par, int idim)
0eea9d4d 269{
40389866 270 // evaluate Chebyshev parameterization derivative over dimd dimention for idim-th output dimension of 3d->DimOut function
271 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
272 return GetChebCalc(idim)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];
0eea9d4d 273 //
274}
275
1cf34ee8 276//__________________________________________________________________________________________
5406439e 277inline Float_t AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t *par, int idim)
1cf34ee8 278{
279 // evaluate Chebyshev parameterization 2ns derivative over dimd1 and dimd2 dimensions for idim-th output dimension of 3d->DimOut function
280 for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
281 return GetChebCalc(idim)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2];
282 //
283}
284
40389866 285//__________________________________________________________________________________________
1d18ebe0 286inline Float_t AliCheb3D::MapToInternal(Float_t x,Int_t d) const
287{
288 // map x to [-1:1]
289#ifdef _BRING_TO_BOUNDARY_
290 T res = (x-fBOffset[d])*fBScale[d];
291 if (res<-1) return -1;
292 if (res> 1) return 1;
293 return res;
294#else
295 return (x-fBOffset[d])*fBScale[d];
296#endif
297}
298
299//__________________________________________________________________________________________
300inline Double_t AliCheb3D::MapToInternal(Double_t x,Int_t d) const
0eea9d4d 301{
40389866 302 // map x to [-1:1]
303#ifdef _BRING_TO_BOUNDARY_
ff66b122 304 T res = (x-fBOffset[d])*fBScale[d];
40389866 305 if (res<-1) return -1;
306 if (res> 1) return 1;
307 return res;
308#else
309 return (x-fBOffset[d])*fBScale[d];
310#endif
0eea9d4d 311}
312
313#endif