Line data Source code
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:
26 : AliCheb2DStack();
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 0 : Int_t GetNSlices() const {return fNSlices;}
32 0 : Int_t GetDimOut() const {return fDimOut;}
33 0 : const Float_t* GetBoundMin() const {return fBMin;}
34 0 : const Float_t* GetBoundMax() const {return fBMax;}
35 : //
36 0 : void SetXRowInv(const float* xi) {fRowXI = xi;}
37 0 : 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 : Bool_t IsInside(const float *par) const;
41 : //
42 : void Print(const Option_t* opt="") const;
43 : virtual void PrintSlice(int isl, const Option_t* opt) const = 0;
44 : //
45 : static float ChebEval1D(float x, const float* array, int ncf);
46 : static float ChebEval1D(float x, const Short_t* array, int ncf);
47 0 : static void SetDefPrecision(float prc=1e-4) {fgkDefPrec = prc>1e-8 ? prc:1e-8;}
48 0 : static float GetDefPrecision() {return fgkDefPrec;}
49 : //
50 : protected:
51 : void MapToInternal(int slice, const float* xy, float *xyint) const;
52 : void MapToInternal(int slice, const float* xy, float &x1, float &x2) const;
53 : float MapToExternal(int slice, float x,int dim) const;
54 : float* DefineGrid(int slice, int dim, const int np[2]) const;
55 : void CheckDimensions(const int *np) const;
56 : Int_t CalcChebCoefs(const float *funval,int np, float *outCoefs, float prec);
57 : //
58 : protected:
59 : //
60 : Int_t fDimOut; // each slice maps 2D to fDimOut values
61 : Int_t fNSlices; // number of slices in the stack
62 : Int_t fNParams; // number of parameterizations = fNSlices*fDimOut
63 : Int_t fNCoefsTot; // dimension of coeffs array for all slices
64 : Int_t fNRowsTot; // total number of 1D Cheb. param rows
65 : Float_t fBMin[2]; // min boundaries in each dimension
66 : Float_t fBMax[2]; // max boundaries in each dimension
67 : Float_t fBScaleZ; // scale for Z boundary mapping to [-1:1] interval
68 : Float_t fBOffsetZ; // offset for Z boundary mapping to [-1:1] interval
69 : Float_t fDead[2]; // dead zone in cm to account if X for each row is provided
70 : const Float_t* fRowXI; //! optional external!!! set 1/X for each row if dead zones to be accounted
71 : //
72 : UChar_t* fNRows; //[fNParams] N of used rows in the 2D coeffs matrix of each param
73 : UChar_t* fNCols; //[fNRowsTot] N of used columns in each row
74 : Int_t* fCoeffsEntry; //[fNSlices] start of the coeffs array in fCoeffs for each slice
75 : Int_t* fColEntry; //[fNSlices] start of the Ncolumns array in fNCols for each slice
76 : //
77 : static Float_t fgkDefPrec; // default precision
78 : static Float_t fWSpace[kMaxPoints]; // workspace
79 :
80 : //
81 : private:
82 : AliCheb2DStack(const AliCheb2DStack& src); // dummy
83 : AliCheb2DStack& operator=(const AliCheb2DStack& rhs); // dummy
84 : //
85 176 : ClassDef(AliCheb2DStack,2) // stack of 2D->fDimOut Chebyshev parameterization slices
86 : };
87 :
88 :
89 : //_________________________________________________________
90 : inline void AliCheb2DStack::MapToInternal(int slice, const float* xy, float args[2]) const
91 : {
92 : // map xy to [-1:1]
93 : // for Z we don't have dead zones
94 : args[kz] = (xy[kz]-fBOffsetZ)*fBScaleZ;
95 : float tmn = fBMin[ktgp], tmx = fBMax[ktgp];
96 : if (fRowXI) {
97 : tmn += fDead[0]*fRowXI[slice];
98 : tmx -= fDead[1]*fRowXI[slice];
99 : }
100 : args[ktgp] = 2.f*(xy[ktgp]-tmn)/(tmx-tmn) - 1.f;
101 : #ifdef _BRING_TO_BOUNDARY2D_
102 : if (args[kz]<-1.0f) args[kz]=-1.0f;
103 : else if (args[kz]> 1.0f) args[kz]=1.0f;
104 : if (args[ktgp]<-1.0f) args[ktgp]=-1.0f;
105 : else if (args[ktgp]> 1.0f) args[ktgp]=1.0f;
106 : #endif
107 : }
108 :
109 : //_________________________________________________________
110 : inline void AliCheb2DStack::MapToInternal(int slice, const float* xy, float &x0, float &x1) const
111 : {
112 : // map xy to [-1:1]
113 0 : x1 = (xy[kz]-fBOffsetZ)*fBScaleZ;
114 0 : float tmn = fBMin[ktgp], tmx = fBMax[ktgp];
115 0 : if (fRowXI) {
116 0 : tmn += fDead[0]*fRowXI[slice];
117 0 : tmx -= fDead[1]*fRowXI[slice];
118 0 : }
119 0 : x0 = 2.f*(xy[ktgp]-tmn)/(tmx-tmn) - 1.f;
120 : //
121 : #ifdef _BRING_TO_BOUNDARY2D_
122 0 : if (x0 < -1.0f) x0 =-1.0f; else if (x0 > 1.0f) x0 = 1.0f;
123 0 : if (x1 < -1.0f) x1 =-1.0f; else if (x1 > 1.0f) x1 = 1.0f;
124 : #endif
125 0 : }
126 :
127 : //__________________________________________________________________________________________
128 : inline float AliCheb2DStack::MapToExternal(int slice, float x,int dim) const
129 : {
130 : // map from [-1:1] to x for dim-th input dimension
131 0 : if (dim==kz) return x/fBScaleZ+fBOffsetZ;
132 : // for y/x need to account for dead zone
133 0 : float tmn = fBMin[ktgp], tmx = fBMax[ktgp];
134 0 : if (fRowXI) {
135 0 : tmn += fDead[0]*fRowXI[slice];
136 0 : tmx -= fDead[1]*fRowXI[slice];
137 0 : }
138 0 : return 0.5*(x+1.0f)*(tmx-tmn)+tmn;
139 0 : }
140 :
141 : //__________________________________________________________________________________________
142 : inline Bool_t AliCheb2DStack::IsInside(const float *par) const
143 : {
144 : // check if the point is inside of the fitted box
145 0 : for (int i=2;i--;) if (fBMin[i]>par[i] || par[i]>fBMax[i]) return kFALSE;
146 0 : return kTRUE;
147 0 : }
148 :
149 : //__________________________________________________________________________________________
150 : inline float AliCheb2DStack::ChebEval1D(float x, const float* array, int ncf)
151 : {
152 : // evaluate 1D Chebyshev parameterization. x is the argument mapped to [-1:1] interval
153 0 : if (!ncf) return 0;
154 0 : float b0(array[--ncf]), b1(0), b2(0), x2(x+x);
155 0 : for (int i=ncf;i--;) {
156 : b2 = b1;
157 : b1 = b0;
158 0 : b0 = array[i] + x2*b1 -b2;
159 : }
160 0 : return b0 - x*b1;
161 : //
162 0 : }
163 :
164 : //__________________________________________________________________________________________
165 : inline float AliCheb2DStack::ChebEval1D(float x, const Short_t* array, int ncf)
166 : {
167 : // evaluate 1D Chebyshev parameterization. x is the argument mapped to [-1:1] interval
168 0 : if (!ncf) return 0;
169 0 : float b0(array[--ncf]), b1(0), b2(0), x2(x+x);
170 0 : for (int i=ncf;i--;) {
171 : b2 = b1;
172 : b1 = b0;
173 0 : b0 = array[i] + x2*b1 -b2;
174 : }
175 0 : return b0 - x*b1;
176 : //
177 0 : }
178 :
179 :
180 :
181 : #endif
|
