]>
Commit | Line | Data |
---|---|---|
5406439e | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
0eea9d4d | 15 | |
16 | #include <TString.h> | |
17 | #include <TSystem.h> | |
0eea9d4d | 18 | #include <TROOT.h> |
40389866 | 19 | #include <TRandom.h> |
5406439e | 20 | #include <stdio.h> |
21 | #include <TMethodCall.h> | |
22 | #include <TMath.h> | |
23 | #include <TH1.h> | |
0eea9d4d | 24 | #include "AliCheb3D.h" |
5406439e | 25 | #include "AliCheb3DCalc.h" |
b8bec44f | 26 | #include "AliLog.h" |
0eea9d4d | 27 | |
0eea9d4d | 28 | ClassImp(AliCheb3D) |
29 | ||
40389866 | 30 | //__________________________________________________________________________________________ |
31 | AliCheb3D::AliCheb3D() : | |
32 | fDimOut(0), | |
33 | fPrec(0), | |
34 | fChebCalc(1), | |
35 | fMaxCoefs(0), | |
36 | fResTmp(0), | |
37 | fGrid(0), | |
38 | fUsrFunName(""), | |
39 | fUsrMacro(0) | |
0eea9d4d | 40 | { |
5e6a3170 | 41 | // Default constructor |
c8f37c50 | 42 | for (int i=3;i--;) { |
d97d27ff | 43 | fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = fArgsTmp[i] = 0; |
c8f37c50 | 44 | fNPoints[i] = 0; |
34473f9e | 45 | fGridOffs[i] = 0; |
c8f37c50 | 46 | } |
0eea9d4d | 47 | } |
48 | ||
40389866 | 49 | //__________________________________________________________________________________________ |
50 | AliCheb3D::AliCheb3D(const AliCheb3D& src) : | |
51 | TNamed(src), | |
52 | fDimOut(src.fDimOut), | |
53 | fPrec(src.fPrec), | |
54 | fChebCalc(1), | |
55 | fMaxCoefs(src.fMaxCoefs), | |
56 | fResTmp(0), | |
57 | fGrid(0), | |
58 | fUsrFunName(src.fUsrFunName), | |
59 | fUsrMacro(0) | |
0eea9d4d | 60 | { |
40389866 | 61 | // read coefs from text file |
62 | for (int i=3;i--;) { | |
63 | fBMin[i] = src.fBMin[i]; | |
64 | fBMax[i] = src.fBMax[i]; | |
65 | fBScale[i] = src.fBScale[i]; | |
66 | fBOffset[i] = src.fBOffset[i]; | |
67 | fNPoints[i] = src.fNPoints[i]; | |
c8f37c50 | 68 | fGridOffs[i] = src.fGridOffs[i]; |
f9efbbf6 | 69 | fArgsTmp[i] = 0; |
40389866 | 70 | } |
71 | for (int i=0;i<fDimOut;i++) { | |
72 | AliCheb3DCalc* cbc = src.GetChebCalc(i); | |
73 | if (cbc) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(*cbc),i); | |
74 | } | |
0eea9d4d | 75 | } |
76 | ||
40389866 | 77 | //__________________________________________________________________________________________ |
78 | AliCheb3D::AliCheb3D(const char* inpFile) : | |
79 | fDimOut(0), | |
80 | fPrec(0), | |
81 | fChebCalc(1), | |
82 | fMaxCoefs(0), | |
83 | fResTmp(0), | |
84 | fGrid(0), | |
85 | fUsrFunName(""), | |
86 | fUsrMacro(0) | |
0bc7b414 | 87 | { |
40389866 | 88 | // read coefs from text file |
c8f37c50 | 89 | for (int i=3;i--;) { |
90 | fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
91 | fNPoints[i] = 0; | |
34473f9e | 92 | fGridOffs[i] = 0; |
f9efbbf6 | 93 | fArgsTmp[i] = 0; |
c8f37c50 | 94 | } |
40389866 | 95 | LoadData(inpFile); |
c437b1a5 | 96 | } |
97 | ||
40389866 | 98 | //__________________________________________________________________________________________ |
99 | AliCheb3D::AliCheb3D(FILE* stream) : | |
100 | fDimOut(0), | |
101 | fPrec(0), | |
102 | fChebCalc(1), | |
103 | fMaxCoefs(0), | |
104 | fResTmp(0), | |
105 | fGrid(0), | |
106 | fUsrFunName(""), | |
107 | fUsrMacro(0) | |
c437b1a5 | 108 | { |
40389866 | 109 | // read coefs from stream |
c8f37c50 | 110 | for (int i=3;i--;) { |
111 | fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
112 | fNPoints[i] = 0; | |
34473f9e | 113 | fGridOffs[i] = 0; |
f9efbbf6 | 114 | fArgsTmp[i] = 0; |
c8f37c50 | 115 | } |
40389866 | 116 | LoadData(stream); |
0bc7b414 | 117 | } |
0eea9d4d | 118 | |
119 | //__________________________________________________________________________________________ | |
120 | #ifdef _INC_CREATION_ALICHEB3D_ | |
5e6a3170 | 121 | AliCheb3D::AliCheb3D(const char* funName, int DimOut, const Float_t *bmin, const Float_t *bmax, Int_t *npoints, Float_t prec) : |
40389866 | 122 | TNamed(funName,funName), |
123 | fDimOut(0), | |
124 | fPrec(TMath::Max(1.E-12f,prec)), | |
125 | fChebCalc(1), | |
126 | fMaxCoefs(0), | |
127 | fResTmp(0), | |
128 | fGrid(0), | |
129 | fUsrFunName("") , | |
130 | fUsrMacro(0) | |
0eea9d4d | 131 | { |
132 | // Construct the parameterization for the function | |
133 | // funName : name of the file containing the function: void funName(Float_t * inp,Float_t * out) | |
134 | // DimOut : dimension of the vector computed by the user function | |
135 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
136 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
137 | // npoints : array of 3 elements with the number of points to compute in each of 3 dimension | |
138 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
139 | // | |
0eea9d4d | 140 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} |
c8f37c50 | 141 | for (int i=3;i--;) { |
142 | fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
143 | fNPoints[i] = 0; | |
144 | fGridOffs[i] = 0.; | |
f9efbbf6 | 145 | fArgsTmp[i] = 0; |
c8f37c50 | 146 | } |
0eea9d4d | 147 | SetDimOut(DimOut); |
148 | PrepareBoundaries(bmin,bmax); | |
149 | DefineGrid(npoints); | |
150 | SetUsrFunction(funName); | |
151 | ChebFit(); | |
152 | // | |
153 | } | |
154 | #endif | |
155 | ||
156 | //__________________________________________________________________________________________ | |
157 | #ifdef _INC_CREATION_ALICHEB3D_ | |
40389866 | 158 | AliCheb3D::AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec) : |
159 | fDimOut(0), | |
160 | fPrec(TMath::Max(1.E-12f,prec)), | |
161 | fChebCalc(1), | |
162 | fMaxCoefs(0), | |
163 | fResTmp(0), | |
164 | fGrid(0), | |
165 | fUsrFunName(""), | |
166 | fUsrMacro(0) | |
0eea9d4d | 167 | { |
168 | // Construct the parameterization for the function | |
169 | // ptr : pointer on the function: void fun(Float_t * inp,Float_t * out) | |
170 | // DimOut : dimension of the vector computed by the user function | |
171 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
172 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
173 | // npoints : array of 3 elements with the number of points to compute in each of 3 dimension | |
174 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
175 | // | |
0eea9d4d | 176 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} |
c8f37c50 | 177 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} |
178 | for (int i=3;i--;) { | |
179 | fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
180 | fNPoints[i] = 0; | |
181 | fGridOffs[i] = 0.; | |
f9efbbf6 | 182 | fArgsTmp[i] = 0; |
c8f37c50 | 183 | } |
0eea9d4d | 184 | SetDimOut(DimOut); |
185 | PrepareBoundaries(bmin,bmax); | |
186 | DefineGrid(npoints); | |
187 | SetUsrFunction(ptr); | |
188 | ChebFit(); | |
189 | // | |
190 | } | |
191 | #endif | |
192 | ||
40389866 | 193 | //__________________________________________________________________________________________ |
194 | #ifdef _INC_CREATION_ALICHEB3D_ | |
195 | AliCheb3D::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) : | |
196 | fDimOut(0), | |
197 | fPrec(TMath::Max(1.E-12f,prec)), | |
198 | fChebCalc(1), | |
199 | fMaxCoefs(0), | |
200 | fResTmp(0), | |
201 | fGrid(0), | |
202 | fUsrFunName(""), | |
203 | fUsrMacro(0) | |
204 | { | |
205 | // Construct very economic parameterization for the function | |
206 | // ptr : pointer on the function: void fun(Float_t * inp,Float_t * out) | |
207 | // DimOut : dimension of the vector computed by the user function | |
208 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
209 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
210 | // npX : array of 3 elements with the number of points to compute in each dimension for 1st component | |
211 | // npY : array of 3 elements with the number of points to compute in each dimension for 2nd component | |
212 | // npZ : array of 3 elements with the number of points to compute in each dimension for 3d component | |
213 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
214 | // | |
215 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} | |
c8f37c50 | 216 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} |
217 | for (int i=3;i--;) { | |
218 | fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
219 | fNPoints[i] = 0; | |
220 | fGridOffs[i] = 0.; | |
f9efbbf6 | 221 | fArgsTmp[i] = 0; |
c8f37c50 | 222 | } |
40389866 | 223 | SetDimOut(DimOut); |
224 | PrepareBoundaries(bmin,bmax); | |
225 | SetUsrFunction(ptr); | |
226 | // | |
227 | DefineGrid(npX); | |
228 | ChebFit(0); | |
229 | DefineGrid(npY); | |
230 | ChebFit(1); | |
231 | DefineGrid(npZ); | |
232 | ChebFit(2); | |
233 | // | |
234 | } | |
235 | #endif | |
236 | ||
237 | ||
238 | //__________________________________________________________________________________________ | |
239 | #ifdef _INC_CREATION_ALICHEB3D_ | |
1d18ebe0 | 240 | AliCheb3D::AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Float_t prec, Bool_t run) : |
40389866 | 241 | fDimOut(0), |
242 | fPrec(TMath::Max(1.E-12f,prec)), | |
243 | fChebCalc(1), | |
244 | fMaxCoefs(0), | |
245 | fResTmp(0), | |
246 | fGrid(0), | |
247 | fUsrFunName(""), | |
248 | fUsrMacro(0) | |
249 | { | |
250 | // Construct very economic parameterization for the function with automatic calculation of the root's grid | |
251 | // ptr : pointer on the function: void fun(Float_t * inp,Float_t * out) | |
252 | // DimOut : dimension of the vector computed by the user function | |
253 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
254 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
255 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
256 | // | |
257 | if (DimOut!=3) {Error("AliCheb3D","This constructor works only for 3D fits, %dD fit was requested\n",fDimOut); exit(1);} | |
c8f37c50 | 258 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} |
259 | for (int i=3;i--;) { | |
260 | fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
261 | fNPoints[i] = 0; | |
262 | fGridOffs[i] = 0.; | |
f9efbbf6 | 263 | fArgsTmp[i] = 0; |
c8f37c50 | 264 | } |
40389866 | 265 | SetDimOut(DimOut); |
266 | PrepareBoundaries(bmin,bmax); | |
267 | SetUsrFunction(ptr); | |
268 | // | |
1d18ebe0 | 269 | if (run) { |
270 | int gridNC[3][3]; | |
271 | EstimateNPoints(prec,gridNC); | |
272 | DefineGrid(gridNC[0]); | |
273 | ChebFit(0); | |
274 | DefineGrid(gridNC[1]); | |
275 | ChebFit(1); | |
276 | DefineGrid(gridNC[2]); | |
277 | ChebFit(2); | |
278 | } | |
40389866 | 279 | // |
280 | } | |
281 | #endif | |
282 | ||
283 | ||
284 | //__________________________________________________________________________________________ | |
285 | AliCheb3D& AliCheb3D::operator=(const AliCheb3D& rhs) | |
286 | { | |
5406439e | 287 | // assignment operator |
288 | // | |
40389866 | 289 | if (this != &rhs) { |
290 | Clear(); | |
291 | fDimOut = rhs.fDimOut; | |
292 | fPrec = rhs.fPrec; | |
293 | fMaxCoefs = rhs.fMaxCoefs; | |
294 | fUsrFunName = rhs.fUsrFunName; | |
295 | fUsrMacro = 0; | |
296 | for (int i=3;i--;) { | |
297 | fBMin[i] = rhs.fBMin[i]; | |
298 | fBMax[i] = rhs.fBMax[i]; | |
299 | fBScale[i] = rhs.fBScale[i]; | |
300 | fBOffset[i] = rhs.fBOffset[i]; | |
301 | fNPoints[i] = rhs.fNPoints[i]; | |
302 | } | |
303 | for (int i=0;i<fDimOut;i++) { | |
304 | AliCheb3DCalc* cbc = rhs.GetChebCalc(i); | |
305 | if (cbc) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(*cbc),i); | |
306 | } | |
307 | } | |
308 | return *this; | |
309 | // | |
310 | } | |
0bc7b414 | 311 | |
0eea9d4d | 312 | //__________________________________________________________________________________________ |
5406439e | 313 | void AliCheb3D::Clear(const Option_t*) |
0eea9d4d | 314 | { |
5406439e | 315 | // clear all dynamic structures |
316 | // | |
0eea9d4d | 317 | if (fResTmp) { delete[] fResTmp; fResTmp = 0; } |
318 | if (fGrid) { delete[] fGrid; fGrid = 0; } | |
319 | if (fUsrMacro) { delete fUsrMacro; fUsrMacro = 0;} | |
d8d3669f | 320 | fChebCalc.SetOwner(kTRUE); |
0eea9d4d | 321 | fChebCalc.Delete(); |
322 | // | |
323 | } | |
324 | ||
325 | //__________________________________________________________________________________________ | |
5406439e | 326 | void AliCheb3D::Print(const Option_t* opt) const |
0eea9d4d | 327 | { |
5406439e | 328 | // print info |
329 | // | |
0eea9d4d | 330 | printf("%s: Chebyshev parameterization for 3D->%dD function. Precision: %e\n",GetName(),fDimOut,fPrec); |
331 | printf("Region of validity: [%+.5e:%+.5e] [%+.5e:%+.5e] [%+.5e:%+.5e]\n",fBMin[0],fBMax[0],fBMin[1],fBMax[1],fBMin[2],fBMax[2]); | |
332 | TString opts = opt; opts.ToLower(); | |
333 | if (opts.Contains("l")) for (int i=0;i<fDimOut;i++) {printf("Output dimension %d:\n",i+1); GetChebCalc(i)->Print();} | |
334 | // | |
335 | } | |
336 | ||
0eea9d4d | 337 | //__________________________________________________________________________________________ |
5406439e | 338 | void AliCheb3D::PrepareBoundaries(const Float_t *bmin, const Float_t *bmax) |
0eea9d4d | 339 | { |
340 | // Set and check boundaries defined by user, prepare coefficients for their conversion to [-1:1] interval | |
341 | // | |
342 | for (int i=3;i--;) { | |
343 | fBMin[i] = bmin[i]; | |
344 | fBMax[i] = bmax[i]; | |
345 | fBScale[i] = bmax[i]-bmin[i]; | |
346 | if (fBScale[i]<=0) { | |
b8bec44f | 347 | AliFatal(Form("Boundaries for %d-th dimension are not increasing: %+.4e %+.4e\nStop\n",i,fBMin[i],fBMax[i])); |
0eea9d4d | 348 | } |
349 | fBOffset[i] = bmin[i] + fBScale[i]/2.0; | |
350 | fBScale[i] = 2./fBScale[i]; | |
351 | } | |
352 | // | |
353 | } | |
354 | ||
40389866 | 355 | |
0eea9d4d | 356 | //__________________________________________________________________________________________ |
357 | #ifdef _INC_CREATION_ALICHEB3D_ | |
40389866 | 358 | |
359 | // Pointer on user function (faster altrnative to TMethodCall) | |
360 | void (*gUsrFunAliCheb3D) (float* ,float* ); | |
361 | ||
362 | void AliCheb3D::EvalUsrFunction() | |
363 | { | |
364 | // call user supplied function | |
365 | if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); | |
366 | else fUsrMacro->Execute(); | |
367 | } | |
368 | ||
0eea9d4d | 369 | void AliCheb3D::SetUsrFunction(const char* name) |
370 | { | |
371 | // load user macro with function definition and compile it | |
372 | gUsrFunAliCheb3D = 0; | |
373 | fUsrFunName = name; | |
374 | gSystem->ExpandPathName(fUsrFunName); | |
375 | if (fUsrMacro) delete fUsrMacro; | |
376 | TString tmpst = fUsrFunName; | |
377 | tmpst += "+"; // prepare filename to compile | |
378 | if (gROOT->LoadMacro(tmpst.Data())) {Error("SetUsrFunction","Failed to load user function from %s\nStop\n",name); exit(1);} | |
379 | fUsrMacro = new TMethodCall(); | |
380 | tmpst = tmpst.Data() + tmpst.Last('/')+1; //Strip away any path preceding the macro file name | |
381 | int dot = tmpst.Last('.'); | |
382 | if (dot>0) tmpst.Resize(dot); | |
383 | fUsrMacro->InitWithPrototype(tmpst.Data(),"Float_t *,Float_t *"); | |
384 | long args[2]; | |
385 | args[0] = (long)fArgsTmp; | |
386 | args[1] = (long)fResTmp; | |
387 | fUsrMacro->SetParamPtrs(args); | |
388 | // | |
389 | } | |
390 | #endif | |
391 | ||
392 | //__________________________________________________________________________________________ | |
393 | #ifdef _INC_CREATION_ALICHEB3D_ | |
394 | void AliCheb3D::SetUsrFunction(void (*ptr)(float*,float*)) | |
395 | { | |
5406439e | 396 | // assign user training function |
397 | // | |
0eea9d4d | 398 | if (fUsrMacro) delete fUsrMacro; |
399 | fUsrMacro = 0; | |
400 | fUsrFunName = ""; | |
401 | gUsrFunAliCheb3D = ptr; | |
402 | } | |
403 | #endif | |
404 | ||
405 | //__________________________________________________________________________________________ | |
406 | #ifdef _INC_CREATION_ALICHEB3D_ | |
da7cd221 | 407 | void AliCheb3D::EvalUsrFunction(const Float_t *x, Float_t *res) |
5406439e | 408 | { |
409 | // evaluate user function value | |
410 | // | |
0eea9d4d | 411 | for (int i=3;i--;) fArgsTmp[i] = x[i]; |
412 | if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); | |
413 | else fUsrMacro->Execute(); | |
414 | for (int i=fDimOut;i--;) res[i] = fResTmp[i]; | |
415 | } | |
416 | #endif | |
417 | ||
418 | //__________________________________________________________________________________________ | |
419 | #ifdef _INC_CREATION_ALICHEB3D_ | |
5406439e | 420 | Int_t AliCheb3D::CalcChebCoefs(const Float_t *funval,int np, Float_t *outCoefs, Float_t prec) |
0eea9d4d | 421 | { |
422 | // Calculate Chebyshev coeffs using precomputed function values at np roots. | |
423 | // If prec>0, estimate the highest coeff number providing the needed precision | |
424 | // | |
425 | double sm; // do summations in double to minimize the roundoff error | |
426 | for (int ic=0;ic<np;ic++) { // compute coeffs | |
427 | sm = 0; | |
428 | for (int ir=0;ir<np;ir++) { | |
429 | float rt = TMath::Cos( ic*(ir+0.5)*TMath::Pi()/np); | |
430 | sm += funval[ir]*rt; | |
431 | } | |
432 | outCoefs[ic] = Float_t( sm * ((ic==0) ? 1./np : 2./np) ); | |
433 | } | |
434 | // | |
435 | if (prec<=0) return np; | |
436 | // | |
437 | sm = 0; | |
438 | int cfMax = 0; | |
439 | for (cfMax=np;cfMax--;) { | |
440 | sm += TMath::Abs(outCoefs[cfMax]); | |
441 | if (sm>=prec) break; | |
442 | } | |
443 | if (++cfMax==0) cfMax=1; | |
444 | return cfMax; | |
445 | // | |
446 | } | |
447 | #endif | |
448 | ||
449 | //__________________________________________________________________________________________ | |
450 | #ifdef _INC_CREATION_ALICHEB3D_ | |
451 | void AliCheb3D::DefineGrid(Int_t* npoints) | |
452 | { | |
453 | // prepare the grid of Chebyshev roots in each dimension | |
454 | const int kMinPoints = 1; | |
455 | int ntot = 0; | |
456 | fMaxCoefs = 1; | |
457 | for (int id=3;id--;) { | |
458 | fNPoints[id] = npoints[id]; | |
459 | if (fNPoints[id]<kMinPoints) { | |
40389866 | 460 | Error("DefineGrid","at %d-th dimension %d point is requested, at least %d is needed\nStop\n",id,fNPoints[id],kMinPoints); |
0eea9d4d | 461 | exit(1); |
462 | } | |
463 | ntot += fNPoints[id]; | |
464 | fMaxCoefs *= fNPoints[id]; | |
465 | } | |
40389866 | 466 | printf("Computing Chebyshev nodes on [%2d/%2d/%2d] grid\n",npoints[0],npoints[1],npoints[2]); |
467 | if (fGrid) delete[] fGrid; | |
0eea9d4d | 468 | fGrid = new Float_t [ntot]; |
469 | // | |
470 | int curp = 0; | |
471 | for (int id=3;id--;) { | |
472 | int np = fNPoints[id]; | |
473 | fGridOffs[id] = curp; | |
474 | for (int ip=0;ip<np;ip++) { | |
475 | Float_t x = TMath::Cos( TMath::Pi()*(ip+0.5)/np ); | |
476 | fGrid[curp++] = MapToExternal(x,id); | |
477 | } | |
478 | } | |
479 | // | |
480 | } | |
481 | #endif | |
482 | ||
483 | //__________________________________________________________________________________________ | |
484 | #ifdef _INC_CREATION_ALICHEB3D_ | |
485 | Int_t AliCheb3D::ChebFit() | |
486 | { | |
487 | // prepare parameterization for all output dimensions | |
488 | int ir=0; | |
489 | for (int i=fDimOut;i--;) ir+=ChebFit(i); | |
490 | return ir; | |
491 | } | |
492 | #endif | |
493 | ||
494 | //__________________________________________________________________________________________ | |
495 | #ifdef _INC_CREATION_ALICHEB3D_ | |
496 | Int_t AliCheb3D::ChebFit(int dmOut) | |
497 | { | |
498 | // prepare paramaterization of 3D function for dmOut-th dimension | |
499 | int maxDim = 0; | |
500 | for (int i=0;i<3;i++) if (maxDim<fNPoints[i]) maxDim = fNPoints[i]; | |
501 | Float_t *fvals = new Float_t [ fNPoints[0] ]; | |
502 | Float_t *tmpCoef3D = new Float_t [ fNPoints[0]*fNPoints[1]*fNPoints[2] ]; | |
503 | Float_t *tmpCoef2D = new Float_t [ fNPoints[0]*fNPoints[1] ]; | |
504 | Float_t *tmpCoef1D = new Float_t [ maxDim ]; | |
505 | // | |
1d18ebe0 | 506 | Float_t rTiny = 0.1*fPrec/Float_t(maxDim); // neglect coefficient below this threshold |
0eea9d4d | 507 | // |
508 | // 1D Cheb.fit for 0-th dimension at current steps of remaining dimensions | |
509 | int ncmax = 0; | |
510 | // | |
40389866 | 511 | printf("Dim%d : 00.00%% Done",dmOut);fflush(stdout); |
0eea9d4d | 512 | AliCheb3DCalc* cheb = GetChebCalc(dmOut); |
513 | // | |
40389866 | 514 | float ncals2count = fNPoints[2]*fNPoints[1]*fNPoints[0]; |
515 | float ncals = 0; | |
516 | float frac = 0; | |
517 | float fracStep = 0.001; | |
518 | // | |
0eea9d4d | 519 | for (int id2=fNPoints[2];id2--;) { |
520 | fArgsTmp[2] = fGrid[ fGridOffs[2]+id2 ]; | |
521 | // | |
522 | for (int id1=fNPoints[1];id1--;) { | |
523 | fArgsTmp[1] = fGrid[ fGridOffs[1]+id1 ]; | |
524 | // | |
525 | for (int id0=fNPoints[0];id0--;) { | |
526 | fArgsTmp[0] = fGrid[ fGridOffs[0]+id0 ]; | |
527 | EvalUsrFunction(); // compute function values at Chebyshev roots of 0-th dimension | |
528 | fvals[id0] = fResTmp[dmOut]; | |
40389866 | 529 | float fr = (++ncals)/ncals2count; |
530 | if (fr-frac>=fracStep) { | |
531 | frac = fr; | |
532 | printf("\b\b\b\b\b\b\b\b\b\b\b"); | |
533 | printf("%05.2f%% Done",fr*100); | |
534 | fflush(stdout); | |
535 | } | |
536 | // | |
0eea9d4d | 537 | } |
538 | int nc = CalcChebCoefs(fvals,fNPoints[0], tmpCoef1D, fPrec); | |
539 | for (int id0=fNPoints[0];id0--;) tmpCoef2D[id1 + id0*fNPoints[1]] = tmpCoef1D[id0]; | |
540 | if (ncmax<nc) ncmax = nc; // max coefs to be kept in dim0 to guarantee needed precision | |
541 | } | |
542 | // | |
543 | // once each 1d slice of given 2d slice is parametrized, parametrize the Cheb.coeffs | |
544 | for (int id0=fNPoints[0];id0--;) { | |
545 | CalcChebCoefs( tmpCoef2D+id0*fNPoints[1], fNPoints[1], tmpCoef1D, -1); | |
546 | for (int id1=fNPoints[1];id1--;) tmpCoef3D[id2 + fNPoints[2]*(id1+id0*fNPoints[1])] = tmpCoef1D[id1]; | |
547 | } | |
548 | } | |
549 | // | |
550 | // now fit the last dimensions Cheb.coefs | |
551 | for (int id0=fNPoints[0];id0--;) { | |
552 | for (int id1=fNPoints[1];id1--;) { | |
553 | CalcChebCoefs( tmpCoef3D+ fNPoints[2]*(id1+id0*fNPoints[1]), fNPoints[2], tmpCoef1D, -1); | |
554 | for (int id2=fNPoints[2];id2--;) tmpCoef3D[id2+ fNPoints[2]*(id1+id0*fNPoints[1])] = tmpCoef1D[id2]; // store on place | |
555 | } | |
556 | } | |
557 | // | |
558 | // now find 2D surface which separates significant coefficients of 3D matrix from nonsignificant ones (up to fPrec) | |
1d18ebe0 | 559 | UShort_t *tmpCoefSurf = new UShort_t[ fNPoints[0]*fNPoints[1] ]; |
0eea9d4d | 560 | for (int id0=fNPoints[0];id0--;) for (int id1=fNPoints[1];id1--;) tmpCoefSurf[id1+id0*fNPoints[1]]=0; |
561 | Double_t resid = 0; | |
562 | for (int id0=fNPoints[0];id0--;) { | |
563 | for (int id1=fNPoints[1];id1--;) { | |
564 | for (int id2=fNPoints[2];id2--;) { | |
565 | int id = id2 + fNPoints[2]*(id1+id0*fNPoints[1]); | |
566 | Float_t cfa = TMath::Abs(tmpCoef3D[id]); | |
5406439e | 567 | if (cfa < rTiny) {tmpCoef3D[id] = 0; continue;} // neglect coefs below the threshold |
0eea9d4d | 568 | resid += cfa; |
569 | if (resid<fPrec) continue; // this coeff is negligible | |
570 | // otherwise go back 1 step | |
571 | resid -= cfa; | |
572 | tmpCoefSurf[id1+id0*fNPoints[1]] = id2+1; // how many coefs to keep | |
573 | break; | |
574 | } | |
575 | } | |
576 | } | |
577 | /* | |
578 | printf("\n\nCoeffs\n"); | |
579 | int cnt = 0; | |
580 | for (int id0=0;id0<fNPoints[0];id0++) { | |
581 | for (int id1=0;id1<fNPoints[1];id1++) { | |
582 | for (int id2=0;id2<fNPoints[2];id2++) { | |
583 | printf("%2d%2d%2d %+.4e |",id0,id1,id2,tmpCoef3D[cnt++]); | |
584 | } | |
585 | printf("\n"); | |
586 | } | |
587 | printf("\n"); | |
588 | } | |
589 | */ | |
590 | // see if there are rows to reject, find max.significant column at each row | |
5406439e | 591 | int nRows = fNPoints[0]; |
1d18ebe0 | 592 | UShort_t *tmpCols = new UShort_t[nRows]; |
0eea9d4d | 593 | for (int id0=fNPoints[0];id0--;) { |
594 | int id1 = fNPoints[1]; | |
595 | while (id1>0 && tmpCoefSurf[(id1-1)+id0*fNPoints[1]]==0) id1--; | |
596 | tmpCols[id0] = id1; | |
597 | } | |
598 | // find max significant row | |
5406439e | 599 | for (int id0=nRows;id0--;) {if (tmpCols[id0]>0) break; nRows--;} |
0eea9d4d | 600 | // find max significant column and fill the permanent storage for the max sigificant column of each row |
5406439e | 601 | cheb->InitRows(nRows); // create needed arrays; |
1d18ebe0 | 602 | UShort_t *nColsAtRow = cheb->GetNColsAtRow(); |
603 | UShort_t *colAtRowBg = cheb->GetColAtRowBg(); | |
5406439e | 604 | int nCols = 0; |
8f499a66 | 605 | int nElemBound2D = 0; |
5406439e | 606 | for (int id0=0;id0<nRows;id0++) { |
607 | nColsAtRow[id0] = tmpCols[id0]; // number of columns to store for this row | |
8f499a66 | 608 | colAtRowBg[id0] = nElemBound2D; // begining of this row in 2D boundary surface |
609 | nElemBound2D += tmpCols[id0]; | |
5406439e | 610 | if (nCols<nColsAtRow[id0]) nCols = nColsAtRow[id0]; |
0eea9d4d | 611 | } |
5406439e | 612 | cheb->InitCols(nCols); |
0eea9d4d | 613 | delete[] tmpCols; |
614 | // | |
615 | // create the 2D matrix defining the boundary of significance for 3D coeffs.matrix | |
616 | // and count the number of siginifacnt coefficients | |
617 | // | |
8f499a66 | 618 | cheb->InitElemBound2D(nElemBound2D); |
1d18ebe0 | 619 | UShort_t *coefBound2D0 = cheb->GetCoefBound2D0(); |
620 | UShort_t *coefBound2D1 = cheb->GetCoefBound2D1(); | |
0eea9d4d | 621 | fMaxCoefs = 0; // redefine number of coeffs |
5406439e | 622 | for (int id0=0;id0<nRows;id0++) { |
623 | int nCLoc = nColsAtRow[id0]; | |
624 | int col0 = colAtRowBg[id0]; | |
0eea9d4d | 625 | for (int id1=0;id1<nCLoc;id1++) { |
5406439e | 626 | coefBound2D0[col0 + id1] = tmpCoefSurf[id1+id0*fNPoints[1]]; // number of coefs to store for 3-d dimension |
627 | coefBound2D1[col0 + id1] = fMaxCoefs; | |
628 | fMaxCoefs += coefBound2D0[col0 + id1]; | |
0eea9d4d | 629 | } |
630 | } | |
631 | // | |
632 | // create final compressed 3D matrix for significant coeffs | |
633 | cheb->InitCoefs(fMaxCoefs); | |
5406439e | 634 | Float_t *coefs = cheb->GetCoefs(); |
0eea9d4d | 635 | int count = 0; |
5406439e | 636 | for (int id0=0;id0<nRows;id0++) { |
637 | int ncLoc = nColsAtRow[id0]; | |
638 | int col0 = colAtRowBg[id0]; | |
0eea9d4d | 639 | for (int id1=0;id1<ncLoc;id1++) { |
5406439e | 640 | int ncf2 = coefBound2D0[col0 + id1]; |
0eea9d4d | 641 | for (int id2=0;id2<ncf2;id2++) { |
5406439e | 642 | coefs[count++] = tmpCoef3D[id2 + fNPoints[2]*(id1+id0*fNPoints[1])]; |
0eea9d4d | 643 | } |
644 | } | |
645 | } | |
646 | /* | |
647 | printf("\n\nNewSurf\n"); | |
648 | for (int id0=0;id0<fNPoints[0];id0++) { | |
649 | for (int id1=0;id1<fNPoints[1];id1++) { | |
650 | printf("(%2d %2d) %2d |",id0,id1,tmpCoefSurf[id1+id0*fNPoints[1]]); | |
651 | } | |
652 | printf("\n"); | |
653 | } | |
654 | */ | |
655 | // | |
656 | delete[] tmpCoefSurf; | |
657 | delete[] tmpCoef1D; | |
658 | delete[] tmpCoef2D; | |
659 | delete[] tmpCoef3D; | |
660 | delete[] fvals; | |
661 | // | |
40389866 | 662 | printf("\b\b\b\b\b\b\b\b\b\b\b\b"); |
663 | printf("100.00%% Done\n"); | |
0eea9d4d | 664 | return 1; |
665 | } | |
666 | #endif | |
667 | ||
668 | //_______________________________________________ | |
669 | #ifdef _INC_CREATION_ALICHEB3D_ | |
670 | void AliCheb3D::SaveData(const char* outfile,Bool_t append) const | |
671 | { | |
672 | // writes coefficients data to output text file, optionallt appending on the end of existing file | |
673 | TString strf = outfile; | |
674 | gSystem->ExpandPathName(strf); | |
675 | FILE* stream = fopen(strf,append ? "a":"w"); | |
676 | SaveData(stream); | |
677 | fclose(stream); | |
678 | // | |
679 | } | |
680 | #endif | |
681 | ||
682 | //_______________________________________________ | |
683 | #ifdef _INC_CREATION_ALICHEB3D_ | |
684 | void AliCheb3D::SaveData(FILE* stream) const | |
685 | { | |
686 | // writes coefficients data to existing output stream | |
687 | // | |
688 | fprintf(stream,"\n# These are automatically generated data for the Chebyshev interpolation of 3D->%dD function\n",fDimOut); | |
689 | fprintf(stream,"#\nSTART %s\n",GetName()); | |
690 | fprintf(stream,"# Dimensionality of the output\n%d\n",fDimOut); | |
691 | fprintf(stream,"# Interpolation abs. precision\n%+.8e\n",fPrec); | |
692 | // | |
693 | fprintf(stream,"# Lower boundaries of interpolation region\n"); | |
694 | for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMin[i]); | |
695 | fprintf(stream,"# Upper boundaries of interpolation region\n"); | |
696 | for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMax[i]); | |
40389866 | 697 | fprintf(stream,"# Parameterization for each output dimension follows:\n"); |
0eea9d4d | 698 | // |
699 | for (int i=0;i<fDimOut;i++) GetChebCalc(i)->SaveData(stream); | |
700 | fprintf(stream,"#\nEND %s\n#\n",GetName()); | |
701 | // | |
702 | } | |
703 | #endif | |
704 | ||
7f42fc51 | 705 | //__________________________________________________________________________________________ |
706 | #ifdef _INC_CREATION_ALICHEB3D_ | |
707 | void AliCheb3D::InvertSign() | |
708 | { | |
709 | // invert the sign of all parameterizations | |
710 | for (int i=fDimOut;i--;) { | |
711 | AliCheb3DCalc* par = GetChebCalc(i); | |
712 | int ncf = par->GetNCoefs(); | |
713 | float *coefs = par->GetCoefs(); | |
714 | for (int j=ncf;j--;) coefs[j] = -coefs[j]; | |
715 | } | |
716 | } | |
717 | #endif | |
718 | ||
719 | ||
0eea9d4d | 720 | //_______________________________________________ |
721 | void AliCheb3D::LoadData(const char* inpFile) | |
722 | { | |
5406439e | 723 | // load coefficients data from txt file |
724 | // | |
0eea9d4d | 725 | TString strf = inpFile; |
726 | gSystem->ExpandPathName(strf); | |
727 | FILE* stream = fopen(strf.Data(),"r"); | |
728 | LoadData(stream); | |
729 | fclose(stream); | |
730 | // | |
731 | } | |
732 | ||
733 | //_______________________________________________ | |
734 | void AliCheb3D::LoadData(FILE* stream) | |
735 | { | |
5406439e | 736 | // load coefficients data from stream |
737 | // | |
b8bec44f | 738 | if (!stream) {AliFatal("No stream provided.\nStop");} |
0eea9d4d | 739 | TString buffs; |
740 | Clear(); | |
741 | AliCheb3DCalc::ReadLine(buffs,stream); | |
b8bec44f | 742 | if (!buffs.BeginsWith("START")) {AliFatal(Form("Expected: \"START <fit_name>\", found \"%s\"\nStop\n",buffs.Data()));} |
0eea9d4d | 743 | SetName(buffs.Data()+buffs.First(' ')+1); |
744 | // | |
745 | AliCheb3DCalc::ReadLine(buffs,stream); // N output dimensions | |
746 | fDimOut = buffs.Atoi(); | |
b8bec44f | 747 | if (fDimOut<1) {AliFatal(Form("Expected: '<number_of_output_dimensions>', found \"%s\"\nStop\n",buffs.Data()));} |
0eea9d4d | 748 | // |
749 | SetDimOut(fDimOut); | |
750 | // | |
751 | AliCheb3DCalc::ReadLine(buffs,stream); // Interpolation abs. precision | |
752 | fPrec = buffs.Atof(); | |
b8bec44f | 753 | if (fPrec<=0) {AliFatal(Form("Expected: '<abs.precision>', found \"%s\"\nStop\n",buffs.Data()));} |
0eea9d4d | 754 | // |
755 | for (int i=0;i<3;i++) { // Lower boundaries of interpolation region | |
756 | AliCheb3DCalc::ReadLine(buffs,stream); | |
757 | fBMin[i] = buffs.Atof(); | |
758 | } | |
759 | for (int i=0;i<3;i++) { // Upper boundaries of interpolation region | |
760 | AliCheb3DCalc::ReadLine(buffs,stream); | |
761 | fBMax[i] = buffs.Atof(); | |
762 | } | |
763 | PrepareBoundaries(fBMin,fBMax); | |
764 | // | |
765 | // data for each output dimension | |
766 | for (int i=0;i<fDimOut;i++) GetChebCalc(i)->LoadData(stream); | |
767 | // | |
768 | // check end_of_data record | |
769 | AliCheb3DCalc::ReadLine(buffs,stream); | |
770 | if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) { | |
b8bec44f | 771 | AliFatal(Form("Expected \"END %s\", found \"%s\".\nStop\n",GetName(),buffs.Data())); |
0eea9d4d | 772 | } |
773 | // | |
774 | } | |
775 | ||
776 | //_______________________________________________ | |
5406439e | 777 | void AliCheb3D::SetDimOut(const int d) |
0eea9d4d | 778 | { |
5406439e | 779 | // init output dimensions |
0eea9d4d | 780 | fDimOut = d; |
781 | if (fResTmp) delete fResTmp; | |
40389866 | 782 | fResTmp = new Float_t[fDimOut]; |
0eea9d4d | 783 | fChebCalc.Delete(); |
784 | for (int i=0;i<d;i++) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(),i); | |
785 | } | |
786 | ||
787 | //_______________________________________________ | |
788 | void AliCheb3D::ShiftBound(int id,float dif) | |
789 | { | |
5406439e | 790 | // modify the bounds of the grid |
791 | // | |
0eea9d4d | 792 | if (id<0||id>2) {printf("Maximum 3 dimensions are supported\n"); return;} |
793 | fBMin[id] += dif; | |
794 | fBMax[id] += dif; | |
795 | fBOffset[id] += dif; | |
796 | } | |
797 | ||
798 | //_______________________________________________ | |
799 | #ifdef _INC_CREATION_ALICHEB3D_ | |
800 | TH1* AliCheb3D::TestRMS(int idim,int npoints,TH1* histo) | |
801 | { | |
802 | // fills the difference between the original function and parameterization (for idim-th component of the output) | |
803 | // to supplied histogram. Calculations are done in npoints random points. | |
804 | // If the hostgram was not supplied, it will be created. It is up to the user to delete it! | |
805 | if (!fUsrMacro) { | |
806 | printf("No user function is set\n"); | |
807 | return 0; | |
808 | } | |
809 | if (!histo) histo = new TH1D(GetName(),"Control: Function - Parametrization",100,-2*fPrec,2*fPrec); | |
810 | for (int ip=npoints;ip--;) { | |
811 | gRandom->RndmArray(3,(Float_t *)fArgsTmp); | |
812 | for (int i=3;i--;) fArgsTmp[i] = fBMin[i] + fArgsTmp[i]*(fBMax[i]-fBMin[i]); | |
813 | EvalUsrFunction(); | |
814 | Float_t valFun = fResTmp[idim]; | |
815 | Eval(fArgsTmp,fResTmp); | |
816 | Float_t valPar = fResTmp[idim]; | |
817 | histo->Fill(valFun - valPar); | |
818 | } | |
819 | return histo; | |
820 | // | |
821 | } | |
822 | #endif | |
40389866 | 823 | |
824 | //_______________________________________________ | |
825 | #ifdef _INC_CREATION_ALICHEB3D_ | |
1d18ebe0 | 826 | |
827 | void AliCheb3D::EstimateNPoints(float Prec, int gridBC[3][3],Int_t npd1,Int_t npd2,Int_t npd3) | |
828 | { | |
829 | // Estimate number of points to generate a training data | |
830 | // | |
831 | const int kScp = 9; | |
832 | const float kScl[9] = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}; | |
833 | // | |
834 | const float sclDim[2] = {0.001,0.999}; | |
835 | const int compDim[3][2] = { {1,2}, {2,0}, {0,1} }; | |
836 | static float xyz[3]; | |
837 | Int_t npdTst[3] = {npd1,npd2,npd3}; | |
838 | // | |
839 | ||
840 | for (int i=3;i--;)for (int j=3;j--;) gridBC[i][j] = -1; | |
841 | // | |
842 | for (int idim=0;idim<3;idim++) { | |
843 | float dimMN = fBMin[idim] + sclDim[0]*(fBMax[idim]-fBMin[idim]); | |
844 | float dimMX = fBMin[idim] + sclDim[1]*(fBMax[idim]-fBMin[idim]); | |
845 | // | |
846 | int id1 = compDim[idim][0]; // 1st fixed dim | |
847 | int id2 = compDim[idim][1]; // 2nd fixed dim | |
848 | for (int i1=0;i1<kScp;i1++) { | |
849 | xyz[ id1 ] = fBMin[id1] + kScl[i1]*( fBMax[id1]-fBMin[id1] ); | |
850 | for (int i2=0;i2<kScp;i2++) { | |
851 | xyz[ id2 ] = fBMin[id2] + kScl[i2]*( fBMax[id2]-fBMin[id2] ); | |
852 | int* npt = GetNCNeeded(xyz,idim, dimMN,dimMX, Prec, npdTst[idim]); // npoints for Bx,By,Bz | |
853 | for (int ib=0;ib<3;ib++) if (npt[ib]>gridBC[ib][idim]) gridBC[ib][idim] = npt[ib]; | |
854 | } | |
855 | } | |
856 | } | |
857 | } | |
858 | ||
859 | /* | |
40389866 | 860 | void AliCheb3D::EstimateNPoints(float Prec, int gridBC[3][3]) |
861 | { | |
5406439e | 862 | // Estimate number of points to generate a training data |
863 | // | |
40389866 | 864 | const float sclA[9] = {0.1, 0.5, 0.9, 0.1, 0.5, 0.9, 0.1, 0.5, 0.9} ; |
865 | const float sclB[9] = {0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9} ; | |
866 | const float sclDim[2] = {0.01,0.99}; | |
867 | const int compDim[3][2] = { {1,2}, {2,0}, {0,1} }; | |
868 | static float xyz[3]; | |
869 | // | |
870 | for (int i=3;i--;)for (int j=3;j--;) gridBC[i][j] = -1; | |
871 | // | |
872 | for (int idim=0;idim<3;idim++) { | |
873 | float dimMN = fBMin[idim] + sclDim[0]*(fBMax[idim]-fBMin[idim]); | |
874 | float dimMX = fBMin[idim] + sclDim[1]*(fBMax[idim]-fBMin[idim]); | |
875 | // | |
876 | for (int it=0;it<9;it++) { // test in 9 points | |
877 | int id1 = compDim[idim][0]; // 1st fixed dim | |
878 | int id2 = compDim[idim][1]; // 2nd fixed dim | |
879 | xyz[ id1 ] = fBMin[id1] + sclA[it]*( fBMax[id1]-fBMin[id1] ); | |
880 | xyz[ id2 ] = fBMin[id2] + sclB[it]*( fBMax[id2]-fBMin[id2] ); | |
881 | // | |
882 | int* npt = GetNCNeeded(xyz,idim, dimMN,dimMX, Prec); // npoints for Bx,By,Bz | |
2572efdf | 883 | for (int ib=0;ib<3;ib++) if (npt[ib]>gridBC[ib][idim]) gridBC[ib][idim] = npt[ib];//+2; |
40389866 | 884 | // |
885 | } | |
886 | } | |
887 | } | |
888 | ||
889 | ||
890 | int* AliCheb3D::GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec) | |
891 | { | |
2572efdf | 892 | // estimate needed number of chebyshev coefs for given function description in DimVar dimension |
40389866 | 893 | // The values for two other dimensions must be set beforehand |
894 | // | |
895 | static int curNC[3]; | |
896 | static int retNC[3]; | |
897 | const int kMaxPoint = 400; | |
898 | float* gridVal = new float[3*kMaxPoint]; | |
899 | float* coefs = new float[3*kMaxPoint]; | |
900 | // | |
901 | float scale = mx-mn; | |
902 | float offs = mn + scale/2.0; | |
903 | scale = 2./scale; | |
904 | // | |
905 | int curNP; | |
906 | int maxNC=-1; | |
907 | int maxNCPrev=-1; | |
908 | for (int i=0;i<3;i++) retNC[i] = -1; | |
909 | for (int i=0;i<3;i++) fArgsTmp[i] = xyz[i]; | |
910 | // | |
2572efdf | 911 | for (curNP=3; curNP<kMaxPoint; curNP+=3) { |
40389866 | 912 | maxNCPrev = maxNC; |
913 | // | |
914 | for (int i=0;i<curNP;i++) { // get function values on Cheb. nodes | |
915 | float x = TMath::Cos( TMath::Pi()*(i+0.5)/curNP ); | |
916 | fArgsTmp[DimVar] = x/scale+offs; // map to requested interval | |
917 | EvalUsrFunction(); | |
918 | for (int ib=3;ib--;) gridVal[ib*kMaxPoint + i] = fResTmp[ib]; | |
919 | } | |
920 | // | |
921 | for (int ib=0;ib<3;ib++) { | |
922 | curNC[ib] = AliCheb3D::CalcChebCoefs(&gridVal[ib*kMaxPoint], curNP, &coefs[ib*kMaxPoint],prec); | |
923 | if (maxNC < curNC[ib]) maxNC = curNC[ib]; | |
924 | if (retNC[ib] < curNC[ib]) retNC[ib] = curNC[ib]; | |
925 | } | |
926 | if ( (curNP-maxNC)>3 && (maxNC-maxNCPrev)<1 ) break; | |
927 | maxNCPrev = maxNC; | |
928 | // | |
929 | } | |
930 | delete[] gridVal; | |
931 | delete[] coefs; | |
932 | return retNC; | |
933 | // | |
934 | } | |
1d18ebe0 | 935 | */ |
936 | ||
937 | ||
938 | int* AliCheb3D::GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec, Int_t npCheck) | |
939 | { | |
940 | // estimate needed number of chebyshev coefs for given function description in DimVar dimension | |
941 | // The values for two other dimensions must be set beforehand | |
942 | // | |
943 | static int retNC[3]; | |
944 | static int npChLast = 0; | |
945 | static float *gridVal=0,*coefs=0; | |
946 | if (npCheck<3) npCheck = 3; | |
947 | if (npChLast<npCheck) { | |
948 | if (gridVal) delete[] gridVal; | |
949 | if (coefs) delete[] coefs; | |
950 | gridVal = new float[3*npCheck]; | |
951 | coefs = new float[3*npCheck]; | |
952 | npChLast = npCheck; | |
953 | } | |
954 | // | |
955 | float scale = mx-mn; | |
956 | float offs = mn + scale/2.0; | |
957 | scale = 2./scale; | |
958 | // | |
959 | for (int i=0;i<3;i++) fArgsTmp[i] = xyz[i]; | |
960 | for (int i=0;i<npCheck;i++) { | |
961 | fArgsTmp[DimVar] = TMath::Cos( TMath::Pi()*(i+0.5)/npCheck)/scale+offs; // map to requested interval | |
962 | EvalUsrFunction(); | |
963 | for (int ib=3;ib--;) gridVal[ib*npCheck + i] = fResTmp[ib]; | |
964 | } | |
965 | // | |
966 | for (int ib=0;ib<3;ib++) retNC[ib] = AliCheb3D::CalcChebCoefs(&gridVal[ib*npCheck], npCheck, &coefs[ib*npCheck],prec); | |
967 | return retNC; | |
968 | // | |
969 | } | |
970 | ||
971 | ||
40389866 | 972 | #endif |