modifications to satisfy the coding conventions
[u/mrichter/AliRoot.git] / STEER / AliCheb3D.h
1 // Author: ruben.shahoyan@cern.ch   09/09/2006
2
3 ////////////////////////////////////////////////////////////////////////////////
4 //                                                                            //
5 // AliCheb3D produces the interpolation of the user 3D->NDimOut arbitrary     //
6 // function supplied in "void (*fcn)(float* inp,float* out)" format           //
7 // either in a separate macro file or as a function pointer.                  //
8 // Only coefficients needed to guarantee the requested precision are kept.    //
9 //                                                                            //
10 // The user-callable methods are:                                             //
11 // To create the interpolation use:                                           //
12 // AliCheb3D(const char* funName,  // name of the file with user function     //
13 //          or                                                                //
14 // AliCheb3D(void (*ptr)(float*,float*),// pointer on the  user function      //
15 //        Int_t     DimOut,     // dimensionality of the function's output    // 
16 //        Float_t  *bmin,       // lower 3D bounds of interpolation domain    // 
17 //        Float_t  *bmax,       // upper 3D bounds of interpolation domain    // 
18 //        Int_t    *npoints,    // number of points in each of 3 input        //
19 //                              // dimension, defining the interpolation grid //
20 //        Float_t   prec=1E-6); // requested max.absolute difference between  //
21 //                              // the interpolation and any point on grid    //
22 //                                                                            //
23 // To test obtained parameterization use the method                           //
24 // TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0);                    // 
25 // it will compare the user output of the user function and interpolation     //
26 // for idim-th output dimension and fill the difference in the supplied       //
27 // histogram. If no histogram is supplied, it will be created.                //
28 //                                                                            //
29 // To save the interpolation data:                                            //
30 // SaveData(const char* filename, Bool_t append )                             //
31 // write text file with data. If append is kTRUE and the output file already  //
32 // exists, data will be added in the end of the file.                         //
33 // Alternatively, SaveData(FILE* stream) will write the data to               //
34 // already existing stream.                                                   //
35 //                                                                            //
36 // To read back already stored interpolation use either the constructor       // 
37 // AliCheb3D(const char* inpFile);                                            //
38 // or the default constructor AliCheb3D() followed by                         //
39 // AliCheb3D::LoadData(const char* inpFile);                                  //
40 //                                                                            //
41 // To compute the interpolation use Eval(float* par,float *res) method, with  //
42 // par being 3D vector of arguments (inside the validity region) and res is   //
43 // the array of DimOut elements for the output.                               //
44 //                                                                            //
45 // If only one component (say, idim-th) of the output is needed, use faster   //
46 // Float_t Eval(Float_t *par,int idim) method.                                //
47 //                                                                            //
48 // void Print(option="") will print the name, the ranges of validity and      //
49 // the absolute precision of the parameterization. Option "l" will also print //
50 // the information about the number of coefficients for each output           //
51 // dimension.                                                                 //
52 //                                                                            //
53 // NOTE: during the evaluation no check is done for parameter vector being    //
54 // outside the interpolation region. If there is such a risk, use             //
55 // Bool_t IsInside(float *par) method. Chebyshev parameterization is not      //
56 // good for extrapolation!                                                    //
57 //                                                                            //
58 // For the properties of Chebyshev parameterization see:                      //
59 // H.Wind, CERN EP Internal Report, 81-12/Rev.                                //
60 //                                                                            //
61 ////////////////////////////////////////////////////////////////////////////////
62
63
64 #ifndef ALICHEB3D_H
65 #define ALICHEB3D_H
66
67 #include <TNamed.h>
68 #include <TObjArray.h>
69 #include "AliCheb3DCalc.h"
70
71 class TString;
72 class TSystem;
73 class TRandom;
74 class TH1;
75 class TMethodCall;
76 class TRandom;
77 class TROOT;
78 class stdio;
79
80
81
82 class AliCheb3D: public TNamed 
83 {
84  public:
85   AliCheb3D();
86   AliCheb3D(const AliCheb3D& src);
87   AliCheb3D(const char* inpFile);
88   AliCheb3D(FILE* stream);
89   //
90 #ifdef _INC_CREATION_ALICHEB3D_
91   AliCheb3D(const char* funName, Int_t DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t  prec=1E-6);
92   AliCheb3D(void (*ptr)(float*,float*), Int_t DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t  prec=1E-6);
93   AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npX,Int_t *npY,Int_t *npZ, Float_t prec=1E-6);
94   AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t  *bmin,Float_t  *bmax, Float_t prec=1E-6);
95 #endif
96   //
97   ~AliCheb3D()                                                                 {Clear();}
98   //
99   AliCheb3D&   operator=(const AliCheb3D& rhs);
100   void         Eval(const Float_t  *par,Float_t  *res);
101   Float_t      Eval(const Float_t  *par,int idim);
102   //
103   void         EvalDeriv(int dimd, const Float_t  *par, Float_t  *res);
104   void         EvalDeriv2(int dimd1, int dimd2, const Float_t  *par,Float_t  *res);
105   Float_t      EvalDeriv(int dimd, const Float_t  *par, int idim);
106   Float_t      EvalDeriv2(int dimd1,int dimd2, const Float_t  *par, int idim);
107   void         EvalDeriv3D(const Float_t *par, Float_t dbdr[3][3]); 
108   void         EvalDeriv3D2(const Float_t *par, Float_t dbdrdr[3][3][3]); 
109   void         Print(const Option_t* opt="")                             const;
110   Bool_t       IsInside(const Float_t  *par)                             const;
111   Bool_t       IsInside(const Double_t *par)                             const;
112   AliCheb3DCalc*  GetChebCalc(int i)                                     const {return (AliCheb3DCalc*)fChebCalc.UncheckedAt(i);}
113   Float_t      GetBoundMin(int i)                                        const {return fBMin[i];}
114   Float_t      GetBoundMax(int i)                                        const {return fBMax[i];}
115   Float_t*     GetBoundMin()                                             const {return (float*)fBMin;}
116   Float_t*     GetBoundMax()                                             const {return (float*)fBMax;}
117   Float_t      GetPrecision()                                            const {return fPrec;}
118   void         ShiftBound(int id,float dif);
119   //
120   void         LoadData(const char* inpFile);
121   void         LoadData(FILE* stream);
122   //
123 #ifdef _INC_CREATION_ALICHEB3D_
124   int*         GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec);
125   void         EstimateNPoints(float Prec, int gridBC[3][3]);
126   void         SaveData(const char* outfile,Bool_t append=kFALSE)        const;
127   void         SaveData(FILE* stream=stdout)                             const;
128   //
129   void         SetUsrFunction(const char* name);
130   void         SetUsrFunction(void (*ptr)(float*,float*));
131   void         EvalUsrFunction(const Float_t  *x, const Float_t  *res);
132   TH1*         TestRMS(int idim,int npoints = 1000,TH1* histo=0);
133   static Int_t CalcChebCoefs(const Float_t  *funval,int np, Float_t  *outCoefs, Float_t  prec=-1);
134 #endif
135   //
136  protected:
137   void         Clear(const Option_t* option = "");
138   void         SetDimOut(const int d);
139   void         PrepareBoundaries(const Float_t  *bmin,const Float_t  *bmax);
140   //
141 #ifdef _INC_CREATION_ALICHEB3D_
142   void         EvalUsrFunction();
143   void         DefineGrid(Int_t* npoints);
144   Int_t        ChebFit();                                                                 // fit all output dimensions
145   Int_t        ChebFit(int dmOut);
146 #endif
147   //
148   Float_t      MapToInternal(Float_t  x,Int_t d) const; // map x to [-1:1]
149   Float_t      MapToExternal(Float_t  x,Int_t d) const {return x/fBScale[d]+fBOffset[d];}   // map from [-1:1] to x
150   //  
151  protected:
152   Int_t        fDimOut;            // dimension of the ouput array
153   Float_t      fPrec;              // requested precision
154   Float_t      fBMin[3];           // min boundaries in each dimension
155   Float_t      fBMax[3];           // max boundaries in each dimension  
156   Float_t      fBScale[3];         // scale for boundary mapping to [-1:1] interval
157   Float_t      fBOffset[3];        // offset for boundary mapping to [-1:1] interval
158   TObjArray    fChebCalc;          // Chebyshev parameterization for each output dimension
159   //
160   Int_t        fMaxCoefs;          //! max possible number of coefs per parameterization
161   Int_t        fNPoints[3];        //! number of used points in each dimension
162   Float_t      fArgsTmp[3];        //! temporary vector for coefs caluclation
163   Float_t      fBuff[6];           //! buffer for coordinate transformations
164   Float_t *    fResTmp;            //! temporary vector for results of user function caluclation
165   Float_t *    fGrid;              //! temporary buffer for Chebyshef roots grid
166   Int_t        fGridOffs[3];       //! start of grid for each dimension
167   TString      fUsrFunName;        //! name of user macro containing the function of  "void (*fcn)(float*,float*)" format
168   TMethodCall* fUsrMacro;          //! Pointer to MethodCall for function from user macro 
169   //
170   ClassDef(AliCheb3D,1)  // Chebyshev parametrization for 3D->N function
171 };
172
173 //__________________________________________________________________________________________
174 inline Bool_t  AliCheb3D::IsInside(const Float_t  *par) const 
175 {
176   // check if the point is inside of the fitted box
177   const float kTol = 1.e-4; 
178   for (int i=3;i--;) if (fBMin[i]-par[i]>kTol || par[i]-fBMax[i]>kTol) return kFALSE;
179   //if(par[i]<fBMin[i]||par[i]>fBMax[i]) return kFALSE;
180   return kTRUE;
181 }
182
183 //__________________________________________________________________________________________
184 inline Bool_t  AliCheb3D::IsInside(const Double_t  *par) const 
185 {
186   // check if the point is inside of the fitted box
187   const float kTol = 1.e-4; 
188   for (int i=3;i--;) if (fBMin[i]-par[i]>kTol || par[i]-fBMax[i]>kTol) return kFALSE;
189   //if(par[i]<fBMin[i]||par[i]>fBMax[i]) return kFALSE;
190   return kTRUE;
191 }
192
193 //__________________________________________________________________________________________
194 inline void AliCheb3D::Eval(const Float_t  *par, Float_t  *res)
195 {
196   // evaluate Chebyshev parameterization for 3d->DimOut function
197   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
198   for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp);
199   //
200 }
201
202 //__________________________________________________________________________________________
203 inline Float_t AliCheb3D::Eval(const Float_t  *par, int idim)
204 {
205   // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
206   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
207   return GetChebCalc(idim)->Eval(fArgsTmp);
208   //
209 }
210
211 //__________________________________________________________________________________________
212 inline void AliCheb3D::EvalDeriv3D(const Float_t *par, Float_t dbdr[3][3])
213 {
214   // return gradient matrix
215   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
216   for (int ib=3;ib--;) for (int id=3;id--;) dbdr[ib][id] = GetChebCalc(ib)->EvalDeriv(id,fArgsTmp)*fBScale[id];
217 }
218
219 //__________________________________________________________________________________________
220 inline void AliCheb3D::EvalDeriv3D2(const Float_t *par, Float_t dbdrdr[3][3][3])
221 {
222   // return gradient matrix
223   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
224   for (int ib=3;ib--;) for (int id=3;id--;)for (int id1=3;id1--;) 
225     dbdrdr[ib][id][id1] = GetChebCalc(ib)->EvalDeriv2(id,id1,fArgsTmp)*fBScale[id]*fBScale[id1];
226 }
227
228 //__________________________________________________________________________________________
229 inline void AliCheb3D::EvalDeriv(int dimd, const Float_t  *par, Float_t  *res)
230 {
231   // evaluate Chebyshev parameterization derivative for 3d->DimOut function
232   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
233   for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];;
234   //
235 }
236
237 //__________________________________________________________________________________________
238 inline void AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t  *par, Float_t  *res)
239 {
240   // evaluate Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function
241   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
242   for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2];
243   //
244 }
245
246 //__________________________________________________________________________________________
247 inline Float_t AliCheb3D::EvalDeriv(int dimd, const Float_t  *par, int idim)
248 {
249   // evaluate Chebyshev parameterization derivative over dimd dimention for idim-th output dimension of 3d->DimOut function
250   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
251   return GetChebCalc(idim)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];
252   //
253 }
254
255 //__________________________________________________________________________________________
256 inline Float_t AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t  *par, int idim)
257 {
258   // evaluate Chebyshev parameterization 2ns derivative over dimd1 and dimd2 dimensions for idim-th output dimension of 3d->DimOut function
259   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
260   return GetChebCalc(idim)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2];
261   //
262 }
263
264 //__________________________________________________________________________________________
265 inline Float_t AliCheb3D::MapToInternal(Float_t  x,Int_t d) const
266 {
267   // map x to [-1:1]
268 #ifdef _BRING_TO_BOUNDARY_
269   float res = (x-fBOffset[d])*fBScale[d];
270   if (res<-1) return -1;
271   if (res> 1) return 1;
272   return res;
273 #else
274   return (x-fBOffset[d])*fBScale[d];
275 #endif
276 }
277
278 #endif