Lines getting the matched track moved to a method in AliCalorimeterUtils. Lines copie...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoYlm.cxx
CommitLineData
17b6dded 1////////////////////////////////////////////////////////////////////////////////
2// //
3// AliFemtoYlm - the class to calculate varous components of spherical //
4// harmonics //
5// //
6// Authors: Adam Kisiel kisiel@mps.ohio-state.edu //
7// //
8////////////////////////////////////////////////////////////////////////////////
9
10#include "AliFemtoYlm.h"
11#include <TMath.h>
12#include <iostream>
13
14double *AliFemtoYlm::fgPrefactors = 0x0;
15int *AliFemtoYlm::fgPrefshift = 0x0;
16int *AliFemtoYlm::fgPlmshift = 0x0;
17
18AliFemtoYlm::AliFemtoYlm() {
19 InitializeYlms();
20}
21
22AliFemtoYlm::~AliFemtoYlm() {}
23
24
25AliFemtoYlm::AliFemtoYlm(const AliFemtoYlm& aYlm){
04328c66 26 fgPrefshift = aYlm.fgPrefshift;
17b6dded 27 InitializeYlms();
28}
29
30AliFemtoYlm& AliFemtoYlm::operator=(const AliFemtoYlm& aYlm){
31 if (this == &aYlm)
32 return *this;
33
34 InitializeYlms();
35
36 return *this;
37}
38
39std::complex<double> AliFemtoYlm::Ceiphi(double phi){
40 return std::complex<double>(cos(phi),sin(phi));
41}
42
43double AliFemtoYlm::Legendre(int ell, int em, double ctheta){
44 // Calculate a single Legendre value
45 // *** Warning - NOT optimal - calculated all Plms up to L !!!
46 double lbuf[36];
47 AliFemtoYlm::LegendreUpToYlm(ell, ctheta, lbuf);
48
49 return lbuf[fgPlmshift[ell]-abs(em)];
50}
51
52std::complex<double> AliFemtoYlm::Ylm(int ell,int m,double theta,double phi){
53 // Calculate Ylm spherical input
54 double ctheta;
55 std::complex<double> answer;
56 std::complex<double> ci(0.0,1.0);
57 ctheta=cos(theta);
58 answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi);
59
60 return answer;
61}
62
63std::complex<double> AliFemtoYlm::Ylm(int ell, int m, double x, double y, double z){
64 // Calculate Ylm cartesian input
65 std::complex<double> answer;
66 double ctheta,phi;
67 double r = sqrt(x*x+y*y+z*z);
68 if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0;
69 else ctheta=z/r;
70 phi=atan2(y,x);
71 answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi);
72
73 return answer;
74}
75
76double AliFemtoYlm::ReYlm(int ell, int m, double theta, double phi){
77 return real(AliFemtoYlm::Ylm(ell,m,theta,phi));
78}
79
80double AliFemtoYlm::ImYlm(int ell, int m, double theta, double phi){
81 return imag(AliFemtoYlm::Ylm(ell,m,theta,phi));
82}
83
84double AliFemtoYlm::ReYlm(int ell, int m, double x,double y,double z){
85 return real(AliFemtoYlm::Ylm(ell,m,x,y,z));
86}
87
88double AliFemtoYlm::ImYlm(int ell, int m, double x,double y,double z){
89 return imag(AliFemtoYlm::Ylm(ell,m,x,y,z));
90}
91
92void AliFemtoYlm::InitializeYlms()
93{
94 // Calculate prefactors for fast Ylm calculation
95
96 double oneoversqrtpi = 1.0/TMath::Sqrt(TMath::Pi());
97
98 fgPrefactors = (double *) malloc(sizeof(double) * 36);
99 fgPrefshift = (int *) malloc(sizeof(int) * 6);
100 fgPlmshift = (int *) malloc(sizeof(int) * 6);
101
102 // l=0 prefactors
103 fgPrefactors[0] = 0.5*oneoversqrtpi;
104
105 // l=1 prefactors
106 fgPrefactors[1] = 0.5*sqrt(3.0/2.0)*oneoversqrtpi;
107 fgPrefactors[2] = 0.5*sqrt(3.0)*oneoversqrtpi;
108 fgPrefactors[3] = -fgPrefactors[1];
109
110 // l=2 prefactors
111 fgPrefactors[4] = 0.25*sqrt(15.0/2.0)*oneoversqrtpi;
112 fgPrefactors[5] = 0.5*sqrt(15.0/2.0)*oneoversqrtpi;
113 fgPrefactors[6] = 0.25*sqrt(5.0)*oneoversqrtpi;
114 fgPrefactors[7] = -fgPrefactors[5];
115 fgPrefactors[8] = fgPrefactors[4];
116
117 // l=3 prefactors
118 fgPrefactors[9] = 0.125*sqrt(35.0)*oneoversqrtpi;
119 fgPrefactors[10] = 0.25*sqrt(105.0/2.0)*oneoversqrtpi;
120 fgPrefactors[11] = 0.125*sqrt(21.0)*oneoversqrtpi;
121 fgPrefactors[12] = 0.25*sqrt(7.0)*oneoversqrtpi;
122 fgPrefactors[13] = -fgPrefactors[11];
123 fgPrefactors[14] = fgPrefactors[10];
124 fgPrefactors[15] = -fgPrefactors[9];
125
126 // l=4 prefactors
127 fgPrefactors[16] = 3.0/16.0*sqrt(35.0/2.0)*oneoversqrtpi;
128 fgPrefactors[17] = 3.0/8.0*sqrt(35.0)*oneoversqrtpi;
129 fgPrefactors[18] = 3.0/8.0*sqrt(5.0/2.0)*oneoversqrtpi;
130 fgPrefactors[19] = 3.0/8.0*sqrt(5.0)*oneoversqrtpi;
131 fgPrefactors[20] = 3.0/16.0*oneoversqrtpi;
132 fgPrefactors[21] = -fgPrefactors[19];
133 fgPrefactors[22] = fgPrefactors[18];
134 fgPrefactors[23] = -fgPrefactors[17];
135 fgPrefactors[24] = fgPrefactors[16];
136
137 // l=5 prefactors
138 fgPrefactors[25] = 3.0/32.0*sqrt(77.0)*oneoversqrtpi;
139 fgPrefactors[26] = 3.0/16.0*sqrt(385.0/2.0)*oneoversqrtpi;
140 fgPrefactors[27] = 1.0/32.0*sqrt(385.0)*oneoversqrtpi;
141 fgPrefactors[28] = 1.0/8.0*sqrt(1155.0/2.0)*oneoversqrtpi;
142 fgPrefactors[29] = 1.0/16.0*sqrt(165.0/2.0)*oneoversqrtpi;
143 fgPrefactors[30] = 1.0/16.0*sqrt(11.0)*oneoversqrtpi;
144 fgPrefactors[31] = -fgPrefactors[29];
145 fgPrefactors[32] = fgPrefactors[28];
146 fgPrefactors[33] = -fgPrefactors[27];
147 fgPrefactors[34] = fgPrefactors[26];
148 fgPrefactors[35] = -fgPrefactors[25];
149
150 fgPrefshift[0] = 0;
151 fgPrefshift[1] = 2;
152 fgPrefshift[2] = 6;
153 fgPrefshift[3] = 12;
154 fgPrefshift[4] = 20;
155 fgPrefshift[5] = 30;
156
157 fgPlmshift[0] = 0;
158 fgPlmshift[1] = 2;
159 fgPlmshift[2] = 5;
160 fgPlmshift[3] = 9;
161 fgPlmshift[4] = 14;
162 fgPlmshift[5] = 20;
163}
164
165void AliFemtoYlm::LegendreUpToYlm(int lmax, double ctheta, double *lbuf)
166{
167 // Calculate a set of legendre polynomials up to a given l
168 // with spherical input
169 double sins[6];
170 double coss[6];
171 sins[0] = 0.0;
172 coss[0] = 1.0;
173 sins[1] = sqrt(1-ctheta*ctheta);
174 coss[1] = ctheta;
175 for (int iter=2; iter<6; iter++) {
176 sins[iter] = sins[iter-1]*sins[1];
177 coss[iter] = coss[iter-1]*coss[1];
178 }
179
180 // Legendre polynomials l=0
181 lbuf[0] = 1.0;
182
183 // Legendre polynomials l=1
184 if (lmax>0) {
185 lbuf[1] = sins[1];
186 lbuf[2] = coss[1];
187 }
188
189 // Legendre polynomials l=2
190 if (lmax>1) {
191 lbuf[3] = sins[2];
192 lbuf[4] = sins[1]*coss[1];
193 lbuf[5] = 3*coss[2]-1;
194 }
195
196 // Legendre polynomials l=3
197 if (lmax>2) {
198 lbuf[6] = sins[3];
199 lbuf[7] = sins[2]*coss[1];
200 lbuf[8] = (5*coss[2]-1)*sins[1];
201 lbuf[9] = 5*coss[3]-3*coss[1];
202 }
203
204 // Legendre polynomials l=4
205 if (lmax>3) {
206 lbuf[10] = sins[4];
207 lbuf[11] = sins[3]*coss[1];
208 lbuf[12] = (7*coss[2]-1)*sins[2];
209 lbuf[13] = (7*coss[3]-3*coss[1])*sins[1];
210 lbuf[14] = 35*coss[4]-30*coss[2]+3;
211 }
212
213 // Legendre polynomials l=5
214 if (lmax>4) {
215 lbuf[15] = sins[5];
216 lbuf[16] = sins[4]*coss[1];
217 lbuf[17] = (9*coss[2]-1)*sins[3];
218 lbuf[18] = (3*coss[3]-1*coss[1])*sins[2];
219 lbuf[19] = (21*coss[4]-14*coss[2]+1)*sins[1];
220 lbuf[20] = 63*coss[5]-70*coss[3]+15*coss[1];
221 }
222}
223
224void AliFemtoYlm::YlmUpToL(int lmax, double x, double y, double z, std::complex<double> *ylms)
225{
226 // Calculate a set of Ylms up to a given l
227 // with cartesian input
228 double ctheta,phi;
229
230 double r = sqrt(x*x+y*y+z*z);
231 if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0;
232 else ctheta=z/r;
233 phi=atan2(y,x);
234
235 YlmUpToL(lmax, ctheta, phi, ylms);
236
237}
238
239void AliFemtoYlm::YlmUpToL(int lmax, double ctheta, double phi, std::complex<double> *ylms)
240{
241 // Calculate a set of Ylms up to a given l
242 // with spherical input
243 int lcur = 0;
244 double lpol;
245
246 double coss[6];
247 double sins[6];
248
249 double lbuf[36];
250 LegendreUpToYlm(lmax, ctheta, lbuf);
251
252 for (int iter=1; iter<=lmax; iter++) {
253 coss[iter-1] = cos(iter*phi);
254 sins[iter-1] = sin(iter*phi);
255 }
256 ylms[lcur++] = fgPrefactors[0]*lbuf[0] * std::complex<double>(1,0);
257
258 for (int il = 1; il<=lmax; il++) {
259 // First im = 0
260 ylms[lcur+il] = fgPrefactors[fgPrefshift[il]]*lbuf[fgPlmshift[il]]*std::complex<double>(1.0,0.0);
261 // Im != 0
262 for (int im=1; im<=il; im++) {
263 lpol = lbuf[fgPlmshift[il]-im];
264 ylms[lcur+il-im] = fgPrefactors[fgPrefshift[il]-im]*lpol*std::complex<double>(coss[im-1],-sins[im-1]);
265 ylms[lcur+il+im] = fgPrefactors[fgPrefshift[il]+im]*lpol*std::complex<double>(coss[im-1],sins[im-1]);
266 }
267 lcur += 2*il + 1;
268 }
269}