]>
Commit | Line | Data |
---|---|---|
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 | 71 | class TString; |
72 | class TSystem; | |
73 | class TRandom; | |
5406439e | 74 | class TH1; |
75 | class TMethodCall; | |
76 | class TRandom; | |
77 | class TROOT; | |
78 | class stdio; | |
79 | ||
40389866 | 80 | |
81 | ||
0eea9d4d | 82 | class 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_ | |
91 | AliCheb3D(const char* funName, Int_t DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec=1E-6); | |
92 | AliCheb3D(void (*ptr)(float*,float*), Int_t DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec=1E-6); | |
40389866 | 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); |
1d18ebe0 | 94 | AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Float_t prec=1E-6, Bool_t run=kTRUE); |
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_ | |
1d18ebe0 | 127 | int* GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec, Int_t npCheck=30); |
128 | void EstimateNPoints(float Prec, int gridBC[3][3],Int_t npd1=30,Int_t npd2=30,Int_t npd3=30); | |
0eea9d4d | 129 | void SaveData(const char* outfile,Bool_t append=kFALSE) const; |
130 | void SaveData(FILE* stream=stdout) const; | |
131 | // | |
132 | void SetUsrFunction(const char* name); | |
133 | void SetUsrFunction(void (*ptr)(float*,float*)); | |
da7cd221 | 134 | void EvalUsrFunction(const Float_t *x, Float_t *res); |
0eea9d4d | 135 | TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0); |
5406439e | 136 | static Int_t CalcChebCoefs(const Float_t *funval,int np, Float_t *outCoefs, Float_t prec=-1); |
0eea9d4d | 137 | #endif |
138 | // | |
139 | protected: | |
5406439e | 140 | void Clear(const Option_t* option = ""); |
141 | void SetDimOut(const int d); | |
142 | void PrepareBoundaries(const Float_t *bmin,const Float_t *bmax); | |
0eea9d4d | 143 | // |
144 | #ifdef _INC_CREATION_ALICHEB3D_ | |
145 | void EvalUsrFunction(); | |
146 | void DefineGrid(Int_t* npoints); | |
147 | Int_t ChebFit(); // fit all output dimensions | |
148 | Int_t ChebFit(int dmOut); | |
2572efdf | 149 | void SetPrecision(float prec) {fPrec = prec;} |
0eea9d4d | 150 | #endif |
151 | // | |
1d18ebe0 | 152 | Float_t MapToInternal(Float_t x,Int_t d) const; // map x to [-1:1] |
153 | Float_t MapToExternal(Float_t x,Int_t d) const {return x/fBScale[d]+fBOffset[d];} // map from [-1:1] to x | |
154 | Double_t MapToInternal(Double_t x,Int_t d) const; // map x to [-1:1] | |
155 | Double_t MapToExternal(Double_t x,Int_t d) const {return x/fBScale[d]+fBOffset[d];} // map from [-1:1] to x | |
40389866 | 156 | // |
0eea9d4d | 157 | protected: |
158 | Int_t fDimOut; // dimension of the ouput array | |
159 | Float_t fPrec; // requested precision | |
160 | Float_t fBMin[3]; // min boundaries in each dimension | |
161 | Float_t fBMax[3]; // max boundaries in each dimension | |
162 | Float_t fBScale[3]; // scale for boundary mapping to [-1:1] interval | |
163 | Float_t fBOffset[3]; // offset for boundary mapping to [-1:1] interval | |
164 | TObjArray fChebCalc; // Chebyshev parameterization for each output dimension | |
165 | // | |
166 | Int_t fMaxCoefs; //! max possible number of coefs per parameterization | |
167 | Int_t fNPoints[3]; //! number of used points in each dimension | |
168 | Float_t fArgsTmp[3]; //! temporary vector for coefs caluclation | |
0eea9d4d | 169 | Float_t * fResTmp; //! temporary vector for results of user function caluclation |
170 | Float_t * fGrid; //! temporary buffer for Chebyshef roots grid | |
171 | Int_t fGridOffs[3]; //! start of grid for each dimension | |
172 | TString fUsrFunName; //! name of user macro containing the function of "void (*fcn)(float*,float*)" format | |
173 | TMethodCall* fUsrMacro; //! Pointer to MethodCall for function from user macro | |
174 | // | |
1d18ebe0 | 175 | ClassDef(AliCheb3D,2) // Chebyshev parametrization for 3D->N function |
0eea9d4d | 176 | }; |
177 | ||
0eea9d4d | 178 | //__________________________________________________________________________________________ |
1d18ebe0 | 179 | inline Bool_t AliCheb3D::IsInside(const Float_t *par) const |
0eea9d4d | 180 | { |
181 | // check if the point is inside of the fitted box | |
1d18ebe0 | 182 | for (int i=3;i--;) if (fBMin[i]>par[i] || par[i]>fBMax[i]) return kFALSE; |
0eea9d4d | 183 | return kTRUE; |
184 | } | |
185 | ||
0eea9d4d | 186 | //__________________________________________________________________________________________ |
1d18ebe0 | 187 | inline Bool_t AliCheb3D::IsInside(const Double_t *par) const |
188 | { | |
189 | // check if the point is inside of the fitted box | |
190 | for (int i=3;i--;) if (fBMin[i]>par[i] || par[i]>fBMax[i]) return kFALSE; | |
191 | return kTRUE; | |
192 | } | |
193 | ||
194 | //__________________________________________________________________________________________ | |
195 | inline void AliCheb3D::Eval(const Float_t *par, Float_t *res) | |
0eea9d4d | 196 | { |
197 | // evaluate Chebyshev parameterization for 3d->DimOut function | |
198 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
199 | for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp); | |
200 | // | |
201 | } | |
1d18ebe0 | 202 | //__________________________________________________________________________________________ |
203 | inline void AliCheb3D::Eval(const Double_t *par, Double_t *res) | |
204 | { | |
205 | // evaluate Chebyshev parameterization for 3d->DimOut function | |
206 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
207 | for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp); | |
208 | // | |
209 | } | |
210 | ||
211 | //__________________________________________________________________________________________ | |
212 | inline Double_t AliCheb3D::Eval(const Double_t *par, int idim) | |
213 | { | |
214 | // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function | |
215 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
216 | return GetChebCalc(idim)->Eval(fArgsTmp); | |
217 | // | |
218 | } | |
0eea9d4d | 219 | |
220 | //__________________________________________________________________________________________ | |
1d18ebe0 | 221 | inline Float_t AliCheb3D::Eval(const Float_t *par, int idim) |
0eea9d4d | 222 | { |
223 | // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function | |
224 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
225 | return GetChebCalc(idim)->Eval(fArgsTmp); | |
226 | // | |
227 | } | |
228 | ||
40389866 | 229 | //__________________________________________________________________________________________ |
5406439e | 230 | inline void AliCheb3D::EvalDeriv3D(const Float_t *par, Float_t dbdr[3][3]) |
40389866 | 231 | { |
232 | // return gradient matrix | |
233 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
234 | for (int ib=3;ib--;) for (int id=3;id--;) dbdr[ib][id] = GetChebCalc(ib)->EvalDeriv(id,fArgsTmp)*fBScale[id]; | |
235 | } | |
236 | ||
1cf34ee8 | 237 | //__________________________________________________________________________________________ |
5406439e | 238 | inline void AliCheb3D::EvalDeriv3D2(const Float_t *par, Float_t dbdrdr[3][3][3]) |
1cf34ee8 | 239 | { |
240 | // return gradient matrix | |
241 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
242 | for (int ib=3;ib--;) for (int id=3;id--;)for (int id1=3;id1--;) | |
243 | dbdrdr[ib][id][id1] = GetChebCalc(ib)->EvalDeriv2(id,id1,fArgsTmp)*fBScale[id]*fBScale[id1]; | |
244 | } | |
245 | ||
40389866 | 246 | //__________________________________________________________________________________________ |
5406439e | 247 | inline void AliCheb3D::EvalDeriv(int dimd, const Float_t *par, Float_t *res) |
40389866 | 248 | { |
249 | // evaluate Chebyshev parameterization derivative for 3d->DimOut function | |
250 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
251 | for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];; | |
252 | // | |
253 | } | |
254 | ||
1cf34ee8 | 255 | //__________________________________________________________________________________________ |
5406439e | 256 | inline void AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t *par, Float_t *res) |
1cf34ee8 | 257 | { |
258 | // evaluate Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function | |
259 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
260 | for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2]; | |
261 | // | |
262 | } | |
263 | ||
40389866 | 264 | //__________________________________________________________________________________________ |
5406439e | 265 | inline Float_t AliCheb3D::EvalDeriv(int dimd, const Float_t *par, int idim) |
0eea9d4d | 266 | { |
40389866 | 267 | // evaluate Chebyshev parameterization derivative over dimd dimention for idim-th output dimension of 3d->DimOut function |
268 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
269 | return GetChebCalc(idim)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd]; | |
0eea9d4d | 270 | // |
271 | } | |
272 | ||
1cf34ee8 | 273 | //__________________________________________________________________________________________ |
5406439e | 274 | inline Float_t AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t *par, int idim) |
1cf34ee8 | 275 | { |
276 | // evaluate Chebyshev parameterization 2ns derivative over dimd1 and dimd2 dimensions for idim-th output dimension of 3d->DimOut function | |
277 | for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); | |
278 | return GetChebCalc(idim)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2]; | |
279 | // | |
280 | } | |
281 | ||
40389866 | 282 | //__________________________________________________________________________________________ |
1d18ebe0 | 283 | inline Float_t AliCheb3D::MapToInternal(Float_t x,Int_t d) const |
284 | { | |
285 | // map x to [-1:1] | |
286 | #ifdef _BRING_TO_BOUNDARY_ | |
287 | T res = (x-fBOffset[d])*fBScale[d]; | |
288 | if (res<-1) return -1; | |
289 | if (res> 1) return 1; | |
290 | return res; | |
291 | #else | |
292 | return (x-fBOffset[d])*fBScale[d]; | |
293 | #endif | |
294 | } | |
295 | ||
296 | //__________________________________________________________________________________________ | |
297 | inline Double_t AliCheb3D::MapToInternal(Double_t x,Int_t d) const | |
0eea9d4d | 298 | { |
40389866 | 299 | // map x to [-1:1] |
300 | #ifdef _BRING_TO_BOUNDARY_ | |
ff66b122 | 301 | T res = (x-fBOffset[d])*fBScale[d]; |
40389866 | 302 | if (res<-1) return -1; |
303 | if (res> 1) return 1; | |
304 | return res; | |
305 | #else | |
306 | return (x-fBOffset[d])*fBScale[d]; | |
307 | #endif | |
0eea9d4d | 308 | } |
309 | ||
310 | #endif |