[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemtoUser / AliFemtoYlm.cxx

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// AliFemtoYlm - the class to calculate varous components of spherical        //
//  harmonics                                                                 //
//                                                                            //
// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                             //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoYlm.h"
#include <TMath.h>
#include <iostream>

double *AliFemtoYlm::fgPrefactors = 0x0;
int    *AliFemtoYlm::fgPrefshift = 0x0;
int    *AliFemtoYlm::fgPlmshift = 0x0;

AliFemtoYlm::AliFemtoYlm() {
  InitializeYlms();
}

AliFemtoYlm::~AliFemtoYlm() {
  free(fgPrefactors);
  free(fgPrefshift);
  free(fgPlmshift);
}


AliFemtoYlm::AliFemtoYlm(const AliFemtoYlm& aYlm){
  fgPrefshift = aYlm.fgPrefshift;
  InitializeYlms();
}

AliFemtoYlm& AliFemtoYlm::operator=(const AliFemtoYlm& aYlm){
  if (this == &aYlm)
    return *this;

  InitializeYlms();

  return *this;
}

std::complex<double> AliFemtoYlm::Ceiphi(double phi){
  return std::complex<double>(cos(phi),sin(phi));
}

double AliFemtoYlm::Legendre(int ell, int em, double ctheta){
  // Calculate a single Legendre value
  // *** Warning - NOT optimal - calculated all Plms up to L !!!
  double lbuf[36];
  AliFemtoYlm::LegendreUpToYlm(ell, ctheta, lbuf);

  return lbuf[fgPlmshift[ell]-abs(em)];
}

std::complex<double> AliFemtoYlm::Ylm(int ell,int m,double theta,double phi){
  // Calculate Ylm spherical input
  double ctheta;
  std::complex<double> answer;
  std::complex<double> ci(0.0,1.0);
  ctheta=cos(theta);
  answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi);

  return answer;
}

std::complex<double> AliFemtoYlm::Ylm(int ell, int m, double x, double y, double z){
  // Calculate Ylm cartesian input
  std::complex<double> answer; 
  double ctheta,phi;
  double r = sqrt(x*x+y*y+z*z);
  if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0;
  else ctheta=z/r;
  phi=atan2(y,x);
  answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi);

  return answer;	
}

double AliFemtoYlm::ReYlm(int ell, int m, double theta, double phi){
	return real(AliFemtoYlm::Ylm(ell,m,theta,phi));
}

double AliFemtoYlm::ImYlm(int ell, int m, double theta, double phi){
	return imag(AliFemtoYlm::Ylm(ell,m,theta,phi));
}

double AliFemtoYlm::ReYlm(int ell, int m, double x,double y,double z){
	return real(AliFemtoYlm::Ylm(ell,m,x,y,z));
}

double AliFemtoYlm::ImYlm(int ell, int m, double x,double y,double z){
	return imag(AliFemtoYlm::Ylm(ell,m,x,y,z));
}

void AliFemtoYlm::InitializeYlms()
{
  // Calculate prefactors for fast Ylm calculation

  double oneoversqrtpi = 1.0/TMath::Sqrt(TMath::Pi());

  if(fgPrefactors!=NULL)
    free(fgPrefactors);
  if(fgPrefshift!=NULL)
    free(fgPrefshift);
  if(fgPlmshift!=NULL)
    free(fgPlmshift);

  fgPrefactors = (double *) malloc(sizeof(double) * 36);
  fgPrefshift  = (int *) malloc(sizeof(int) * 6);
  fgPlmshift   = (int *) malloc(sizeof(int) * 6);

  // l=0 prefactors
  fgPrefactors[0] = 0.5*oneoversqrtpi;

  // l=1 prefactors
  fgPrefactors[1] = 0.5*sqrt(3.0/2.0)*oneoversqrtpi;
  fgPrefactors[2] = 0.5*sqrt(3.0)*oneoversqrtpi;
  fgPrefactors[3] = -fgPrefactors[1];

  // l=2 prefactors
  fgPrefactors[4] = 0.25*sqrt(15.0/2.0)*oneoversqrtpi;
  fgPrefactors[5] = 0.5*sqrt(15.0/2.0)*oneoversqrtpi;
  fgPrefactors[6] = 0.25*sqrt(5.0)*oneoversqrtpi;
  fgPrefactors[7] = -fgPrefactors[5];
  fgPrefactors[8] = fgPrefactors[4];

  // l=3 prefactors
  fgPrefactors[9]  = 0.125*sqrt(35.0)*oneoversqrtpi;
  fgPrefactors[10] = 0.25*sqrt(105.0/2.0)*oneoversqrtpi;
  fgPrefactors[11] = 0.125*sqrt(21.0)*oneoversqrtpi;
  fgPrefactors[12] = 0.25*sqrt(7.0)*oneoversqrtpi;
  fgPrefactors[13] = -fgPrefactors[11];
  fgPrefactors[14] = fgPrefactors[10];
  fgPrefactors[15] = -fgPrefactors[9];

  // l=4 prefactors
  fgPrefactors[16] = 3.0/16.0*sqrt(35.0/2.0)*oneoversqrtpi;
  fgPrefactors[17] = 3.0/8.0*sqrt(35.0)*oneoversqrtpi;
  fgPrefactors[18] = 3.0/8.0*sqrt(5.0/2.0)*oneoversqrtpi;
  fgPrefactors[19] = 3.0/8.0*sqrt(5.0)*oneoversqrtpi;
  fgPrefactors[20] = 3.0/16.0*oneoversqrtpi;
  fgPrefactors[21] = -fgPrefactors[19];
  fgPrefactors[22] = fgPrefactors[18];
  fgPrefactors[23] = -fgPrefactors[17];
  fgPrefactors[24] = fgPrefactors[16];

  // l=5 prefactors
  fgPrefactors[25] = 3.0/32.0*sqrt(77.0)*oneoversqrtpi;
  fgPrefactors[26] = 3.0/16.0*sqrt(385.0/2.0)*oneoversqrtpi;
  fgPrefactors[27] = 1.0/32.0*sqrt(385.0)*oneoversqrtpi;
  fgPrefactors[28] = 1.0/8.0*sqrt(1155.0/2.0)*oneoversqrtpi;
  fgPrefactors[29] = 1.0/16.0*sqrt(165.0/2.0)*oneoversqrtpi;
  fgPrefactors[30] = 1.0/16.0*sqrt(11.0)*oneoversqrtpi;
  fgPrefactors[31] = -fgPrefactors[29];
  fgPrefactors[32] = fgPrefactors[28];
  fgPrefactors[33] = -fgPrefactors[27];
  fgPrefactors[34] = fgPrefactors[26];
  fgPrefactors[35] = -fgPrefactors[25];

  fgPrefshift[0] = 0;
  fgPrefshift[1] = 2;
  fgPrefshift[2] = 6;
  fgPrefshift[3] = 12;
  fgPrefshift[4] = 20;
  fgPrefshift[5] = 30;

  fgPlmshift[0] = 0;
  fgPlmshift[1] = 2;
  fgPlmshift[2] = 5;
  fgPlmshift[3] = 9;
  fgPlmshift[4] = 14;
  fgPlmshift[5] = 20;
}

