PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDirectYlm.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoCorrFctnDirectYlm - Correlation function that is binned in Ylms    //
   4 // directly. Provides a way to store the numerator and denominator            //
   5 // in Ylms directly and correctly calculate the correlation                   //
   6 // function from them.                                                        //
   7 //                                                                            //
   8 // Authors: Adam Kisiel kisiel@mps.ohio-state.edu                             //
   9 //                                                                            //
  10 ////////////////////////////////////////////////////////////////////////////////
  11 #include "AliFemtoCorrFctnDirectYlm.h"
  12 #include <TMath.h>
  13 #include <iostream>
  14
  15 using namespace std;
  16
  17 AliFemtoCorrFctnDirectYlm::AliFemtoCorrFctnDirectYlm(const char *name, int maxl, int ibin=30, double vmin=0.0, double vmax=0.3):
  18   fnumsreal(0),
  19   fnumsimag(0),
  20   fdensreal(0),
  21   fdensimag(0),
  22   fbinctn(0),
  23   fbinctd(0),
  24   fcovnum(0),
  25   fcovden(0),
  26   fcovmnum(0),
  27   fcovmden(0),
  28   fMaxL(0),
  29   fMaxJM(0),
  30   fels(0),
  31   fems(0),
  32   felsi(0),
  33   femsi(0),
  34   fYlmBuffer(0),
  35   factorials(0),
  36   fSout(0.0),
  37   fSside(0.0),
  38   fSlong(0.0)
  39 {
  40   // Main constructor
  41   fMaxL = maxl;
  42   fMaxJM = (maxl+1)*(maxl+1);
  43
  44   cout <<  "Size is " << sizeof(double) << " " << sizeof(complex<double>) << endl;
  45
  46   // Fill in factorials table
  47   factorials = (double *) malloc(sizeof(double) * (4 * (maxl + 1)));
  48   int fac = 1;
  49   factorials[0] = 1;
  50   for (int iter=1; iter<4*(maxl+1); iter++)
  51     {
  52       fac *= iter;
  53       factorials[iter] = fac;
  54     }
  55
  56   // Fill in els and ems table
  57   int el = 0;
  58   int em = 0;
  59   int il = 0;
  60   fels = (double *) malloc(sizeof(double) * (fMaxJM));
  61   fems = (double *) malloc(sizeof(double) * (fMaxJM));
  62   felsi = (int *) malloc(sizeof(int) * (fMaxJM));
  63   femsi = (int *) malloc(sizeof(int) * (fMaxJM));
  64   do {
  65     fels[il] = el;
  66     fems[il] = em;
  67     felsi[il] = (int) el;
  68     femsi[il] = (int) em;
  69
  70     cout << "il el em " << il << " " << felsi[il] << " " << femsi[il] << endl;
  71     em++;
  72     il++;
  73     if (em > el) {
  74       el++;
  75       em = -el;
  76     }
  77   }
  78   while (el <= maxl);
  79
  80   for (il=0; il<fMaxJM; il++)
  81     cout << "il el em " << il << " " << felsi[il] << " " << femsi[il] << endl;
  82
  83   // Create numerator and denominator historgrams
  84   //  int sthp = sizeof(TH1D *);
  85   //  fnumsreal = (TH1D **) malloc(sthp * fMaxJM);
  86 //   fnumsreal = new TH1D * [fMaxJM];
  87 //   fnumsimag = new TH1D * [fMaxJM];
  88 //   fdensreal = new TH1D * [fMaxJM];
  89 //   fdensimag = new TH1D * [fMaxJM];
  90   fnumsreal = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
  91   fnumsimag = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
  92   fdensreal = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
  93   fdensimag = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
  94
  95   char bufname[200];
  96   for (int ihist=0; ihist<fMaxJM; ihist++) {
  97     sprintf(bufname, "NumReYlm%i%i%s", felsi[ihist], femsi[ihist]<0 ? felsi[ihist]-femsi[ihist] : femsi[ihist], name);
  98     fnumsreal[ihist] = new TH1D(bufname, bufname, ibin, vmin, vmax);
  99     sprintf(bufname, "NumImYlm%i%i%s", felsi[ihist], femsi[ihist]<0 ? felsi[ihist]-femsi[ihist] : femsi[ihist], name);
 100     fnumsimag[ihist] = new TH1D(bufname, bufname, ibin, vmin, vmax);
 101     sprintf(bufname, "DenReYlm%i%i%s", felsi[ihist], femsi[ihist]<0 ? felsi[ihist]-femsi[ihist] : femsi[ihist], name);
 102     fdensreal[ihist] = new TH1D(bufname, bufname, ibin, vmin, vmax);
 103     sprintf(bufname, "DenImYlm%i%i%s", felsi[ihist], femsi[ihist]<0 ? felsi[ihist]-femsi[ihist] : femsi[ihist], name);
 104     fdensimag[ihist] = new TH1D(bufname, bufname, ibin, vmin, vmax);
 105
 106     fnumsreal[ihist]->Sumw2();
 107     fnumsimag[ihist]->Sumw2();
 108     fdensreal[ihist]->Sumw2();
 109     fdensimag[ihist]->Sumw2();
 110   }
 111
 112   sprintf(bufname, "BinCountNum%s", name);
 113   fbinctn = new TH1D(bufname, bufname, ibin, vmin, vmax);
 114
 115   sprintf(bufname, "BinCountDen%s", name);
 116   fbinctd = new TH1D(bufname, bufname, ibin, vmin, vmax);
 117
 118   fYlmBuffer = (complex<double> *) malloc(sizeof(complex<double>) * fMaxJM);
 119
 120   // Covariance matrices
 121   fcovmnum = (double *) malloc(sizeof(double) * fMaxJM * fMaxJM * 4 * ibin);
 122   fcovmden = (double *) malloc(sizeof(double) * fMaxJM * fMaxJM * 4 * ibin);
 123
 124   fcovnum = 0;
 125   fcovden = 0;
 126
 127   AliFemtoYlm::InitializeYlms();
 128 }
 129
 130
 131 AliFemtoCorrFctnDirectYlm::AliFemtoCorrFctnDirectYlm():
 132   fnumsreal(0),
 133   fnumsimag(0),
 134   fdensreal(0),
 135   fdensimag(0),
 136   fbinctn(0),
 137   fbinctd(0),
 138   fcovnum(0),
 139   fcovden(0),
 140   fcovmnum(0),
 141   fcovmden(0),
 142   fMaxL(0),
 143   fMaxJM(0),
 144   fels(0),
 145   fems(0),
 146   felsi(0),
 147   femsi(0),
 148   fYlmBuffer(0),
 149   factorials(0),
 150   fSout(0.0),
 151   fSside(0.0),
 152   fSlong(0.0)
 153 {
 154   // Default constructor
 155   AliFemtoCorrFctnDirectYlm("AliFemtoCorrFctnDirectYlm",2);
 156 }
 157
 158 AliFemtoCorrFctnDirectYlm::AliFemtoCorrFctnDirectYlm(const AliFemtoCorrFctnDirectYlm& aCorrFctn):
 159   AliFemtoCorrFctn(),
 160   fnumsreal(0),
 161   fnumsimag(0),
 162   fdensreal(0),
 163   fdensimag(0),
 164   fbinctn(0),
 165   fbinctd(0),
 166   fcovnum(0),
 167   fcovden(0),
 168   fcovmnum(0),
 169   fcovmden(0),
 170   fMaxL(0),
 171   fMaxJM(0),
 172   fels(0),
 173   fems(0),
 174   felsi(0),
 175   femsi(0),
 176   fYlmBuffer(0),
 177   factorials(0),
 178   fSout(0.0),
 179   fSside(0.0),
 180   fSlong(0.0)
 181 {
 182   // Copy constructor
 183   int ibin = aCorrFctn.fbinctn->GetNbinsX();
 184
 185   fMaxL = aCorrFctn.fMaxL;
 186   fMaxJM = (fMaxL+1)*(fMaxL+1);
 187
 188   // Fill in factorials table
 189   factorials = (double *) malloc(sizeof(double) * (4 * (fMaxL + 1)));
 190   for (int iter=1; iter<4*(fMaxL+1); iter++)
 191     {
 192       factorials[iter] = aCorrFctn.factorials[iter];
 193     }
 194
 195   // Fill in els and ems table
 196   int el = 0;
 197   int em = 0;
 198   int il = 0;
 199   fels = (double *) malloc(sizeof(double) * (fMaxJM));
 200   fems = (double *) malloc(sizeof(double) * (fMaxJM));
 201   felsi = (int *) malloc(sizeof(int) * (fMaxJM));
 202   femsi = (int *) malloc(sizeof(int) * (fMaxJM));
 203   do {
 204     fels[il] = el;
 205     fems[il] = em;
 206     felsi[il] = (int) el;
 207     femsi[il] = (int) em;
 208
 209     em++;
 210     il++;
 211     if (em > el) {
 212       el++;
 213       em = -el;
 214     }
 215   }
 216   while (el <= fMaxL);
 217
 218   fnumsreal = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
 219   fnumsimag = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
 220   fdensreal = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
 221   fdensimag = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
 222
 223   for (int ihist=0; ihist<fMaxJM; ihist++) {
 224     if (aCorrFctn.fnumsreal[ihist])
 225       fnumsreal[ihist] = new TH1D(*aCorrFctn.fnumsreal[ihist]);
 226     else
 227       fnumsreal[ihist] = 0;
 228     if (aCorrFctn.fnumsimag[ihist])
 229       fnumsimag[ihist] = new TH1D(*aCorrFctn.fnumsimag[ihist]);
 230     else
 231       fnumsimag[ihist] = 0;
 232     if (aCorrFctn.fdensreal[ihist])
 233       fdensreal[ihist] = new TH1D(*aCorrFctn.fdensreal[ihist]);
 234     else
 235       fdensreal[ihist] = 0;
 236     if (aCorrFctn.fdensimag[ihist])
 237       fdensimag[ihist] = new TH1D(*aCorrFctn.fdensimag[ihist]);
 238     else
 239       fdensimag[ihist] = 0;
 240   }
 241
 242   if (aCorrFctn.fbinctn)
 243     fbinctn = new TH1D(*aCorrFctn.fbinctn);
 244   else
 245     fbinctn = 0;
 246   if (aCorrFctn.fbinctd)
 247     fbinctd = new TH1D(*aCorrFctn.fbinctd);
 248   else
 249     fbinctd = 0;
 250
 251   fYlmBuffer = (complex<double> *) malloc(sizeof(complex<double>) * fMaxJM);
 252
 253   // Covariance matrices
 254   fcovmnum = (double *) malloc(sizeof(double) * fMaxJM * fMaxJM * 4 * ibin);
 255   fcovmden = (double *) malloc(sizeof(double) * fMaxJM * fMaxJM * 4 * ibin);
 256
 257   for (int iter=0; iter<fMaxJM * fMaxJM * 4 * ibin; iter++) {
 258     fcovmnum[iter] = aCorrFctn.fcovmnum[iter];
 259     fcovmden[iter] = aCorrFctn.fcovmden[iter];
 260   }
 261
 262   if (aCorrFctn.fcovnum)
 263     fcovnum = new TH3D(*aCorrFctn.fcovnum);
 264   else
 265     fcovnum = 0;
 266   if (aCorrFctn.fcovden)
 267     fcovden = new TH3D(*aCorrFctn.fcovden);
 268   else
 269     fcovden = 0;
 270
 271   fSout = aCorrFctn.fSout;
 272   fSside = aCorrFctn.fSside;
 273   fSlong = aCorrFctn.fSlong;
 274 }
 275
 276 AliFemtoCorrFctnDirectYlm& AliFemtoCorrFctnDirectYlm::operator=(const AliFemtoCorrFctnDirectYlm& aCorrFctn)
 277 {
 278   // assignment operator
 279   if (this == &aCorrFctn)
 280     return *this;
 281
 282   int ibin = aCorrFctn.fbinctn->GetNbinsX();
 283
 284   fMaxL = aCorrFctn.fMaxL;
 285   fMaxJM = (fMaxL+1)*(fMaxL+1);
 286
 287   // Fill in factorials table
 288   factorials = (double *) malloc(sizeof(double) * (4 * (fMaxL + 1)));
 289   for (int iter=1; iter<4*(fMaxL+1); iter++)
 290     {
 291       factorials[iter] = aCorrFctn.factorials[iter];
 292     }
 293
 294   // Fill in els and ems table
 295   int el = 0;
 296   int em = 0;
 297   int il = 0;
 298   fels = (double *) malloc(sizeof(double) * (fMaxJM));
 299   fems = (double *) malloc(sizeof(double) * (fMaxJM));
 300   felsi = (int *) malloc(sizeof(int) * (fMaxJM));
 301   femsi = (int *) malloc(sizeof(int) * (fMaxJM));
 302   do {
 303     fels[il] = el;
 304     fems[il] = em;
 305     felsi[il] = (int) el;
 306     femsi[il] = (int) em;
 307
 308     em++;
 309     il++;
 310     if (em > el) {
 311       el++;
 312       em = -el;
 313     }
 314   }
 315   while (el <= fMaxL);
 316
 317   fnumsreal = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
 318   fnumsimag = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
 319   fdensreal = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
 320   fdensimag = (TH1D **) malloc(sizeof(TH1D *) * fMaxJM);
 321
 322   for (int ihist=0; ihist<fMaxJM; ihist++) {
 323     if (aCorrFctn.fnumsreal[ihist])
 324       fnumsreal[ihist] = new TH1D(*aCorrFctn.fnumsreal[ihist]);
 325     else
 326       fnumsreal[ihist] = 0;
 327     if (aCorrFctn.fnumsimag[ihist])
 328       fnumsimag[ihist] = new TH1D(*aCorrFctn.fnumsimag[ihist]);
 329     else
 330       fnumsimag[ihist] = 0;
 331     if (aCorrFctn.fdensreal[ihist])
 332       fdensreal[ihist] = new TH1D(*aCorrFctn.fdensreal[ihist]);
 333     else
 334       fdensreal[ihist] = 0;
 335     if (aCorrFctn.fdensimag[ihist])
 336       fdensimag[ihist] = new TH1D(*aCorrFctn.fdensimag[ihist]);
 337     else
 338       fdensimag[ihist] = 0;
 339   }
 340
 341   if (aCorrFctn.fbinctn)
 342     fbinctn = new TH1D(*aCorrFctn.fbinctn);
 343   else
 344     fbinctn = 0;
 345   if (aCorrFctn.fbinctd)
 346     fbinctd = new TH1D(*aCorrFctn.fbinctd);
 347   else
 348     fbinctd = 0;
 349
 350   fYlmBuffer = (complex<double> *) malloc(sizeof(complex<double>) * fMaxJM);
 351
 352   // Covariance matrices
 353   fcovmnum = (double *) malloc(sizeof(double) * fMaxJM * fMaxJM * 4 * ibin);
 354   fcovmden = (double *) malloc(sizeof(double) * fMaxJM * fMaxJM * 4 * ibin);
 355
 356   for (int iter=0; iter<fMaxJM * fMaxJM * 4 * ibin; iter++) {
 357     fcovmnum[iter] = aCorrFctn.fcovmnum[iter];
 358     fcovmden[iter] = aCorrFctn.fcovmden[iter];
 359   }
 360
 361   if (aCorrFctn.fcovnum)
 362     fcovnum = new TH3D(*aCorrFctn.fcovnum);
 363   else
 364     fcovnum = 0;
 365   if (aCorrFctn.fcovden)
 366     fcovden = new TH3D(*aCorrFctn.fcovden);
 367   else
 368     fcovden = 0;
 369
 370   fSout = aCorrFctn.fSout;
 371   fSside = aCorrFctn.fSside;
 372   fSlong = aCorrFctn.fSlong;
 373
 374   return *this;
 375 }
 376
 377 AliFemtoCorrFctnDirectYlm::~AliFemtoCorrFctnDirectYlm()
 378 {
 379   // Destructor
 380   for (int ihist=0; ihist<fMaxJM; ihist++) {
 381     delete fnumsreal[ihist];
 382     delete fnumsimag[ihist];
 383     delete fdensreal[ihist];
 384     delete fdensimag[ihist];
 385   }
 386
 387   delete fbinctn;
 388   delete fbinctd;
 389
 390   //  delete fnumsreal;
 391   //  delete fnumsimag;
 392   //  delete fdensreal;
 393   //  delete fdensimag;
 394
 395   free( fnumsreal);
 396   free( fnumsimag);
 397   free( fdensreal);
 398   free( fdensimag);
 399
 400   free(factorials);
 401   free(fels);
 402   free(fems);
 403   free(felsi);
 404   free(femsi);
 405   free(fYlmBuffer);
 406
 407   free(fcovmnum);
 408   free(fcovmden);
 409
 410   if (fcovnum) delete fcovnum;
 411   if (fcovden) delete fcovden;
 412 }
 413
 414 double AliFemtoCorrFctnDirectYlm::ClebschGordan(double aJot1, double aEm1, double aJot2, double aEm2, double aJot, double aEm)
 415 {
 416   // Calculate Clebsh-Gordan coefficient
 417   int mint, maxt;
 418   double cgc = 0.0;
 419   int titer;
 420   double coef;
 421
 422   maxt = lrint(aJot1 + aJot2 - aJot);
 423   mint = 0;
 424   if (lrint(aJot1 - aEm1) < maxt) maxt = lrint(aJot1 - aEm1);
 425   if (lrint(aJot2 + aEm2) < maxt) maxt = lrint(aJot2 + aEm2);
 426   if (lrint(-(aJot-aJot2+aEm1)) > mint) mint = lrint(-(aJot-aJot2+aEm1));
 427   if (lrint(-(aJot-aJot1-aEm2)) > mint) mint = lrint(-(aJot-aJot1-aEm2));
 428
 429   for (titer = mint; titer<=maxt; titer ++)
 430     {
 431       coef = TMath::Power(-1, titer);
 432       coef *= TMath::Sqrt((2*aJot+1)*
 433                           factorials[lrint(aJot1+aEm1)] *
 434                           factorials[lrint(aJot1-aEm1)] *
 435                           factorials[lrint(aJot2+aEm2)] *
 436                           factorials[lrint(aJot2-aEm2)] *
 437                           factorials[lrint(aJot+aEm)] *
 438                           factorials[lrint(aJot-aEm)]);
 439       coef /= (factorials[titer] *
 440                factorials[lrint(aJot1+aJot2-aJot-titer)] *
 441                factorials[lrint(aJot1-aEm1-titer)] *
 442                factorials[lrint(aJot2+aEm2-titer)] *
 443                factorials[lrint(aJot-aJot2+aEm1+titer)] *
 444                factorials[lrint(aJot-aJot1-aEm2+titer)]);
 445
 446       cgc += coef;
 447     }
 448
 449   cgc *= DeltaJ(aJot1, aJot2, aJot);
 450
 451   return cgc;
 452 }
 453
 454 double AliFemtoCorrFctnDirectYlm::DeltaJ(double aJot1, double aJot2, double aJot)
 455 {
 456   // Calculate J for the Clebsh-Gordan coefficient
 457   if ((aJot1+aJot2-aJot) < 0) {
 458     //    cout << "J1+J2-J3 < 0 !!!" << " " << aJot1 << " " << aJot2 << " " << aJot << endl;
 459     return 0;
 460   }
 461   if ((aJot1-aJot2+aJot) < 0) {
 462     //    cout << "J1-J2+J3 < 0 !!!" << " " << aJot1 << " " << aJot2 << " " << aJot << endl;
 463     return 0;
 464   }
 465   if ((-aJot1+aJot2+aJot) < 0) {
 466     //    cout << "-J1+J2+J3 < 0 !!!" << " " << aJot1 << " " << aJot2 << " " << aJot << endl;
 467     return 0;
 468   }
 469   if ((aJot1+aJot2+aJot+1) < 0) {
 470     //    cout << "J1+J2+J3+1 < 0 !!!" << " " << aJot1 << " " << aJot2 << " " << aJot << endl;
 471     return 0;
 472   }
 473   double res = TMath::Sqrt(1.0 *
 474                            factorials[lrint(aJot1+aJot2-aJot)] *
 475                            factorials[lrint(aJot1-aJot2+aJot)] *
 476                            factorials[lrint(-aJot1+aJot2+aJot)] /
 477                            factorials[lrint(aJot1+aJot2+aJot+1)]);
 478
 479   return res;
 480 }
 481
 482 double AliFemtoCorrFctnDirectYlm::WignerSymbol(double aJot1, double aEm1, double aJot2, double aEm2, double aJot, double aEm)
 483 {
 484   // Get Wigner symbol
 485   if (lrint(aEm1+aEm2+aEm) != 0.0)
 486     return 0.0;
 487   double cge = ClebschGordan(aJot1, aEm1, aJot2, aEm2, aJot, -aEm);
 488   if (lrint(abs(aJot1 - aJot2 - aEm)) % 2)
 489     cge *= -1.0;
 490   cge /= sqrt(2*aJot + 1);
 491
 492   if (cge == -0.0) cge = 0.0;
 493
 494   return cge;
 495 }
 496
 497
 498 void AliFemtoCorrFctnDirectYlm::GetMtilde(complex<double> *aMat, double *aMTilde)
 499 {
 500   // Create the Mtilde for a given q bin
 501   double lzero, mzero;
 502   double lprim, mprim;
 503   double lbis, mbis;
 504
 505   int lzeroi, mzeroi;
 506   int lprimi, mprimi;
 507   int lbisi, mbisi;
 508
 509   complex<double> mcomp;
 510
 511   for (int izero = 0; izero<GetMaxJM(); izero++) {
 512     GetElEmForIndex(izero, &lzero, &mzero);
 513     GetElEmForIndex(izero, &lzeroi, &mzeroi);
 514     for (int ibis = 0; ibis<GetMaxJM(); ibis++) {
 515       GetElEmForIndex(ibis, &lbis, &mbis);
 516       GetElEmForIndex(ibis, &lbisi, &mbisi);
 517       complex<double> val = complex<double>(0.0, 0.0);
 518       for (int iprim = 0; iprim<GetMaxJM(); iprim++) {
 519
 520         GetElEmForIndex(iprim, &lprim, &mprim);
 521         GetElEmForIndex(iprim, &lprimi, &mprimi);
 522
 523         if (abs(mzeroi) % 2) mcomp = complex<double>(-1.0, 0.0); // (-1)^m
 524         else mcomp = complex<double>(1.0, 0.0);
 525
 526         mcomp *= sqrt((2*lzero+1)*(2*lprim+1)*(2*lbis+1));   // P1
 527         mcomp *= WignerSymbol(lzero, 0, lprim, 0, lbis, 0); // W1
 528         mcomp *= WignerSymbol(lzero, -mzero, lprim, mprim, lbis, mbis); // W2
 529         mcomp *= aMat[iprim];
 530         val += mcomp;
 531       }
 532       aMTilde[(izero*2)*(2*GetMaxJM()) + (ibis*2)]     =  real(val);
 533       aMTilde[(izero*2+1)*(2*GetMaxJM()) + (ibis*2)]   =  imag(val);
 534       if (imag(val) != 0.0)
 535         aMTilde[(izero*2)*(2*GetMaxJM()) + (ibis*2+1)]   = -imag(val);
 536       else
 537         aMTilde[(izero*2)*(2*GetMaxJM()) + (ibis*2+1)]   = 0.0;
 538       aMTilde[(izero*2+1)*(2*GetMaxJM()) + (ibis*2+1)] =  real(val);
 539
 540     }
 541   }
 542 }
 543
 544 int  AliFemtoCorrFctnDirectYlm::GetMaxJM() const
 545 { return fMaxJM; }
 546
 547 void AliFemtoCorrFctnDirectYlm::GetElEmForIndex(int aIndex, double *aEl, double *aEm) const
 548 {
 549   // Get l,m for a given index
 550   *aEl = fels[aIndex];
 551   *aEm = fems[aIndex];
 552 }
 553
 554 void AliFemtoCorrFctnDirectYlm::GetElEmForIndex(int aIndex, int *aEl, int *aEm) const
 555 {
 556   // Get l,m for a given index
 557   *aEl = felsi[aIndex];
 558   *aEm = femsi[aIndex];
 559 }
 560
 561 int AliFemtoCorrFctnDirectYlm::GetBin(int qbin, int ilmzero, int zeroimag, int ilmprim, int primimag)
 562 {
 563   return (qbin*GetMaxJM()*GetMaxJM()*4 +
 564           (ilmprim*2 + primimag) * GetMaxJM()*2 +
 565           ilmzero*2 + zeroimag);
 566 }
 567
 568 void AliFemtoCorrFctnDirectYlm::AddRealPair(double qout, double qside, double qlong, double weight)
 569 {
 570   // Fill numerator
 571   double kv = sqrt(qout*qout + qside*qside + qlong*qlong);
 572   int nqbin = fbinctn->GetXaxis()->FindFixBin(kv) - 1;
 573
 574   // Use saved ylm values for same qout, qside, qlong
 575   if ((qout != fSout) || (qside != fSside) || (qlong != fSlong)) {
 576     AliFemtoYlm::YlmUpToL(fMaxL, qout, qside, qlong, fYlmBuffer);
 577     fSout = qout; fSside = qside; fSlong = qlong;
 578   }
 579   for (int ilm=0; ilm<GetMaxJM(); ilm++) {
 580     //    fYlmBuffer[ilm] = AliFemtoYlm::Ylm(elsi[ilm], emsi[ilm], qout, qside, qlong);
 581
 582     fnumsreal[ilm]->Fill(kv, real(fYlmBuffer[ilm])*weight);
 583     fnumsimag[ilm]->Fill(kv, -imag(fYlmBuffer[ilm])*weight);
 584
 585     fbinctn->Fill(kv, 1.0);
 586   }
 587
 588   // Fill in the error matrix
 589   //  int tabshift = nqbin*GetMaxJM()*GetMaxJM()*4;
 590   if (nqbin < fbinctn->GetNbinsX())
 591     for (int ilmzero=0; ilmzero<GetMaxJM(); ilmzero++)
 592       for (int ilmprim=0; ilmprim<GetMaxJM(); ilmprim++) {
 593         fcovmnum[GetBin(nqbin, ilmzero, 0, ilmprim, 0)] += real(fYlmBuffer[ilmzero])*real(fYlmBuffer[ilmprim])*weight*weight;
 594         fcovmnum[GetBin(nqbin, ilmzero, 0, ilmprim, 1)] += real(fYlmBuffer[ilmzero])*-imag(fYlmBuffer[ilmprim])*weight*weight;
 595         fcovmnum[GetBin(nqbin, ilmzero, 1, ilmprim, 0)] += -imag(fYlmBuffer[ilmzero])*real(fYlmBuffer[ilmprim])*weight*weight;
 596         fcovmnum[GetBin(nqbin, ilmzero, 1, ilmprim, 1)] += -imag(fYlmBuffer[ilmzero])*-imag(fYlmBuffer[ilmprim])*weight*weight;
 597
 598       }
 599
 600 }
 601
 602 void AliFemtoCorrFctnDirectYlm::AddMixedPair(double qout, double qside, double qlong, double weight)
 603 {
 604   // Fill denominator
 605   double kv = sqrt(qout*qout + qside*qside + qlong*qlong);
 606
 607   // Use saved ylm values for same qout, qside, qlong
 608   if ((qout != fSout) || (qside != fSside) || (qlong != fSlong)) {
 609     AliFemtoYlm::YlmUpToL(fMaxL, qout, qside, qlong, fYlmBuffer);
 610     fSout = qout; fSside = qside; fSlong = qlong;
 611   }
 612   for (int ilm=0; ilm<GetMaxJM(); ilm++) {
 613     //    fYlmBuffer[ilm] = AliFemtoYlm::Ylm(elsi[ilm], emsi[ilm], qout, qside, qlong);
 614
 615     fdensreal[ilm]->Fill(kv, real(fYlmBuffer[ilm])*weight);
 616     fdensimag[ilm]->Fill(kv, -imag(fYlmBuffer[ilm])*weight);
 617
 618     fbinctd->Fill(kv, 1.0);
 619   }
 620
 621   // Fill in the error matrix
 622   int nqbin = fbinctn->GetXaxis()->FindFixBin(kv) - 1;
 623   //  int tabshift = nqbin*GetMaxJM()*GetMaxJM()*4;
 624   if (nqbin < fbinctn->GetNbinsX())
 625     for (int ilmzero=0; ilmzero<GetMaxJM(); ilmzero++)
 626       for (int ilmprim=0; ilmprim<GetMaxJM(); ilmprim++) {
 627         fcovmden[GetBin(nqbin, ilmzero, 0, ilmprim, 0)] += real(fYlmBuffer[ilmzero])*real(fYlmBuffer[ilmprim]);
 628         fcovmden[GetBin(nqbin, ilmzero, 0, ilmprim, 1)] += real(fYlmBuffer[ilmzero])*-imag(fYlmBuffer[ilmprim]);
 629         fcovmden[GetBin(nqbin, ilmzero, 1, ilmprim, 0)] += -imag(fYlmBuffer[ilmzero])*real(fYlmBuffer[ilmprim]);
 630         fcovmden[GetBin(nqbin, ilmzero, 1, ilmprim, 1)] += -imag(fYlmBuffer[ilmzero])*-imag(fYlmBuffer[ilmprim]);
 631
 632     }
 633 }
 634
 635 void AliFemtoCorrFctnDirectYlm::AddRealPair(double *qvec, double weight) {
 636   AddRealPair(qvec[0], qvec[1], qvec[2], weight);
 637 }
 638
 639 void AliFemtoCorrFctnDirectYlm::AddMixedPair(double *qvec, double weight) {
 640   AddMixedPair(qvec[0], qvec[1], qvec[2], weight);
 641 }
 642
 643 void AliFemtoCorrFctnDirectYlm::Finish()
 644 {
 645   PackCovariances();
 646 }
 647
 648 void AliFemtoCorrFctnDirectYlm::Write()
 649 {
 650   // Write out output histograms
 651   for (int ilm=0; ilm<fMaxJM; ilm++) {
 652     fnumsreal[ilm]->Write();
 653     fdensreal[ilm]->Write();
 654     fnumsimag[ilm]->Write();
 655     fdensimag[ilm]->Write();
 656   }
 657   if (fcovnum) fcovnum->Write();
 658   if (fcovden) fcovden->Write();
 659 }
 660
 661 TList* AliFemtoCorrFctnDirectYlm::GetOutputList()
 662 {
 663   // Prepare the list of objects to be written to the output
 664   TList *tOutputList = new TList();
 665
 666   for (int ilm=0; ilm<fMaxJM; ilm++) {
 667     tOutputList->Add(fnumsreal[ilm]);
 668     tOutputList->Add(fdensreal[ilm]);
 669     tOutputList->Add(fnumsimag[ilm]);
 670     tOutputList->Add(fdensimag[ilm]);
 671   }
 672   if (fcovnum) tOutputList->Add(fcovnum);
 673   if (fcovden) tOutputList->Add(fcovden);
 674
 675   return tOutputList;
 676 }
 677
 678
 679 void AliFemtoCorrFctnDirectYlm::ReadFromFile(TFile *infile, const char *name, int maxl)
 680 {
 681   // Raad in the numerator and denominator from file
 682   if (maxl != fMaxL) {
 683     cout << "Cannot read function for L " << maxl << " into a container with L "<< fMaxL << endl;
 684     return;
 685   }
 686   cout << "Reading in numerators and denominators" << endl;
 687
 688   char bufname[200];
 689   for (int ihist=0; ihist<fMaxJM; ihist++) {
 690     sprintf(bufname, "NumReYlm%i%i%s", felsi[ihist], femsi[ihist]<0 ? felsi[ihist]-femsi[ihist] : femsi[ihist], name);
 691     if (fnumsreal[ihist]) delete fnumsreal[ihist];
 692     fnumsreal[ihist] = new TH1D(*((TH1D *) infile->Get(bufname)));
 693
 694     sprintf(bufname, "NumImYlm%i%i%s", felsi[ihist], femsi[ihist]<0 ? felsi[ihist]-femsi[ihist] : femsi[ihist], name);
 695     if (fnumsimag[ihist]) delete fnumsimag[ihist];
 696     fnumsimag[ihist] = new TH1D(*((TH1D *) infile->Get(bufname)));
 697
 698     sprintf(bufname, "DenReYlm%i%i%s", felsi[ihist], femsi[ihist]<0 ? felsi[ihist]-femsi[ihist] : femsi[ihist], name);
 699     if (fdensreal[ihist]) delete fdensreal[ihist];
 700     fdensreal[ihist] = new TH1D(*((TH1D *) infile->Get(bufname)));
 701
 702     sprintf(bufname, "DenImYlm%i%i%s", felsi[ihist], femsi[ihist]<0 ? felsi[ihist]-femsi[ihist] : femsi[ihist], name);
 703     if (fdensimag[ihist]) delete fdensimag[ihist];
 704     fdensimag[ihist] = new TH1D(*((TH1D *) infile->Get(bufname)));
 705   }
 706
 707   if (fcovnum) delete fcovnum;
 708   sprintf(bufname, "covNum%s", name);
 709   fcovnum = new TH3D (*((TH3D *) infile->Get(bufname)));
 710
 711   if (fcovden) delete fcovden;
 712   sprintf(bufname, "CovDen%s", name);
 713   fcovden = new TH3D (*((TH3D *) infile->Get(bufname)));
 714
 715   if ((fcovnum) && (fcovden)) {
 716     cout << "Unpacking covariance matrices from file " << endl;
 717     UnpackCovariances();
 718   }
 719   else {
 720
 721     cout << "Creating fake covariance matrices" << endl;
 722
 723     for (int ibin=1; ibin<=fnumsreal[0]->GetNbinsX(); ibin++) {
 724       double nent = fnumsreal[0]->GetEntries();
 725       double nentd = fdensreal[0]->GetEntries();
 726       for (int ilmx=0; ilmx<GetMaxJM(); ilmx++) {
 727         for (int ilmy=0; ilmy<GetMaxJM(); ilmy++) {
 728           double t1t2rr = fnumsreal[ilmx]->GetBinContent(ibin)*fnumsreal[ilmy]->GetBinContent(ibin)/nent/nent;
 729           double t1t2ri = fnumsreal[ilmx]->GetBinContent(ibin)*fnumsimag[ilmy]->GetBinContent(ibin)/nent/nent;
 730           double t1t2ir = fnumsimag[ilmx]->GetBinContent(ibin)*fnumsreal[ilmy]->GetBinContent(ibin)/nent/nent;
 731           double t1t2ii = fnumsimag[ilmx]->GetBinContent(ibin)*fnumsimag[ilmy]->GetBinContent(ibin)/nent/nent;
 732           if (ilmx == ilmy) {
 733             fcovmnum[GetBin(ibin-1, ilmx, 0, ilmy, 0)] = nent*(TMath::Power(fnumsreal[ilmx]->GetBinError(ibin)/nent,2)*(nent-1) + t1t2rr);
 734             fcovmnum[GetBin(ibin-1, ilmx, 0, ilmy, 1)] = nent*t1t2ri;
 735             fcovmnum[GetBin(ibin-1, ilmx, 1, ilmy, 0)] = nent*t1t2ir;
 736             fcovmnum[GetBin(ibin-1, ilmx, 1, ilmy, 1)] = nent*(TMath::Power(fnumsimag[ilmx]->GetBinError(ibin)/nent,2)*(nent-1) + t1t2rr);
 737           }
 738           else {
 739             fcovmnum[GetBin(ibin-1, ilmx, 0, ilmy, 0)] = nent*t1t2rr;
 740             fcovmnum[GetBin(ibin-1, ilmx, 0, ilmy, 1)] = nent*t1t2ri;
 741             fcovmnum[GetBin(ibin-1, ilmx, 1, ilmy, 0)] = nent*t1t2ir;
 742             fcovmnum[GetBin(ibin-1, ilmx, 1, ilmy, 1)] = nent*t1t2ii;
 743           }
 744           t1t2rr = fdensreal[ilmx]->GetBinContent(ibin)*fdensreal[ilmy]->GetBinContent(ibin)/nentd/nentd;
 745           t1t2ri = fdensreal[ilmx]->GetBinContent(ibin)*fdensimag[ilmy]->GetBinContent(ibin)/nentd/nentd;
 746           t1t2ir = fdensimag[ilmx]->GetBinContent(ibin)*fdensreal[ilmy]->GetBinContent(ibin)/nentd/nentd;
 747           t1t2ii = fdensimag[ilmx]->GetBinContent(ibin)*fdensimag[ilmy]->GetBinContent(ibin)/nentd/nentd;
 748
 749           fcovmden[GetBin(ibin-1, ilmx, 0, ilmy, 0)] = nentd*t1t2rr;
 750           fcovmden[GetBin(ibin-1, ilmx, 0, ilmy, 1)] = nentd*t1t2ri;
 751           fcovmden[GetBin(ibin-1, ilmx, 1, ilmy, 0)] = nentd*t1t2ir;
 752           fcovmden[GetBin(ibin-1, ilmx, 1, ilmy, 1)] = nentd*t1t2ii;
 753         }
 754       }
 755     }
 756   }
 757
 758   // Recalculating the correlation functions
 759   Finish();
 760 }
 761
 762 int AliFemtoCorrFctnDirectYlm::PackYlmVector(const double *invec, double *outvec)
 763 {
 764   // Pack a vector in l,m into an array using
 765   // only independent components
 766   int ioutcount = 0;
 767   int em, el;
 768   for (int ilm=0; ilm<GetMaxJM(); ilm++) {
 769     GetElEmForIndex(ilm, &el, &em);
 770     outvec[ioutcount++] = invec[ilm*2];
 771     if (em == 0)
 772       continue;
 773     outvec[ioutcount++] = invec[ilm*2 + 1];
 774   }
 775
 776   return ioutcount;
 777 }
 778
 779 int AliFemtoCorrFctnDirectYlm::PackYlmMatrix(const double *inmat, double *outmat)
 780 {
 781   // Pack a matrix in l,m x l,m into an array using
 782   // only independent components
 783   int ioutcountz = 0;
 784   int ioutcountp = 0;
 785   int emz, elz;
 786   int emp, elp;
 787   int finalsize = 0;
 788
 789   for (int ilm=0; ilm<GetMaxJM(); ilm++) {
 790     GetElEmForIndex(ilm, &elz, &emz);
 791     finalsize++;
 792     if (emz == 0) continue;
 793     finalsize++;
 794   }
 795
 796   for (int ilmz=0; ilmz<GetMaxJM(); ilmz++) {
 797     GetElEmForIndex(ilmz, &elz, &emz);
 798     ioutcountp = 0;
 799     for (int ilmp=0; ilmp<GetMaxJM(); ilmp++) {
 800       GetElEmForIndex(ilmp, &elp, &emp);
 801       outmat[ioutcountz*finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 0, ilmp, 0)];
 802       ioutcountp++;
 803       if (emp == 0) continue;
 804       outmat[ioutcountz*finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 0, ilmp, 1)];
 805       ioutcountp++;
 806     }
 807     ioutcountz++;
 808
 809     if (emz == 0) continue;
 810     ioutcountp = 0;
 811     for (int ilmp=0; ilmp<GetMaxJM(); ilmp++) {
 812       GetElEmForIndex(ilmp, &elp, &emp);
 813       outmat[ioutcountz*finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 1, ilmp, 0)];
 814       ioutcountp++;
 815       if (emp == 0) continue;
 816       outmat[ioutcountz*finalsize + ioutcountp] = inmat[GetBin(0, ilmz, 1, ilmp, 1)];
 817       ioutcountp++;
 818     }
 819     ioutcountz++;
 820   }
 821
 822   return ioutcountz;
 823 }
 824
 825 void AliFemtoCorrFctnDirectYlm::PackCovariances()
 826 {
 827   // Migrate the covariance matrix into a 3D histogram for storage
 828   char bufname[200];
 829   sprintf(bufname, "CovNum%s", fnumsreal[0]->GetName()+10);
 830
 831   if (fcovnum) delete fcovnum;
 832   fcovnum = new TH3D(bufname,bufname,
 833                     fnumsreal[0]->GetNbinsX(), fnumsreal[0]->GetXaxis()->GetXmin(), fnumsreal[0]->GetXaxis()->GetXmax(),
 834                     GetMaxJM()*2, -0.5, GetMaxJM()*2 - 0.5,
 835                     GetMaxJM()*2, -0.5, GetMaxJM()*2 - 0.5);
 836
 837   for (int ibin=1; ibin<=fcovnum->GetNbinsX(); ibin++)
 838     for (int ilmz=0; ilmz<GetMaxJM()*2; ilmz++)
 839       for (int ilmp=0; ilmp<GetMaxJM()*2; ilmp++)
 840         fcovnum->SetBinContent(ibin, ilmz+1, ilmp+1, fcovmnum[GetBin(ibin-1, ilmz/2, ilmz%2, ilmp/2, ilmp%2)]);
 841
 842   if (fcovden) delete fcovden;
 843   sprintf(bufname, "CovDen%s", fnumsreal[0]->GetName()+10);
 844   fcovden  = new TH3D(bufname,bufname,
 845                      fdensreal[0]->GetNbinsX(), fdensreal[0]->GetXaxis()->GetXmin(), fdensreal[0]->GetXaxis()->GetXmax(),
 846                      GetMaxJM()*2, -0.5, GetMaxJM()*2 - 0.5,
 847                      GetMaxJM()*2, -0.5, GetMaxJM()*2 - 0.5);
 848
 849   for (int ibin=1; ibin<=fcovden->GetNbinsX(); ibin++)
 850     for (int ilmz=0; ilmz<GetMaxJM()*2; ilmz++)
 851       for (int ilmp=0; ilmp<GetMaxJM()*2; ilmp++)
 852         fcovden->SetBinContent(ibin, ilmz+1, ilmp+1, fcovmden[GetBin(ibin-1, ilmz/2, ilmz%2, ilmp/2, ilmp%2)]);
 853
 854 }
 855
 856 void AliFemtoCorrFctnDirectYlm::UnpackCovariances()
 857 {
 858   // Extract the covariance matrices from storage
 859   if (fcovnum) {
 860     for (int ibin=1; ibin<=fcovnum->GetNbinsX(); ibin++)
 861       for (int ilmz=0; ilmz<GetMaxJM()*2; ilmz++)
 862         for (int ilmp=0; ilmp<GetMaxJM()*2; ilmp++)
 863           fcovmnum[GetBin(ibin-1, ilmz/2, ilmz%2, ilmp/2, ilmp%2)] = fcovnum->GetBinContent(ibin, ilmz+1, ilmp+1);
 864
 865   }
 866   if (fcovden) {
 867     for (int ibin=1; ibin<=fcovden->GetNbinsX(); ibin++)
 868       for (int ilmz=0; ilmz<GetMaxJM()*2; ilmz++)
 869         for (int ilmp=0; ilmp<GetMaxJM()*2; ilmp++)
 870           fcovmden[GetBin(ibin-1, ilmz/2, ilmz%2, ilmp/2, ilmp%2)] = fcovden->GetBinContent(ibin, ilmz+1, ilmp+1);
 871   }
 872 }
 873
 874 int AliFemtoCorrFctnDirectYlm::GetIndexForLM(int el, int em) const
 875 {
 876   // Get array index for a given l,m
 877   for (int iter=0; iter<fMaxJM; iter++)
 878     if ((el == felsi[iter]) && (em == femsi[iter]))
 879       return iter;
 880   return -1;
 881 }
 882
 883 TH1D *AliFemtoCorrFctnDirectYlm::GetNumRealHist(int el, int em)
 884 {
 885   // Get numerator hist for a given l,m
 886   if (GetIndexForLM(el, em)>=0)
 887     return fnumsreal[GetIndexForLM(el, em)];
 888   else
 889     return 0;
 890 }
 891 TH1D *AliFemtoCorrFctnDirectYlm::GetNumImagHist(int el, int em)
 892 {
 893   // Get numerator hist for a given l,m
 894   if (GetIndexForLM(el, em)>=0)
 895     return fnumsimag[GetIndexForLM(el, em)];
 896   else
 897     return 0;
 898 }
 899
 900 TH1D *AliFemtoCorrFctnDirectYlm::GetDenRealHist(int el, int em)
 901 {
 902   // Get denominator hist for a given l,m
 903   if (GetIndexForLM(el, em)>=0)
 904     return fdensreal[GetIndexForLM(el, em)];
 905   else
 906     return 0;
 907 }
 908 TH1D *AliFemtoCorrFctnDirectYlm::GetDenImagHist(int el, int em)
 909 {
 910   // Get denominator hist for a given l,m
 911   if (GetIndexForLM(el, em)>=0)
 912     return fdensimag[GetIndexForLM(el, em)];
 913   else
 914     return 0;
 915 }
 916
 917 AliFemtoString AliFemtoCorrFctnDirectYlm::Report()
 918 {
 919   return "AliFemtoCorrFctnDirectYlm::Finish";
 920 }
 921
 922 void AliFemtoCorrFctnDirectYlm::AddRealPair(AliFemtoPair* aPair)
 923 {
 924   AddRealPair(aPair->QOutPf(), aPair->QSidePf(), aPair->QLongPf(), 1.0);
 925 }
 926 void AliFemtoCorrFctnDirectYlm::AddMixedPair(AliFemtoPair* aPair)
 927 {
 928   AddMixedPair(aPair->QOutPf(), aPair->QSidePf(), aPair->QLongPf(), 1.0);
 929 }
 930