AliRoot Core  3dc7879 (3dc7879)
AliCheb2DStack.h
Go to the documentation of this file.
1 #ifndef ALICHEB2DSTACK_H
2 #define ALICHEB2DSTACK_H
3 
4 #include <TObject.h>
5 
6 /*****************************************************************************
7  * Base class for stack of 2D->ND Chebishev parameterizations *
8  * *
9  * Author: ruben.shahoyan@cern.ch *
10  *****************************************************************************/
11 
12 // when _BRING_TO_BOUNDARY2D_ is defined, the point outside of the fitted folume is assumed
13 // to be on the surface
14 #define _BRING_TO_BOUNDARY2D_
15 //
16 
17 typedef void (*stFun_t)(int,float*,float*); // original distribution provided via such function
18 
19 class AliCheb2DStack : public TObject
20 {
21  public:
22  enum {kMaxPoints=255}; // max number of points in input dimension
23  enum {ktgp,kz};
24  //
25  public:
27  virtual ~AliCheb2DStack();
28  //
29  AliCheb2DStack(int nSlices, int dimOut, const float bmin[2], const float bmax[2],const float* dead=0, const float *rowXI=0);
30  //
31  Int_t GetNSlices() const {return fNSlices;}
32  Int_t GetDimOut() const {return fDimOut;}
33  const Float_t* GetBoundMin() const {return fBMin;}
34  const Float_t* GetBoundMax() const {return fBMax;}
35  //
36  void SetXRowInv(const float* xi) {fRowXI = xi;}
37  const float* GetXRowInv() const {return fRowXI;}
38  virtual void Eval(int sliceID, const float *par, float *res) const = 0;
39  virtual Float_t Eval(int sliceID, int dimOut, const float *par) const = 0;
40  virtual void EvalDeriv(int sliceID, int dim, const Float_t *par, float* res) const = 0;
41 
42  Bool_t IsInside(const float *par) const;
43  //
44  void Print(const Option_t* opt="") const;
45  virtual void PrintSlice(int isl, const Option_t* opt) const = 0;
46  //
47  static float ChebEval1D(float x, const float* array, int ncf);
48  static float ChebEval1D(float x, const Short_t* array, int ncf);
49  static float ChebEval1Deriv(float x, const float* array, int ncf);
50  static float ChebEval1Deriv(float x, const Short_t* array, int ncf);
51  static void SetDefPrecision(float prc=1e-4) {fgkDefPrec = prc>1e-8 ? prc:1e-8;}
52  static float GetDefPrecision() {return fgkDefPrec;}
53  //
54  protected:
55  void MapToInternal(int slice, const float* xy, float *xyint) const;
56  void MapToInternal(int slice, const float* xy, float &x1, float &x2) const;
57  float MapToExternal(int slice, float x,int dim) const;
58  float* DefineGrid(int slice, int dim, const int np[2]) const;
59  void CheckDimensions(const int *np) const;
60  Int_t CalcChebCoefs(const float *funval,int np, float *outCoefs, float prec);
61  //
62  protected:
63  //
64  Int_t fDimOut; // each slice maps 2D to fDimOut values
65  Int_t fNSlices; // number of slices in the stack
66  Int_t fNParams; // number of parameterizations = fNSlices*fDimOut
67  Int_t fNCoefsTot; // dimension of coeffs array for all slices
68  Int_t fNRowsTot; // total number of 1D Cheb. param rows
69  Float_t fBMin[2]; // min boundaries in each dimension
70  Float_t fBMax[2]; // max boundaries in each dimension
71  Float_t fBScaleZ; // scale for Z boundary mapping to [-1:1] interval
72  Float_t fBOffsetZ; // offset for Z boundary mapping to [-1:1] interval
73  Float_t fDead[2]; // dead zone in cm to account if X for each row is provided
74  const Float_t* fRowXI;
75  //
76  UChar_t* fNRows; //[fNParams] N of used rows in the 2D coeffs matrix of each param
77  UChar_t* fNCols; //[fNRowsTot] N of used columns in each row
78  Int_t* fCoeffsEntry; //[fNSlices] start of the coeffs array in fCoeffs for each slice
79  Int_t* fColEntry; //[fNSlices] start of the Ncolumns array in fNCols for each slice
80  //
81  static Float_t fgkDefPrec; // default precision
82  static Float_t fWSpace[kMaxPoints]; // workspace
83 
84  //
85  private:
86  AliCheb2DStack(const AliCheb2DStack& src); // dummy
87  AliCheb2DStack& operator=(const AliCheb2DStack& rhs); // dummy
88  //
89  ClassDef(AliCheb2DStack,2) // stack of 2D->fDimOut Chebyshev parameterization slices
90 };
91 
92 
93 //_________________________________________________________
94 inline void AliCheb2DStack::MapToInternal(int slice, const float* xy, float args[2]) const
95 {
96  // map xy to [-1:1]
97  // for Z we don't have dead zones
98  args[kz] = (xy[kz]-fBOffsetZ)*fBScaleZ;
99  float tmn = fBMin[ktgp], tmx = fBMax[ktgp];
100  if (fRowXI) {
101  tmn += fDead[0]*fRowXI[slice];
102  tmx -= fDead[1]*fRowXI[slice];
103  }
104  args[ktgp] = 2.f*(xy[ktgp]-tmn)/(tmx-tmn) - 1.f;
105 #ifdef _BRING_TO_BOUNDARY2D_
106  if (args[kz]<-1.0f) args[kz]=-1.0f;
107  else if (args[kz]> 1.0f) args[kz]=1.0f;
108  if (args[ktgp]<-1.0f) args[ktgp]=-1.0f;
109  else if (args[ktgp]> 1.0f) args[ktgp]=1.0f;
110 #endif
111 }
112 
113 //_________________________________________________________
114 inline void AliCheb2DStack::MapToInternal(int slice, const float* xy, float &x0, float &x1) const
115 {
116  // map xy to [-1:1]
117  x1 = (xy[kz]-fBOffsetZ)*fBScaleZ;
118  float tmn = fBMin[ktgp], tmx = fBMax[ktgp];
119  if (fRowXI) {
120  tmn += fDead[0]*fRowXI[slice];
121  tmx -= fDead[1]*fRowXI[slice];
122  }
123  x0 = 2.f*(xy[ktgp]-tmn)/(tmx-tmn) - 1.f;
124  //
125 #ifdef _BRING_TO_BOUNDARY2D_
126  if (x0 < -1.0f) x0 =-1.0f; else if (x0 > 1.0f) x0 = 1.0f;
127  if (x1 < -1.0f) x1 =-1.0f; else if (x1 > 1.0f) x1 = 1.0f;
128 #endif
129 }
130 
131 //__________________________________________________________________________________________
132 inline float AliCheb2DStack::MapToExternal(int slice, float x,int dim) const
133 {
134  // map from [-1:1] to x for dim-th input dimension
135  if (dim==kz) return x/fBScaleZ+fBOffsetZ;
136  // for y/x need to account for dead zone
137  float tmn = fBMin[ktgp], tmx = fBMax[ktgp];
138  if (fRowXI) {
139  tmn += fDead[0]*fRowXI[slice];
140  tmx -= fDead[1]*fRowXI[slice];
141  }
142  return 0.5*(x+1.0f)*(tmx-tmn)+tmn;
143 }
144 
145 //__________________________________________________________________________________________
146 inline Bool_t AliCheb2DStack::IsInside(const float *par) const
147 {
148  // check if the point is inside of the fitted box
149  for (int i=2;i--;) if (fBMin[i]>par[i] || par[i]>fBMax[i]) return kFALSE;
150  return kTRUE;
151 }
152 
153 //__________________________________________________________________________________________
154 inline float AliCheb2DStack::ChebEval1D(float x, const float* array, int ncf)
155 {
156  // evaluate 1D Chebyshev parameterization. x is the argument mapped to [-1:1] interval
157  if (!ncf) return 0;
158  float b0(array[--ncf]), b1(0), b2(0), x2(x+x);
159  for (int i=ncf;i--;) {
160  b2 = b1;
161  b1 = b0;
162  b0 = array[i] + x2*b1 -b2;
163  }
164  return b0 - x*b1;
165  //
166 }
167 
168 //__________________________________________________________________________________________
169 inline float AliCheb2DStack::ChebEval1D(float x, const Short_t* array, int ncf)
170 {
171  // evaluate 1D Chebyshev parameterization. x is the argument mapped to [-1:1] interval
172  if (!ncf) return 0;
173  float b0(array[--ncf]), b1(0), b2(0), x2(x+x);
174  for (int i=ncf;i--;) {
175  b2 = b1;
176  b1 = b0;
177  b0 = array[i] + x2*b1 -b2;
178  }
179  return b0 - x*b1;
180  //
181 }
182 
183 
184 //__________________________________________________________________________________________
185 inline float AliCheb2DStack::ChebEval1Deriv(float x, const float* array, int ncf)
186 {
187  // evaluate 1D Chebyshev parameterization. x is the argument mapped to [-1:1] interval
188  if (--ncf<1) return 0;
189  Float_t b0, b1(0), b2(0), x2(x+x);
190  float dcf0(0),dcf1,dcf2(0);
191  b0 = dcf1 = 2*ncf*array[ncf];
192  if (!(--ncf)) return b0/2;
193 
194  for (int i=ncf;i--;) {
195  b2 = b1;
196  b1 = b0;
197  dcf0 = dcf2 + 2*(i+1)*array[i+1];
198  b0 = dcf0 + x2*b1 -b2;
199  dcf2 = dcf1;
200  dcf1 = dcf0;
201  }
202  return b0 - x*b1 - dcf0/2;
203  //
204 }
205 
206 //__________________________________________________________________________________________
207 inline float AliCheb2DStack::ChebEval1Deriv(float x, const Short_t* array, int ncf)
208 {
209  // evaluate 1D Chebyshev parameterization. x is the argument mapped to [-1:1] interval
210  if (--ncf<1) return 0;
211  Float_t b0, b1(0), b2(0), x2(x+x);
212  float dcf0(0),dcf1,dcf2(0);
213  b0 = dcf1 = 2*ncf*array[ncf];
214  if (!(--ncf)) return b0/2;
215 
216  for (int i=ncf;i--;) {
217  b2 = b1;
218  b1 = b0;
219  dcf0 = dcf2 + 2*(i+1)*array[i+1];
220  b0 = dcf0 + x2*b1 -b2;
221  dcf2 = dcf1;
222  dcf1 = dcf0;
223  }
224  return b0 - x*b1 - dcf0/2;
225  //
226 }
227 
228 
229 #endif
void MapToInternal(int slice, const float *xy, float *xyint) const
const Float_t * fRowXI
Int_t CalcChebCoefs(const float *funval, int np, float *outCoefs, float prec)
UChar_t * fNRows
optional external!!! set 1/X for each row if dead zones to be accounted
Int_t GetDimOut() const
static float GetDefPrecision()
virtual void PrintSlice(int isl, const Option_t *opt) const =0
void SetXRowInv(const float *xi)
AliCheb2DStack & operator=(const AliCheb2DStack &rhs)
Int_t * fCoeffsEntry
Float_t fDead[2]
Float_t fBMax[2]
const Float_t * GetBoundMin() const
const Float_t * GetBoundMax() const
TObjArray * array
Definition: AnalyzeLaser.C:12
UChar_t * fNCols
Bool_t IsInside(const float *par) const
static Float_t fWSpace[kMaxPoints]
Float_t fBMin[2]
static float ChebEval1D(float x, const float *array, int ncf)
float MapToExternal(int slice, float x, int dim) const
Int_t GetNSlices() const
void Print(const Option_t *opt="") const
virtual ~AliCheb2DStack()
virtual void EvalDeriv(int sliceID, int dim, const Float_t *par, float *res) const =0
static float ChebEval1Deriv(float x, const float *array, int ncf)
const float * GetXRowInv() const
void(* stFun_t)(int, float *, float *)
void CheckDimensions(const int *np) const
TF1 * f
Definition: interpolTest.C:21
float * DefineGrid(int slice, int dim, const int np[2]) const
static Float_t fgkDefPrec
virtual void Eval(int sliceID, const float *par, float *res) const =0
static void SetDefPrecision(float prc=1e-4)
void res(Char_t i)
Definition: Resolution.C:2