]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDirectYlm.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Herwig event header.
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoCorrFctnDirectYlm.cxx
CommitLineData
17b6dded 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>
17b6dded 13#include <iostream>
14
15using namespace std;
16
963893bf 17AliFemtoCorrFctnDirectYlm::AliFemtoCorrFctnDirectYlm(const char *name, int maxl, int ibin=30, double vmin=0.0, double vmax=0.3, int aUseLCMS=0):
17b6dded 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),
963893bf 38 fSlong(0.0),
39 fUseLCMS(aUseLCMS)
17b6dded 40{
41 // Main constructor
42 fMaxL = maxl;
43 fMaxJM = (maxl+1)*(maxl+1);
44
e82b5d68 45 // *DEB* cout << "Size is " << sizeof(double) << " " << sizeof(complex<double>) << endl;
17b6dded 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
e82b5d68 71 // *DEB* cout << "il el em " << il << " " << felsi[il] << " " << femsi[il] << endl;
17b6dded 72 em++;
73 il++;
74 if (em > el) {
75 el++;
76 em = -el;
77 }
78 }
79 while (el <= maxl);
80
e82b5d68 81 // *DEB* for (il=0; il<fMaxJM; il++)
82 // *DEB* cout << "il el em " << il << " " << felsi[il] << " " << femsi[il] << endl;
17b6dded 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
132AliFemtoCorrFctnDirectYlm::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),
963893bf 153 fSlong(0.0),
154 fUseLCMS(0)
17b6dded 155{
156 // Default constructor
157 AliFemtoCorrFctnDirectYlm("AliFemtoCorrFctnDirectYlm",2);
158}
159
160AliFemtoCorrFctnDirectYlm::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),
963893bf 182 fSlong(0.0),
183 fUseLCMS(0)
17b6dded 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;
44c6b6dc 277
278 if (aCorrFctn.fPairCut)
279 fPairCut = aCorrFctn.fPairCut;
963893bf 280
281 fUseLCMS = aCorrFctn.fUseLCMS;
17b6dded 282}
283
284AliFemtoCorrFctnDirectYlm& 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
44c6b6dc 382 if (aCorrFctn.fPairCut)
383 fPairCut = aCorrFctn.fPairCut;
384
963893bf 385 fUseLCMS = aCorrFctn.fUseLCMS;
386
17b6dded 387 return *this;
388}
389
390AliFemtoCorrFctnDirectYlm::~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;
44c6b6dc 425
426 if (fPairCut) delete fPairCut;
17b6dded 427}
428
429double 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
469double 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
497double 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
513void 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
559int AliFemtoCorrFctnDirectYlm::GetMaxJM() const
560{ return fMaxJM; }
561
562void 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
569void 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
576int 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
583void 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
617void 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
650void AliFemtoCorrFctnDirectYlm::AddRealPair(double *qvec, double weight) {
651 AddRealPair(qvec[0], qvec[1], qvec[2], weight);
652}
653
654void AliFemtoCorrFctnDirectYlm::AddMixedPair(double *qvec, double weight) {
655 AddMixedPair(qvec[0], qvec[1], qvec[2], weight);
656}
657
658void AliFemtoCorrFctnDirectYlm::Finish()
659{
660 PackCovariances();
661}
662
663void AliFemtoCorrFctnDirectYlm::Write()
664{
665 // Write out output histograms
1fb889b2 666 if ((!fcovnum) || (!fcovden))
667 PackCovariances();
668
17b6dded 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
679TList* 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
697void 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
780int 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
797int 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
843void 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
874void 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
892int 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
901TH1D *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}
909TH1D *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
918TH1D *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}
926TH1D *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
935AliFemtoString AliFemtoCorrFctnDirectYlm::Report()
936{
937 return "AliFemtoCorrFctnDirectYlm::Finish";
938}
939
940void AliFemtoCorrFctnDirectYlm::AddRealPair(AliFemtoPair* aPair)
941{
963893bf 942 // Fill in the numerator
44c6b6dc 943 if (fPairCut)
944 if (!fPairCut->Pass(aPair)) return;
945
963893bf 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);
17b6dded 950}
963893bf 951
17b6dded 952void AliFemtoCorrFctnDirectYlm::AddMixedPair(AliFemtoPair* aPair)
953{
963893bf 954 // Fill in the denominator
44c6b6dc 955 if (fPairCut)
956 if (!fPairCut->Pass(aPair)) return;
957
963893bf 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
964void AliFemtoCorrFctnDirectYlm::SetUseLCMS(int aUseLCMS)
965{
966 fUseLCMS = aUseLCMS;
967}
968
969int AliFemtoCorrFctnDirectYlm::GetUseLCMS()
970{
971 return fUseLCMS;
17b6dded 972}
973