STEER/AliCheb3DCalc.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 "AliCheb3DCalc.h"
  81
  82 ClassImp(AliCheb3DCalc)
  83
  84
  85 AliCheb3DCalc::AliCheb3DCalc():
  86     TNamed("", ""),
  87     fNCoefs(0),
  88     fNRows(0),
  89     fNCols(0),
  90     fNElemBound2D(0),
  91     fNColsAtRow(0),
  92     fColAtRowBg(0),
  93     fCoefBound2D0(0),
  94     fCoefBound2D1(0),
  95     fCoefs(0),
  96     fTmpCf1(0),
  97     fTmpCf0(0)
  98 {
  99     // Default constructor
 100     Init0();
 101 }
 102
 103 AliCheb3DCalc::AliCheb3DCalc(FILE* stream):
 104     TNamed("", ""),
 105     fNCoefs(0),
 106     fNRows(0),
 107     fNCols(0),
 108     fNElemBound2D(0),
 109     fNColsAtRow(0),
 110     fColAtRowBg(0),
 111     fCoefBound2D0(0),
 112     fCoefBound2D1(0),
 113     fCoefs(0),
 114     fTmpCf1(0),
 115     fTmpCf0(0)
 116 {
 117     // Default constructor
 118     Init0();
 119     LoadData(stream);
 120 }
 121
 122
 123 AliCheb3DCalc::AliCheb3DCalc(const AliCheb3DCalc& src) :
 124     TNamed(src),
 125     fNCoefs(src.fNCoefs),
 126     fNRows(src.fNRows),
 127     fNCols(src.fNCols),
 128     fNElemBound2D(src.fNElemBound2D),
 129     fNColsAtRow(0),
 130     fColAtRowBg(0),
 131     fCoefBound2D0(0),
 132     fCoefBound2D1(0),
 133     fCoefs(0),
 134     fTmpCf1(0),
 135     fTmpCf0(0)
 136 {
 137     // Copy constructor
 138     if (src.fNColsAtRow) {
 139         fNColsAtRow = new Int_t[fNRows];
 140         for (int i=fNRows;i--;) fNColsAtRow[i] = src.fNColsAtRow[i];
 141     }
 142     if (src.fColAtRowBg) {
 143         fColAtRowBg = new Int_t[fNRows];
 144         for (int i=fNRows;i--;) fColAtRowBg[i] = src.fColAtRowBg[i];
 145     }
 146     if (src.fCoefBound2D0) {
 147         fCoefBound2D0 = new Int_t[fNElemBound2D];
 148         for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = src.fCoefBound2D0[i];
 149     }
 150     if (src.fCoefBound2D1) {
 151         fCoefBound2D1 = new Int_t[fNElemBound2D];
 152         for (int i=fNElemBound2D;i--;) fCoefBound2D1[i] = src.fCoefBound2D1[i];
 153     }
 154     if (src.fCoefs) {
 155         fCoefs = new Float_t[fNCoefs];
 156         for (int i=fNCoefs;i--;) fCoefs[i] = src.fCoefs[i];
 157     }
 158     if (src.fTmpCf1) fTmpCf1 = new Float_t[fNCols];
 159     if (src.fTmpCf0) fTmpCf0 = new Float_t[fNRows];
 160 }
 161
 162 AliCheb3DCalc& AliCheb3DCalc::operator=(const AliCheb3DCalc& rhs)
 163 {
 164     // Assignment operator
 165     if (this != &rhs) {
 166         Clear();
 167         SetName(rhs.GetName());
 168         SetTitle(rhs.GetTitle());
 169         fNCoefs = rhs.fNCoefs;
 170         fNRows  = rhs.fNRows;
 171         fNCols  = rhs.fNCols;
 172         if (rhs.fNColsAtRow) {
 173             fNColsAtRow = new Int_t[fNRows];
 174             for (int i=fNRows;i--;) fNColsAtRow[i] = rhs.fNColsAtRow[i];
 175         }
 176         if (rhs.fColAtRowBg) {
 177             fColAtRowBg = new Int_t[fNRows];
 178             for (int i=fNRows;i--;) fColAtRowBg[i] = rhs.fColAtRowBg[i];
 179         }
 180         if (rhs.fCoefBound2D0) {
 181             fCoefBound2D0 = new Int_t[fNElemBound2D];
 182             for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = rhs.fCoefBound2D0[i];
 183         }
 184         if (rhs.fCoefBound2D1) {
 185             fCoefBound2D1 = new Int_t[fNElemBound2D];
 186             for (int i=fNElemBound2D;i--;) fCoefBound2D1[i] = rhs.fCoefBound2D1[i];
 187         }
 188         if (rhs.fCoefs) {
 189             fCoefs = new Float_t[fNCoefs];
 190             for (int i=fNCoefs;i--;) fCoefs[i] = rhs.fCoefs[i];
 191         }
 192         if (rhs.fTmpCf1) fTmpCf1 = new Float_t[fNCols];
 193         if (rhs.fTmpCf0) fTmpCf0 = new Float_t[fNRows];
 194     }
 195     return *this;
 196 }
 197
 198 //__________________________________________________________________________________________
 199 void AliCheb3DCalc::Clear(Option_t*)
 200 {
 201   // delete all dynamycally allocated structures
 202   if (fTmpCf1)       { delete[] fTmpCf1;  fTmpCf1 = 0;}
 203   if (fTmpCf0)       { delete[] fTmpCf0;  fTmpCf0 = 0;}
 204   if (fCoefs)        { delete[] fCoefs;   fCoefs  = 0;}
 205   if (fCoefBound2D0) { delete[] fCoefBound2D0; fCoefBound2D0 = 0; }
 206   if (fCoefBound2D1) { delete[] fCoefBound2D1; fCoefBound2D1 = 0; }
 207   if (fNColsAtRow)   { delete[] fNColsAtRow;   fNColsAtRow = 0; }
 208   if (fColAtRowBg)   { delete[] fColAtRowBg;   fColAtRowBg = 0; }
 209   //
 210 }
 211
 212 //__________________________________________________________________________________________
 213 void AliCheb3DCalc::Init0()
 214 {
 215   // reset everything to 0
 216   fNCoefs = fNRows = fNCols = fNElemBound2D = 0;
 217   fCoefs = 0;
 218   fCoefBound2D0 = fCoefBound2D1 = 0;
 219   fNColsAtRow = fColAtRowBg = 0;
 220   fTmpCf0 = fTmpCf1 = 0;
 221 }
 222
 223 //__________________________________________________________________________________________
 224 void AliCheb3DCalc::Print(Option_t* ) const
 225 {
 226   // Print the Chebychev paramterization data
 227   printf("Chebyshev parameterization data %s for 3D->1 function.\n",GetName());
 228   int nmax3d = 0;
 229   for (int i=fNElemBound2D;i--;) if (fCoefBound2D0[i]>nmax3d) nmax3d = fCoefBound2D0[i];
 230   printf("%d coefficients in %dx%dx%d matrix\n",fNCoefs,fNRows,fNCols,nmax3d);
 231   //
 232 }
 233
 234 //__________________________________________________________________________________________
 235 Float_t  AliCheb3DCalc::Eval(Float_t  *par) const
 236 {
 237   // evaluate Chebyshev parameterization for 3D function.
 238   // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
 239   Float_t  &z = par[2];
 240   Float_t  &y = par[1];
 241   Float_t  &x = par[0];
 242   //
 243   int ncfRC;
 244   for (int id0=fNRows;id0--;) {
 245     int nCLoc = fNColsAtRow[id0];                   // number of significant coefs on this row
 246     int col0  = fColAtRowBg[id0];                   // beginning of local column in the 2D boundary matrix
 247     for (int id1=nCLoc;id1--;) {
 248       int id = id1+col0;
 249       fTmpCf1[id1] = (ncfRC=fCoefBound2D0[id]) ? ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC) : 0.0;
 250     }
 251     fTmpCf0[id0] = nCLoc>0 ? ChebEval1D(y,fTmpCf1,nCLoc):0.0;
 252   }
 253   return ChebEval1D(x,fTmpCf0,fNRows);
 254   //
 255 }
 256
 257 //_______________________________________________
 258 #ifdef _INC_CREATION_ALICHEB3D_
 259 void AliCheb3DCalc::SaveData(const char* outfile,Bool_t append) const
 260 {
 261   // writes coefficients data to output text file, optionallt appending on the end of existing file
 262   TString strf = outfile;
 263   gSystem->ExpandPathName(strf);
 264   FILE* stream = fopen(strf,append ? "a":"w");
 265   SaveData(stream);
 266   fclose(stream);
 267   //
 268 }
 269 #endif
 270
 271 //_______________________________________________
 272 #ifdef _INC_CREATION_ALICHEB3D_
 273 void AliCheb3DCalc::SaveData(FILE* stream) const
 274 {
 275   // writes coefficients data to existing output stream
 276   // Note: fNCols, fNElemBound2D and fColAtRowBg is not stored, will be computed on fly during the loading of this file
 277   fprintf(stream,"#\nSTART %s\n",GetName());
 278   fprintf(stream,"# Number of rows\n%d\n",fNRows);
 279   //
 280   fprintf(stream,"# Number of columns per row\n");
 281   for (int i=0;i<fNRows;i++) fprintf(stream,"%d\n",fNColsAtRow[i]);
 282   //
 283   fprintf(stream,"# Number of Coefs in each significant block of third dimension\n");
 284   for (int i=0;i<fNElemBound2D;i++) fprintf(stream,"%d\n",fCoefBound2D0[i]);
 285   //
 286   fprintf(stream,"# Coefficients\n");
 287   for (int i=0;i<fNCoefs;i++) fprintf(stream,"%+.8e\n",fCoefs[i]);
 288   fprintf(stream,"END %s\n",GetName());
 289   //
 290 }
 291 #endif
 292
 293 //_______________________________________________
 294 void AliCheb3DCalc::LoadData(FILE* stream)
 295 {
 296   // Load coefs. from the stream
 297   if (!stream) {Error("LoadData","No stream provided.\nStop"); exit(1);}
 298   TString buffs;
 299   Clear();
 300   ReadLine(buffs,stream);
 301   if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <fit_name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
 302   if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
 303   //
 304   ReadLine(buffs,stream); // NRows
 305   fNRows = buffs.Atoi();
 306   if (fNRows<1) {Error("LoadData","Expected: '<number_of_rows>', found \"%s\"\nStop\n",buffs.Data());exit(1);}
 307   //
 308   fNCols = 0;
 309   fNElemBound2D = 0;
 310   InitRows(fNRows);
 311   //
 312   for (int id0=0;id0<fNRows;id0++) {
 313     ReadLine(buffs,stream);               // n.cols at this row
 314     fNColsAtRow[id0] = buffs.Atoi();
 315     fColAtRowBg[id0] = fNElemBound2D;     // begining of this row in 2D boundary surface
 316     fNElemBound2D   += fNColsAtRow[id0];
 317     if (fNCols<fNColsAtRow[id0]) fNCols = fNColsAtRow[id0];
 318   }
 319   InitCols(fNCols);
 320   //
 321   fNCoefs = 0;
 322   InitElemBound2D(fNElemBound2D);
 323   //
 324   for (int i=0;i<fNElemBound2D;i++) {
 325     ReadLine(buffs,stream);               // n.coeffs at 3-d dimension for the given column/row
 326     fCoefBound2D0[i] = buffs.Atoi();
 327     fCoefBound2D1[i] = fNCoefs;
 328     fNCoefs += fCoefBound2D0[i];
 329   }
 330   //
 331   if (fNCoefs<=0) {Error("LoadData","Negtive (%d) number of Chebychef coeffs. is obtained.\nStop\n",fNCoefs);exit(1);}
 332   //
 333   InitCoefs(fNCoefs);
 334   for (int i=0;i<fNCoefs;i++) {
 335     ReadLine(buffs,stream);
 336     fCoefs[i] = buffs.Atof();
 337   }
 338   // check end_of_data record
 339   ReadLine(buffs,stream);
 340   if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {
 341     Error("LoadData","Expected \"END %s\", found \"%s\".\nStop\n",GetName(),buffs.Data());
 342     exit(1);
 343   }
 344   //
 345 }
 346
 347 //_______________________________________________
 348 void AliCheb3DCalc::ReadLine(TString& str,FILE* stream)
 349 {
 350   // read single line from the stream, skipping empty and commented lines. EOF is not expected
 351   while (str.Gets(stream)) {
 352     str = str.Strip(TString::kBoth,' ');
 353     if (str.IsNull()||str.BeginsWith("#")) continue;
 354     return;
 355   }
 356   fprintf(stderr,"AliCheb3D::ReadLine: Failed to read from stream.\nStop");exit(1); // normally, should not reach here
 357 }
 358
 359 //_______________________________________________
 360 void AliCheb3DCalc::InitCols(int nc)
 361 {
 362   // Set max.number of significant columns in the coefs matrix
 363   fNCols = nc;
 364   if (fTmpCf1) delete[] fTmpCf1;
 365   fTmpCf1 = new Float_t [fNCols];
 366 }
 367
 368 //_______________________________________________
 369 void AliCheb3DCalc::InitRows(int nr)
 370 {
 371   // Set max.number of significant rows in the coefs matrix
 372   if (fNColsAtRow) delete[] fNColsAtRow;
 373   if (fColAtRowBg) delete[] fColAtRowBg;
 374   if (fTmpCf0)     delete[] fTmpCf0;
 375   fNRows = nr;
 376   fNColsAtRow = new Int_t[fNRows];
 377   fTmpCf0     = new Float_t [fNRows];
 378   fColAtRowBg = new Int_t[fNRows];
 379   for (int i=fNRows;i--;) fNColsAtRow[i] = fColAtRowBg[i] = 0;
 380 }
 381
 382 //_______________________________________________
 383 void AliCheb3DCalc::InitElemBound2D(int ne)
 384 {
 385   // Set max number of significant coefs for given row/column of coefs 3D matrix
 386   if (fCoefBound2D0) delete[] fCoefBound2D0;
 387   if (fCoefBound2D1) delete[] fCoefBound2D1;
 388   fNElemBound2D = ne;
 389   fCoefBound2D0 = new Int_t[fNElemBound2D];
 390   fCoefBound2D1 = new Int_t[fNElemBound2D];
 391   for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = fCoefBound2D1[i] = 0;
 392 }
 393
 394 //_______________________________________________
 395 void AliCheb3DCalc::InitCoefs(int nc)
 396 {
 397   // Set total number of significant coefs
 398   if (fCoefs) delete[] fCoefs;
 399   fNCoefs = nc;
 400   fCoefs = new Float_t [fNCoefs];
 401   for (int i=fNCoefs;i--;) fCoefs[i] = 0.0;
 402 }
 403
 404