Constant extrapolation for gain calibration
[u/mrichter/AliRoot.git] / TPC / AliTPCMonitorFFT.cxx
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$
18 Revision 1.2  2007/10/12 13:36:27  cvetan
19 Coding convention fixes from Stefan
20
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.
23
24 */ 
25
26 ////////////////////////////////////////////////////////////////////////
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 ////
37 /////////////////////////////////////////////////////////////////////////
38
39
40 #include "AliTPCMonitorFFT.h"
41
42
43 ClassImp(AliTPCMonitorFFT)
44   
45 //__________________________________________________________________
46   AliTPCMonitorFFT::AliTPCMonitorFFT()
47 {
48   // Constuctor
49 }
50  
51 //__________________________________________________________________
52 AliTPCMonitorFFT::~AliTPCMonitorFFT()
53 {
54   // Destructor
55 }
56
57
58 //__________________________________________________________________
59 Int_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 //__________________________________________________________________
69 Int_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 //__________________________________________________________________
78 Int_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   
89   const Double_t knorm = 1.0 / n;
90   size_t i;
91   for (i = 0; i < n; i++)
92     {
93       REAL(data,stride,i) *= knorm;
94       IMAG(data,stride,i) *= knorm;
95     }
96   
97   return ret;
98
99
100 //__________________________________________________________________
101 Int_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     {
133       Double_t wreal = 1.0;
134       Double_t wimag = 0.0;
135
136       const Double_t ktheta = 2.0 * ((Int_t) sign) * M_PI / ((Double_t) (2 * dual));
137
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;
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             {
148               const size_t ki = b + a;
149               const size_t kj = b + a + dual;
150               
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;
160             }
161
162           /* trignometric recurrence for w-> exp(i ktheta) w */
163
164           {
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;
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 //__________________________________________________________________
183 Int_t AliTPCMonitorFFT::ComplexBitReverseOrderWrap(Double_t* data, Int_t stride, size_t n, Int_t logn) const
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   
191   logn = 0 ; /* not needed for this algorithm */
192   
193   for (i = 0; i < n - 1; i++)
194     {
195       size_t k = n / 2 ;
196       
197       if (i < j)
198         {
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;
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 //__________________________________________________________________
221 Int_t AliTPCMonitorFFT::FFTBinaryLogn(size_t n) const
222 {
223
224   // Return log on base 2
225   size_t ntest ;
226   size_t binarylogn = 0 ;
227   size_t k = 1;
228   
229   while (k < n)
230     {
231       k *= 2;
232       binarylogn++;
233     }
234   
235   ntest = (1 << binarylogn) ;
236   
237   if (n != ntest )       
238     {
239       return -1 ; /* n is not a power of 2 */
240     } 
241   
242   return binarylogn;
243 }
244