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