3 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * See cxx source for full Copyright notice */
8 // Author: ruben.shahoyan@cern.ch 09/09/2006
10 ////////////////////////////////////////////////////////////////////////////////
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. //
17 // The user-callable methods are: //
18 // To create the interpolation use: //
19 // AliCheb3D(const char* funName, // name of the file with user function //
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 //
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. //
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. //
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); //
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. //
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. //
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 //
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! //
65 // For the properties of Chebyshev parameterization see: //
66 // H.Wind, CERN EP Internal Report, 81-12/Rev. //
68 ////////////////////////////////////////////////////////////////////////////////
73 #include <TMethodCall.h>
76 #include <TObjArray.h>
78 #include "AliCheb3DCalc.h"
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
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);
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);
101 ~AliCheb3D() {Clear();}
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);
113 void LoadData(const char* inpFile);
114 void LoadData(FILE* stream);
116 #ifdef _INC_CREATION_ALICHEB3D_
117 void SaveData(const char* outfile,Bool_t append=kFALSE) const;
118 void SaveData(FILE* stream=stdout) const;
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);
128 void Clear(Option_t* option = "");
129 void SetDimOut(int d);
130 void PrepareBoundaries(Float_t *bmin,Float_t *bmax);
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);
140 void Cyl2CartCyl(float *rphiz, float *b) const;
141 void Cart2Cyl(float *xyz,float *rphiz) const;
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
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
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
165 ClassDef(AliCheb3D,1) // Chebyshev parametrization for 3D->N function
168 // Pointer on user function (faster altrnative to TMethodCall)
169 #ifdef _INC_CREATION_ALICHEB3D_
170 void (*gUsrFunAliCheb3D) (float* ,float* );
173 //__________________________________________________________________________________________
174 #ifdef _INC_CREATION_ALICHEB3D_
175 inline void AliCheb3D::EvalUsrFunction()
177 // call user supplied function
178 if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp);
179 else fUsrMacro->Execute();
183 //__________________________________________________________________________________________
184 inline Bool_t AliCheb3D::IsInside(Float_t *par) const
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;
191 //__________________________________________________________________________________________
192 inline void AliCheb3D::Eval(Float_t *par, Float_t *res)
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);
200 //__________________________________________________________________________________________
201 inline Float_t AliCheb3D::Eval(Float_t *par, int idim)
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);
209 //__________________________________________________________________________________________________
210 inline void AliCheb3D::Cyl2CartCyl(float *rphiz, float *b) const
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);
220 //__________________________________________________________________________________________________
221 inline void AliCheb3D::Cart2Cyl(float *xyz,float *rphiz) const
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]);