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