void AliFemtoYlm::LegendreUpToYlm(int lmax, double ctheta, double *lbuf)
{
  // Calculate a set of legendre polynomials up to a given l
  // with spherical input
  double sins[6];
  double coss[6];
  sins[0] = 0.0;
  coss[0] = 1.0;
  sins[1] = sqrt(1-ctheta*ctheta);
  coss[1] = ctheta;
  for (int iter=2; iter<6; iter++) {
    sins[iter] = sins[iter-1]*sins[1];
    coss[iter] = coss[iter-1]*coss[1];
  }

  // Legendre polynomials l=0
  lbuf[0] = 1.0;

  // Legendre polynomials l=1
  if (lmax>0) {
    lbuf[1] = sins[1];
    lbuf[2] = coss[1];
  }

  // Legendre polynomials l=2
  if (lmax>1) {
    lbuf[3] = sins[2];
    lbuf[4] = sins[1]*coss[1];
    lbuf[5] = 3*coss[2]-1;
  }

  // Legendre polynomials l=3
  if (lmax>2) {
    lbuf[6] = sins[3];
    lbuf[7] = sins[2]*coss[1];
    lbuf[8] = (5*coss[2]-1)*sins[1];
    lbuf[9] = 5*coss[3]-3*coss[1];
  }

  // Legendre polynomials l=4
  if (lmax>3) {
    lbuf[10] = sins[4];
    lbuf[11] = sins[3]*coss[1];
    lbuf[12] = (7*coss[2]-1)*sins[2];
    lbuf[13] = (7*coss[3]-3*coss[1])*sins[1];
    lbuf[14] = 35*coss[4]-30*coss[2]+3;
  }

  // Legendre polynomials l=5
  if (lmax>4) {
    lbuf[15] = sins[5];
    lbuf[16] = sins[4]*coss[1];
    lbuf[17] = (9*coss[2]-1)*sins[3];
    lbuf[18] = (3*coss[3]-1*coss[1])*sins[2];
    lbuf[19] = (21*coss[4]-14*coss[2]+1)*sins[1];
    lbuf[20] = 63*coss[5]-70*coss[3]+15*coss[1];
  }
}

void AliFemtoYlm::YlmUpToL(int lmax, double x, double y, double z, std::complex<double> *ylms)
{
  // Calculate a set of Ylms up to a given l
  // with cartesian input
  double ctheta,phi;

  double r = sqrt(x*x+y*y+z*z);
  if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0;
  else ctheta=z/r;
  phi=atan2(y,x);
  
  YlmUpToL(lmax, ctheta, phi, ylms);

}

void AliFemtoYlm::YlmUpToL(int lmax, double ctheta, double phi, std::complex<double> *ylms)
{
  // Calculate a set of Ylms up to a given l
  // with spherical input
  int lcur = 0;  
  double lpol;
  
  double coss[6];
  double sins[6];

  double lbuf[36];
  LegendreUpToYlm(lmax, ctheta, lbuf);

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