STEER/AliCheb3D.h

   1 #ifndef ALICHEB3D_H
   2 #define ALICHEB3D_H
   3 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  * See cxx source for full Copyright notice                               */
   5
   6 /* $Id$ */
   7
   8 // Author: ruben.shahoyan@cern.ch   09/09/2006
   9 //
  10 ////////////////////////////////////////////////////////////////////////////////
  11 //                                                                            //
  12 // AliCheb3D produces the interpolation of the user 3D->NDimOut arbitrary     //
  13 // function supplied in "void (*fcn)(float* inp,float* out)" format           //
  14 // either in a separate macro file or as a function pointer.                  //
  15 // Only coefficients needed to guarantee the requested precision are kept.    //
  16 //                                                                            //
  17 // The user-callable methods are:                                             //
  18 // To create the interpolation use:                                           //
  19 // AliCheb3D(const char* funName,  // name of the file with user function     //
  20 //          or                                                                //
  21 // AliCheb3D(void (*ptr)(float*,float*),// pointer on the  user function      //
  22 //        Int_t     DimOut,     // dimensionality of the function's output    //
  23 //        Float_t  *bmin,       // lower 3D bounds of interpolation domain    //
  24 //        Float_t  *bmax,       // upper 3D bounds of interpolation domain    //
  25 //        Int_t    *npoints,    // number of points in each of 3 input        //
  26 //                              // dimension, defining the interpolation grid //
  27 //        Float_t   prec=1E-6); // requested max.absolute difference between  //
  28 //                              // the interpolation and any point on grid    //
  29 //                                                                            //
  30 // To test obtained parameterization use the method                           //
  31 // TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0);                    //
  32 // it will compare the user output of the user function and interpolation     //
  33 // for idim-th output dimension and fill the difference in the supplied       //
  34 // histogram. If no histogram is supplied, it will be created.                //
  35 //                                                                            //
  36 // To save the interpolation data:                                            //
  37 // SaveData(const char* filename, Bool_t append )                             //
  38 // write text file with data. If append is kTRUE and the output file already  //
  39 // exists, data will be added in the end of the file.                         //
  40 // Alternatively, SaveData(FILE* stream) will write the data to               //
  41 // already existing stream.                                                   //
  42 //                                                                            //
  43 // To read back already stored interpolation use either the constructor       //
  44 // AliCheb3D(const char* inpFile);                                            //
  45 // or the default constructor AliCheb3D() followed by                         //
  46 // AliCheb3D::LoadData(const char* inpFile);                                  //
  47 //                                                                            //
  48 // To compute the interpolation use Eval(float* par,float *res) method, with  //
  49 // par being 3D vector of arguments (inside the validity region) and res is   //
  50 // the array of DimOut elements for the output.                               //
  51 //                                                                            //
  52 // If only one component (say, idim-th) of the output is needed, use faster   //
  53 // Float_t Eval(Float_t *par,int idim) method.                                //
  54 //                                                                            //
  55 // void Print(option="") will print the name, the ranges of validity and      //
  56 // the absolute precision of the parameterization. Option "l" will also print //
  57 // the information about the number of coefficients for each output           //
  58 // dimension.                                                                 //
  59 //                                                                            //
  60 // NOTE: during the evaluation no check is done for parameter vector being    //
  61 // outside the interpolation region. If there is such a risk, use             //
  62 // Bool_t IsInside(float *par) method. Chebyshev parameterization is not      //
  63 // good for extrapolation!                                                    //
  64 //                                                                            //
  65 // For the properties of Chebyshev parameterization see:                      //
  66 // H.Wind, CERN EP Internal Report, 81-12/Rev.                                //
  67 //                                                                            //
  68 ////////////////////////////////////////////////////////////////////////////////
  69
  70
  71 #include <stdio.h>
  72 #include <TNamed.h>
  73 #include <TMethodCall.h>
  74 #include <TMath.h>
  75 #include <TH1.h>
  76 #include <TObjArray.h>
  77
  78 class TString;
  79 class TSystem;
  80 class TRandom;
  81
  82
  83 // to decrease the compilable code size comment this define. This will exclude the routines
  84 // used for the calculation and saving of the coefficients.
  85 // #define _INC_CREATION_ALICHEB3D_
  86
  87 class AliCheb3DCalc: public TNamed
  88 {
  89  public:
  90   AliCheb3DCalc();
  91   AliCheb3DCalc(FILE* stream); // read coefs from text file
  92   AliCheb3DCalc(const AliCheb3DCalc &src);
  93   AliCheb3DCalc& operator= (const AliCheb3DCalc &rhs);
  94   ~AliCheb3DCalc() {Clear();}
  95   //
  96   void       Print(Option_t* opt="")                                   const;
  97   void       LoadData(FILE* stream);
  98   Float_t    Eval(Float_t  *par)                                       const;
  99   //
 100 #ifdef _INC_CREATION_ALICHEB3D_
 101   void       SaveData(const char* outfile,Bool_t append=kFALSE)        const;
 102   void       SaveData(FILE* stream=stdout)                             const;
 103 #endif
 104   //
 105   static void ReadLine(TString& str,FILE* stream);
 106   //
 107  protected:
 108   //
 109   void       Clear(Option_t* option = "");
 110   void       Init0();
 111   Float_t    ChebEval1D(Float_t  x, const Float_t * array, int ncf)     const;
 112   void       InitRows(int nr);
 113   void       InitCols(int nc);
 114   void       InitElemBound2D(int ne);
 115   void       InitCoefs(int nc);
 116   Int_t*     GetNColsAtRow()                                            const {return fNColsAtRow;}
 117   Int_t*     GetColAtRowBg()                                            const {return fColAtRowBg;}
 118   Int_t*     GetCoefBound2D0()                                          const {return fCoefBound2D0;}
 119   Int_t*     GetCoefBound2D1()                                          const {return fCoefBound2D1;}
 120   Float_t *  GetCoefs()                                                 const {return fCoefs;}
 121   //
 122  protected:
 123   Int_t      fNCoefs;            // total number of coeeficients
 124   Int_t      fNRows;             // number of significant rows in the 3D coeffs matrix
 125   Int_t      fNCols;             // max number of significant cols in the 3D coeffs matrix
 126   Int_t      fNElemBound2D;      // number of elements (fNRows*fNCols) to store for the 2D boundary of significant coeffs
 127   Int_t*     fNColsAtRow;        //[fNRows] number of sighificant columns (2nd dim) at each row of 3D coefs matrix
 128   Int_t*     fColAtRowBg;        //[fNRows] beginnig of significant columns (2nd dim) for row in the 2D boundary matrix
 129   Int_t*     fCoefBound2D0;      //[fNElemBound2D] 2D matrix defining the boundary of significance for 3D coeffs.matrix (Ncoefs for col/row)
 130   Int_t*     fCoefBound2D1;      //[fNElemBound2D] 2D matrix defining the start beginnig of significant coeffs for col/row
 131   Float_t *  fCoefs;             //[fNCoefs] array of Chebyshev coefficients
 132   //
 133   Float_t *  fTmpCf1;            //[fNCols] temp. coeffs for 2d summation
 134   Float_t *  fTmpCf0;            //[fNRows] temp. coeffs for 1d summation
 135   //
 136   ClassDef(AliCheb3DCalc,1)  // Class for interpolation of 3D->1 function by Chebyshev parametrization
 137 };
 138
 139
 140 class AliCheb3D: public TNamed
 141 {
 142  public:
 143     AliCheb3D();
 144     AliCheb3D(const char* inpFile); // read coefs from text file
 145     AliCheb3D(FILE*);               // read coefs from stream
 146     AliCheb3D(const AliCheb3D &src);
 147     AliCheb3D& operator= (const AliCheb3D &rhs);
 148
 149   //
 150 #ifdef _INC_CREATION_ALICHEB3D_
 151   AliCheb3D(const char* funName, Int_t DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t  prec=1E-6);
 152   AliCheb3D(void (*ptr)(float*,float*), Int_t DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t  prec=1E-6);
 153 #endif
 154   //
 155   ~AliCheb3D()                                                                 {Clear();}
 156   //
 157   void         Eval(Float_t  *par,Float_t  *res);
 158   Float_t      Eval(Float_t  *par,int idim);
 159   void         Print(Option_t* opt="")                                   const;
 160   Bool_t       IsInside(Float_t  *par)                                   const;
 161   AliCheb3DCalc*  GetChebCalc(int i)                                     const {return (AliCheb3DCalc*)fChebCalc.UncheckedAt(i);}
 162   Float_t      GetBoundMin(int i)                                        const {return fBMin[i];}
 163   Float_t      GetBoundMax(int i)                                        const {return fBMax[i];}
 164   Float_t      GetPrecision()                                            const {return fPrec;}
 165   void         ShiftBound(int id,float dif);
 166   //
 167   void         LoadData(const char* inpFile);
 168   void         LoadData(FILE* stream);
 169   //
 170 #ifdef _INC_CREATION_ALICHEB3D_
 171   void         SaveData(const char* outfile,Bool_t append=kFALSE)        const;
 172   void         SaveData(FILE* stream=stdout)                             const;
 173   //
 174   void         SetUsrFunction(const char* name);
 175   void         SetUsrFunction(void (*ptr)(float*,float*));
 176   void         EvalUsrFunction(Float_t  *x, Float_t  *res);
 177   TH1*         TestRMS(int idim,int npoints = 1000,TH1* histo=0);
 178 #endif
 179   //
 180  protected:
 181   void         Init0();
 182   void         Clear(Option_t* option = "");
 183   void         SetDimOut(int d);
 184   void         PrepareBoundaries(Float_t  *bmin,Float_t  *bmax);
 185   //
 186 #ifdef _INC_CREATION_ALICHEB3D_
 187   void         EvalUsrFunction();
 188   void         DefineGrid(Int_t* npoints);
 189   Int_t        ChebFit();                                                                 // fit all output dimensions
 190   Int_t        ChebFit(int dmOut);
 191   Int_t        CalcChebCoefs(Float_t  *funval,int np, Float_t  *outCoefs, Float_t  prec=-1);
 192 #endif
 193   //
 194   void         Cyl2CartCyl(float *rphiz, float *b) const;
 195   void         Cart2Cyl(float *xyz,float *rphiz) const;
 196   //
 197   Float_t      MapToInternal(Float_t  x,Int_t d) const {return (x-fBOffset[d])*fBScale[d];} // map x to [-1:1]
 198   Float_t      MapToExternal(Float_t  x,Int_t d) const {return x/fBScale[d]+fBOffset[d];}   // map from [-1:1] to x
 199   //
 200  protected:
 201   Int_t        fDimOut;            // dimension of the ouput array
 202   Float_t      fPrec;              // requested precision
 203   Float_t      fBMin[3];           // min boundaries in each dimension
 204   Float_t      fBMax[3];           // max boundaries in each dimension
 205   Float_t      fBScale[3];         // scale for boundary mapping to [-1:1] interval
 206   Float_t      fBOffset[3];        // offset for boundary mapping to [-1:1] interval
 207   TObjArray    fChebCalc;          // Chebyshev parameterization for each output dimension
 208   //
 209   Int_t        fMaxCoefs;          //! max possible number of coefs per parameterization
 210   Int_t        fNPoints[3];        //! number of used points in each dimension
 211   Float_t      fArgsTmp[3];        //! temporary vector for coefs caluclation
 212   Float_t      fBuff[6];           //! buffer for coordinate transformations
 213   Float_t *    fResTmp;            //! temporary vector for results of user function caluclation
 214   Float_t *    fGrid;              //! temporary buffer for Chebyshef roots grid
 215   Int_t        fGridOffs[3];       //! start of grid for each dimension
 216   TString      fUsrFunName;        //! name of user macro containing the function of  "void (*fcn)(float*,float*)" format
 217   TMethodCall* fUsrMacro;          //! Pointer to MethodCall for function from user macro
 218   //
 219   ClassDef(AliCheb3D,1)  // Chebyshev parametrization for 3D->N function
 220 };
 221
 222 // Pointer on user function (faster altrnative to TMethodCall)
 223 #ifdef _INC_CREATION_ALICHEB3D_
 224 void (*gUsrFunAliCheb3D) (float* ,float* );
 225 #endif
 226
 227 //__________________________________________________________________________________________
 228 #ifdef _INC_CREATION_ALICHEB3D_
 229 inline void AliCheb3D::EvalUsrFunction()
 230 {
 231   // call user supplied function
 232   if   (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp);
 233   else fUsrMacro->Execute();
 234 }
 235 #endif
 236
 237 //__________________________________________________________________________________________
 238 inline Bool_t  AliCheb3D::IsInside(Float_t  *par) const
 239 {
 240   // check if the point is inside of the fitted box
 241   for (int i=3;i--;) if(par[i]<fBMin[i]||par[i]>fBMax[i]) return kFALSE;
 242   return kTRUE;
 243 }
 244
 245 //__________________________________________________________________________________________
 246 inline Float_t AliCheb3DCalc::ChebEval1D(Float_t  x, const Float_t * array, int ncf ) const
 247 {
 248   // evaluate 1D Chebyshev parameterization. x is the argument mapped to [-1:1] interval
 249   Float_t b0, b1, b2;
 250   Float_t x2 = x+x;
 251   b0 = array[--ncf];
 252   b1 = b2 = 0;
 253   for (int i=ncf;i--;) {
 254     b2 = b1;
 255     b1 = b0;
 256     b0 = array[i] + x2*b1 -b2;
 257   }
 258   return b0 - x*b1;
 259 }
 260
 261 //__________________________________________________________________________________________
 262 inline void AliCheb3D::Eval(Float_t  *par, Float_t  *res)
 263 {
 264   // evaluate Chebyshev parameterization for 3d->DimOut function
 265   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 266   for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp);
 267   //
 268 }
 269
 270 //__________________________________________________________________________________________
 271 inline Float_t AliCheb3D::Eval(Float_t  *par, int idim)
 272 {
 273   // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
 274   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 275   return GetChebCalc(idim)->Eval(fArgsTmp);
 276   //
 277 }
 278
 279 //__________________________________________________________________________________________________
 280 inline void AliCheb3D::Cyl2CartCyl(float *rphiz, float *b) const
 281 {
 282   // convert field in cylindrical coordinates to cartesian system, point is in cyl.system
 283   float btr = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
 284   float ang = TMath::ATan2(b[1],b[0]) + rphiz[1];
 285   b[0] = btr*TMath::Cos(ang);
 286   b[1] = btr*TMath::Sin(ang);
 287   //
 288 }
 289
 290 //__________________________________________________________________________________________________
 291 inline void AliCheb3D::Cart2Cyl(float *xyz,float *rphiz) const
 292 {
 293   // convert cartesian coordinate to cylindrical one
 294   rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 295   rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
 296   rphiz[2] = xyz[2];
 297 }
 298
 299 #endif