Patch for the tracker
[u/mrichter/AliRoot.git] / TPC / AliTPCMonitorFFT.cxx
CommitLineData
48265b32 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/*
17$Log$
fb3305d1 18Revision 1.2 2007/10/12 13:36:27 cvetan
19Coding convention fixes from Stefan
20
ca7b8371 21Revision 1.1 2007/09/17 10:23:31 cvetan
22New TPC monitoring package from Stefan Kniege. The monitoring package can be started by running TPCMonitor.C macro located in macros folder.
23
48265b32 24*/
25
ca7b8371 26////////////////////////////////////////////////////////////////////////
fb3305d1 27////
28//// AliTPCMonitorFFT class
29////
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.
33////
34//// Author: Stefan Kniege, IKF, Frankfurt
35////
36////
ca7b8371 37/////////////////////////////////////////////////////////////////////////
38
39
48265b32 40#include "AliTPCMonitorFFT.h"
41
42
43ClassImp(AliTPCMonitorFFT)
44
45//__________________________________________________________________
46 AliTPCMonitorFFT::AliTPCMonitorFFT()
47{
48 // Constuctor
49}
50
51//__________________________________________________________________
52AliTPCMonitorFFT::~AliTPCMonitorFFT()
53{
54 // Destructor
55}
56
57
58//__________________________________________________________________
59Int_t AliTPCMonitorFFT::ComplexRadix2ForwardWrap(Double_t* data, Int_t stride, size_t n)
60{
61
62 // Wrapper function for forward fourier transformation
63 Int_t direction = 1;
64 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
65 return ret ;
66}
67
68//__________________________________________________________________
69Int_t AliTPCMonitorFFT::ComplexRadix2BackwardWrap(Double_t* data,Int_t stride, size_t n)
70{
71 // Wrapper function for backward fourier transformation
72 Int_t direction = -1;
73 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
74 return ret ;
75}
76
77//__________________________________________________________________
78Int_t AliTPCMonitorFFT::ComplexRadix2InverseWrap(Double_t* data,Int_t stride, size_t n)
79{
80 // Wrapper function for inverse fourier transformation
81 Int_t direction = -1;
82 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
83
84 if (ret)
85 {
86 return ret;
87 }
88
fb3305d1 89 const Double_t knorm = 1.0 / n;
48265b32 90 size_t i;
91 for (i = 0; i < n; i++)
92 {
fb3305d1 93 REAL(data,stride,i) *= knorm;
94 IMAG(data,stride,i) *= knorm;
48265b32 95 }
96
97 return ret;
98}
99
100//__________________________________________________________________
101Int_t AliTPCMonitorFFT::ComplexRadix2TransformWrap(Double_t* data, Int_t stride, size_t n, Int_t sign)
102{
103 // Wrapper function for fourier transformation depending on sign ( forward 1;backward -1 )
104 Int_t result ;
105 size_t dual;
106 size_t bit;
107 size_t logn = 0;
108 Int_t status;
109
110 if (n == 1) /* identity operation */
111 {
112 return 0 ;
113 }
114 /* make sure that n is a power of 2 */
115
116 result = FFTBinaryLogn(n) ;
117
118 if (result == -1)
119 {
120 cout << " not a power of 2 " << endl;
121 }
122 else
123 {
124 logn = result ;
125 }
126
127 /* apply fft recursion */
128
129 dual = n / 2;
130
131 for (bit = 0; bit < logn; bit++)
132 {
fb3305d1 133 Double_t wreal = 1.0;
134 Double_t wimag = 0.0;
48265b32 135
fb3305d1 136 const Double_t ktheta = 2.0 * ((Int_t) sign) * M_PI / ((Double_t) (2 * dual));
48265b32 137
fb3305d1 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;
48265b32 141
142 size_t a, b;
143
144 for (b = 0; b < dual; b++)
145 {
146 for (a = 0; a < n; a+= 2 * dual)
147 {
fb3305d1 148 const size_t ki = b + a;
149 const size_t kj = b + a + dual;
48265b32 150
fb3305d1 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);
155
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;
48265b32 160 }
161
fb3305d1 162 /* trignometric recurrence for w-> exp(i ktheta) w */
48265b32 163
164 {
fb3305d1 165 const Double_t ktmpreal = wreal - ks * wimag - ks2 * wreal;
166 const Double_t ktmpimag = wimag + ks * wreal - ks2 * wimag;
167 wreal = ktmpreal;
168 wimag = ktmpimag;
48265b32 169 }
170 }
171 dual /= 2;
172 }
173
174 /* bit reverse the ordering of output data for decimation in
175 frequency algorithm */
176
177 //status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
178 status = ComplexBitReverseOrderWrap(data, stride, n, logn) ;
179 return 0;
180}
181
182//__________________________________________________________________
5dbad769 183Int_t AliTPCMonitorFFT::ComplexBitReverseOrderWrap(Double_t* data, Int_t stride, size_t n, Int_t /*logn*/) const
48265b32 184{
185 // Wrapper function from gnu scientific library
186 /* This is the Goldrader bit-reversal algorithm */
187
188 size_t i;
189 size_t j = 0;
190
5dbad769 191// logn = 0 ; /* not needed for this algorithm */
48265b32 192
193 for (i = 0; i < n - 1; i++)
194 {
195 size_t k = n / 2 ;
196
197 if (i < j)
198 {
fb3305d1 199 const Double_t ktmpreal = REAL(data,stride,i);
200 const Double_t ktmpimag = IMAG(data,stride,i);
48265b32 201 REAL(data,stride,i) = REAL(data,stride,j);
202 IMAG(data,stride,i) = IMAG(data,stride,j);
fb3305d1 203 REAL(data,stride,j) = ktmpreal;
204 IMAG(data,stride,j) = ktmpimag;
48265b32 205 }
206
207 while (k <= j)
208 {
209 j = j - k ;
210 k = k / 2 ;
211 }
212
213 j += k ;
214 }
215
216 return 0;
217
218}
219
220//__________________________________________________________________
fb3305d1 221Int_t AliTPCMonitorFFT::FFTBinaryLogn(size_t n) const
48265b32 222{
223
224 // Return log on base 2
225 size_t ntest ;
fb3305d1 226 size_t binarylogn = 0 ;
48265b32 227 size_t k = 1;
228
229 while (k < n)
230 {
231 k *= 2;
fb3305d1 232 binarylogn++;
48265b32 233 }
234
fb3305d1 235 ntest = (1 << binarylogn) ;
48265b32 236
237 if (n != ntest )
238 {
239 return -1 ; /* n is not a power of 2 */
240 }
241
fb3305d1 242 return binarylogn;
48265b32 243}
244