]>
Commit | Line | Data |
---|---|---|
0eea9d4d | 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 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
18 | // Author: ruben.shahoyan@cern.ch 09/09/2006 | |
19 | // | |
20 | //////////////////////////////////////////////////////////////////////////////// | |
21 | // // | |
22 | // AliCheb3D produces the interpolation of the user 3D->NDimOut arbitrary // | |
23 | // function supplied in "void (*fcn)(float* inp,float* out)" format // | |
24 | // either in a separate macro file or as a function pointer. // | |
25 | // Only coefficients needed to guarantee the requested precision are kept. // | |
26 | // // | |
27 | // The user-callable methods are: // | |
28 | // To create the interpolation use: // | |
29 | // AliCheb3D(const char* funName, // name of the file with user function // | |
30 | // or // | |
31 | // AliCheb3D(void (*ptr)(float*,float*),// pointer on the user function // | |
32 | // Int_t DimOut, // dimensionality of the function's output // | |
33 | // Float_t *bmin, // lower 3D bounds of interpolation domain // | |
34 | // Float_t *bmax, // upper 3D bounds of interpolation domain // | |
35 | // Int_t *npoints, // number of points in each of 3 input // | |
36 | // // dimension, defining the interpolation grid // | |
37 | // Float_t prec=1E-6); // requested max.absolute difference between // | |
38 | // // the interpolation and any point on grid // | |
39 | // // | |
40 | // To test obtained parameterization use the method // | |
41 | // TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0); // | |
42 | // it will compare the user output of the user function and interpolation // | |
43 | // for idim-th output dimension and fill the difference in the supplied // | |
44 | // histogram. If no histogram is supplied, it will be created. // | |
45 | // // | |
46 | // To save the interpolation data: // | |
47 | // SaveData(const char* filename, Bool_t append ) // | |
48 | // write text file with data. If append is kTRUE and the output file already // | |
49 | // exists, data will be added in the end of the file. // | |
50 | // Alternatively, SaveData(FILE* stream) will write the data to // | |
51 | // already existing stream. // | |
52 | // // | |
53 | // To read back already stored interpolation use either the constructor // | |
54 | // AliCheb3D(const char* inpFile); // | |
55 | // or the default constructor AliCheb3D() followed by // | |
56 | // AliCheb3D::LoadData(const char* inpFile); // | |
57 | // // | |
58 | // To compute the interpolation use Eval(float* par,float *res) method, with // | |
59 | // par being 3D vector of arguments (inside the validity region) and res is // | |
60 | // the array of DimOut elements for the output. // | |
61 | // // | |
62 | // If only one component (say, idim-th) of the output is needed, use faster // | |
63 | // Float_t Eval(Float_t *par,int idim) method. // | |
64 | // // | |
65 | // void Print(option="") will print the name, the ranges of validity and // | |
66 | // the absolute precision of the parameterization. Option "l" will also print // | |
67 | // the information about the number of coefficients for each output // | |
68 | // dimension. // | |
69 | // // | |
70 | // NOTE: during the evaluation no check is done for parameter vector being // | |
71 | // outside the interpolation region. If there is such a risk, use // | |
72 | // Bool_t IsInside(float *par) method. Chebyshev parameterization is not // | |
73 | // good for extrapolation! // | |
74 | // // | |
75 | // For the properties of Chebyshev parameterization see: // | |
76 | // H.Wind, CERN EP Internal Report, 81-12/Rev. // | |
77 | // // | |
78 | //////////////////////////////////////////////////////////////////////////////// | |
79 | ||
80 | #include <TString.h> | |
81 | #include <TSystem.h> | |
82 | #include <TRandom.h> | |
83 | #include <TROOT.h> | |
84 | #include "AliCheb3D.h" | |
0bc7b414 | 85 | #include "AliLog.h" |
0eea9d4d | 86 | |
87 | ||
88 | ||
89 | ClassImp(AliCheb3DCalc) | |
90 | ||
91 | ||
92 | AliCheb3DCalc::AliCheb3DCalc(): | |
93 | TNamed("", ""), | |
94 | fNCoefs(0), | |
95 | fNRows(0), | |
96 | fNCols(0), | |
97 | fNElemBound2D(0), | |
98 | fNColsAtRow(0), | |
99 | fColAtRowBg(0), | |
100 | fCoefBound2D0(0), | |
101 | fCoefBound2D1(0), | |
102 | fCoefs(0), | |
103 | fTmpCf1(0), | |
104 | fTmpCf0(0) | |
105 | { | |
106 | // Default constructor | |
107 | Init0(); | |
108 | } | |
109 | ||
110 | AliCheb3DCalc::AliCheb3DCalc(FILE* stream): | |
111 | TNamed("", ""), | |
112 | fNCoefs(0), | |
113 | fNRows(0), | |
114 | fNCols(0), | |
115 | fNElemBound2D(0), | |
116 | fNColsAtRow(0), | |
117 | fColAtRowBg(0), | |
118 | fCoefBound2D0(0), | |
119 | fCoefBound2D1(0), | |
120 | fCoefs(0), | |
121 | fTmpCf1(0), | |
122 | fTmpCf0(0) | |
123 | { | |
124 | // Default constructor | |
125 | Init0(); | |
126 | LoadData(stream); | |
127 | } | |
128 | ||
c437b1a5 | 129 | |
130 | AliCheb3DCalc::AliCheb3DCalc(const AliCheb3DCalc& src) : | |
131 | TNamed(src), | |
132 | fNCoefs(src.fNCoefs), | |
133 | fNRows(src.fNRows), | |
134 | fNCols(src.fNCols), | |
135 | fNElemBound2D(src.fNElemBound2D), | |
136 | fNColsAtRow(0), | |
137 | fColAtRowBg(0), | |
138 | fCoefBound2D0(0), | |
139 | fCoefBound2D1(0), | |
140 | fCoefs(0), | |
141 | fTmpCf1(0), | |
0bc7b414 | 142 | fTmpCf0(0) |
143 | { | |
144 | // Copy constructor | |
c437b1a5 | 145 | if (src.fNColsAtRow) { |
146 | fNColsAtRow = new Int_t[fNRows]; | |
147 | for (int i=fNRows;i--;) fNColsAtRow[i] = src.fNColsAtRow[i]; | |
148 | } | |
149 | if (src.fColAtRowBg) { | |
150 | fColAtRowBg = new Int_t[fNRows]; | |
151 | for (int i=fNRows;i--;) fColAtRowBg[i] = src.fColAtRowBg[i]; | |
152 | } | |
153 | if (src.fCoefBound2D0) { | |
154 | fCoefBound2D0 = new Int_t[fNElemBound2D]; | |
155 | for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = src.fCoefBound2D0[i]; | |
156 | } | |
157 | if (src.fCoefBound2D1) { | |
158 | fCoefBound2D1 = new Int_t[fNElemBound2D]; | |
159 | for (int i=fNElemBound2D;i--;) fCoefBound2D1[i] = src.fCoefBound2D1[i]; | |
160 | } | |
161 | if (src.fCoefs) { | |
162 | fCoefs = new Float_t[fNCoefs]; | |
163 | for (int i=fNCoefs;i--;) fCoefs[i] = src.fCoefs[i]; | |
164 | } | |
165 | if (src.fTmpCf1) fTmpCf1 = new Float_t[fNCols]; | |
166 | if (src.fTmpCf0) fTmpCf0 = new Float_t[fNRows]; | |
0bc7b414 | 167 | } |
168 | ||
c437b1a5 | 169 | AliCheb3DCalc& AliCheb3DCalc::operator=(const AliCheb3DCalc& rhs) |
0bc7b414 | 170 | { |
c437b1a5 | 171 | // Assignment operator |
172 | if (this != &rhs) { | |
173 | Clear(); | |
174 | SetName(rhs.GetName()); | |
175 | SetTitle(rhs.GetTitle()); | |
176 | fNCoefs = rhs.fNCoefs; | |
177 | fNRows = rhs.fNRows; | |
178 | fNCols = rhs.fNCols; | |
179 | if (rhs.fNColsAtRow) { | |
180 | fNColsAtRow = new Int_t[fNRows]; | |
181 | for (int i=fNRows;i--;) fNColsAtRow[i] = rhs.fNColsAtRow[i]; | |
182 | } | |
183 | if (rhs.fColAtRowBg) { | |
184 | fColAtRowBg = new Int_t[fNRows]; | |
185 | for (int i=fNRows;i--;) fColAtRowBg[i] = rhs.fColAtRowBg[i]; | |
186 | } | |
187 | if (rhs.fCoefBound2D0) { | |
188 | fCoefBound2D0 = new Int_t[fNElemBound2D]; | |
189 | for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = rhs.fCoefBound2D0[i]; | |
190 | } | |
191 | if (rhs.fCoefBound2D1) { | |
192 | fCoefBound2D1 = new Int_t[fNElemBound2D]; | |
193 | for (int i=fNElemBound2D;i--;) fCoefBound2D1[i] = rhs.fCoefBound2D1[i]; | |
194 | } | |
195 | if (rhs.fCoefs) { | |
196 | fCoefs = new Float_t[fNCoefs]; | |
197 | for (int i=fNCoefs;i--;) fCoefs[i] = rhs.fCoefs[i]; | |
198 | } | |
199 | if (rhs.fTmpCf1) fTmpCf1 = new Float_t[fNCols]; | |
200 | if (rhs.fTmpCf0) fTmpCf0 = new Float_t[fNRows]; | |
201 | } | |
202 | return *this; | |
0bc7b414 | 203 | } |
204 | ||
0eea9d4d | 205 | //__________________________________________________________________________________________ |
206 | void AliCheb3DCalc::Clear(Option_t*) | |
207 | { | |
208 | // delete all dynamycally allocated structures | |
209 | if (fTmpCf1) { delete[] fTmpCf1; fTmpCf1 = 0;} | |
210 | if (fTmpCf0) { delete[] fTmpCf0; fTmpCf0 = 0;} | |
211 | if (fCoefs) { delete[] fCoefs; fCoefs = 0;} | |
212 | if (fCoefBound2D0) { delete[] fCoefBound2D0; fCoefBound2D0 = 0; } | |
213 | if (fCoefBound2D1) { delete[] fCoefBound2D1; fCoefBound2D1 = 0; } | |
214 | if (fNColsAtRow) { delete[] fNColsAtRow; fNColsAtRow = 0; } | |
215 | if (fColAtRowBg) { delete[] fColAtRowBg; fColAtRowBg = 0; } | |
216 | // | |
217 | } | |
218 | ||
219 | //__________________________________________________________________________________________ | |
220 | void AliCheb3DCalc::Init0() | |
221 | { | |
222 | // reset everything to 0 | |
223 | fNCoefs = fNRows = fNCols = fNElemBound2D = 0; | |
224 | fCoefs = 0; | |
225 | fCoefBound2D0 = fCoefBound2D1 = 0; | |
226 | fNColsAtRow = fColAtRowBg = 0; | |
227 | fTmpCf0 = fTmpCf1 = 0; | |
228 | } | |
229 | ||
230 | //__________________________________________________________________________________________ | |
231 | void AliCheb3DCalc::Print(Option_t* ) const | |
232 | { | |
233 | printf("Chebyshev parameterization data %s for 3D->1 function.\n",GetName()); | |
234 | int nmax3d = 0; | |
235 | for (int i=fNElemBound2D;i--;) if (fCoefBound2D0[i]>nmax3d) nmax3d = fCoefBound2D0[i]; | |
236 | printf("%d coefficients in %dx%dx%d matrix\n",fNCoefs,fNRows,fNCols,nmax3d); | |
237 | // | |
238 | } | |
239 | ||
240 | //__________________________________________________________________________________________ | |
241 | Float_t AliCheb3DCalc::Eval(Float_t *par) const | |
242 | { | |
243 | // evaluate Chebyshev parameterization for 3D function. | |
244 | // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval | |
245 | Float_t &z = par[2]; | |
246 | Float_t &y = par[1]; | |
247 | Float_t &x = par[0]; | |
248 | // | |
249 | int ncfRC; | |
250 | for (int id0=fNRows;id0--;) { | |
251 | int nCLoc = fNColsAtRow[id0]; // number of significant coefs on this row | |
252 | int Col0 = fColAtRowBg[id0]; // beginning of local column in the 2D boundary matrix | |
253 | for (int id1=nCLoc;id1--;) { | |
254 | int id = id1+Col0; | |
255 | fTmpCf1[id1] = (ncfRC=fCoefBound2D0[id]) ? ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC) : 0.0; | |
256 | } | |
257 | fTmpCf0[id0] = nCLoc>0 ? ChebEval1D(y,fTmpCf1,nCLoc):0.0; | |
258 | } | |
259 | return ChebEval1D(x,fTmpCf0,fNRows); | |
260 | // | |
261 | } | |
262 | ||
263 | //_______________________________________________ | |
264 | #ifdef _INC_CREATION_ALICHEB3D_ | |
265 | void AliCheb3DCalc::SaveData(const char* outfile,Bool_t append) const | |
266 | { | |
267 | // writes coefficients data to output text file, optionallt appending on the end of existing file | |
268 | TString strf = outfile; | |
269 | gSystem->ExpandPathName(strf); | |
270 | FILE* stream = fopen(strf,append ? "a":"w"); | |
271 | SaveData(stream); | |
272 | fclose(stream); | |
273 | // | |
274 | } | |
275 | #endif | |
276 | ||
277 | //_______________________________________________ | |
278 | #ifdef _INC_CREATION_ALICHEB3D_ | |
279 | void AliCheb3DCalc::SaveData(FILE* stream) const | |
280 | { | |
281 | // writes coefficients data to existing output stream | |
282 | // Note: fNCols, fNElemBound2D and fColAtRowBg is not stored, will be computed on fly during the loading of this file | |
283 | fprintf(stream,"#\nSTART %s\n",GetName()); | |
284 | fprintf(stream,"# Number of rows\n%d\n",fNRows); | |
285 | // | |
286 | fprintf(stream,"# Number of columns per row\n"); | |
287 | for (int i=0;i<fNRows;i++) fprintf(stream,"%d\n",fNColsAtRow[i]); | |
288 | // | |
289 | fprintf(stream,"# Number of Coefs in each significant block of third dimension\n"); | |
290 | for (int i=0;i<fNElemBound2D;i++) fprintf(stream,"%d\n",fCoefBound2D0[i]); | |
291 | // | |
292 | fprintf(stream,"# Coefficients\n"); | |
293 | for (int i=0;i<fNCoefs;i++) fprintf(stream,"%+.8e\n",fCoefs[i]); | |
294 | fprintf(stream,"END %s\n",GetName()); | |
295 | // | |
296 | } | |
297 | #endif | |
298 | ||
299 | //_______________________________________________ | |
300 | void AliCheb3DCalc::LoadData(FILE* stream) | |
301 | { | |
302 | // Load coefs. from the stream | |
303 | if (!stream) {Error("LoadData","No stream provided.\nStop"); exit(1);} | |
304 | TString buffs; | |
305 | Clear(); | |
306 | ReadLine(buffs,stream); | |
307 | if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <fit_name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);} | |
308 | if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1); | |
309 | // | |
310 | ReadLine(buffs,stream); // NRows | |
311 | fNRows = buffs.Atoi(); | |
312 | if (fNRows<1) {Error("LoadData","Expected: '<number_of_rows>', found \"%s\"\nStop\n",buffs.Data());exit(1);} | |
313 | // | |
314 | fNCols = 0; | |
315 | fNElemBound2D = 0; | |
316 | InitRows(fNRows); | |
317 | // | |
318 | for (int id0=0;id0<fNRows;id0++) { | |
319 | ReadLine(buffs,stream); // n.cols at this row | |
320 | fNColsAtRow[id0] = buffs.Atoi(); | |
321 | fColAtRowBg[id0] = fNElemBound2D; // begining of this row in 2D boundary surface | |
322 | fNElemBound2D += fNColsAtRow[id0]; | |
323 | if (fNCols<fNColsAtRow[id0]) fNCols = fNColsAtRow[id0]; | |
324 | } | |
325 | InitCols(fNCols); | |
326 | // | |
327 | fNCoefs = 0; | |
328 | InitElemBound2D(fNElemBound2D); | |
329 | // | |
330 | for (int i=0;i<fNElemBound2D;i++) { | |
331 | ReadLine(buffs,stream); // n.coeffs at 3-d dimension for the given column/row | |
332 | fCoefBound2D0[i] = buffs.Atoi(); | |
333 | fCoefBound2D1[i] = fNCoefs; | |
334 | fNCoefs += fCoefBound2D0[i]; | |
335 | } | |
336 | // | |
337 | if (fNCoefs<=0) {Error("LoadData","Negtive (%d) number of Chebychef coeffs. is obtained.\nStop\n",fNCoefs);exit(1);} | |
338 | // | |
339 | InitCoefs(fNCoefs); | |
340 | for (int i=0;i<fNCoefs;i++) { | |
341 | ReadLine(buffs,stream); | |
342 | fCoefs[i] = buffs.Atof(); | |
343 | } | |
344 | // check end_of_data record | |
345 | ReadLine(buffs,stream); | |
346 | if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) { | |
347 | Error("LoadData","Expected \"END %s\", found \"%s\".\nStop\n",GetName(),buffs.Data()); | |
348 | exit(1); | |
349 | } | |
350 | // | |
351 | } | |
352 | ||
353 | //_______________________________________________ | |
354 | void AliCheb3DCalc::ReadLine(TString& str,FILE* stream) | |
355 | { | |
356 | // read single line from the stream, skipping empty and commented lines. EOF is not expected | |
357 | while (str.Gets(stream)) { | |
358 | str = str.Strip(TString::kBoth,' '); | |
359 | if (str.IsNull()||str.BeginsWith("#")) continue; | |
360 | return; | |
361 | } | |
362 | fprintf(stderr,"AliCheb3D::ReadLine: Failed to read from stream.\nStop");exit(1); // normally, should not reach here | |
363 | } | |
364 | ||
365 | //_______________________________________________ | |
366 | void AliCheb3DCalc::InitCols(int nc) | |
367 | { | |
368 | // Set max.number of significant columns in the coefs matrix | |
369 | fNCols = nc; | |
370 | if (fTmpCf1) delete[] fTmpCf1; | |
371 | fTmpCf1 = new Float_t [fNCols]; | |
372 | } | |
373 | ||
374 | //_______________________________________________ | |
375 | void AliCheb3DCalc::InitRows(int nr) | |
376 | { | |
377 | // Set max.number of significant rows in the coefs matrix | |
378 | if (fNColsAtRow) delete[] fNColsAtRow; | |
379 | if (fColAtRowBg) delete[] fColAtRowBg; | |
380 | if (fTmpCf0) delete[] fTmpCf0; | |
381 | fNRows = nr; | |
382 | fNColsAtRow = new Int_t[fNRows]; | |
383 | fTmpCf0 = new Float_t [fNRows]; | |
384 | fColAtRowBg = new Int_t[fNRows]; | |
385 | for (int i=fNRows;i--;) fNColsAtRow[i] = fColAtRowBg[i] = 0; | |
386 | } | |
387 | ||
388 | //_______________________________________________ | |
389 | void AliCheb3DCalc::InitElemBound2D(int ne) | |
390 | { | |
391 | // Set max number of significant coefs for given row/column of coefs 3D matrix | |
392 | if (fCoefBound2D0) delete[] fCoefBound2D0; | |
393 | if (fCoefBound2D1) delete[] fCoefBound2D1; | |
394 | fNElemBound2D = ne; | |
395 | fCoefBound2D0 = new Int_t[fNElemBound2D]; | |
396 | fCoefBound2D1 = new Int_t[fNElemBound2D]; | |
397 | for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = fCoefBound2D1[i] = 0; | |
398 | } | |
399 | ||
400 | //_______________________________________________ | |
401 | void AliCheb3DCalc::InitCoefs(int nc) | |
402 | { | |
403 | // Set total number of significant coefs | |
404 | if (fCoefs) delete[] fCoefs; | |
405 | fNCoefs = nc; | |
406 | fCoefs = new Float_t [fNCoefs]; | |
407 | for (int i=fNCoefs;i--;) fCoefs[i] = 0.0; | |
408 | } | |
409 | ||
410 | ||
411 | ||
412 | ||
413 | ClassImp(AliCheb3D) | |
414 | ||
415 | AliCheb3D::AliCheb3D(): | |
416 | TNamed("", ""), | |
417 | fDimOut(0), | |
418 | fPrec(0.), | |
0bc7b414 | 419 | fChebCalc(), |
0eea9d4d | 420 | fMaxCoefs(0), |
421 | fResTmp(0), | |
422 | fGrid(0), | |
0bc7b414 | 423 | fUsrFunName(), |
0eea9d4d | 424 | fUsrMacro(0) |
425 | { | |
426 | // Default constructor | |
427 | Init0(); | |
428 | } | |
429 | ||
430 | AliCheb3D::AliCheb3D(const char* inputFile): | |
431 | TNamed("", ""), | |
432 | fDimOut(0), | |
433 | fPrec(0.), | |
0bc7b414 | 434 | fChebCalc(), |
0eea9d4d | 435 | fMaxCoefs(0), |
436 | fResTmp(0), | |
437 | fGrid(0), | |
0bc7b414 | 438 | fUsrFunName(), |
0eea9d4d | 439 | fUsrMacro(0) |
440 | { | |
441 | // Default constructor | |
442 | Init0(); | |
443 | LoadData(inputFile); | |
444 | } | |
445 | ||
446 | ||
447 | ||
448 | AliCheb3D::AliCheb3D(FILE* stream): | |
449 | TNamed("", ""), | |
450 | fDimOut(0), | |
451 | fPrec(0.), | |
0bc7b414 | 452 | fChebCalc(), |
0eea9d4d | 453 | fMaxCoefs(0), |
454 | fResTmp(0), | |
455 | fGrid(0), | |
0bc7b414 | 456 | fUsrFunName(), |
0eea9d4d | 457 | fUsrMacro(0) |
458 | { | |
459 | // Default constructor | |
460 | Init0(); | |
461 | LoadData(stream); | |
462 | } | |
463 | ||
c437b1a5 | 464 | AliCheb3D::AliCheb3D(const AliCheb3D& src) : |
465 | TNamed(src), | |
466 | fDimOut(src.fDimOut), | |
467 | fPrec(src.fPrec), | |
468 | fChebCalc(1), | |
469 | fMaxCoefs(src.fMaxCoefs), | |
470 | fResTmp(0), | |
471 | fGrid(0), | |
472 | fUsrFunName(src.fUsrFunName), | |
473 | fUsrMacro(0) | |
0bc7b414 | 474 | { |
475 | // Copy constructor | |
c437b1a5 | 476 | // read coefs from text file |
477 | for (int i=3;i--;) { | |
478 | fBMin[i] = src.fBMin[i]; | |
479 | fBMax[i] = src.fBMax[i]; | |
480 | fBScale[i] = src.fBScale[i]; | |
481 | fBOffset[i] = src.fBOffset[i]; | |
482 | fNPoints[i] = src.fNPoints[i]; | |
483 | } | |
484 | for (int i=0;i<fDimOut;i++) { | |
485 | AliCheb3DCalc* cbc = src.GetChebCalc(i); | |
486 | if (cbc) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(*cbc),i); | |
487 | } | |
488 | } | |
489 | ||
490 | AliCheb3D& AliCheb3D::operator=(const AliCheb3D& rhs) | |
491 | { | |
492 | // Assignment operator | |
493 | if (this != &rhs) { | |
494 | Clear(); | |
495 | fDimOut = rhs.fDimOut; | |
496 | fPrec = rhs.fPrec; | |
497 | fMaxCoefs = rhs.fMaxCoefs; | |
498 | fUsrFunName = rhs.fUsrFunName; | |
499 | fUsrMacro = 0; | |
500 | for (int i=3;i--;) { | |
501 | fBMin[i] = rhs.fBMin[i]; | |
502 | fBMax[i] = rhs.fBMax[i]; | |
503 | fBScale[i] = rhs.fBScale[i]; | |
504 | fBOffset[i] = rhs.fBOffset[i]; | |
505 | fNPoints[i] = rhs.fNPoints[i]; | |
506 | } | |
507 | for (int i=0;i<fDimOut;i++) { | |
508 | AliCheb3DCalc* cbc = rhs.GetChebCalc(i); | |
509 | if (cbc) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(*cbc),i); | |
510 | } | |
511 | } | |
512 | return *this; | |
513 | // | |
0bc7b414 | 514 | } |
0eea9d4d | 515 | |
c437b1a5 | 516 | |
0eea9d4d | 517 | //__________________________________________________________________________________________ |
518 | #ifdef _INC_CREATION_ALICHEB3D_ | |
519 | AliCheb3D::AliCheb3D(const char* funName, int DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec) : TNamed(funName,funName) | |
520 | { | |
521 | // Construct the parameterization for the function | |
522 | // funName : name of the file containing the function: void funName(Float_t * inp,Float_t * out) | |
523 | // DimOut : dimension of the vector computed by the user function | |
524 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
525 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
526 | // npoints : array of 3 elements with the number of points to compute in each of 3 dimension | |
527 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
528 | // | |
529 | Init0(); | |
530 | fPrec = TMath::Max(1.E-12f,prec); | |
531 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} | |
532 | SetDimOut(DimOut); | |
533 | PrepareBoundaries(bmin,bmax); | |
534 | DefineGrid(npoints); | |
535 | SetUsrFunction(funName); | |
536 | ChebFit(); | |
537 | // | |
538 | } | |
539 | #endif | |
540 | ||
541 | //__________________________________________________________________________________________ | |
542 | #ifdef _INC_CREATION_ALICHEB3D_ | |
543 | AliCheb3D::AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec) : TNamed("AliCheb3D","AliCheb3D") | |
544 | { | |
545 | // Construct the parameterization for the function | |
546 | // ptr : pointer on the function: void fun(Float_t * inp,Float_t * out) | |
547 | // DimOut : dimension of the vector computed by the user function | |
548 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
549 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
550 | // npoints : array of 3 elements with the number of points to compute in each of 3 dimension | |
551 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
552 | // | |
553 | Init0(); | |
554 | fPrec = TMath::Max(1.E-12f,prec); | |
555 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} | |
556 | SetDimOut(DimOut); | |
557 | PrepareBoundaries(bmin,bmax); | |
558 | DefineGrid(npoints); | |
559 | SetUsrFunction(ptr); | |
560 | ChebFit(); | |
561 | // | |
562 | } | |
563 | #endif | |
564 | ||
0bc7b414 | 565 | |
0eea9d4d | 566 | //__________________________________________________________________________________________ |
567 | void AliCheb3D::Clear(Option_t*) | |
568 | { | |
569 | if (fResTmp) { delete[] fResTmp; fResTmp = 0; } | |
570 | if (fGrid) { delete[] fGrid; fGrid = 0; } | |
571 | if (fUsrMacro) { delete fUsrMacro; fUsrMacro = 0;} | |
572 | fChebCalc.Delete(); | |
573 | // | |
574 | } | |
575 | ||
576 | //__________________________________________________________________________________________ | |
577 | void AliCheb3D::Print(Option_t* opt) const | |
578 | { | |
579 | printf("%s: Chebyshev parameterization for 3D->%dD function. Precision: %e\n",GetName(),fDimOut,fPrec); | |
580 | printf("Region of validity: [%+.5e:%+.5e] [%+.5e:%+.5e] [%+.5e:%+.5e]\n",fBMin[0],fBMax[0],fBMin[1],fBMax[1],fBMin[2],fBMax[2]); | |
581 | TString opts = opt; opts.ToLower(); | |
582 | if (opts.Contains("l")) for (int i=0;i<fDimOut;i++) {printf("Output dimension %d:\n",i+1); GetChebCalc(i)->Print();} | |
583 | // | |
584 | } | |
585 | ||
586 | //__________________________________________________________________________________________ | |
587 | void AliCheb3D::Init0() | |
588 | { | |
589 | for (int i=3;i--;) fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
590 | fMaxCoefs = 0; | |
591 | fGrid = 0; | |
592 | fResTmp = 0; | |
593 | fUsrFunName = ""; | |
594 | fUsrMacro = 0; | |
595 | #ifdef _INC_CREATION_ALICHEB3D_ | |
596 | gUsrFunAliCheb3D = 0; | |
597 | #endif | |
598 | } | |
599 | ||
600 | //__________________________________________________________________________________________ | |
601 | void AliCheb3D::PrepareBoundaries(Float_t *bmin,Float_t *bmax) | |
602 | { | |
603 | // Set and check boundaries defined by user, prepare coefficients for their conversion to [-1:1] interval | |
604 | // | |
605 | for (int i=3;i--;) { | |
606 | fBMin[i] = bmin[i]; | |
607 | fBMax[i] = bmax[i]; | |
608 | fBScale[i] = bmax[i]-bmin[i]; | |
609 | if (fBScale[i]<=0) { | |
610 | Error("PrepareBoundaries","Boundaries for %d-th dimension are not increasing: %+.4e %+.4e\nStop\n",i,fBMin[i],fBMax[i]); | |
611 | exit(1); | |
612 | } | |
613 | fBOffset[i] = bmin[i] + fBScale[i]/2.0; | |
614 | fBScale[i] = 2./fBScale[i]; | |
615 | } | |
616 | // | |
617 | } | |
618 | ||
619 | //__________________________________________________________________________________________ | |
620 | #ifdef _INC_CREATION_ALICHEB3D_ | |
621 | void AliCheb3D::SetUsrFunction(const char* name) | |
622 | { | |
623 | // load user macro with function definition and compile it | |
624 | gUsrFunAliCheb3D = 0; | |
625 | fUsrFunName = name; | |
626 | gSystem->ExpandPathName(fUsrFunName); | |
627 | if (fUsrMacro) delete fUsrMacro; | |
628 | TString tmpst = fUsrFunName; | |
629 | tmpst += "+"; // prepare filename to compile | |
630 | if (gROOT->LoadMacro(tmpst.Data())) {Error("SetUsrFunction","Failed to load user function from %s\nStop\n",name); exit(1);} | |
631 | fUsrMacro = new TMethodCall(); | |
632 | tmpst = tmpst.Data() + tmpst.Last('/')+1; //Strip away any path preceding the macro file name | |
633 | int dot = tmpst.Last('.'); | |
634 | if (dot>0) tmpst.Resize(dot); | |
635 | fUsrMacro->InitWithPrototype(tmpst.Data(),"Float_t *,Float_t *"); | |
636 | long args[2]; | |
637 | args[0] = (long)fArgsTmp; | |
638 | args[1] = (long)fResTmp; | |
639 | fUsrMacro->SetParamPtrs(args); | |
640 | // | |
641 | } | |
642 | #endif | |
643 | ||
644 | //__________________________________________________________________________________________ | |
645 | #ifdef _INC_CREATION_ALICHEB3D_ | |
646 | void AliCheb3D::SetUsrFunction(void (*ptr)(float*,float*)) | |
647 | { | |
648 | if (fUsrMacro) delete fUsrMacro; | |
649 | fUsrMacro = 0; | |
650 | fUsrFunName = ""; | |
651 | gUsrFunAliCheb3D = ptr; | |
652 | } | |
653 | #endif | |
654 | ||
655 | //__________________________________________________________________________________________ | |
656 | #ifdef _INC_CREATION_ALICHEB3D_ | |
657 | void AliCheb3D::EvalUsrFunction(Float_t *x, Float_t *res) { | |
658 | for (int i=3;i--;) fArgsTmp[i] = x[i]; | |
659 | if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); | |
660 | else fUsrMacro->Execute(); | |
661 | for (int i=fDimOut;i--;) res[i] = fResTmp[i]; | |
662 | } | |
663 | #endif | |
664 | ||
665 | //__________________________________________________________________________________________ | |
666 | #ifdef _INC_CREATION_ALICHEB3D_ | |
667 | Int_t AliCheb3D::CalcChebCoefs(Float_t *funval,int np, Float_t *outCoefs, Float_t prec) | |
668 | { | |
669 | // Calculate Chebyshev coeffs using precomputed function values at np roots. | |
670 | // If prec>0, estimate the highest coeff number providing the needed precision | |
671 | // | |
672 | double sm; // do summations in double to minimize the roundoff error | |
673 | for (int ic=0;ic<np;ic++) { // compute coeffs | |
674 | sm = 0; | |
675 | for (int ir=0;ir<np;ir++) { | |
676 | float rt = TMath::Cos( ic*(ir+0.5)*TMath::Pi()/np); | |
677 | sm += funval[ir]*rt; | |
678 | } | |
679 | outCoefs[ic] = Float_t( sm * ((ic==0) ? 1./np : 2./np) ); | |
680 | } | |
681 | // | |
682 | if (prec<=0) return np; | |
683 | // | |
684 | sm = 0; | |
685 | int cfMax = 0; | |
686 | for (cfMax=np;cfMax--;) { | |
687 | sm += TMath::Abs(outCoefs[cfMax]); | |
688 | if (sm>=prec) break; | |
689 | } | |
690 | if (++cfMax==0) cfMax=1; | |
691 | return cfMax; | |
692 | // | |
693 | } | |
694 | #endif | |
695 | ||
696 | //__________________________________________________________________________________________ | |
697 | #ifdef _INC_CREATION_ALICHEB3D_ | |
698 | void AliCheb3D::DefineGrid(Int_t* npoints) | |
699 | { | |
700 | // prepare the grid of Chebyshev roots in each dimension | |
701 | const int kMinPoints = 1; | |
702 | int ntot = 0; | |
703 | fMaxCoefs = 1; | |
704 | for (int id=3;id--;) { | |
705 | fNPoints[id] = npoints[id]; | |
706 | if (fNPoints[id]<kMinPoints) { | |
707 | Error("DefineGrid","at %d-th dimension %d point is requested, at least %d is needed\nStop\n",fNPoints[id],kMinPoints); | |
708 | exit(1); | |
709 | } | |
710 | ntot += fNPoints[id]; | |
711 | fMaxCoefs *= fNPoints[id]; | |
712 | } | |
713 | fGrid = new Float_t [ntot]; | |
714 | // | |
715 | int curp = 0; | |
716 | for (int id=3;id--;) { | |
717 | int np = fNPoints[id]; | |
718 | fGridOffs[id] = curp; | |
719 | for (int ip=0;ip<np;ip++) { | |
720 | Float_t x = TMath::Cos( TMath::Pi()*(ip+0.5)/np ); | |
721 | fGrid[curp++] = MapToExternal(x,id); | |
722 | } | |
723 | } | |
724 | // | |
725 | } | |
726 | #endif | |
727 | ||
728 | //__________________________________________________________________________________________ | |
729 | #ifdef _INC_CREATION_ALICHEB3D_ | |
730 | Int_t AliCheb3D::ChebFit() | |
731 | { | |
732 | // prepare parameterization for all output dimensions | |
733 | int ir=0; | |
734 | for (int i=fDimOut;i--;) ir+=ChebFit(i); | |
735 | return ir; | |
736 | } | |
737 | #endif | |
738 | ||
739 | //__________________________________________________________________________________________ | |
740 | #ifdef _INC_CREATION_ALICHEB3D_ | |
741 | Int_t AliCheb3D::ChebFit(int dmOut) | |
742 | { | |
743 | // prepare paramaterization of 3D function for dmOut-th dimension | |
744 | int maxDim = 0; | |
745 | for (int i=0;i<3;i++) if (maxDim<fNPoints[i]) maxDim = fNPoints[i]; | |
746 | Float_t *fvals = new Float_t [ fNPoints[0] ]; | |
747 | Float_t *tmpCoef3D = new Float_t [ fNPoints[0]*fNPoints[1]*fNPoints[2] ]; | |
748 | Float_t *tmpCoef2D = new Float_t [ fNPoints[0]*fNPoints[1] ]; | |
749 | Float_t *tmpCoef1D = new Float_t [ maxDim ]; | |
750 | // | |
751 | Float_t RTiny = fPrec/Float_t(maxDim); // neglect coefficient below this threshold | |
752 | // | |
753 | // 1D Cheb.fit for 0-th dimension at current steps of remaining dimensions | |
754 | int ncmax = 0; | |
755 | // | |
756 | AliCheb3DCalc* cheb = GetChebCalc(dmOut); | |
757 | // | |
758 | for (int id2=fNPoints[2];id2--;) { | |
759 | fArgsTmp[2] = fGrid[ fGridOffs[2]+id2 ]; | |
760 | // | |
761 | for (int id1=fNPoints[1];id1--;) { | |
762 | fArgsTmp[1] = fGrid[ fGridOffs[1]+id1 ]; | |
763 | // | |
764 | for (int id0=fNPoints[0];id0--;) { | |
765 | fArgsTmp[0] = fGrid[ fGridOffs[0]+id0 ]; | |
766 | EvalUsrFunction(); // compute function values at Chebyshev roots of 0-th dimension | |
767 | fvals[id0] = fResTmp[dmOut]; | |
768 | } | |
769 | int nc = CalcChebCoefs(fvals,fNPoints[0], tmpCoef1D, fPrec); | |
770 | for (int id0=fNPoints[0];id0--;) tmpCoef2D[id1 + id0*fNPoints[1]] = tmpCoef1D[id0]; | |
771 | if (ncmax<nc) ncmax = nc; // max coefs to be kept in dim0 to guarantee needed precision | |
772 | } | |
773 | // | |
774 | // once each 1d slice of given 2d slice is parametrized, parametrize the Cheb.coeffs | |
775 | for (int id0=fNPoints[0];id0--;) { | |
776 | CalcChebCoefs( tmpCoef2D+id0*fNPoints[1], fNPoints[1], tmpCoef1D, -1); | |
777 | for (int id1=fNPoints[1];id1--;) tmpCoef3D[id2 + fNPoints[2]*(id1+id0*fNPoints[1])] = tmpCoef1D[id1]; | |
778 | } | |
779 | } | |
780 | // | |
781 | // now fit the last dimensions Cheb.coefs | |
782 | for (int id0=fNPoints[0];id0--;) { | |
783 | for (int id1=fNPoints[1];id1--;) { | |
784 | CalcChebCoefs( tmpCoef3D+ fNPoints[2]*(id1+id0*fNPoints[1]), fNPoints[2], tmpCoef1D, -1); | |
785 | for (int id2=fNPoints[2];id2--;) tmpCoef3D[id2+ fNPoints[2]*(id1+id0*fNPoints[1])] = tmpCoef1D[id2]; // store on place | |
786 | } | |
787 | } | |
788 | // | |
789 | // now find 2D surface which separates significant coefficients of 3D matrix from nonsignificant ones (up to fPrec) | |
790 | int *tmpCoefSurf = new Int_t[ fNPoints[0]*fNPoints[1] ]; | |
791 | for (int id0=fNPoints[0];id0--;) for (int id1=fNPoints[1];id1--;) tmpCoefSurf[id1+id0*fNPoints[1]]=0; | |
792 | Double_t resid = 0; | |
793 | for (int id0=fNPoints[0];id0--;) { | |
794 | for (int id1=fNPoints[1];id1--;) { | |
795 | for (int id2=fNPoints[2];id2--;) { | |
796 | int id = id2 + fNPoints[2]*(id1+id0*fNPoints[1]); | |
797 | Float_t cfa = TMath::Abs(tmpCoef3D[id]); | |
798 | if (cfa < RTiny) {tmpCoef3D[id] = 0; continue;} // neglect coeefs below the threshold | |
799 | ||
800 | resid += cfa; | |
801 | if (resid<fPrec) continue; // this coeff is negligible | |
802 | // otherwise go back 1 step | |
803 | resid -= cfa; | |
804 | tmpCoefSurf[id1+id0*fNPoints[1]] = id2+1; // how many coefs to keep | |
805 | break; | |
806 | } | |
807 | } | |
808 | } | |
809 | /* | |
810 | printf("\n\nCoeffs\n"); | |
811 | int cnt = 0; | |
812 | for (int id0=0;id0<fNPoints[0];id0++) { | |
813 | for (int id1=0;id1<fNPoints[1];id1++) { | |
814 | for (int id2=0;id2<fNPoints[2];id2++) { | |
815 | printf("%2d%2d%2d %+.4e |",id0,id1,id2,tmpCoef3D[cnt++]); | |
816 | } | |
817 | printf("\n"); | |
818 | } | |
819 | printf("\n"); | |
820 | } | |
821 | */ | |
822 | // see if there are rows to reject, find max.significant column at each row | |
823 | int NRows = fNPoints[0]; | |
824 | int *tmpCols = new int[NRows]; | |
825 | for (int id0=fNPoints[0];id0--;) { | |
826 | int id1 = fNPoints[1]; | |
827 | while (id1>0 && tmpCoefSurf[(id1-1)+id0*fNPoints[1]]==0) id1--; | |
828 | tmpCols[id0] = id1; | |
829 | } | |
830 | // find max significant row | |
831 | for (int id0=NRows;id0--;) {if (tmpCols[id0]>0) break; NRows--;} | |
832 | // find max significant column and fill the permanent storage for the max sigificant column of each row | |
833 | cheb->InitRows(NRows); // create needed arrays; | |
834 | int *NColsAtRow = cheb->GetNColsAtRow(); | |
835 | int *ColAtRowBg = cheb->GetColAtRowBg(); | |
836 | int NCols = 0; | |
837 | int NElemBound2D = 0; | |
838 | for (int id0=0;id0<NRows;id0++) { | |
839 | NColsAtRow[id0] = tmpCols[id0]; // number of columns to store for this row | |
840 | ColAtRowBg[id0] = NElemBound2D; // begining of this row in 2D boundary surface | |
841 | NElemBound2D += tmpCols[id0]; | |
842 | if (NCols<NColsAtRow[id0]) NCols = NColsAtRow[id0]; | |
843 | } | |
844 | cheb->InitCols(NCols); | |
845 | delete[] tmpCols; | |
846 | // | |
847 | // create the 2D matrix defining the boundary of significance for 3D coeffs.matrix | |
848 | // and count the number of siginifacnt coefficients | |
849 | // | |
850 | cheb->InitElemBound2D(NElemBound2D); | |
851 | int *CoefBound2D0 = cheb->GetCoefBound2D0(); | |
852 | int *CoefBound2D1 = cheb->GetCoefBound2D1(); | |
853 | fMaxCoefs = 0; // redefine number of coeffs | |
854 | for (int id0=0;id0<NRows;id0++) { | |
855 | int nCLoc = NColsAtRow[id0]; | |
856 | int Col0 = ColAtRowBg[id0]; | |
857 | for (int id1=0;id1<nCLoc;id1++) { | |
858 | CoefBound2D0[Col0 + id1] = tmpCoefSurf[id1+id0*fNPoints[1]]; // number of coefs to store for 3-d dimension | |
859 | CoefBound2D1[Col0 + id1] = fMaxCoefs; | |
860 | fMaxCoefs += CoefBound2D0[Col0 + id1]; | |
861 | } | |
862 | } | |
863 | // | |
864 | // create final compressed 3D matrix for significant coeffs | |
865 | cheb->InitCoefs(fMaxCoefs); | |
866 | Float_t *Coefs = cheb->GetCoefs(); | |
867 | int count = 0; | |
868 | for (int id0=0;id0<NRows;id0++) { | |
869 | int ncLoc = NColsAtRow[id0]; | |
870 | int Col0 = ColAtRowBg[id0]; | |
871 | for (int id1=0;id1<ncLoc;id1++) { | |
872 | int ncf2 = CoefBound2D0[Col0 + id1]; | |
873 | for (int id2=0;id2<ncf2;id2++) { | |
874 | Coefs[count++] = tmpCoef3D[id2 + fNPoints[2]*(id1+id0*fNPoints[1])]; | |
875 | } | |
876 | } | |
877 | } | |
878 | /* | |
879 | printf("\n\nNewSurf\n"); | |
880 | for (int id0=0;id0<fNPoints[0];id0++) { | |
881 | for (int id1=0;id1<fNPoints[1];id1++) { | |
882 | printf("(%2d %2d) %2d |",id0,id1,tmpCoefSurf[id1+id0*fNPoints[1]]); | |
883 | } | |
884 | printf("\n"); | |
885 | } | |
886 | */ | |
887 | // | |
888 | delete[] tmpCoefSurf; | |
889 | delete[] tmpCoef1D; | |
890 | delete[] tmpCoef2D; | |
891 | delete[] tmpCoef3D; | |
892 | delete[] fvals; | |
893 | // | |
894 | return 1; | |
895 | } | |
896 | #endif | |
897 | ||
898 | //_______________________________________________ | |
899 | #ifdef _INC_CREATION_ALICHEB3D_ | |
900 | void AliCheb3D::SaveData(const char* outfile,Bool_t append) const | |
901 | { | |
902 | // writes coefficients data to output text file, optionallt appending on the end of existing file | |
903 | TString strf = outfile; | |
904 | gSystem->ExpandPathName(strf); | |
905 | FILE* stream = fopen(strf,append ? "a":"w"); | |
906 | SaveData(stream); | |
907 | fclose(stream); | |
908 | // | |
909 | } | |
910 | #endif | |
911 | ||
912 | //_______________________________________________ | |
913 | #ifdef _INC_CREATION_ALICHEB3D_ | |
914 | void AliCheb3D::SaveData(FILE* stream) const | |
915 | { | |
916 | // writes coefficients data to existing output stream | |
917 | // | |
918 | fprintf(stream,"\n# These are automatically generated data for the Chebyshev interpolation of 3D->%dD function\n",fDimOut); | |
919 | fprintf(stream,"#\nSTART %s\n",GetName()); | |
920 | fprintf(stream,"# Dimensionality of the output\n%d\n",fDimOut); | |
921 | fprintf(stream,"# Interpolation abs. precision\n%+.8e\n",fPrec); | |
922 | // | |
923 | fprintf(stream,"# Lower boundaries of interpolation region\n"); | |
924 | for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMin[i]); | |
925 | fprintf(stream,"# Upper boundaries of interpolation region\n"); | |
926 | for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMax[i]); | |
927 | fprintf(stream,"# Parameterization for each output dimension follows:\n",GetName()); | |
928 | // | |
929 | for (int i=0;i<fDimOut;i++) GetChebCalc(i)->SaveData(stream); | |
930 | fprintf(stream,"#\nEND %s\n#\n",GetName()); | |
931 | // | |
932 | } | |
933 | #endif | |
934 | ||
935 | //_______________________________________________ | |
936 | void AliCheb3D::LoadData(const char* inpFile) | |
937 | { | |
938 | TString strf = inpFile; | |
939 | gSystem->ExpandPathName(strf); | |
940 | FILE* stream = fopen(strf.Data(),"r"); | |
941 | LoadData(stream); | |
942 | fclose(stream); | |
943 | // | |
944 | } | |
945 | ||
946 | //_______________________________________________ | |
947 | void AliCheb3D::LoadData(FILE* stream) | |
948 | { | |
949 | if (!stream) {Error("LoadData","No stream provided.\nStop"); exit(1);} | |
950 | TString buffs; | |
951 | Clear(); | |
952 | AliCheb3DCalc::ReadLine(buffs,stream); | |
953 | if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <fit_name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);} | |
954 | SetName(buffs.Data()+buffs.First(' ')+1); | |
955 | // | |
956 | AliCheb3DCalc::ReadLine(buffs,stream); // N output dimensions | |
957 | fDimOut = buffs.Atoi(); | |
958 | if (fDimOut<1) {Error("LoadData","Expected: '<number_of_output_dimensions>', found \"%s\"\nStop\n",buffs.Data());exit(1);} | |
959 | // | |
960 | SetDimOut(fDimOut); | |
961 | // | |
962 | AliCheb3DCalc::ReadLine(buffs,stream); // Interpolation abs. precision | |
963 | fPrec = buffs.Atof(); | |
964 | if (fPrec<=0) {Error("LoadData","Expected: '<abs.precision>', found \"%s\"\nStop\n",buffs.Data());exit(1);} | |
965 | // | |
966 | for (int i=0;i<3;i++) { // Lower boundaries of interpolation region | |
967 | AliCheb3DCalc::ReadLine(buffs,stream); | |
968 | fBMin[i] = buffs.Atof(); | |
969 | } | |
970 | for (int i=0;i<3;i++) { // Upper boundaries of interpolation region | |
971 | AliCheb3DCalc::ReadLine(buffs,stream); | |
972 | fBMax[i] = buffs.Atof(); | |
973 | } | |
974 | PrepareBoundaries(fBMin,fBMax); | |
975 | // | |
976 | // data for each output dimension | |
977 | for (int i=0;i<fDimOut;i++) GetChebCalc(i)->LoadData(stream); | |
978 | // | |
979 | // check end_of_data record | |
980 | AliCheb3DCalc::ReadLine(buffs,stream); | |
981 | if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) { | |
982 | Error("LoadData","Expected \"END %s\", found \"%s\".\nStop\n",GetName(),buffs.Data()); | |
983 | exit(1); | |
984 | } | |
985 | // | |
986 | } | |
987 | ||
988 | //_______________________________________________ | |
989 | void AliCheb3D::SetDimOut(int d) | |
990 | { | |
991 | fDimOut = d; | |
992 | if (fResTmp) delete fResTmp; | |
993 | fResTmp = new Float_t[fDimOut]; // RRR | |
994 | fChebCalc.Delete(); | |
995 | for (int i=0;i<d;i++) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(),i); | |
996 | } | |
997 | ||
998 | //_______________________________________________ | |
999 | void AliCheb3D::ShiftBound(int id,float dif) | |
1000 | { | |
1001 | if (id<0||id>2) {printf("Maximum 3 dimensions are supported\n"); return;} | |
1002 | fBMin[id] += dif; | |
1003 | fBMax[id] += dif; | |
1004 | fBOffset[id] += dif; | |
1005 | } | |
1006 | ||
1007 | //_______________________________________________ | |
1008 | #ifdef _INC_CREATION_ALICHEB3D_ | |
1009 | TH1* AliCheb3D::TestRMS(int idim,int npoints,TH1* histo) | |
1010 | { | |
1011 | // fills the difference between the original function and parameterization (for idim-th component of the output) | |
1012 | // to supplied histogram. Calculations are done in npoints random points. | |
1013 | // If the hostgram was not supplied, it will be created. It is up to the user to delete it! | |
1014 | if (!fUsrMacro) { | |
1015 | printf("No user function is set\n"); | |
1016 | return 0; | |
1017 | } | |
1018 | if (!histo) histo = new TH1D(GetName(),"Control: Function - Parametrization",100,-2*fPrec,2*fPrec); | |
1019 | for (int ip=npoints;ip--;) { | |
1020 | gRandom->RndmArray(3,(Float_t *)fArgsTmp); | |
1021 | for (int i=3;i--;) fArgsTmp[i] = fBMin[i] + fArgsTmp[i]*(fBMax[i]-fBMin[i]); | |
1022 | EvalUsrFunction(); | |
1023 | Float_t valFun = fResTmp[idim]; | |
1024 | Eval(fArgsTmp,fResTmp); | |
1025 | Float_t valPar = fResTmp[idim]; | |
1026 | histo->Fill(valFun - valPar); | |
1027 | } | |
1028 | return histo; | |
1029 | // | |
1030 | } | |
1031 | #endif |