source: branches/removing_cpp11_features/mhft/MHexagonFFT.cc@ 19608

Last change on this file since 19608 was 5691, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 15.7 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Christoph Kolodziejski, 11/2004 <mailto:>
19! Author(s): Thomas Bretz, 11/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2004
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28//////////////////////////////////////////////////////////////////////////////
29
30#include "MHexagonFFT.h"
31
32#include <TMath.h>
33
34#include "MLog.h"
35#include "MLogManip.h"
36
37ClassImp(MHexagonFFT);
38
39using namespace std;
40
41// ---------------------------------------------------------------------------
42//
43// Default Constructor
44// Initializes random number generator and default variables
45//
46MHexagonFFT::MHexagonFFT()
47{
48}
49
50/*
51void MHexagonFFT::Prepare(Int_t num, Float_t scale)
52{
53 //Int_t num = 34;
54 //Double_t scale = dist_y/dist_x;
55
56 Int_t cnt = 108*num*num;
57
58 psire.Set(num*num);
59 psiim.Set(num*num);
60
61// for(int j=0; j<num; j++)
62// {
63// for(int n=0; n<num; n++)
64// {
65// if (arr_k_m_id[j][n]<0)
66// continue;
67//
68// Double_t sumre=0;
69// Double_t sumim=0;
70
71 for(int k=0; k<num; k++)
72 {
73 for(int m=0; m<34-k; m++)
74 {
75 //Int_t new_ID=arr_k_m_id[k][m];
76 //if (new_ID<0)
77 // continue;
78
79 Double_t dx = 0.5*(m-k)/num;
80 Double_t dy = 0.5*(m+k)/num*scale;
81
82 dx *= TMath::TwoPi()/3;
83 dy *= TMath::TwoPi()/fgSqrt3;
84
85 const Double_t cos1 = cos(dy*(j+n));
86 const Double_t cos2 = cos(dy*j);
87 const Double_t cos3 = cos(dy*n);
88
89 //Alternatie nach Paper:
90 psire[m*num+k] = 0.5*(
91 +cos1*cos(dx*(j-n))
92 +cos2*cos(dx*(j+2*n))
93 +cos3*cos(dx*(2*j+n)));
94
95 psiim[m*num+k] = 0.5*(
96 +cos1*sin(dx*(j-n))
97 +cos2*sin(dx*(j+2*n))
98 -cos3*sin(dx*(2*j+n)));
99
100
101// psi_im *= i_inv;
102//
103// Double_t factor = (i_inv==1?1.:P_j_n(k,m));
104//
105// sumre += factor * (inre[new_ID]*psi_re - inim[new_ID]*psi_im);
106// sumim += factor * (inre[new_ID]*psi_im + inim[new_ID]*psi_re);
107}
108 }
109
110// Double_t factor = (i_inv==1?1.:P_j_n(j,n)/cnt);
111//
112// outre[arr_k_m_id[j][n]] = factor * sumre;
113// outim[arr_k_m_id[j][n]] = factor * sumim;
114// }
115// }
116}
117*/
118
119void MHexagonFFT::Prepare(Float_t scale, Int_t num)
120{
121 static const Double_t fgSqrt3 = TMath::Sqrt(3.);
122
123 fNum = num;
124
125 MArrayD fCosX(num*num*num*3);
126 MArrayD fCosY(num*num*num*3);
127 MArrayD fSin(num*num*num*3);
128
129 for(int j=0; j<3*num; j++)
130 {
131 for(int k=0; k<num; k++)
132 {
133 for(int m=0; m<num-k; m++)
134 {
135 Double_t dx = 0.5*(m-k)/num;
136 Double_t dy = 0.5*(m+k)/num*scale;
137
138 dx *= TMath::TwoPi()/3;
139 dy *= TMath::TwoPi()/fgSqrt3;
140
141 const Int_t idx = (m*num + k)*3*num;
142
143 fCosX[idx+j] = TMath::Cos(dx*j);
144 fCosY[idx+j] = TMath::Cos(dy*j);
145 fSin [idx+j] = TMath::Sin(dx*j);
146 }
147 }
148 }
149
150 //fPsiRe.Set(num*num*num*num);
151 //fPsiIm.Set(num*num*num*num);
152 fPsi.Set(num*num*num*num*2);
153
154 Int_t lim = num*(num+1)/2;
155
156 fM.Set(lim);
157 fK.Set(lim);
158 fP.Set(lim);
159 fIdx.Set(lim);
160
161 for(int j=0; j<num; j++)
162 {
163 for(int n=0; n<num-j; n++)
164 {
165 int idx0 = num-n-1;
166 int idx1 = idx0*(idx0+1)/2 + j;
167
168 fM[idx1]=n;
169 fK[idx1]=j;
170
171 fP[idx1]=P(j,n);
172
173 for(int k=0; k<num; k++)
174 {
175 for(int m=0; m<num-k; m++)
176 {
177 const Int_t idx = (m*num + k)*3*num;
178
179 const Float_t cos1 = fCosY[idx+j+n];
180 const Float_t cos2 = fCosY[idx+j];
181 const Float_t cos3 = fCosY[idx+n];
182
183 int idx2 = num-m-1;
184 int idx3 = idx2*(idx2+1)/2 + k;
185 const Int_t id1 = idx1*lim + idx3;
186
187 //fPsiRe[id1]
188 Double_t fPsiRe
189 = 2*(
190 +cos1*fCosX[idx+TMath::Abs(j-n)]
191 +cos2*fCosX[idx+j+2*n]
192 +cos3*fCosX[idx+2*j+n]);
193
194 //fPsiIm[id1] = 2*(
195 Double_t fPsiIm = 2*(
196 +cos1*fSin[idx+TMath::Abs(j-n)]*TMath::Sign(1, j-n)
197 +cos2*fSin[idx+j+2*n]
198 -cos3*fSin[idx+2*j+n]);
199
200 fPsi[id1*2] = fPsiRe;//fPsiRe[id1];
201 fPsi[id1*2+1] = fPsiIm;//fPsiIm[id1];
202 }
203 }
204
205 }
206 }
207
208 for (int idx1=0; idx1<lim; idx1++)
209 {
210 int n = fM[idx1];
211 int j = fK[idx1];
212
213 int idx0;
214 for (idx0=0; idx0<lim; idx0++)
215 if (fM[idx0]==j && fK[idx0]==n)
216 break;
217
218 fIdx[idx1]=idx0;
219 }
220
221}
222/*
223void MHexagonFFT::DGT4(const MArrayD &inre,
224 const MArrayD &inim,
225 MArrayD &outre,
226 MArrayD &outim,
227 Float_t scale,
228 Bool_t fwd)
229{
230 Int_t num = (Int_t)TMath::Sqrt((Float_t)inim.GetSize());
231 Int_t cnt = 108*num*num;
232 Int_t lim = num*(num+1)/2;
233
234 Float_t *endp = fP.GetArray()+lim;
235
236 for (int idx1=0; idx1<lim; idx1++)
237 {
238 if (fK[idx1]>fM[idx1] && fwd)
239 continue;
240
241 Double_t sumre=0;
242 Double_t sumim=0;
243
244 Float_t *psi = fPsi.GetArray() + idx1*lim*2;
245 Float_t *p = fP.GetArray();
246 Double_t *im = inim.GetArray();
247 Double_t *re = inre.GetArray();
248
249 while (p<endp)
250 {
251 const Float_t factor1 = (fwd?*p:1.);
252
253 const Float_t psire = *psi++;
254 const Float_t psiim = *psi++;
255
256 sumre += factor1 * (*re * psire - *im * psiim);
257 if (fwd)
258 sumim += factor1 * (*re * psiim + *im * psire);
259
260 im++;
261 re++;
262 p++;
263 }
264
265 const Double_t factor2 = fwd?fP[idx1]/cnt:1.;
266
267 outre[idx1] = factor2 * sumre;
268 outim[idx1] = -factor2 * sumim;
269 }
270
271 if (!fwd)
272 return;
273
274 for (int idx1=0; idx1<lim; idx1++)
275 {
276 if (fK[idx1]<fM[idx1])
277 continue;
278
279 outre[idx1] = outre[fIdx[idx1]];
280 outim[idx1] = -outim[fIdx[idx1]];
281 }
282 }*/
283
284void MHexagonFFT::TransformFastFWD(const MArrayD &inre,
285 MArrayD &outre,
286 MArrayD &outim) const
287{
288 const UInt_t num = fP.GetSize();
289
290 if (inre.GetSize()!=num)
291 {
292 cout << "ERROR - MHexagonFFT prepared for different size." << endl;
293 return;
294 }
295
296 outre.Set(num);
297 outim.Set(num);
298
299 const Int_t cnt = 108*fNum*fNum;
300
301 const Float_t *endp = fP.GetArray()+num;
302
303 for (UInt_t idx1=0; idx1<num; idx1++)
304 {
305 if (fK[idx1]>fM[idx1])
306 continue;
307
308 Double_t sumre=0;
309 Double_t sumim=0;
310
311 Float_t *psi = fPsi.GetArray() + idx1*num*2;
312 Float_t *p = fP.GetArray();
313 Double_t *re = inre.GetArray();
314
315 // 1st access to psi: const Float_t psire = *psi++;
316 // 2nd access to psi: const Float_t psiim = *psi++;
317 // sumre += f * *psire;
318 // sumim += f * *psiim;
319 while (p<endp)
320 {
321 const Double_t f = *p++ * *re++;
322
323 sumre += f * *psi++;
324 sumim += f * *psi++;
325 }
326
327 const Double_t factor2 = fP[idx1]/cnt;
328
329 outre[fIdx[idx1]] = outre[idx1] = factor2 * sumre;
330 outim[fIdx[idx1]] = -(outim[idx1] = -factor2 * sumim);
331
332 /*
333 outre[idx1] = factor2 * sumre;
334 outim[idx1] = -factor2 * sumim;
335
336 outre[fIdx[idx1]] = outre[idx1];
337 outim[fIdx[idx1]] = -outim[idx1];
338 */
339 }
340 /*
341 for (UInt_t idx1=0; idx1<num; idx1++)
342 {
343 if (fK[idx1]<fM[idx1])
344 continue;
345
346 outre[idx1] = outre[fIdx[idx1]];
347 outim[idx1] = -outim[fIdx[idx1]];
348 }
349 */
350}
351
352void MHexagonFFT::TransformFastBWD(const MArrayD &inre,
353 const MArrayD &inim,
354 MArrayD &outre) const
355{
356 const UInt_t num = fP.GetSize();
357
358 if (inre.GetSize()!=num)
359 {
360 cout << "ERROR - MHexagonFFT prepared for different size." << endl;
361 return;
362 }
363 if (inim.GetSize()!=num)
364 {
365 cout << "ERROR - MHexagonFFT prepared for different size." << endl;
366 return;
367 }
368 outre.Set(num);
369
370
371 const Double_t *endre = inre.GetArray()+num;
372
373 for (UInt_t idx1=0; idx1<num; idx1++)
374 {
375 Float_t *psi = fPsi.GetArray() + idx1*num*2;
376 Double_t *im = inim.GetArray();
377 Double_t *re = inre.GetArray();
378
379 Double_t sumre=0;
380 while (re<endre)
381 {
382 const Float_t psire = *psi++;
383 const Float_t psiim = *psi++;
384
385 sumre += *re++ * psire - *im++ * psiim;
386 }
387
388 outre[idx1] = sumre;
389 }
390}
391/*
392void MHexagonFFT::DGT3(const MArrayD &inre,
393 const MArrayD &inim,
394 MArrayD &outre,
395 MArrayD &outim,
396 Float_t scale,
397 Bool_t fwd)
398{
399 Int_t num = (Int_t)TMath::Sqrt((Float_t)inim.GetSize());
400 Int_t cnt = 108*num*num;
401
402 for(int j=0; j<num; j++)
403 {
404 for(int n=0; n<num-j; n++)
405 {
406 if (j-n>0 && fwd)
407 continue;
408
409 Double_t sumre=0;
410 Double_t sumim=0;
411
412 Int_t lim = num*(num+1)/2;
413 for (int idx0=0; idx0<lim; idx0++)
414
415// for(int k=0; k<num; k++)
416 {
417 int m = fM[idx0];
418 int k = fK[idx0];
419 // for(int m=0; m<num-k; m++)
420
421 {
422 const Int_t id = k*num + m;
423 const Int_t id1 = (((id*num)+n)*num+j)*2;
424
425 //Alternatie nach Paper:
426 const Float_t psire = fPsi[id1]; //fPsiRe[(id*num+n)*num+j];
427 const Float_t psiim = fPsi[id1+1]; //fPsiIm[(id*num+n)*num+j]*inv;
428
429 const Float_t factor1 = fwd==1?P(k,m):1.;
430
431 sumre += factor1 * (inre[id]*psire - inim[id]*psiim);
432 if (fwd)
433 sumim += factor1 * (inre[id]*psiim + inim[id]*psire);
434 }
435 }
436
437 const Double_t factor2 = fwd==1?P(j,n)/cnt:1.;
438
439 outre[j*num+n] = factor2 * sumre;
440 outim[j*num+n] = -factor2 * sumim;
441
442 if (!fwd)
443 continue;
444
445 outre[n*num+j] = factor2 * sumre;
446 outim[n*num+j] = factor2 * sumim;
447 }
448 }
449}
450
451void MHexagonFFT::DGT2(const MArrayD &inre,
452 const MArrayD &inim,
453 MArrayD &outre,
454 MArrayD &outim,
455 Float_t scale,
456 Bool_t fwd)
457{
458 Int_t num = (Int_t)TMath::Sqrt((Float_t)inim.GetSize());
459 Int_t cnt = 108*num*num;
460 Int_t inv = fwd?-1:1;
461
462 for(int j=0; j<num; j++)
463 {
464 for(int n=0; n<num-j; n++)
465 {
466 if (j-n>0 && fwd)
467 continue;
468
469 Double_t sumre=0;
470 Double_t sumim=0;
471
472 for(int k=0; k<num; k++)
473 {
474 for(int m=0; m<num-k; m++)
475 {
476 const Int_t idx = (m*num + k)*3*num;
477
478 const Float_t cos1 = fCosY[idx+j+n];
479 const Float_t cos2 = fCosY[idx+j];
480 const Float_t cos3 = fCosY[idx+n];
481
482 //Alternatie nach Paper:
483 const Float_t psire = 2*(
484 +cos1*fCosX[idx+TMath::Abs(j-n)]
485 +cos2*fCosX[idx+j+2*n]
486 +cos3*fCosX[idx+2*j+n]);
487
488 const Float_t psiim = 2*inv*(
489 +cos1*fSin[idx+TMath::Abs(j-n)]*TMath::Sign(1, j-n)
490 +cos2*fSin[idx+j+2*n]
491 -cos3*fSin[idx+2*j+n]);
492
493 const Float_t factor = (fwd==1?P(k,m):1.);
494
495 sumre += factor * (inre[k*num+m]*psire - inim[k*num+m]*psiim);
496 sumim += factor * (inre[k*num+m]*psiim + inim[k*num+m]*psire);
497 }
498 }
499
500 const Double_t factor = (fwd==1?P(j,n)/cnt:1.);
501
502 outre[j*num+n] = factor * sumre;
503 outim[j*num+n] = factor * sumim;
504
505 if (fwd)
506 {
507 outre[n*num+j] = factor * sumre;
508 outim[n*num+j] = -factor * sumim;
509 }
510 }
511 }
512}
513*/
514void MHexagonFFT::TransformSlow(const MArrayD &inre, const MArrayD &inim,
515 MArrayD &outre, MArrayD &outim,
516 Float_t scale, Bool_t fwd)
517{
518 static const Double_t fgSqrt3 = TMath::Sqrt(3.);
519
520 Int_t num = (Int_t)TMath::Sqrt((Float_t)inim.GetSize());
521 Int_t cnt = 108*num*num;
522 Int_t inv = fwd?-1:1;
523
524 for(int j=0; j<num; j++)
525 {
526 for(int n=0; n<num-j; n++)
527 {
528 if (j-n>0 && fwd)
529 continue;
530
531 Double_t sumre=0;
532 Double_t sumim=0;
533
534 for(int k=0; k<num; k++)
535 {
536 for(int m=0; m<num-k; m++)
537 {
538 Double_t dx = 0.5*(m-k)/num;
539 Double_t dy = 0.5*(m+k)/num*scale;
540
541 dx *= TMath::TwoPi()/3;
542 dy *= TMath::TwoPi()/fgSqrt3;
543
544 const Double_t cos1 = TMath::Cos(dy*(j+n));
545 const Double_t cos2 = TMath::Cos(dy*j);
546 const Double_t cos3 = TMath::Cos(dy*n);
547
548 //Alternatie nach Paper:
549 const Double_t psire = 2*(
550 +cos1*TMath::Cos(dx*(j-n))
551 +cos2*TMath::Cos(dx*(j+2*n))
552 +cos3*TMath::Cos(dx*(2*j+n)));
553
554 const Double_t psiim = 2*inv*(
555 +cos1*TMath::Sin(dx*(j-n))
556 +cos2*TMath::Sin(dx*(j+2*n))
557 -cos3*TMath::Sin(dx*(2*j+n)));
558
559 const Double_t factor = (fwd==1?P(k,m):1.);
560
561 sumre += factor * (inre[k*num+m]*psire - inim[k*num+m]*psiim);
562 sumim += factor * (inre[k*num+m]*psiim + inim[k*num+m]*psire);
563 }
564 }
565
566 const Double_t factor = (fwd==1?P(j,n)/cnt:1.);
567
568 outre[j*num+n] = factor * sumre;
569 outim[j*num+n] = factor * sumim;
570
571 if (fwd)
572 {
573 outre[n*num+j] = factor * sumre;
574 outim[n*num+j] = -factor * sumim;
575 }
576 }
577 }
578}
Note: See TracBrowser for help on using the repository browser.