Completely Updated (Mario Sitta)
[u/mrichter/AliRoot.git] / 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 #include "AliCheb3DCalc.h"
79
80 class TString;
81 class TSystem;
82 class TRandom;
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 class AliCheb3D: public TNamed 
87 {
88  public:
89     AliCheb3D();
90     AliCheb3D(const char* inpFile);      // read coefs from text file
91     AliCheb3D(FILE* file);               // read coefs from stream
92     AliCheb3D(const AliCheb3D &src);
93     AliCheb3D& operator= (const AliCheb3D &rhs);
94   
95   //
96 #ifdef _INC_CREATION_ALICHEB3D_
97   AliCheb3D(const char* funName, Int_t DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t  prec=1E-6);
98   AliCheb3D(void (*ptr)(float*,float*), Int_t DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t  prec=1E-6);
99 #endif
100   //
101   ~AliCheb3D()                                                                 {Clear();}
102   //
103   void         Eval(Float_t  *par,Float_t  *res);
104   Float_t      Eval(Float_t  *par,int idim);
105   void         Print(Option_t* opt="")                                   const;
106   Bool_t       IsInside(Float_t  *par)                                   const;
107   AliCheb3DCalc*  GetChebCalc(int i)                                     const {return (AliCheb3DCalc*)fChebCalc.UncheckedAt(i);}
108   Float_t      GetBoundMin(int i)                                        const {return fBMin[i];}
109   Float_t      GetBoundMax(int i)                                        const {return fBMax[i];}
110   Float_t      GetPrecision()                                            const {return fPrec;}
111   void         ShiftBound(int id,float dif);
112   //
113   void         LoadData(const char* inpFile);
114   void         LoadData(FILE* stream);
115   //
116 #ifdef _INC_CREATION_ALICHEB3D_
117   void         SaveData(const char* outfile,Bool_t append=kFALSE)        const;
118   void         SaveData(FILE* stream=stdout)                             const;
119   //
120   void         SetUsrFunction(const char* name);
121   void         SetUsrFunction(void (*ptr)(float*,float*));
122   void         EvalUsrFunction(Float_t  *x, Float_t  *res);
123   TH1*         TestRMS(int idim,int npoints = 1000,TH1* histo=0);
124 #endif
125   //
126  protected:
127   void         Init0();
128   void         Clear(Option_t* option = "");
129   void         SetDimOut(int d);
130   void         PrepareBoundaries(Float_t  *bmin,Float_t  *bmax);
131   //
132 #ifdef _INC_CREATION_ALICHEB3D_
133   void         EvalUsrFunction();
134   void         DefineGrid(Int_t* npoints);
135   Int_t        ChebFit();                                                                 // fit all output dimensions
136   Int_t        ChebFit(int dmOut);
137   Int_t        CalcChebCoefs(Float_t  *funval,int np, Float_t  *outCoefs, Float_t  prec=-1);
138 #endif
139   //
140   void         Cyl2CartCyl(float *rphiz, float *b) const;
141   void         Cart2Cyl(float *xyz,float *rphiz) const;
142   //
143   Float_t      MapToInternal(Float_t  x,Int_t d) const {return (x-fBOffset[d])*fBScale[d];} // map x to [-1:1]
144   Float_t      MapToExternal(Float_t  x,Int_t d) const {return x/fBScale[d]+fBOffset[d];}   // map from [-1:1] to x
145   //
146  protected:
147   Int_t        fDimOut;            // dimension of the ouput array
148   Float_t      fPrec;              // requested precision
149   Float_t      fBMin[3];           // min boundaries in each dimension
150   Float_t      fBMax[3];           // max boundaries in each dimension  
151   Float_t      fBScale[3];         // scale for boundary mapping to [-1:1] interval
152   Float_t      fBOffset[3];        // offset for boundary mapping to [-1:1] interval
153   TObjArray    fChebCalc;          // Chebyshev parameterization for each output dimension
154   //
155   Int_t        fMaxCoefs;          //! max possible number of coefs per parameterization
156   Int_t        fNPoints[3];        //! number of used points in each dimension
157   Float_t      fArgsTmp[3];        //! temporary vector for coefs caluclation
158   Float_t      fBuff[6];           //! buffer for coordinate transformations
159   Float_t *    fResTmp;            //! temporary vector for results of user function caluclation
160   Float_t *    fGrid;              //! temporary buffer for Chebyshef roots grid
161   Int_t        fGridOffs[3];       //! start of grid for each dimension
162   TString      fUsrFunName;        //! name of user macro containing the function of  "void (*fcn)(float*,float*)" format
163   TMethodCall* fUsrMacro;          //! Pointer to MethodCall for function from user macro 
164   //
165   ClassDef(AliCheb3D,1)  // Chebyshev parametrization for 3D->N function
166 };
167
168 // Pointer on user function (faster altrnative to TMethodCall)
169 #ifdef _INC_CREATION_ALICHEB3D_
170 void (*gUsrFunAliCheb3D) (float* ,float* );
171 #endif
172
173 //__________________________________________________________________________________________
174 #ifdef _INC_CREATION_ALICHEB3D_
175 inline void AliCheb3D::EvalUsrFunction() 
176 {
177   // call user supplied function
178   if   (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp);
179   else fUsrMacro->Execute(); 
180 }
181 #endif
182
183 //__________________________________________________________________________________________
184 inline Bool_t  AliCheb3D::IsInside(Float_t  *par) const 
185 {
186   // check if the point is inside of the fitted box
187   for (int i=3;i--;) if(par[i]<fBMin[i]||par[i]>fBMax[i]) return kFALSE;
188   return kTRUE;
189 }
190
191 //__________________________________________________________________________________________
192 inline void AliCheb3D::Eval(Float_t  *par, Float_t  *res)
193 {
194   // evaluate Chebyshev parameterization for 3d->DimOut function
195   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
196   for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp);
197   //
198 }
199
200 //__________________________________________________________________________________________
201 inline Float_t AliCheb3D::Eval(Float_t  *par, int idim)
202 {
203   // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
204   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
205   return GetChebCalc(idim)->Eval(fArgsTmp);
206   //
207 }
208
209 //__________________________________________________________________________________________________
210 inline void AliCheb3D::Cyl2CartCyl(float *rphiz, float *b) const
211 {
212   // convert field in cylindrical coordinates to cartesian system, point is in cyl.system
213   float btr = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
214   float ang = TMath::ATan2(b[1],b[0]) + rphiz[1];
215   b[0] = btr*TMath::Cos(ang);
216   b[1] = btr*TMath::Sin(ang);
217   //
218 }
219
220 //__________________________________________________________________________________________________
221 inline void AliCheb3D::Cart2Cyl(float *xyz,float *rphiz) const
222 {
223   // convert cartesian coordinate to cylindrical one
224   rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
225   rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
226   rphiz[2] = xyz[2];
227 }
228
229 #endif