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