]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/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
In AOD Event Reader: TOF mismatch with new method + macros for proton femtoscopy...
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemtoUser / AliFemtoCorrFctnDirectYlm.cxx
CommitLineData
76ce4b5b 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
15using namespace std;
16
17AliFemtoCorrFctnDirectYlm::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
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),
153 fSlong(0.0),
154 fUseLCMS(0)
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),
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
286AliFemtoCorrFctnDirectYlm& 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
412AliFemtoCorrFctnDirectYlm::~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
451double 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
491double 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
519double 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
535void 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
581int AliFemtoCorrFctnDirectYlm::GetMaxJM() const
582{ return fMaxJM; }
583
584void 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
591void 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
598int 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
605void 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
639void 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
672void AliFemtoCorrFctnDirectYlm::AddRealPair(double *qvec, double weight) {
673 AddRealPair(qvec[0], qvec[1], qvec[2], weight);
674}
675
676void AliFemtoCorrFctnDirectYlm::AddMixedPair(double *qvec, double weight) {
677 AddMixedPair(qvec[0], qvec[1], qvec[2], weight);
678}
679
680void AliFemtoCorrFctnDirectYlm::Finish()
681{
682 PackCovariances();
683}
684
685void 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
701TList* 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
722void 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
805int 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
822int 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
868void 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
902void 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
920int 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
929TH1D *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}
937TH1D *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
946TH1D *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}
954TH1D *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
963AliFemtoString AliFemtoCorrFctnDirectYlm::Report()
964{
965 return "AliFemtoCorrFctnDirectYlm::Finish";
966}
967
968void 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
980void 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
992void AliFemtoCorrFctnDirectYlm::SetUseLCMS(int aUseLCMS)
993{
994 fUseLCMS = aUseLCMS;
995}
996
997int AliFemtoCorrFctnDirectYlm::GetUseLCMS()
998{
999 return fUseLCMS;
1000}
1001