Copy constructors and asignment operators implemented (R. Shahoyan)
[u/mrichter/AliRoot.git] / STEER / AliCheb3D.cxx
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"
85 #include "AliLog.h"
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
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), 
142     fTmpCf0(0)
143 {
144     // Copy constructor
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];
167 }
168
169 AliCheb3DCalc& AliCheb3DCalc::operator=(const AliCheb3DCalc& rhs)
170 {
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;
203 }
204
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.),
419     fChebCalc(),
420     fMaxCoefs(0),
421     fResTmp(0),
422     fGrid(0),
423     fUsrFunName(),
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.),
434     fChebCalc(),
435     fMaxCoefs(0),
436     fResTmp(0),
437     fGrid(0),
438     fUsrFunName(),
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.),
452     fChebCalc(),
453     fMaxCoefs(0),
454     fResTmp(0),
455     fGrid(0),
456     fUsrFunName(),
457     fUsrMacro(0)             
458 {
459     // Default constructor
460     Init0();
461     LoadData(stream);
462 }
463
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)
474 {
475     // Copy constructor
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     //
514 }
515
516
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
565
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