1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.2 2007/10/12 13:36:27 cvetan
19 Coding convention fixes from Stefan
21 Revision 1.1 2007/09/17 10:23:31 cvetan
22 New TPC monitoring package from Stefan Kniege. The monitoring package can be started by running TPCMonitor.C macro located in macros folder.
26 ////////////////////////////////////////////////////////////////////////
28 //// AliTPCMonitorFFT class
30 //// Wrapper class to perform Fast Fourier Transformations.
31 //// The code is based on the Gnu Scientific Library.
32 //// See documentation of gsl for further details.
34 //// Author: Stefan Kniege, IKF, Frankfurt
37 /////////////////////////////////////////////////////////////////////////
40 #include "AliTPCMonitorFFT.h"
43 ClassImp(AliTPCMonitorFFT)
45 //__________________________________________________________________
46 AliTPCMonitorFFT::AliTPCMonitorFFT()
51 //__________________________________________________________________
52 AliTPCMonitorFFT::~AliTPCMonitorFFT()
58 //__________________________________________________________________
59 Int_t AliTPCMonitorFFT::ComplexRadix2ForwardWrap(Double_t* data, Int_t stride, size_t n)
62 // Wrapper function for forward fourier transformation
64 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
68 //__________________________________________________________________
69 Int_t AliTPCMonitorFFT::ComplexRadix2BackwardWrap(Double_t* data,Int_t stride, size_t n)
71 // Wrapper function for backward fourier transformation
73 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
77 //__________________________________________________________________
78 Int_t AliTPCMonitorFFT::ComplexRadix2InverseWrap(Double_t* data,Int_t stride, size_t n)
80 // Wrapper function for inverse fourier transformation
82 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
89 const Double_t knorm = 1.0 / n;
91 for (i = 0; i < n; i++)
93 REAL(data,stride,i) *= knorm;
94 IMAG(data,stride,i) *= knorm;
100 //__________________________________________________________________
101 Int_t AliTPCMonitorFFT::ComplexRadix2TransformWrap(Double_t* data, Int_t stride, size_t n, Int_t sign)
103 // Wrapper function for fourier transformation depending on sign ( forward 1;backward -1 )
110 if (n == 1) /* identity operation */
114 /* make sure that n is a power of 2 */
116 result = FFTBinaryLogn(n) ;
120 cout << " not a power of 2 " << endl;
127 /* apply fft recursion */
131 for (bit = 0; bit < logn; bit++)
133 Double_t wreal = 1.0;
134 Double_t wimag = 0.0;
136 const Double_t ktheta = 2.0 * ((Int_t) sign) * M_PI / ((Double_t) (2 * dual));
138 const Double_t ks = sin (ktheta);
139 const Double_t kt = sin (ktheta / 2.0);
140 const Double_t ks2 = 2.0 * kt * kt;
144 for (b = 0; b < dual; b++)
146 for (a = 0; a < n; a+= 2 * dual)
148 const size_t ki = b + a;
149 const size_t kj = b + a + dual;
151 const Double_t kt1real = REAL(data,stride,ki) + REAL(data,stride,kj);
152 const Double_t kt1imag = IMAG(data,stride,ki) + IMAG(data,stride,kj);
153 const Double_t kt2real = REAL(data,stride,ki) - REAL(data,stride,kj);
154 const Double_t kt2imag = IMAG(data,stride,ki) - IMAG(data,stride,kj);
156 REAL(data,stride,ki) = kt1real;
157 IMAG(data,stride,ki) = kt1imag;
158 REAL(data,stride,kj) = wreal*kt2real - wimag * kt2imag;
159 IMAG(data,stride,kj) = wreal*kt2imag + wimag * kt2real;
162 /* trignometric recurrence for w-> exp(i ktheta) w */
165 const Double_t ktmpreal = wreal - ks * wimag - ks2 * wreal;
166 const Double_t ktmpimag = wimag + ks * wreal - ks2 * wimag;
174 /* bit reverse the ordering of output data for decimation in
175 frequency algorithm */
177 //status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
178 status = ComplexBitReverseOrderWrap(data, stride, n, logn) ;
182 //__________________________________________________________________
183 Int_t AliTPCMonitorFFT::ComplexBitReverseOrderWrap(Double_t* data, Int_t stride, size_t n, Int_t logn) const
185 // Wrapper function from gnu scientific library
186 /* This is the Goldrader bit-reversal algorithm */
191 logn = 0 ; /* not needed for this algorithm */
193 for (i = 0; i < n - 1; i++)
199 const Double_t ktmpreal = REAL(data,stride,i);
200 const Double_t ktmpimag = IMAG(data,stride,i);
201 REAL(data,stride,i) = REAL(data,stride,j);
202 IMAG(data,stride,i) = IMAG(data,stride,j);
203 REAL(data,stride,j) = ktmpreal;
204 IMAG(data,stride,j) = ktmpimag;
220 //__________________________________________________________________
221 Int_t AliTPCMonitorFFT::FFTBinaryLogn(size_t n) const
224 // Return log on base 2
226 size_t binarylogn = 0 ;
235 ntest = (1 << binarylogn) ;
239 return -1 ; /* n is not a power of 2 */