Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : #include <cstdlib>
17 : #include <TSystem.h>
18 : #include "AliCheb3DCalc.h"
19 : #include "AliLog.h"
20 :
21 176 : ClassImp(AliCheb3DCalc)
22 :
23 : //__________________________________________________________________________________________
24 68760 : AliCheb3DCalc::AliCheb3DCalc() :
25 68760 : fNCoefs(0),
26 68760 : fNRows(0),
27 68760 : fNCols(0),
28 68760 : fNElemBound2D(0),
29 68760 : fNColsAtRow(0),
30 68760 : fColAtRowBg(0),
31 68760 : fCoefBound2D0(0),
32 68760 : fCoefBound2D1(0),
33 68760 : fCoefs(0),
34 68760 : fTmpCf1(0),
35 68760 : fTmpCf0(0),
36 68760 : fPrec(0)
37 343800 : {
38 : // default constructor
39 137520 : }
40 :
41 : //__________________________________________________________________________________________
42 : AliCheb3DCalc::AliCheb3DCalc(const AliCheb3DCalc& src) :
43 0 : TNamed(src),
44 0 : fNCoefs(src.fNCoefs),
45 0 : fNRows(src.fNRows),
46 0 : fNCols(src.fNCols),
47 0 : fNElemBound2D(src.fNElemBound2D),
48 0 : fNColsAtRow(0),
49 0 : fColAtRowBg(0),
50 0 : fCoefBound2D0(0),
51 0 : fCoefBound2D1(0),
52 0 : fCoefs(0),
53 0 : fTmpCf1(0),
54 0 : fTmpCf0(0),
55 0 : fPrec(src.fPrec)
56 0 : {
57 : // copy constructor
58 : //
59 0 : if (src.fNColsAtRow) {
60 0 : fNColsAtRow = new UShort_t[fNRows];
61 0 : for (int i=fNRows;i--;) fNColsAtRow[i] = src.fNColsAtRow[i];
62 0 : }
63 0 : if (src.fColAtRowBg) {
64 0 : fColAtRowBg = new UShort_t[fNRows];
65 0 : for (int i=fNRows;i--;) fColAtRowBg[i] = src.fColAtRowBg[i];
66 0 : }
67 0 : if (src.fCoefBound2D0) {
68 0 : fCoefBound2D0 = new UShort_t[fNElemBound2D];
69 0 : for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = src.fCoefBound2D0[i];
70 0 : }
71 0 : if (src.fCoefBound2D1) {
72 0 : fCoefBound2D1 = new UShort_t[fNElemBound2D];
73 0 : for (int i=fNElemBound2D;i--;) fCoefBound2D1[i] = src.fCoefBound2D1[i];
74 0 : }
75 0 : if (src.fCoefs) {
76 0 : fCoefs = new Float_t[fNCoefs];
77 0 : for (int i=fNCoefs;i--;) fCoefs[i] = src.fCoefs[i];
78 0 : }
79 0 : if (src.fTmpCf1) fTmpCf1 = new Float_t[fNCols];
80 0 : if (src.fTmpCf0) fTmpCf0 = new Float_t[fNRows];
81 0 : }
82 :
83 : //__________________________________________________________________________________________
84 0 : AliCheb3DCalc::AliCheb3DCalc(FILE* stream) :
85 0 : fNCoefs(0),
86 0 : fNRows(0),
87 0 : fNCols(0),
88 0 : fNElemBound2D(0),
89 0 : fNColsAtRow(0),
90 0 : fColAtRowBg(0),
91 0 : fCoefBound2D0(0),
92 0 : fCoefBound2D1(0),
93 0 : fCoefs(0),
94 0 : fTmpCf1(0),
95 0 : fTmpCf0(0),
96 0 : fPrec(0)
97 0 : {
98 : // constructor from coeffs. streem
99 0 : LoadData(stream);
100 0 : }
101 :
102 : //__________________________________________________________________________________________
103 : AliCheb3DCalc& AliCheb3DCalc::operator=(const AliCheb3DCalc& rhs)
104 : {
105 : // assignment operator
106 0 : if (this != &rhs) {
107 0 : Clear();
108 0 : SetName(rhs.GetName());
109 0 : SetTitle(rhs.GetTitle());
110 0 : fNCoefs = rhs.fNCoefs;
111 0 : fNRows = rhs.fNRows;
112 0 : fNCols = rhs.fNCols;
113 0 : fPrec = rhs.fPrec;
114 0 : if (rhs.fNColsAtRow) {
115 0 : fNColsAtRow = new UShort_t[fNRows];
116 0 : for (int i=fNRows;i--;) fNColsAtRow[i] = rhs.fNColsAtRow[i];
117 0 : }
118 0 : if (rhs.fColAtRowBg) {
119 0 : fColAtRowBg = new UShort_t[fNRows];
120 0 : for (int i=fNRows;i--;) fColAtRowBg[i] = rhs.fColAtRowBg[i];
121 0 : }
122 0 : if (rhs.fCoefBound2D0) {
123 0 : fCoefBound2D0 = new UShort_t[fNElemBound2D];
124 0 : for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = rhs.fCoefBound2D0[i];
125 0 : }
126 0 : if (rhs.fCoefBound2D1) {
127 0 : fCoefBound2D1 = new UShort_t[fNElemBound2D];
128 0 : for (int i=fNElemBound2D;i--;) fCoefBound2D1[i] = rhs.fCoefBound2D1[i];
129 0 : }
130 0 : if (rhs.fCoefs) {
131 0 : fCoefs = new Float_t[fNCoefs];
132 0 : for (int i=fNCoefs;i--;) fCoefs[i] = rhs.fCoefs[i];
133 0 : }
134 0 : if (rhs.fTmpCf1) fTmpCf1 = new Float_t[fNCols];
135 0 : if (rhs.fTmpCf0) fTmpCf0 = new Float_t[fNRows];
136 : }
137 0 : return *this;
138 : }
139 :
140 : //__________________________________________________________________________________________
141 : void AliCheb3DCalc::Clear(const Option_t*)
142 : {
143 : // delete all dynamycally allocated structures
144 343560 : if (fTmpCf1) { delete[] fTmpCf1; fTmpCf1 = 0;}
145 274800 : if (fTmpCf0) { delete[] fTmpCf0; fTmpCf0 = 0;}
146 274800 : if (fCoefs) { delete[] fCoefs; fCoefs = 0;}
147 274800 : if (fCoefBound2D0) { delete[] fCoefBound2D0; fCoefBound2D0 = 0; }
148 274800 : if (fCoefBound2D1) { delete[] fCoefBound2D1; fCoefBound2D1 = 0; }
149 274800 : if (fNColsAtRow) { delete[] fNColsAtRow; fNColsAtRow = 0; }
150 274800 : if (fColAtRowBg) { delete[] fColAtRowBg; fColAtRowBg = 0; }
151 : //
152 68760 : }
153 :
154 : //__________________________________________________________________________________________
155 : void AliCheb3DCalc::Print(const Option_t* ) const
156 : {
157 : // print info
158 0 : printf("Chebyshev parameterization data %s for 3D->1 function. Prec:%e\n",GetName(),fPrec);
159 : int nmax3d = 0;
160 0 : for (int i=fNElemBound2D;i--;) if (fCoefBound2D0[i]>nmax3d) nmax3d = fCoefBound2D0[i];
161 0 : printf("%d coefficients in %dx%dx%d matrix\n",fNCoefs,fNRows,fNCols,nmax3d);
162 : //
163 0 : }
164 :
165 : //__________________________________________________________________________________________
166 : Float_t AliCheb3DCalc::EvalDeriv(int dim, const Float_t *par) const
167 : {
168 : // evaluate Chebyshev parameterization derivative in given dimension for 3D function.
169 : // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
170 : //
171 : int ncfRC;
172 0 : for (int id0=fNRows;id0--;) {
173 0 : int nCLoc = fNColsAtRow[id0]; // number of significant coefs on this row
174 0 : if (!nCLoc) {fTmpCf0[id0]=0; continue;}
175 : //
176 0 : int col0 = fColAtRowBg[id0]; // beginning of local column in the 2D boundary matrix
177 0 : for (int id1=nCLoc;id1--;) {
178 0 : int id = id1+col0;
179 0 : if (!(ncfRC=fCoefBound2D0[id])) { fTmpCf1[id1]=0; continue;}
180 0 : if (dim==2) fTmpCf1[id1] = ChebEval1Deriv(par[2],fCoefs + fCoefBound2D1[id], ncfRC);
181 0 : else fTmpCf1[id1] = ChebEval1D(par[2],fCoefs + fCoefBound2D1[id], ncfRC);
182 0 : }
183 0 : if (dim==1) fTmpCf0[id0] = ChebEval1Deriv(par[1],fTmpCf1,nCLoc);
184 0 : else fTmpCf0[id0] = ChebEval1D(par[1],fTmpCf1,nCLoc);
185 0 : }
186 0 : return (dim==0) ? ChebEval1Deriv(par[0],fTmpCf0,fNRows) : ChebEval1D(par[0],fTmpCf0,fNRows);
187 : //
188 : }
189 :
190 : //__________________________________________________________________________________________
191 : Float_t AliCheb3DCalc::EvalDeriv2(int dim1,int dim2, const Float_t *par) const
192 : {
193 : // evaluate Chebyshev parameterization 2n derivative in given dimensions for 3D function.
194 : // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
195 : //
196 0 : Bool_t same = dim1==dim2;
197 : int ncfRC;
198 0 : for (int id0=fNRows;id0--;) {
199 0 : int nCLoc = fNColsAtRow[id0]; // number of significant coefs on this row
200 0 : if (!nCLoc) {fTmpCf0[id0]=0; continue;}
201 : //
202 0 : int col0 = fColAtRowBg[id0]; // beginning of local column in the 2D boundary matrix
203 0 : for (int id1=nCLoc;id1--;) {
204 0 : int id = id1+col0;
205 0 : if (!(ncfRC=fCoefBound2D0[id])) { fTmpCf1[id1]=0; continue;}
206 0 : if (dim1==2||dim2==2) fTmpCf1[id1] = same ? ChebEval1Deriv2(par[2],fCoefs + fCoefBound2D1[id], ncfRC)
207 0 : : ChebEval1Deriv(par[2],fCoefs + fCoefBound2D1[id], ncfRC);
208 0 : else fTmpCf1[id1] = ChebEval1D(par[2],fCoefs + fCoefBound2D1[id], ncfRC);
209 0 : }
210 0 : if (dim1==1||dim2==1) fTmpCf0[id0] = same ? ChebEval1Deriv2(par[1],fTmpCf1,nCLoc):ChebEval1Deriv(par[1],fTmpCf1,nCLoc);
211 0 : else fTmpCf0[id0] = ChebEval1D(par[1],fTmpCf1,nCLoc);
212 0 : }
213 0 : return (dim1==0||dim2==0) ? (same ? ChebEval1Deriv2(par[0],fTmpCf0,fNRows):ChebEval1Deriv(par[0],fTmpCf0,fNRows)) :
214 0 : ChebEval1D(par[0],fTmpCf0,fNRows);
215 : //
216 : }
217 :
218 : //_______________________________________________
219 : #ifdef _INC_CREATION_ALICHEB3D_
220 : void AliCheb3DCalc::SaveData(const char* outfile,Bool_t append) const
221 : {
222 : // writes coefficients data to output text file, optionallt appending on the end of existing file
223 0 : TString strf = outfile;
224 0 : gSystem->ExpandPathName(strf);
225 0 : FILE* stream = fopen(strf,append ? "a":"w");
226 0 : SaveData(stream);
227 0 : fclose(stream);
228 : //
229 0 : }
230 : #endif
231 :
232 : //_______________________________________________
233 : #ifdef _INC_CREATION_ALICHEB3D_
234 : void AliCheb3DCalc::SaveData(FILE* stream) const
235 : {
236 : // writes coefficients data to existing output stream
237 : // Note: fNCols, fNElemBound2D and fColAtRowBg is not stored, will be computed on fly during the loading of this file
238 0 : fprintf(stream,"#\nSTART %s\n",GetName());
239 0 : fprintf(stream,"# Number of rows\n%d\n",fNRows);
240 : //
241 0 : fprintf(stream,"# Number of columns per row\n");
242 0 : for (int i=0;i<fNRows;i++) fprintf(stream,"%d\n",fNColsAtRow[i]);
243 : //
244 0 : fprintf(stream,"# Number of Coefs in each significant block of third dimension\n");
245 0 : for (int i=0;i<fNElemBound2D;i++) fprintf(stream,"%d\n",fCoefBound2D0[i]);
246 : //
247 0 : fprintf(stream,"# Coefficients\n");
248 0 : for (int i=0;i<fNCoefs;i++) fprintf(stream,"%+.8e\n",fCoefs[i]);
249 : //
250 0 : fprintf(stream,"# Precision\n");
251 0 : fprintf(stream,"%+.8e\n",fPrec);
252 : //
253 0 : fprintf(stream,"END %s\n",GetName());
254 : //
255 0 : }
256 : #endif
257 :
258 : //_______________________________________________
259 : void AliCheb3DCalc::LoadData(FILE* stream)
260 : {
261 : // Load coefs. from the stream
262 0 : if (!stream) AliFatal("No stream provided");
263 0 : TString buffs;
264 0 : Clear();
265 0 : ReadLine(buffs,stream);
266 0 : if (!buffs.BeginsWith("START")) AliFatalF("Expected: \"START <fit_name>\", found \"%s\"",buffs.Data());
267 0 : if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
268 : //
269 0 : ReadLine(buffs,stream); // NRows
270 0 : fNRows = buffs.Atoi();
271 0 : if (fNRows<0 && !buffs.IsDigit()) AliFatalF("Expected: '<number_of_rows>', found \"%s\"",buffs.Data());
272 : //
273 0 : fNCols = 0;
274 0 : fNElemBound2D = 0;
275 0 : InitRows(fNRows);
276 : //
277 0 : for (int id0=0;id0<fNRows;id0++) {
278 0 : ReadLine(buffs,stream); // n.cols at this row
279 0 : fNColsAtRow[id0] = buffs.Atoi();
280 0 : fColAtRowBg[id0] = fNElemBound2D; // begining of this row in 2D boundary surface
281 0 : fNElemBound2D += fNColsAtRow[id0];
282 0 : if (fNCols<fNColsAtRow[id0]) fNCols = fNColsAtRow[id0];
283 : }
284 0 : InitCols(fNCols);
285 : //
286 0 : fNCoefs = 0;
287 0 : InitElemBound2D(fNElemBound2D);
288 : //
289 0 : for (int i=0;i<fNElemBound2D;i++) {
290 0 : ReadLine(buffs,stream); // n.coeffs at 3-d dimension for the given column/row
291 0 : fCoefBound2D0[i] = buffs.Atoi();
292 0 : fCoefBound2D1[i] = fNCoefs;
293 0 : fNCoefs += fCoefBound2D0[i];
294 : }
295 : //
296 0 : InitCoefs(fNCoefs);
297 0 : for (int i=0;i<fNCoefs;i++) {
298 0 : ReadLine(buffs,stream);
299 0 : fCoefs[i] = buffs.Atof();
300 : }
301 : //
302 : // read precision
303 0 : ReadLine(buffs,stream);
304 0 : fPrec = buffs.Atof();
305 : //
306 : // check end_of_data record
307 0 : ReadLine(buffs,stream);
308 0 : if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) AliFatalF("Expected \"END %s\", found \"%s\"",GetName(),buffs.Data());
309 : //
310 0 : }
311 :
312 : //_______________________________________________
313 : void AliCheb3DCalc::ReadLine(TString& str,FILE* stream)
314 : {
315 : // read single line from the stream, skipping empty and commented lines. EOF is not expected
316 0 : while (str.Gets(stream)) {
317 0 : str = str.Strip(TString::kBoth,' ');
318 0 : if (str.IsNull()||str.BeginsWith("#")) continue;
319 : return;
320 : }
321 0 : AliFatalGeneral("ReadLine","Failed to read from stream"); // normally, should not reach here
322 0 : }
323 :
324 : //_______________________________________________
325 : void AliCheb3DCalc::InitCols(int nc)
326 : {
327 : // Set max.number of significant columns in the coefs matrix
328 0 : fNCols = nc;
329 0 : if (fTmpCf1) {delete[] fTmpCf1; fTmpCf1 = 0;}
330 0 : if (fNCols>0) fTmpCf1 = new Float_t [fNCols];
331 0 : }
332 :
333 : //_______________________________________________
334 : void AliCheb3DCalc::InitRows(int nr)
335 : {
336 : // Set max.number of significant rows in the coefs matrix
337 0 : if (fNColsAtRow) {delete[] fNColsAtRow; fNColsAtRow = 0;}
338 0 : if (fColAtRowBg) {delete[] fColAtRowBg; fColAtRowBg = 0;}
339 0 : if (fTmpCf0) {delete[] fTmpCf0; fTmpCf0 = 0;}
340 0 : fNRows = nr;
341 0 : if (fNRows>0) {
342 0 : fNColsAtRow = new UShort_t[fNRows];
343 0 : fTmpCf0 = new Float_t [fNRows];
344 0 : fColAtRowBg = new UShort_t[fNRows];
345 0 : for (int i=fNRows;i--;) fNColsAtRow[i] = fColAtRowBg[i] = 0;
346 0 : }
347 0 : }
348 :
349 : //_______________________________________________
350 : void AliCheb3DCalc::InitElemBound2D(int ne)
351 : {
352 : // Set max number of significant coefs for given row/column of coefs 3D matrix
353 0 : if (fCoefBound2D0) {delete[] fCoefBound2D0; fCoefBound2D0 = 0;}
354 0 : if (fCoefBound2D1) {delete[] fCoefBound2D1; fCoefBound2D1 = 0;}
355 0 : fNElemBound2D = ne;
356 0 : if (fNElemBound2D>0) {
357 0 : fCoefBound2D0 = new UShort_t[fNElemBound2D];
358 0 : fCoefBound2D1 = new UShort_t[fNElemBound2D];
359 0 : for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = fCoefBound2D1[i] = 0;
360 0 : }
361 0 : }
362 :
363 : //_______________________________________________
364 : void AliCheb3DCalc::InitCoefs(int nc)
365 : {
366 : // Set total number of significant coefs
367 0 : if (fCoefs) {delete[] fCoefs; fCoefs = 0;}
368 0 : fNCoefs = nc;
369 0 : if (fNCoefs>0) {
370 0 : fCoefs = new Float_t [fNCoefs];
371 0 : for (int i=fNCoefs;i--;) fCoefs[i] = 0.0;
372 0 : }
373 0 : }
374 :
375 : //__________________________________________________________________________________________
376 : Float_t AliCheb3DCalc::ChebEval1Deriv(Float_t x, const Float_t * array, int ncf )
377 : {
378 : // evaluate 1D Chebyshev parameterization's derivative. x is the argument mapped to [-1:1] interval
379 0 : if (--ncf<1) return 0;
380 : Float_t b0, b1, b2;
381 0 : Float_t x2 = x+x;
382 : b1 = b2 = 0;
383 : float dcf0=0,dcf1,dcf2=0;
384 0 : b0 = dcf1 = 2*ncf*array[ncf];
385 0 : if (!(--ncf)) return b0/2;
386 : //
387 0 : for (int i=ncf;i--;) {
388 : b2 = b1;
389 : b1 = b0;
390 0 : dcf0 = dcf2 + 2*(i+1)*array[i+1];
391 0 : b0 = dcf0 + x2*b1 -b2;
392 : dcf2 = dcf1;
393 : dcf1 = dcf0;
394 : }
395 : //
396 0 : return b0 - x*b1 - dcf0/2;
397 0 : }
398 :
399 : //__________________________________________________________________________________________
400 : Float_t AliCheb3DCalc::ChebEval1Deriv2(Float_t x, const Float_t * array, int ncf )
401 : {
402 : // evaluate 1D Chebyshev parameterization's 2nd derivative. x is the argument mapped to [-1:1] interval
403 0 : if (--ncf<2) return 0;
404 : Float_t b0, b1, b2;
405 0 : Float_t x2 = x+x;
406 : b1 = b2 = 0;
407 : float dcf0=0,dcf1=0,dcf2=0;
408 : float ddcf0=0,ddcf1,ddcf2=0;
409 : //
410 0 : dcf2 = 2*ncf*array[ncf];
411 0 : --ncf;
412 :
413 0 : dcf1 = 2*ncf*array[ncf];
414 0 : b0 = ddcf1 = 2*ncf*dcf2;
415 : //
416 0 : if (!(--ncf)) return b0/2;
417 : //
418 0 : for (int i=ncf;i--;) {
419 : b2 = b1;
420 : b1 = b0;
421 0 : dcf0 = dcf2 + 2*(i+1)*array[i+1];
422 0 : ddcf0 = ddcf2 + 2*(i+1)*dcf1;
423 0 : b0 = ddcf0 + x2*b1 -b2;
424 : //
425 : ddcf2 = ddcf1;
426 : ddcf1 = ddcf0;
427 : //
428 : dcf2 = dcf1;
429 : dcf1 = dcf0;
430 : //
431 : }
432 : //
433 0 : return b0 - x*b1 - ddcf0/2;
434 0 : }
435 :
436 : //__________________________________________________________________________________________
437 : Int_t AliCheb3DCalc::GetMaxColsAtRow() const
438 : {
439 : int nmax3d = 0;
440 0 : for (int i=fNElemBound2D;i--;) if (fCoefBound2D0[i]>nmax3d) nmax3d = fCoefBound2D0[i];
441 0 : return nmax3d;
442 : }
|
