/
ThreeDeePolySpace.cxx
276 lines (208 loc) · 8.02 KB
/
ThreeDeePolySpace.cxx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
/* -------------------------------------------------------------------------- */
/* ThreeDeePolySpace.cc */
/* */
/* Collection of fitting routines with orthogonal polynomials. */
/* */
/* A.Kisselev, PNPI, St.Petersburg, Russia. */
/* e-mail: kisselev@hermes.desy.de */
/* -------------------------------------------------------------------------- */
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cassert>
#include <ayk.h>
#include <ThreeDeePolySpace.h>
/* ========================================================================== */
ThreeDeePolySpace::ThreeDeePolySpace(int max_power[3], int _point_num, unsigned char _parity[3])
{
int odd, even;
// Check later, whether this is really needed;
memset(this, 0x00, sizeof(ThreeDeePolySpace));
// Start off "1"; the actual calculation is below;
dim = 1;
// Figure out overall dimension;
for(int ik=0; ik<3; ik++)
{
_max_power[ik] = max_power[ik];
parity[ik] = _parity ? _parity[ik] : _BOTH_;
// Figure out actual number of degrees of freedom;
degrees_of_freedom[ik] = 0;
odd = even = max_power[ik]/2;
if (max_power[ik]%2) even++;
if (parity[ik] & _ODD_) degrees_of_freedom[ik] += odd;
if (parity[ik] & _EVEN_) degrees_of_freedom[ik] += even;
dim *= degrees_of_freedom[ik];
} /*for*/
// Allocate memory for basis and derivatives;
basis = new ThreeDeePolynomial* [dim];
for(int ik=0; ik<dim; ik++)
basis[ik] = new ThreeDeePolynomial(this);
dbasis = new ThreeDeePolynomial** [dim];
for(int ii=0; ii<dim; ii++)
{
dbasis[ii] = new ThreeDeePolynomial* [3];
for(int ik=0; ik<3; ik++)
dbasis[ii][ik] = new ThreeDeePolynomial(this);
} /*for*/
// Working space;
buffer = new ThreeDeePolynomial(this);
// Allocate 3D points;
point_num = _point_num;
points = new ThreeDeePolyPoint [point_num];
// A reasonable default for weight, please;
for(int ik=0; ik<point_num; ik++)
points[ik].weight = 1.;
} /* ThreeDeePolySpace::ThreeDeePolySpace */
/* ========================================================================== */
/* Does not take into account possibility that some of the directions may */
/* have too few points; fix later; */
int ThreeDeePolySpace::getNDF()
{
return point_num - dim;
} /* ThreeDeePolySpace::getNDF */
/* ========================================================================== */
//
// -> do not bother to create operators; just convert functions to methods;
// since at least few of these are logically related to ThreeDeePolySpace rather
// than ThreeDeePolynomial class, put them all into ThreeDeePolySpace;
//
double ThreeDeePolySpace::polyProduct(ThreeDeePolynomial *p1, ThreeDeePolynomial *p2)
{
double ret = 0.0;
for(int pt=0; pt<point_num; pt++)
{
ThreeDeePolyPoint *point = points + pt;
assert(0);
//@@@ ret += point->weight*p1->value(point->xx)*p2->value(point->xx);
} /*for*/
return ret;
} /* ThreeDeePolySpace::polyProduct */
/* -------------------------------------------------------------------------- */
double ThreeDeePolySpace::polyProjection(ThreeDeePolynomial *poly)
{
double ret = 0.0;
if (poly->off) return 0.0;
for(int pt=0; pt<point_num; pt++)
{
ThreeDeePolyPoint *point = points + pt;
assert(0);
//@@@ ret += point->weight*poly->value(point->xx)*point->f;
} /*for*/
return ret;
} /* ThreeDeePolySpace::polyProjection */
/* -------------------------------------------------------------------------- */
void ThreeDeePolySpace::polyCopy(ThreeDeePolynomial *dest, ThreeDeePolynomial *source)
{
int ix, iy, iz;
for(ix=0; ix<_max_power[_X_]; ix++)
{
if (!PARITY_CHECK(this, _X_, ix)) continue;
for(iy=0; iy<_max_power[_Y_]; iy++)
{
if (!PARITY_CHECK(this, _Y_, iy)) continue;
for(iz=0; iz<_max_power[_Z_]; iz++)
{
if (!PARITY_CHECK(this, _Z_, iz)) continue;
dest->cff[ix][iy][iz] = source->cff[ix][iy][iz];
} /*for*/
} /*for*/
} /*for*/
} /* ThreeDeePolySpace::polyCopy */
/* -------------------------------------------------------------------------- */
void ThreeDeePolySpace::buildOrthogonalPolynomials()
{
// Do I actually need this "off_counter"?;
int id = 0, off_counter = 0;
// Start with 1, x, x^2, ...;
for(int ix=0; ix<_max_power[_X_]; ix++)
for(int iy=0; iy<_max_power[_Y_]; iy++)
for(int iz=0; iz<_max_power[_Z_]; iz++)
{
if (!PARITY_CHECK(this, _X_, ix) ||
!PARITY_CHECK(this, _Y_, iy) ||
!PARITY_CHECK(this, _Z_, iz))
continue;
{
ThreeDeePolynomial *poly = basis[id++];
for(int qx=0; qx<_max_power[_X_]; qx++)
for(int qy=0; qy<_max_power[_Y_]; qy++)
for(int qz=0; qz<_max_power[_Z_]; qz++)
poly->cff[qx][qy][qz] = 0.;
poly->cff[ix][iy][iz] = 1.;
}
} /*for..for*/
for(int ii=0; ii<dim; ii++)
{
for(int ik=0; ik<ii; ik++)
{
double projection = polyProduct(basis[ii], basis[ik]);
polyCopy(buffer, basis[ik]);
buffer->multiply(-projection);
basis[ii]->increment(buffer);
} /*for*/
if (basis[ii]->normalize()) off_counter++;
} /*for*/
} /* ThreeDeePolySpace::buildOrthogonalPolynomials */
/* ========================================================================== */
int ThreeDeePolySpace::buildBasisGradients()
{
for(int ii=0; ii<dim; ii++)
if (basis[ii]->calculateGradient(dbasis[ii]))
return -1;
return 0;
} /* ThreeDeePolySpace::buildBasisGradients */
/* ========================================================================== */
//
// -> yes, at least here I should use ThreeDeePolySpace-based routine (not to
// blindly reuse fit->space pointer);
//
void ThreeDeePolySpace::calculateFittingPolynomial(ThreeDeePolynomial *fit)
{
for(int ii=0; ii<dim; ii++)
fit->linear[ii] = polyProjection(basis[ii]);
fit->convertLinearToCff();
} /* ThreeDeePolySpace::calculateFittingPolynomial */
/* ========================================================================== */
//
// -> yes, as any point-coupled routine this one should be ThreeDeePolySpace-related;
//
double ThreeDeePolySpace::getPolyFitChiSquare(ThreeDeePolynomial *fit)
{
double ret = 0.0;
for(int pt=0; pt<point_num; pt++)
{
ThreeDeePolyPoint *point = points + pt;
assert(0);
//@@@ ret += point->weight*SQR(fit->value(point->xx) - point->f);
} /*for*/
return ret;
} /* ThreeDeePolySpace::getPolyFitChiSquare */
/* ========================================================================== */
/* This calculation assumes that chi^2 is correct; in this case errors on */
/* 'linear' coefficients are all equal '1.' and below calculation is Ok; */
double ThreeDeePolySpace::getNaivePolyFitError(double xx[3])
{
double err = 0.0;
for(int ii=0; ii<dim; ii++)
err += SQR(basis[ii]->value(xx));
return sqrt(err);
} /* ThreeDeePolySpace::getNaivePolyFitError */
/* ========================================================================== */
/* This calculation takes chi^2 into account and renormalizes naive error */
/* estimate by sqrt(chi^2/ndf); this way error estimate will not depend on the*/
/* weight scaling factor; so for instance if errors in all points are more or */
/* less the same, they can be safely set to 1. and this function will provide */
/* best guess error estimate; */
//
// -> not sure this has ever been checked;
//
double ThreeDeePolySpace::getRealisticPolyFitError(ThreeDeePolynomial *fit, double xx[3])
{
int ndf = getNDF();
double err = 0.0, chi2 = getPolyFitChiSquare(fit);
for(int ii=0; ii<dim; ii++)
err += SQR(basis[ii]->value(xx));
return sqrt(err*chi2/ndf);
} /* ThreeDeePolySpace::getRealisticPolyFitError */
/* ========================================================================== */