Coding violations corrected.
[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(AliCheb3D)
90
91 AliCheb3D::AliCheb3D():
92     TNamed("", ""),
93     fDimOut(0),
94     fPrec(0.),
95     fChebCalc(),
96     fMaxCoefs(0),
97     fResTmp(0),
98     fGrid(0),
99     fUsrFunName(),
100     fUsrMacro(0)             
101 {
102     // Default constructor
103     Init0();
104 }
105
106 AliCheb3D::AliCheb3D(const char* inputFile):
107     TNamed("", ""),
108     fDimOut(0),
109     fPrec(0.),
110     fChebCalc(),
111     fMaxCoefs(0),
112     fResTmp(0),
113     fGrid(0),
114     fUsrFunName(),
115     fUsrMacro(0)             
116 {
117     // Default constructor
118     Init0();
119     LoadData(inputFile);
120 }
121
122
123
124 AliCheb3D::AliCheb3D(FILE* stream):
125     TNamed("", ""),
126     fDimOut(0),
127     fPrec(0.),
128     fChebCalc(),
129     fMaxCoefs(0),
130     fResTmp(0),
131     fGrid(0),
132     fUsrFunName(),
133     fUsrMacro(0)             
134 {
135     // Default constructor
136     Init0();
137     LoadData(stream);
138 }
139
140 AliCheb3D::AliCheb3D(const AliCheb3D& src) : 
141     TNamed(src),
142     fDimOut(src.fDimOut), 
143     fPrec(src.fPrec), 
144     fChebCalc(1), 
145     fMaxCoefs(src.fMaxCoefs), 
146                                            fResTmp(0),
147     fGrid(0), 
148     fUsrFunName(src.fUsrFunName), 
149     fUsrMacro(0)
150 {
151     // Copy constructor
152     // read coefs from text file
153     for (int i=3;i--;) {
154         fBMin[i]    = src.fBMin[i];
155         fBMax[i]    = src.fBMax[i];
156         fBScale[i]  = src.fBScale[i];
157         fBOffset[i] = src.fBOffset[i];
158         fNPoints[i] = src.fNPoints[i];
159     }
160     for (int i=0;i<fDimOut;i++) {
161         AliCheb3DCalc* cbc = src.GetChebCalc(i);
162         if (cbc) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(*cbc),i);
163     }
164 }
165
166 AliCheb3D& AliCheb3D::operator=(const AliCheb3D& rhs)
167 {
168     // Assignment operator
169     if (this != &rhs) {
170         Clear();
171         fDimOut   = rhs.fDimOut;
172         fPrec     = rhs.fPrec;
173         fMaxCoefs = rhs.fMaxCoefs;
174         fUsrFunName = rhs.fUsrFunName;
175         fUsrMacro   = 0;
176         for (int i=3;i--;) {
177             fBMin[i]    = rhs.fBMin[i];
178             fBMax[i]    = rhs.fBMax[i];
179             fBScale[i]  = rhs.fBScale[i];
180             fBOffset[i] = rhs.fBOffset[i];
181             fNPoints[i] = rhs.fNPoints[i];
182         } 
183         for (int i=0;i<fDimOut;i++) {
184             AliCheb3DCalc* cbc = rhs.GetChebCalc(i);
185             if (cbc) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(*cbc),i);
186         }    
187     }
188     return *this;
189     //
190 }
191
192
193 //__________________________________________________________________________________________
194 #ifdef _INC_CREATION_ALICHEB3D_
195 AliCheb3D::AliCheb3D(const char* funName, int DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t prec) : TNamed(funName,funName)
196 {
197   // Construct the parameterization for the function
198   // funName : name of the file containing the function: void funName(Float_t * inp,Float_t * out)
199   // DimOut  : dimension of the vector computed by the user function
200   // bmin    : array of 3 elements with the lower boundaries of the region where the function is defined
201   // bmax    : array of 3 elements with the upper boundaries of the region where the function is defined
202   // npoints : array of 3 elements with the number of points to compute in each of 3 dimension
203   // prec    : max allowed absolute difference between the user function and computed parameterization on the requested grid
204   //
205   Init0();
206   fPrec = TMath::Max(1.E-12f,prec);
207   if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);}
208   SetDimOut(DimOut);
209   PrepareBoundaries(bmin,bmax);
210   DefineGrid(npoints);
211   SetUsrFunction(funName);
212   ChebFit();
213   //
214 }
215 #endif
216
217 //__________________________________________________________________________________________
218 #ifdef _INC_CREATION_ALICHEB3D_
219 AliCheb3D::AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t prec) : TNamed("AliCheb3D","AliCheb3D")
220 {
221   // Construct the parameterization for the function
222   // ptr     : pointer on the function: void fun(Float_t * inp,Float_t * out)
223   // DimOut  : dimension of the vector computed by the user function
224   // bmin    : array of 3 elements with the lower boundaries of the region where the function is defined
225   // bmax    : array of 3 elements with the upper boundaries of the region where the function is defined
226   // npoints : array of 3 elements with the number of points to compute in each of 3 dimension
227   // prec    : max allowed absolute difference between the user function and computed parameterization on the requested grid
228   //
229   Init0();
230   fPrec = TMath::Max(1.E-12f,prec);
231   if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);}
232   SetDimOut(DimOut);
233   PrepareBoundaries(bmin,bmax);
234   DefineGrid(npoints);
235   SetUsrFunction(ptr);
236   ChebFit();
237   //
238 }
239 #endif
240
241
242 //__________________________________________________________________________________________
243 void AliCheb3D::Clear(Option_t*)
244 {
245 // Clean-up
246   if (fResTmp)        { delete[] fResTmp; fResTmp = 0; }
247   if (fGrid)          { delete[] fGrid;   fGrid   = 0; }
248   if (fUsrMacro)      { delete fUsrMacro; fUsrMacro = 0;}
249   fChebCalc.Delete();
250   //
251 }
252
253 //__________________________________________________________________________________________
254 void AliCheb3D::Print(Option_t* opt) const
255 {
256     // Print Chebyshev parameterisation data
257   printf("%s: Chebyshev parameterization for 3D->%dD function. Precision: %e\n",GetName(),fDimOut,fPrec);
258   printf("Region of validity: [%+.5e:%+.5e] [%+.5e:%+.5e] [%+.5e:%+.5e]\n",fBMin[0],fBMax[0],fBMin[1],fBMax[1],fBMin[2],fBMax[2]);
259   TString opts = opt; opts.ToLower();
260   if (opts.Contains("l")) for (int i=0;i<fDimOut;i++) {printf("Output dimension %d:\n",i+1); GetChebCalc(i)->Print();}
261   //
262 }
263
264 //__________________________________________________________________________________________
265 void AliCheb3D::Init0()
266 {
267   // Initialisation  
268   for (int i=3;i--;) fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0;
269   fMaxCoefs = 0;
270   fGrid = 0;
271   fResTmp = 0;
272   fUsrFunName = "";
273   fUsrMacro = 0;
274 #ifdef _INC_CREATION_ALICHEB3D_
275   gUsrFunAliCheb3D = 0;
276 #endif
277 }
278
279 //__________________________________________________________________________________________
280 void AliCheb3D::PrepareBoundaries(Float_t  *bmin,Float_t  *bmax)
281 {
282   // Set and check boundaries defined by user, prepare coefficients for their conversion to [-1:1] interval
283   //
284   for (int i=3;i--;) {
285     fBMin[i]   = bmin[i];
286     fBMax[i]   = bmax[i];
287     fBScale[i] = bmax[i]-bmin[i];
288     if (fBScale[i]<=0) { 
289       Error("PrepareBoundaries","Boundaries for %d-th dimension are not increasing: %+.4e %+.4e\nStop\n",i,fBMin[i],fBMax[i]);
290       exit(1);
291     }
292     fBOffset[i] = bmin[i] + fBScale[i]/2.0;
293     fBScale[i] = 2./fBScale[i];
294   }
295   //
296 }
297
298 //__________________________________________________________________________________________
299 #ifdef _INC_CREATION_ALICHEB3D_
300 void AliCheb3D::SetUsrFunction(const char* name)
301 {
302   // load user macro with function definition and compile it
303   gUsrFunAliCheb3D = 0; 
304   fUsrFunName = name;
305   gSystem->ExpandPathName(fUsrFunName);
306   if (fUsrMacro) delete fUsrMacro;
307   TString tmpst = fUsrFunName;
308   tmpst += "+"; // prepare filename to compile
309   if (gROOT->LoadMacro(tmpst.Data())) {Error("SetUsrFunction","Failed to load user function from %s\nStop\n",name); exit(1);}
310   fUsrMacro = new TMethodCall();        
311   tmpst = tmpst.Data() + tmpst.Last('/')+1; //Strip away any path preceding the macro file name
312   int dot = tmpst.Last('.');
313   if (dot>0) tmpst.Resize(dot);
314   fUsrMacro->InitWithPrototype(tmpst.Data(),"Float_t *,Float_t *");
315   long args[2];
316   args[0] = (long)fArgsTmp;
317   args[1] = (long)fResTmp;
318   fUsrMacro->SetParamPtrs(args); 
319   //
320 }
321 #endif
322
323 //__________________________________________________________________________________________
324 #ifdef _INC_CREATION_ALICHEB3D_
325 void AliCheb3D::SetUsrFunction(void (*ptr)(float*,float*))
326 {
327   if (fUsrMacro) delete fUsrMacro;
328   fUsrMacro = 0;
329   fUsrFunName = "";
330   gUsrFunAliCheb3D = ptr;
331 }
332 #endif
333
334 //__________________________________________________________________________________________
335 #ifdef _INC_CREATION_ALICHEB3D_
336 void AliCheb3D::EvalUsrFunction(Float_t  *x, Float_t  *res) {
337   for (int i=3;i--;) fArgsTmp[i] = x[i];
338   if   (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp);
339   else fUsrMacro->Execute(); 
340   for (int i=fDimOut;i--;) res[i] = fResTmp[i];
341 }
342 #endif
343
344 //__________________________________________________________________________________________
345 #ifdef _INC_CREATION_ALICHEB3D_
346 Int_t AliCheb3D::CalcChebCoefs(Float_t  *funval,int np, Float_t  *outCoefs, Float_t  prec)
347 {
348   // Calculate Chebyshev coeffs using precomputed function values at np roots.
349   // If prec>0, estimate the highest coeff number providing the needed precision
350   //
351   double sm;                 // do summations in double to minimize the roundoff error
352   for (int ic=0;ic<np;ic++) { // compute coeffs
353     sm = 0;          
354     for (int ir=0;ir<np;ir++) {
355       float  rt = TMath::Cos( ic*(ir+0.5)*TMath::Pi()/np);
356       sm += funval[ir]*rt;
357     }
358     outCoefs[ic] = Float_t( sm * ((ic==0) ? 1./np : 2./np) );
359   }
360   //
361   if (prec<=0) return np;
362   //
363   sm = 0;
364   int cfMax = 0;
365   for (cfMax=np;cfMax--;) {
366     sm += TMath::Abs(outCoefs[cfMax]);
367     if (sm>=prec) break;
368   }
369   if (++cfMax==0) cfMax=1;
370   return cfMax;
371   //
372 }
373 #endif
374
375 //__________________________________________________________________________________________
376 #ifdef _INC_CREATION_ALICHEB3D_
377 void AliCheb3D::DefineGrid(Int_t* npoints)
378 {
379   // prepare the grid of Chebyshev roots in each dimension
380   const int kMinPoints = 1;
381   int ntot = 0;
382   fMaxCoefs = 1;
383   for (int id=3;id--;) { 
384     fNPoints[id] = npoints[id];
385     if (fNPoints[id]<kMinPoints) {
386       Error("DefineGrid","at %d-th dimension %d point is requested, at least %d is needed\nStop\n",fNPoints[id],kMinPoints);
387       exit(1);
388     }
389     ntot += fNPoints[id];
390     fMaxCoefs *= fNPoints[id];
391   }
392   fGrid = new Float_t [ntot];
393   //
394   int curp = 0;
395   for (int id=3;id--;) { 
396     int np = fNPoints[id];
397     fGridOffs[id] = curp;
398     for (int ip=0;ip<np;ip++) {
399       Float_t x = TMath::Cos( TMath::Pi()*(ip+0.5)/np );
400       fGrid[curp++] = MapToExternal(x,id);
401     }
402   }
403   //
404 }
405 #endif
406
407 //__________________________________________________________________________________________
408 #ifdef _INC_CREATION_ALICHEB3D_
409 Int_t AliCheb3D::ChebFit()
410 {
411   // prepare parameterization for all output dimensions
412   int ir=0; 
413   for (int i=fDimOut;i--;) ir+=ChebFit(i); 
414   return ir;
415 }
416 #endif
417
418 //__________________________________________________________________________________________
419 #ifdef _INC_CREATION_ALICHEB3D_
420 Int_t AliCheb3D::ChebFit(int dmOut)
421 {
422   // prepare paramaterization of 3D function for dmOut-th dimension 
423   int maxDim = 0;
424   for (int i=0;i<3;i++) if (maxDim<fNPoints[i]) maxDim = fNPoints[i];
425   Float_t  *fvals      = new Float_t [ fNPoints[0] ];
426   Float_t  *tmpCoef3D  = new Float_t [ fNPoints[0]*fNPoints[1]*fNPoints[2] ]; 
427   Float_t  *tmpCoef2D  = new Float_t [ fNPoints[0]*fNPoints[1] ]; 
428   Float_t  *tmpCoef1D  = new Float_t [ maxDim ];
429   //
430   Float_t RTiny = fPrec/Float_t(maxDim); // neglect coefficient below this threshold
431   //
432   // 1D Cheb.fit for 0-th dimension at current steps of remaining dimensions
433   int ncmax = 0;
434   //
435   AliCheb3DCalc* cheb =  GetChebCalc(dmOut);
436   //
437   for (int id2=fNPoints[2];id2--;) {
438     fArgsTmp[2] = fGrid[ fGridOffs[2]+id2 ];
439     //
440     for (int id1=fNPoints[1];id1--;) {
441       fArgsTmp[1] = fGrid[ fGridOffs[1]+id1 ];
442       //
443       for (int id0=fNPoints[0];id0--;) {
444         fArgsTmp[0] = fGrid[ fGridOffs[0]+id0 ];
445         EvalUsrFunction();     // compute function values at Chebyshev roots of 0-th dimension
446         fvals[id0] =  fResTmp[dmOut];
447       }
448       int nc = CalcChebCoefs(fvals,fNPoints[0], tmpCoef1D, fPrec);
449       for (int id0=fNPoints[0];id0--;) tmpCoef2D[id1 + id0*fNPoints[1]] = tmpCoef1D[id0];
450       if (ncmax<nc) ncmax = nc;              // max coefs to be kept in dim0 to guarantee needed precision
451     }
452     //
453     // once each 1d slice of given 2d slice is parametrized, parametrize the Cheb.coeffs
454     for (int id0=fNPoints[0];id0--;) {
455       CalcChebCoefs( tmpCoef2D+id0*fNPoints[1], fNPoints[1], tmpCoef1D, -1);
456       for (int id1=fNPoints[1];id1--;) tmpCoef3D[id2 + fNPoints[2]*(id1+id0*fNPoints[1])] = tmpCoef1D[id1];
457     }
458   }
459   //
460   // now fit the last dimensions Cheb.coefs
461   for (int id0=fNPoints[0];id0--;) {
462     for (int id1=fNPoints[1];id1--;) {
463       CalcChebCoefs( tmpCoef3D+ fNPoints[2]*(id1+id0*fNPoints[1]), fNPoints[2], tmpCoef1D, -1);
464       for (int id2=fNPoints[2];id2--;) tmpCoef3D[id2+ fNPoints[2]*(id1+id0*fNPoints[1])] = tmpCoef1D[id2]; // store on place
465     }
466   }
467   //
468   // now find 2D surface which separates significant coefficients of 3D matrix from nonsignificant ones (up to fPrec)
469   int *tmpCoefSurf = new Int_t[ fNPoints[0]*fNPoints[1] ];
470   for (int id0=fNPoints[0];id0--;) for (int id1=fNPoints[1];id1--;) tmpCoefSurf[id1+id0*fNPoints[1]]=0;  
471   Double_t resid = 0;
472   for (int id0=fNPoints[0];id0--;) {
473     for (int id1=fNPoints[1];id1--;) {
474       for (int id2=fNPoints[2];id2--;) {
475         int id = id2 + fNPoints[2]*(id1+id0*fNPoints[1]);
476         Float_t  cfa = TMath::Abs(tmpCoef3D[id]);
477         if (cfa < RTiny) {tmpCoef3D[id] = 0; continue;} // neglect coeefs below the threshold
478
479         resid += cfa;
480         if (resid<fPrec) continue; // this coeff is negligible
481         // otherwise go back 1 step
482         resid -= cfa;
483         tmpCoefSurf[id1+id0*fNPoints[1]] = id2+1; // how many coefs to keep
484         break;
485       }
486     }
487   }
488   /*
489   printf("\n\nCoeffs\n");  
490   int cnt = 0;
491   for (int id0=0;id0<fNPoints[0];id0++) {
492     for (int id1=0;id1<fNPoints[1];id1++) {
493       for (int id2=0;id2<fNPoints[2];id2++) {
494         printf("%2d%2d%2d %+.4e |",id0,id1,id2,tmpCoef3D[cnt++]);
495       }
496       printf("\n");
497     }
498     printf("\n");
499   }
500   */
501   // see if there are rows to reject, find max.significant column at each row
502   int NRows = fNPoints[0];
503   int *tmpCols = new int[NRows]; 
504   for (int id0=fNPoints[0];id0--;) {
505     int id1 = fNPoints[1];
506     while (id1>0 && tmpCoefSurf[(id1-1)+id0*fNPoints[1]]==0) id1--;
507     tmpCols[id0] = id1;
508   }
509   // find max significant row
510   for (int id0=NRows;id0--;) {if (tmpCols[id0]>0) break; NRows--;}
511   // find max significant column and fill the permanent storage for the max sigificant column of each row
512   cheb->InitRows(NRows);                  // create needed arrays;
513   int *NColsAtRow = cheb->GetNColsAtRow();
514   int *ColAtRowBg = cheb->GetColAtRowBg();
515   int NCols = 0;
516   int NElemBound2D = 0;
517   for (int id0=0;id0<NRows;id0++) {
518     NColsAtRow[id0] = tmpCols[id0];     // number of columns to store for this row
519     ColAtRowBg[id0] = NElemBound2D;     // begining of this row in 2D boundary surface
520     NElemBound2D += tmpCols[id0];
521     if (NCols<NColsAtRow[id0]) NCols = NColsAtRow[id0];
522   }
523   cheb->InitCols(NCols);
524   delete[] tmpCols;
525   //  
526   // create the 2D matrix defining the boundary of significance for 3D coeffs.matrix 
527   // and count the number of siginifacnt coefficients
528   //
529   cheb->InitElemBound2D(NElemBound2D);
530   int *CoefBound2D0 = cheb->GetCoefBound2D0();
531   int *CoefBound2D1 = cheb->GetCoefBound2D1();
532   fMaxCoefs = 0; // redefine number of coeffs
533   for (int id0=0;id0<NRows;id0++) {
534     int nCLoc = NColsAtRow[id0];
535     int col0  = ColAtRowBg[id0];
536     for (int id1=0;id1<nCLoc;id1++) {
537       CoefBound2D0[col0 + id1] = tmpCoefSurf[id1+id0*fNPoints[1]];  // number of coefs to store for 3-d dimension
538       CoefBound2D1[col0 + id1] = fMaxCoefs;
539       fMaxCoefs += CoefBound2D0[col0 + id1];
540     }
541   }
542   //
543   // create final compressed 3D matrix for significant coeffs
544   cheb->InitCoefs(fMaxCoefs);
545   Float_t  *Coefs = cheb->GetCoefs();
546   int count = 0;
547   for (int id0=0;id0<NRows;id0++) {
548     int ncLoc = NColsAtRow[id0];
549     int col0  = ColAtRowBg[id0];
550     for (int id1=0;id1<ncLoc;id1++) {
551       int ncf2 = CoefBound2D0[col0 + id1];
552       for (int id2=0;id2<ncf2;id2++) {
553         Coefs[count++] = tmpCoef3D[id2 + fNPoints[2]*(id1+id0*fNPoints[1])];
554       }
555     }
556   }
557   /*
558   printf("\n\nNewSurf\n");
559   for (int id0=0;id0<fNPoints[0];id0++) {
560     for (int id1=0;id1<fNPoints[1];id1++) {
561       printf("(%2d %2d) %2d |",id0,id1,tmpCoefSurf[id1+id0*fNPoints[1]]);  
562     }
563     printf("\n");
564   }
565   */
566   //
567   delete[] tmpCoefSurf;
568   delete[] tmpCoef1D;
569   delete[] tmpCoef2D;
570   delete[] tmpCoef3D;
571   delete[] fvals;
572   //
573   return 1;
574 }
575 #endif
576
577 //_______________________________________________
578 #ifdef _INC_CREATION_ALICHEB3D_
579 void AliCheb3D::SaveData(const char* outfile,Bool_t append) const
580 {
581   // writes coefficients data to output text file, optionallt appending on the end of existing file
582   TString strf = outfile;
583   gSystem->ExpandPathName(strf);
584   FILE* stream = fopen(strf,append ? "a":"w");
585   SaveData(stream);
586   fclose(stream);
587   //
588 }
589 #endif
590
591 //_______________________________________________
592 #ifdef _INC_CREATION_ALICHEB3D_
593 void AliCheb3D::SaveData(FILE* stream) const
594 {
595   // writes coefficients data to existing output stream
596   //
597   fprintf(stream,"\n# These are automatically generated data for the Chebyshev interpolation of 3D->%dD function\n",fDimOut); 
598   fprintf(stream,"#\nSTART %s\n",GetName());
599   fprintf(stream,"# Dimensionality of the output\n%d\n",fDimOut);
600   fprintf(stream,"# Interpolation abs. precision\n%+.8e\n",fPrec);
601   //
602   fprintf(stream,"# Lower boundaries of interpolation region\n");
603   for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMin[i]);
604   fprintf(stream,"# Upper boundaries of interpolation region\n");
605   for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMax[i]);
606   fprintf(stream,"# Parameterization for each output dimension follows:\n",GetName());
607   //
608   for (int i=0;i<fDimOut;i++) GetChebCalc(i)->SaveData(stream);
609   fprintf(stream,"#\nEND %s\n#\n",GetName());
610   //
611 }
612 #endif
613
614 //_______________________________________________
615 void AliCheb3D::LoadData(const char* inpFile)
616 {
617   // Load data from input file  
618   TString strf = inpFile;
619   gSystem->ExpandPathName(strf);
620   FILE* stream = fopen(strf.Data(),"r");
621   LoadData(stream);
622   fclose(stream);
623   //
624 }
625
626 //_______________________________________________
627 void AliCheb3D::LoadData(FILE* stream)
628 {
629   // Load data from input stream stream  
630   if (!stream) {Error("LoadData","No stream provided.\nStop"); exit(1);}
631   TString buffs;
632   Clear();
633   AliCheb3DCalc::ReadLine(buffs,stream);
634   if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <fit_name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
635   SetName(buffs.Data()+buffs.First(' ')+1);
636   //
637   AliCheb3DCalc::ReadLine(buffs,stream); // N output dimensions
638   fDimOut = buffs.Atoi(); 
639   if (fDimOut<1) {Error("LoadData","Expected: '<number_of_output_dimensions>', found \"%s\"\nStop\n",buffs.Data());exit(1);}
640   //
641   SetDimOut(fDimOut);
642   //
643   AliCheb3DCalc::ReadLine(buffs,stream); // Interpolation abs. precision
644   fPrec = buffs.Atof();
645   if (fPrec<=0) {Error("LoadData","Expected: '<abs.precision>', found \"%s\"\nStop\n",buffs.Data());exit(1);}
646   //
647   for (int i=0;i<3;i++) { // Lower boundaries of interpolation region
648     AliCheb3DCalc::ReadLine(buffs,stream);
649     fBMin[i] = buffs.Atof(); 
650   }
651   for (int i=0;i<3;i++) { // Upper boundaries of interpolation region
652     AliCheb3DCalc::ReadLine(buffs,stream);
653     fBMax[i] = buffs.Atof(); 
654   }
655   PrepareBoundaries(fBMin,fBMax);
656   //
657   // data for each output dimension
658   for (int i=0;i<fDimOut;i++) GetChebCalc(i)->LoadData(stream);
659   //
660   // check end_of_data record
661   AliCheb3DCalc::ReadLine(buffs,stream);
662   if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {
663     Error("LoadData","Expected \"END %s\", found \"%s\".\nStop\n",GetName(),buffs.Data());
664     exit(1);
665   }
666   //
667 }
668
669 //_______________________________________________
670 void AliCheb3D::SetDimOut(int d)
671 {
672   // Set the dimension of the output array 
673   fDimOut = d;
674   if (fResTmp) delete fResTmp;
675   fResTmp = new Float_t[fDimOut]; // RRR
676   fChebCalc.Delete();
677   for (int i=0;i<d;i++) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(),i);
678 }
679
680 //_______________________________________________
681 void AliCheb3D::ShiftBound(int id,float dif)
682 {
683   //Shift the boundary of dimension id
684   if (id<0||id>2) {printf("Maximum 3 dimensions are supported\n"); return;}
685   fBMin[id] += dif;
686   fBMax[id] += dif;
687   fBOffset[id] += dif;
688 }
689
690 //_______________________________________________
691 #ifdef _INC_CREATION_ALICHEB3D_
692 TH1* AliCheb3D::TestRMS(int idim,int npoints,TH1* histo)
693 {
694   // fills the difference between the original function and parameterization (for idim-th component of the output)
695   // to supplied histogram. Calculations are done in npoints random points. 
696   // If the hostgram was not supplied, it will be created. It is up to the user to delete it! 
697   if (!fUsrMacro) {
698     printf("No user function is set\n");
699     return 0;
700   }
701   if (!histo) histo = new TH1D(GetName(),"Control: Function - Parametrization",100,-2*fPrec,2*fPrec);
702   for (int ip=npoints;ip--;) {
703     gRandom->RndmArray(3,(Float_t *)fArgsTmp);
704     for (int i=3;i--;) fArgsTmp[i] = fBMin[i] + fArgsTmp[i]*(fBMax[i]-fBMin[i]);
705     EvalUsrFunction();
706     Float_t valFun = fResTmp[idim];
707     Eval(fArgsTmp,fResTmp);
708     Float_t valPar = fResTmp[idim];
709     histo->Fill(valFun - valPar);
710   }
711   return histo;
712   //
713 }
714 #endif