Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
FemtoYlmSolver.cxx
1/*
2 * FemtoYlmSolver.cxx
3 *
4 * Created on: 24 sie 2022
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9#include "FemtoYlmSolver.h"
10
11#include "Cout.h"
12#include "FemtoSHCF.h"
13#include "FemtoSHSlice.h"
14#include <TH1.h>
15#include <TMath.h>
16#include <gsl/gsl_blas.h>
17#include <gsl/gsl_cblas.h>
18#include <gsl/gsl_matrix_double.h>
19#include <gsl/gsl_vector_double.h>
20#include <iostream>
21
22//#define _FINISH_DEBUG_
23namespace Hal {
24
25 FemtoYlmSolver::FemtoYlmSolver() {}
26
27 FemtoYlmSolver::~FemtoYlmSolver() {}
28
29 void FemtoYlmSolver::GetMtilde(std::complex<double>* aMat, double* aMTilde) {
30 // Create the Mtilde for a given q bin
31 double lzero, mzero;
32 double lprim, mprim;
33 double lbis, mbis;
34
35 int lzeroi, mzeroi;
36 int lprimi, mprimi;
37 int lbisi, mbisi;
38
39 for (int iz = 0; iz < GetMaxJM() * 2; iz++)
40 for (int ip = 0; ip < GetMaxJM() * 2; ip++)
41 aMTilde[iz * GetMaxJM() * 2 + ip] = 0.0;
42
43 for (int izero = 0; izero < GetMaxJM(); izero++) {
44 GetElEmForIndex(izero, lzero, mzero);
45 GetElEmForIndex(izero, lzeroi, mzeroi);
46 // if (mzero < 0)
47 // continue;
48 for (int ibis = 0; ibis < GetMaxJM(); ibis++) {
49 GetElEmForIndex(ibis, lbis, mbis);
50 GetElEmForIndex(ibis, lbisi, mbisi);
51
52 // if (mbis<0) continue;
53
54 std::complex<double> val = std::complex<double>(0.0, 0.0);
55 std::unique_ptr<std::complex<double>[]> mcomp(new std::complex<double>[fMaxJM]());
56 for (int iprim = 0; iprim < GetMaxJM(); iprim++) {
57 GetElEmForIndex(iprim, lprim, mprim);
58 GetElEmForIndex(iprim, lprimi, mprimi);
59
60 // if (mprim < 0 ) continue;
61
62 if (abs(mzeroi) % 2)
63 mcomp[iprim] = std::complex<double>(-1.0, 0.0); // (-1)^m
64 else
65 mcomp[iprim] = std::complex<double>(1.0, 0.0);
66
67 mcomp[iprim] *= sqrt((2 * lzero + 1) * (2 * lprim + 1) * (2 * lbis + 1)); // P1
68 mcomp[iprim] *= fLmMath.WignerSymbol(lzero, 0, lprim, 0, lbis, 0); // W1
69 mcomp[iprim] *= fLmMath.WignerSymbol(lzero, -mzero, lprim, mprim, lbis,
70 mbis); // W2
71 mcomp[iprim] *= aMat[iprim];
72 // if (
73 val += mcomp[iprim];
74 }
75
76 aMTilde[(izero * 2) * (2 * GetMaxJM()) + (ibis * 2)] = real(val);
77 aMTilde[(izero * 2 + 1) * (2 * GetMaxJM()) + (ibis * 2)] = imag(val);
78 if (imag(val) != 0.0)
79 aMTilde[(izero * 2) * (2 * GetMaxJM()) + (ibis * 2 + 1)] = -imag(val);
80 else
81 aMTilde[(izero * 2) * (2 * GetMaxJM()) + (ibis * 2 + 1)] = 0.0;
82 aMTilde[(izero * 2 + 1) * (2 * GetMaxJM()) + (ibis * 2 + 1)] = real(val);
83 }
84 }
85 }
86
87 void FemtoYlmSolver::InvertYlmIndependentMatrix(double* inmat, double* outmat) const {
88 // Invert the Ylm matrix by inverting only the matrix
89 // with independent elements and filling in the rest
90 // according to sign rules
91
92 std::unique_ptr<double[]> mU(new double[GetMaxJM() * GetMaxJM() * 4]);
93 int isize = PackYlmMatrixIndependentOnly(inmat, mU.get());
94 // cout << "Independent count " << isize << std::endl;
95#ifndef DISABLE_GSL
96 gsl_matrix_view matU = gsl_matrix_view_array(mU.get(), isize, isize);
97#endif
98 // Identity matrix helper for inversion
99 std::unique_ptr<double[]> mI(new double[GetMaxJM() * GetMaxJM() * 4]);
100 for (int iterm = 0; iterm < isize; iterm++)
101 for (int iterp = 0; iterp < isize; iterp++)
102 if (iterm == iterp)
103 mI[iterm * isize + iterp] = 1.0;
104 else
105 mI[iterm * isize + iterp] = 0.0;
106#ifndef DISABLE_GSL
107 gsl_matrix_view matI = gsl_matrix_view_array(mI.get(), isize, isize);
108 // Invert the matrix
109 gsl_blas_dtrsm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, &matU.matrix, &matI.matrix);
110#else
111 std::cout << "GLS NOT ENABLED!" << std::endl;
112#endif
113 UnPackYlmMatrixIndependentOnly(mI.get(), outmat, isize);
114 }
115
116 void FemtoYlmSolver::UnPackYlmMatrixIndependentOnly(double* inmat, double* outmat, int insize) const {
117 int lmax = sqrt(insize) - 1;
118 // cout << "lmax is " << lmax << std::endl;
119 if (0) { lmax *= 2; }
120 int tmax = (lmax + 1) * (lmax + 1) * 2;
121 std::unique_ptr<int[]> indexfrom(new int[tmax]);
122 std::unique_ptr<int[]> multfrom(new int[tmax]);
123
124 int el, em;
125 for (int iter = 0; iter < tmax; iter++) {
126 int im = iter % 2;
127 GetElEmForIndex(iter / 2, el, em);
128 if (em == 0) {
129 if (im == 1) {
130 indexfrom[iter] = 0;
131 multfrom[iter] = 0;
132 } else {
133 indexfrom[iter] = el * el;
134 multfrom[iter] = 1;
135 }
136 } else if (em < 0) {
137 indexfrom[iter] = (el * el) + (-em) * 2 - 1;
138 if (im) indexfrom[iter]++;
139 if ((-em) % 2)
140 if (im)
141 multfrom[iter] = 1;
142 else
143 multfrom[iter] = -1;
144 else if (im)
145 multfrom[iter] = -1;
146 else
147 multfrom[iter] = 1;
148 } else if (em > 0) {
149 indexfrom[iter] = (el * el) + (em) *2 - 1;
150 if (im) indexfrom[iter]++;
151 multfrom[iter] = 1;
152 }
153 }
154
155 // cout << "From Mult " << std::endl;
156 // for (int iter=0; iter<tmax; iter++)
157 // cout << indexfrom[iter] << " ";
158 // cout << std::endl;
159 // for (int iter=0; iter<tmax; iter++)
160 // cout << multfrom[iter] << " ";
161 // cout << std::endl;
162
163 for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++)
164 for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++)
165 outmat[ilmz * GetMaxJM() * 2 + ilmp] =
166 inmat[(indexfrom[ilmz] * insize) + indexfrom[ilmp]] * multfrom[ilmz] * multfrom[ilmp];
167 }
168
169 int FemtoYlmSolver::PackYlmVectorIndependentOnly(double* invec, double* outvec) const {
170 int ioutcount = 0;
171 int em, el;
172 for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
173 GetElEmForIndex(ilm, el, em);
174 if (em < 0) continue; // KURWA
175 outvec[ioutcount++] = invec[ilm * 2];
176 if (em == 0) continue;
177 outvec[ioutcount++] = invec[ilm * 2 + 1];
178 }
179 return ioutcount;
180 }
181
182 int FemtoYlmSolver::PackYlmMatrixIndependentOnly(double* inmat, double* outmat) const {
183 int ioutcountz = 0;
184 int ioutcountp = 0;
185 int emz, elz;
186 int emp, elp;
187 int finalsize = 0;
188
189 for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
190 GetElEmForIndex(ilm, elz, emz);
191 if (emz < 0) continue;
192 finalsize++;
193 if (emz == 0) continue;
194 finalsize++;
195 }
196
197 // cout << "Final size " << finalsize << std::endl;
198 auto GetBin = [](int maxJm, int ilmzero, int zeroimag, int ilmprim, int primimag) -> int {
199 return ((ilmprim * 2 + primimag) * maxJm * 2 + ilmzero * 2 + zeroimag);
200 };
201
202 for (int ilmz = 0; ilmz < GetMaxJM(); ilmz++) {
203 GetElEmForIndex(ilmz, elz, emz);
204 ioutcountp = 0;
205
206 if (emz < 0) continue;
207 for (int ilmp = 0; ilmp < GetMaxJM(); ilmp++) {
208 GetElEmForIndex(ilmp, elp, emp);
209 if (emp < 0) continue;
210 outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(fMaxJM, ilmz, 0, ilmp, 0)];
211 ioutcountp++;
212 if (emp == 0) continue;
213 outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(fMaxJM, ilmz, 0, ilmp, 1)];
214 ioutcountp++;
215 }
216 ioutcountz++;
217
218 if (emz == 0) continue;
219 ioutcountp = 0;
220 for (int ilmp = 0; ilmp < GetMaxJM(); ilmp++) {
221 GetElEmForIndex(ilmp, elp, emp);
222 if (emp < 0) continue;
223 outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(fMaxJM, ilmz, 1, ilmp, 0)];
224 ioutcountp++;
225 if (emp == 0) continue;
226 outmat[ioutcountz * finalsize + ioutcountp] = inmat[GetBin(fMaxJM, ilmz, 1, ilmp, 1)];
227 ioutcountp++;
228 }
229 ioutcountz++;
230 }
231
232 return ioutcountz;
233 }
234
235 FemtoYlmSolver::FemtoYlmSolver(Int_t maxL, FemtoSHCF* cf, Bool_t sumw) :
236 fMaxJM((maxL + 1) * (maxL + 1)),
237 fFactorialsSize((maxL + 1) * 4),
238 fMaxL(maxL),
239 fSumw(sumw),
240 fSlice(FemtoSHSlice(maxL)),
241 fCF(cf) {
242 fLmVals.Resize(maxL);
243 fFactorials.resize(fFactorialsSize);
244 fMaxJM2_4 = 4.0 * fMaxJM * fMaxJM;
245
246 if (fSumw) {
247 fNumEnt = cf->GetNumRe(0, 0)->GetEffectiveEntries(); // KURWA
248 fDenEnt = cf->GetDenRe(0, 0)->GetEffectiveEntries(); // KURWA
249 } else {
250 fNumEnt = cf->GetNumRe(0, 0)->GetEntries();
251 fDenEnt = cf->GetDenRe(0, 0)->GetEntries();
252 }
253
254
255 Double_t fac = 1;
256 fFactorials[0] = 1;
257 for (int iter = 1; iter < fFactorialsSize; iter++) {
258 fac *= iter;
259 fFactorials[iter] = fac;
260 }
261 }
262
263 void FemtoYlmSolver::DoMath(Bool_t recalc) {
264
265 std::unique_ptr<std::complex<double>[]> tMq0(new std::complex<double>[fMaxJM]);
266 std::unique_ptr<std::complex<double>[]> tTq0(new std::complex<double>[fMaxJM]);
267 std::unique_ptr<double[]> tMTilde(new double[fMaxJM * fMaxJM * 4]);
268
269 for (int i = 0; i < fMaxJM; i++) {
270 fSlice.fCFImag[i] = 0;
271 fSlice.fCFReal[i] = 0;
272 fSlice.fCFImagE[i] = 0;
273 fSlice.fCFRealE[i] = 0;
274 }
275
276 for (int i = 0; i < fMaxJM * 2; i++)
277 for (int j = 0; j < fMaxJM * 2; j++)
278 fSlice.fCovCF[i][j] = 0;
279 if (fSlice.fCovNum[0][0] == 0) MakeFakeCovMatrix();
280 if (fDebug) {
281 std::cout << "NORM " << fNormFactor << std::endl;
282 std::cout << "RECALC " << recalc << std::endl;
283
284 std::cout << "Nu ent " << fCF->fNumReal[0]->GetEntries() << " " << fCF->fNumReal[0]->GetBinContent(1) << std::endl;
285 }
286 for (int ilm = 0; ilm < fMaxJM; ilm++) {
287 // cout << fNumImag[ilm]->GetBinContent(ibin) << std::endl;
288 if (recalc) {
289 // fNumEnt = 5.76589e+08; // KURWA
290 tMq0[ilm] =
291 std::complex<double>(fSlice.fDenReal[ilm] / (fDenEnt / fNormFactor), fSlice.fDenImag[ilm] / (fDenEnt / fNormFactor));
292 // std::cout << "tMQ" << ilm << "\t" << fSlice.fDenReal[ilm] << " " << fDenEnt << " " << fSlice.fDenImag[ilm] << " "
293 // << fDenEnt << std::endl;
294 tTq0[ilm] = std::complex<double>(fSlice.fNumReal[ilm] / fNumEnt, fSlice.fNumImag[ilm] / fNumEnt);
295 } else {
296 tMq0[ilm] = std::complex<double>(fSlice.fDenReal[ilm] / fNormFactor, fSlice.fDenImag[ilm] / fNormFactor);
297 tTq0[ilm] = std::complex<double>(fSlice.fNumReal[ilm], fSlice.fNumImag[ilm]);
298 }
299 // cout << imag(tTq0[ilm]) << std::endl;
300 }
301
302 if (fDebug) {
303 std::cout << "TT ";
304 for (int ilm = 0; ilm < fMaxJM; ilm++) {
305 std::cout << tTq0[ilm] << " " << std::endl;
306 }
307 std::cout << std::endl;
308 std::cout << "NN ";
309 for (int ilm = 0; ilm < fMaxJM; ilm++) {
310 std::cout << tMq0[ilm] << " " << std::endl;
311 }
312 std::cout << std::endl;
313 }
314 // Calculate the proper error matrix for T
315 // from the temporary covariance matrices
316 // int tabshift = (ibin-1)*GetMaxJM()*GetMaxJM()*4;
317 if (recalc) {
318 if (fDebug) { /*
319 std::cout << "COVNUM" << std::endl;
320 for (int ilmzero = 0; ilmzero < GetMaxJM() * 2; ilmzero++) {
321 for (int ilmprim = 0; ilmprim < GetMaxJM() * 2; ilmprim++) {
322 std::cout << fSlice.fCovNum[ilmzero][ilmprim] << " ";
323 }
324 std::cout << std::endl;
325 }
326 std::cout << std::endl;*/
327 }
328
329 for (int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++) {
330 const int ilmzero2 = ilmzero * 2;
331 for (int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
332 const int ilmprim2 = ilmprim * 2;
333 if (std::isnan(fSlice.fCovNum[ilmzero2][ilmprim2])) {}
334 if (std::isnan(fSlice.fCovNum[ilmzero2][ilmprim2 + 1])) {}
335 if (std::isnan(fSlice.fCovNum[ilmzero2 + 1][ilmprim2])) {}
336 if (std::isnan(fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1])) {}
337 double a, b, c, d;
338 a = fSlice.fCovNum[ilmzero2][ilmprim2];
339 b = fSlice.fCovNum[ilmzero2][ilmprim2 + 1];
340 c = fSlice.fCovNum[ilmzero2 + 1][ilmprim2];
341 d = fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1];
342 if (fDebug) std::cout << ">" << a << " " << b << " " << c << " " << d << " " << fNumEnt << std::endl;
343
344 fSlice.fCovNum[ilmzero2][ilmprim2] /= fNumEnt;
345 fSlice.fCovNum[ilmzero2][ilmprim2 + 1] /= fNumEnt;
346 fSlice.fCovNum[ilmzero2 + 1][ilmprim2] /= fNumEnt;
347 fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1] /= fNumEnt;
348
349 a = fSlice.fCovNum[ilmzero2][ilmprim2];
350 b = fSlice.fCovNum[ilmzero2][ilmprim2 + 1];
351 c = fSlice.fCovNum[ilmzero2 + 1][ilmprim2];
352 d = fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1];
353 if (fDebug) std::cout << "->" << a << " " << b << " " << c << " " << d << std::endl;
354
355
356 fSlice.fCovNum[ilmzero2][ilmprim2] -= real(tTq0[ilmzero]) * real(tTq0[ilmprim]);
357 fSlice.fCovNum[ilmzero2][ilmprim2 + 1] -= real(tTq0[ilmzero]) * imag(tTq0[ilmprim]);
358 fSlice.fCovNum[ilmzero2 + 1][ilmprim2] -= imag(tTq0[ilmzero]) * real(tTq0[ilmprim]);
359 fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1] -= imag(tTq0[ilmzero]) * imag(tTq0[ilmprim]);
360 a = fSlice.fCovNum[ilmzero2][ilmprim2];
361 b = fSlice.fCovNum[ilmzero2][ilmprim2 + 1];
362 c = fSlice.fCovNum[ilmzero2 + 1][ilmprim2];
363 d = fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1];
364 if (fDebug) std::cout << "-->" << a << " " << b << " " << c << " " << d << std::endl;
365 fSlice.fCovNum[ilmzero2][ilmprim2] /= ((fNumEnt) -1);
366 fSlice.fCovNum[ilmzero2][ilmprim2 + 1] /= ((fNumEnt) -1);
367 fSlice.fCovNum[ilmzero2 + 1][ilmprim2] /= ((fNumEnt) -1);
368 fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1] /= ((fNumEnt) -1);
369 a = fSlice.fCovNum[ilmzero2][ilmprim2];
370 b = fSlice.fCovNum[ilmzero2][ilmprim2 + 1];
371 c = fSlice.fCovNum[ilmzero2 + 1][ilmprim2];
372 d = fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1];
373 if (fDebug) std::cout << "--->" << a << " " << b << " " << c << " " << d << std::endl;
374 }
375 }
376 }
377
378 if (fDebug) {
379 std::cout << "COVNUM" << std::endl;
380 for (int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++) {
381 const int ilmzero2 = ilmzero * 2;
382 for (int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
383 const int ilmprim2 = ilmprim * 2;
384 std::cout << fSlice.fCovNum[ilmzero2 + 0][ilmprim2 + 0] << " ";
385 std::cout << fSlice.fCovNum[ilmzero2 + 0][ilmprim2 + 1] << " ";
386 std::cout << fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 0] << " ";
387 std::cout << fSlice.fCovNum[ilmzero2 + 1][ilmprim2 + 1] << " ";
388 std::cout << std::endl;
389 }
390 std::cout << std::endl;
391 }
392 }
393
394
395 GetMtilde(tMq0.get(), tMTilde.get());
396
397 // Perform the solution for the correlation function itself and the errors
398 // cout << "=============================" << std::endl;
399 // cout << "C calculation for bin " << (ibin-1) << std::endl;
400 // cout << std::endl;
401 // cout << "Input: " << std::endl;
402 // cout << "T vector " << std::endl;
403 // for (int ilm=0; ilm<GetMaxJM(); ilm++)
404 // cout << real(tTq0[ilm]) << " " << imag(tTq0[ilm]) << " ";
405 // cout << std::endl << "M vector " << std::endl;
406 // for (int ilm=0; ilm<GetMaxJM(); ilm++)
407 // cout << real(tMq0[ilm]) << " " << imag(tMq0[ilm]) << " ";
408 // cout << std::endl;
409
410 if (fSlice.fNumReal[0] <= 0) {
411 for (int ilm = 0; ilm < GetMaxJM(); ilm++) {
412 fSlice.fCFReal[ilm] = 0;
413 fSlice.fCFRealE[ilm] = 0;
414 fSlice.fCFImag[ilm] = 0;
415 fSlice.fCFImagE[ilm] = 0;
416 }
417
418 for (int ilmz = 0; ilmz < GetMaxJM(); ilmz++) { // bylo *2
419 for (int ilmp = 0; ilmp < GetMaxJM(); ilmp++) { // bylo *2
420 fSlice.fCovNum[2 * ilmz][2 * ilmp] = 0;
421 fSlice.fCovNum[2 * ilmz][2 * ilmp + 1] = 0;
422 }
423 }
424 return;
425 }
426
427 // Rewrite the new way to use the solving wherever there is inversion
428
429 std::unique_ptr<double[]> mDeltaT(new double[fMaxJM * fMaxJM * 4]);
430 for (int ilmzero = 0; ilmzero < GetMaxJM() * 2; ilmzero++)
431 for (int ilmprim = 0; ilmprim < GetMaxJM() * 2; ilmprim++)
432 mDeltaT[(ilmzero * fMaxJM * 2) + ilmprim] = fSlice.fCovNum[ilmzero][ilmprim];
433
434 if (fDebug) {
435 std::cout << "Delta T matrix " << std::endl;
436 for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++) {
437 for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
438 std::cout.precision(3);
439 std::cout.width(10);
440 std::cout << mDeltaT[ilmz * GetMaxJM() * 2 + ilmp];
441 }
442 std::cout << std::endl;
443 }
444 }
445
446 std::unique_ptr<double[]> mDeltaTPacked(new double[fMaxJM * fMaxJM * 4]);
447 int msize = PackYlmMatrixIndependentOnly(mDeltaT.get(), mDeltaTPacked.get());
448
449 if (fDebug) {
450 std::cout << "Delta T matrix packed " << std::endl;
451 for (int ilmz = 0; ilmz < msize; ilmz++) {
452 for (int ilmp = 0; ilmp < msize; ilmp++) {
453 std::cout.precision(3);
454 std::cout.width(10);
455 std::cout << mDeltaTPacked[ilmz * msize + ilmp];
456 }
457 std::cout << std::endl;
458 }
459 }
460
461 // (1) Solve (DeltaT)^1 Mtilde = Q
462
463 // Prepare halper matrices
464
465 std::unique_ptr<double[]> mM(new double[fMaxJM * fMaxJM * 4]);
466 std::unique_ptr<double[]> mMPacked(new double[fMaxJM * fMaxJM * 4]);
467 for (int iter = 0; iter < fMaxJM * fMaxJM * 4; iter++)
468 mM[iter] = tMTilde[iter];
469 PackYlmMatrixIndependentOnly(mM.get(), mMPacked.get());
470
471 gsl_matrix_view matM = gsl_matrix_view_array(mMPacked.get(), msize, msize);
472 if (fDebug) {
473 std::cout << "Mtilde matrix " << std::endl;
474 for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++) {
475 for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
476 std::cout.precision(3);
477 std::cout.width(10);
478 std::cout << mM[ilmz * GetMaxJM() * 2 + ilmp];
479 }
480 std::cout << std::endl;
481 }
482
483 std::cout << "Mtilde matrix packed " << std::endl;
484 for (int ilmz = 0; ilmz < msize; ilmz++) {
485 for (int ilmp = 0; ilmp < msize; ilmp++) {
486 std::cout.precision(3);
487 std::cout.width(10);
488 std::cout << mMPacked[ilmz * msize + ilmp];
489 }
490 std::cout << std::endl;
491 }
492 }
493
494 // Inverting matrix DeltaT.
495
496 std::unique_ptr<double[]> mU(new double[fMaxJM * fMaxJM * 4]);
497 InvertYlmIndependentMatrix(mDeltaT.get(), mU.get());
498
499 std::unique_ptr<double[]> mDTInvertedPacked(new double[fMaxJM * fMaxJM * 4]);
500 PackYlmMatrixIndependentOnly(mU.get(), mDTInvertedPacked.get());
501
502 gsl_matrix_view matDTI = gsl_matrix_view_array(mDTInvertedPacked.get(), msize, msize);
503
504 if (fDebug) {
505 std::cout << "Delta T matrix inverted packed " << std::endl;
506 for (int ilmz = 0; ilmz < msize; ilmz++) {
507 for (int ilmp = 0; ilmp < msize; ilmp++) {
508 std::cout.precision(3);
509 std::cout.width(10);
510 std::cout << mDTInvertedPacked[ilmz * msize + ilmp];
511 }
512 std::cout << std::endl;
513 }
514 }
515
516 // (2) Multiply DeltaT^1 M = Q
517 std::unique_ptr<double[]> mQ(new double[fMaxJM * fMaxJM * 4]);
518 for (int iter = 0; iter < msize * msize; iter++)
519 mQ[iter] = 0.0;
520 gsl_matrix_view matQ = gsl_matrix_view_array(mQ.get(), msize, msize);
521
522 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &matDTI.matrix, &matM.matrix, 0.0, &matQ.matrix);
523
524
525 std::unique_ptr<double[]> mTest(new double[fMaxJM * fMaxJM * 4]);
526 gsl_matrix_view matTest = gsl_matrix_view_array(mTest.get(), msize, msize);
527
528 std::unique_ptr<double[]> mF(new double[fMaxJM * fMaxJM * 4]);
529 for (int iter = 0; iter < fMaxJM * fMaxJM * 4; iter++)
530 mF[iter] = mDeltaTPacked[iter];
531 gsl_matrix_view matF = gsl_matrix_view_array(mF.get(), msize, msize);
532
533 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &matF.matrix, &matQ.matrix, 0.0, &matTest.matrix);
534
535 if (fDebug) {
536 std::cout << "Test matrix packed - compare to Mtilde" << std::endl;
537 for (int ilmz = 0; ilmz < msize; ilmz++) {
538 for (int ilmp = 0; ilmp < msize; ilmp++) {
539 std::cout.precision(3);
540 std::cout.width(10);
541 std::cout << mTest[ilmz * msize + ilmp];
542 }
543 std::cout << std::endl;
544 }
545 }
546
547 // (2) Multiply Mtilde^T Q = P
548
549 std::unique_ptr<double[]> mP(new double[fMaxJM * fMaxJM * 4]);
550 for (int iter = 0; iter < fMaxJM * fMaxJM * 4; iter++)
551 mP[iter] = 0;
552
553 gsl_matrix_view matP = gsl_matrix_view_array(mP.get(), msize, msize);
554
555 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &matM.matrix, &matQ.matrix, 0.0, &matP.matrix);
556
557 if (fDebug) {
558 std::cout << "P matrix packed " << std::endl;
559 for (int ilmz = 0; ilmz < msize; ilmz++) {
560 for (int ilmp = 0; ilmp < msize; ilmp++) {
561 std::cout.precision(3);
562 std::cout.width(10);
563 std::cout << mP[ilmz * msize + ilmp];
564 }
565 std::cout << std::endl;
566 }
567 }
568
569 // (3) Solve P^-1 Mtilde^T = R
570 std::unique_ptr<double[]> mPUnpacked(new double[fMaxJM * fMaxJM * 4]);
571 UnPackYlmMatrixIndependentOnly(mP.get(), mPUnpacked.get(), msize);
572
573 if (fDebug) {
574 std::cout << "P matrix unpacked " << std::endl;
575 for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++) {
576 for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
577 std::cout.precision(3);
578 std::cout.width(10);
579 std::cout << mPUnpacked[ilmz * GetMaxJM() * 2 + ilmp];
580 }
581 std::cout << std::endl;
582 }
583 }
584
585 // Invert the P matrix
586
587 std::unique_ptr<double[]> mPInverted(new double[fMaxJM * fMaxJM * 4]);
588 InvertYlmIndependentMatrix(mPUnpacked.get(), mPInverted.get());
589
590 std::unique_ptr<double[]> mPInvertedPacked(new double[fMaxJM * fMaxJM * 4]);
591 PackYlmMatrixIndependentOnly(mPInverted.get(), mPInvertedPacked.get());
592
593 gsl_matrix_view matPI = gsl_matrix_view_array(mPInvertedPacked.get(), msize, msize);
594
595 if (fDebug) {
596 std::cout << "P matrix inverted packed " << std::endl;
597 for (int ilmz = 0; ilmz < msize; ilmz++) {
598 for (int ilmp = 0; ilmp < msize; ilmp++) {
599 std::cout.precision(3);
600 std::cout.width(10);
601 std::cout << mPInvertedPacked[ilmz * msize + ilmp];
602 }
603 std::cout << std::endl;
604 }
605 }
606
607 // // gsl_matrix_view matR = gsl_matrix_view_array(mR, msize,
608 // msize);
609
610 // double mG[fMaxJM*fMaxJM*4];
611 // for (int iter=0; iter<fMaxJM*fMaxJM*4; iter++)
612 // mG[iter] = mP[iter];
613 // gsl_matrix_view matG = gsl_matrix_view_array(mG, msize, msize);
614
615 // // Decomposing the M matrix
616 // gsl_linalg_SV_decomp(&matG.matrix, &matS.matrix, &vecST.vector,
617 // &vecWT.vector);
618
619 // for (int itert=0; itert<msize; itert++) {
620 // for (int iterm=0; iterm<msize; iterm++)
621 // vCT[iterm] = mMPacked[iterm*msize + itert];
622 // // Transvere !!! ^^^^^ ^^^^^
623
624 // // Solving the problem
625 // gsl_linalg_SV_solve(&matG.matrix, &matS.matrix, &vecST.vector,
626 // &vecCT.vector, &vecXT.vector);
627
628 // for (int iterm=0; iterm<msize; iterm++)
629 // mR[itert*msize + iterm] = vXT[iterm];
630 // }
631
632
633 std::unique_ptr<double[]> mR(new double[fMaxJM * fMaxJM * 4]);
634 for (int ir = 0; ir < fMaxJM * fMaxJM * 4; ir++)
635 mR[ir] = 0.0;
636 gsl_matrix_view matR = gsl_matrix_view_array(mR.get(), msize, msize);
637
638 // (2) Multiply P^-1 M (Trans) = R
639
640 if (fDebug) {
641 std::cout << "Matrix M Packed " << std::endl;
642 for (int ilmz = 0; ilmz < msize; ilmz++) {
643 for (int ilmp = 0; ilmp < msize; ilmp++) {
644 std::cout.precision(3);
645 std::cout.width(10);
646 std::cout << mMPacked[ilmz * msize + ilmp];
647 }
648 std::cout << std::endl;
649 }
650 }
651
652 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &matPI.matrix, &matM.matrix, 1.0, &matR.matrix);
653
654 if (fDebug) {
655 std::cout << "R matrix packed " << std::endl;
656
657 for (int ilmz = 0; ilmz < msize; ilmz++) {
658 for (int ilmp = 0; ilmp < msize; ilmp++) {
659 std::cout.precision(3);
660 std::cout.width(10);
661 std::cout << mR[ilmz * msize + ilmp];
662 }
663 std::cout << std::endl;
664 }
665 }
666
667 // (4) Solve DeltaT^-1 T = L
668 std::unique_ptr<double[]> vL(new double[fMaxJM * 2]);
669 for (int i = 0; i < fMaxJM * 2; i++)
670 vL[i] = 0;
671 gsl_vector_view vecL = gsl_vector_view_array(vL.get(), msize);
672
673 // // Decomposing the M matrix
674 // gsl_linalg_SV_decomp(&matF.matrix, &matS.matrix, &vecST.vector,
675 // &vecWT.vector);
676
677 std::unique_ptr<double[]> vB(new double[fMaxJM * 2]);
678 for (int iter = 0; iter < GetMaxJM(); iter++) {
679 vB[iter * 2] = real(tTq0[iter]);
680 vB[iter * 2 + 1] = imag(tTq0[iter]);
681 }
682
683 std::unique_ptr<double[]> vBPacked(new double[fMaxJM * 2]);
684 PackYlmVectorIndependentOnly(vB.get(), vBPacked.get());
685
686 gsl_vector_view vecB = gsl_vector_view_array(vBPacked.get(), msize);
687
688 // // Solving the problem
689 // gsl_linalg_SV_solve(&matF.matrix, &matS.matrix, &vecST.vector,
690 // &vecB.vector, &vecL.vector);
691
692 if (fDebug) {
693 std::cout << "L vector packed " << std::endl;
694 for (int ilmp = 0; ilmp < msize; ilmp++) {
695 std::cout.precision(3);
696 std::cout.width(10);
697 std::cout << vL[ilmp];
698 }
699 std::cout << std::endl;
700 }
701
702 // Multiply DeltaT^-1 T = L
703
704 gsl_blas_dgemv(CblasNoTrans, 1.0, &matDTI.matrix, &vecB.vector, 0.0, &vecL.vector);
705
706 // (5) Multiply R L = C
707
708 std::unique_ptr<double[]> vY(new double[fMaxJM * 2]);
709 for (int iter = 0; iter < GetMaxJM() * 2; iter++) {
710 vY[iter] = 0.0;
711 }
712
713 // Prepare inputs for solving the problem
714 gsl_vector_view vecY = gsl_vector_view_array(vY.get(), msize);
715
716 gsl_blas_dgemv(CblasNoTrans, 1.0, &matR.matrix, &vecL.vector, 0.0, &vecY.vector);
717
718 if (fDebug) {
719 std::cout << "C vector packed " << std::endl;
720 for (int ilmp = 0; ilmp < msize; ilmp++) {
721 std::cout.precision(3);
722 std::cout.width(10);
723 std::cout << vY[ilmp];
724 }
725 std::cout << std::endl;
726 }
727 int mpack = 0;
728 int el, em;
729 for (int ilm = 0; ilm < fMaxJM; ilm++) {
730 // fCFReal[ilm]->SetBinContent(ibin, vC[mpack++]);
731 GetElEmForIndex(ilm, el, em);
732 if ((el % 2) == 1) {
733 fSlice.fCFReal[ilm] = 0.0;
734 fSlice.fCFImag[ilm] = 0.0;
735 }
736
737#ifdef FULL_CALC
738 fCFReal[ilm]->SetBinContent(ibin, fNumReal[ilm]->GetBinContent(ibin) / fDenReal[ilm]->GetBinContent(ibin));
739 fCFImag[ilm]->SetBinContent(ibin, fNumImag[ilm]->GetBinContent(ibin) / fDenImag[ilm]->GetBinContent(ibin));
740
741#else
742 if (em < 0) {
743 fSlice.fCFReal[ilm] = 0.0;
744 fSlice.fCFImag[ilm] = 0.0;
745 } else {
746 fSlice.fCFReal[ilm] = vY[mpack++];
747 if (em == 0) {
748 fSlice.fCFImag[ilm] = 0;
749 } else {
750 // fSlice.fCFImag[ilm] = vC[mpack++];
751 fSlice.fCFImag[ilm] = vY[mpack++];
752 }
753 }
754
755#endif
756 }
757
758 // invert the P matrix to get C errors
759 // double mS[fMaxJM*fMaxJM*4];
760
761 // for (int iterz=0; iterz<msize; iterz++)
762 // for (int iterp=0; iterp<msize; iterp++)
763 // if (iterp == iterz)
764 // mS[iterz*msize + iterp] = 1.0;
765 // else
766 // mS[iterz*msize + iterp] = 0.0;
767
768 // gsl_matrix_view matS = gsl_matrix_view_array(mS, msize, msize);
769
770 // Invert V
771
772 // gsl_blas_dtrsm(CblasLeft, CblasUpper, CblasNoTrans,
773 // CblasNonUnit, 1.0, &matP.matrix, &matS.matrix);
774
775 mpack = 0;
776 for (int ilm = 0; ilm < fMaxJM; ilm++) {
777 GetElEmForIndex(ilm, el, em);
778 if (em < 0) {
779 fSlice.fCFRealE[ilm] = 0;
780 fSlice.fCFImagE[ilm] = 0;
781 } else {
782 fSlice.fCFRealE[ilm] = sqrt(fabs(mPInvertedPacked[mpack * msize + mpack]));
783 mpack++;
784 if (em == 0)
785 fSlice.fCFImag[ilm] = 0;
786 else {
787 fSlice.fCFImagE[ilm] = sqrt(fabs(mPInvertedPacked[mpack * msize + mpack]));
788 mpack++;
789 }
790 }
791 }
792
793 for (int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++) {
794 for (int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
795 if (ilmp > ilmz)
796 fSlice.fCovCF[ilmz][ilmp] = mPInverted[ilmz * GetMaxJM() * 2 + ilmp];
797 else
798 fSlice.fCovCF[ilmz][ilmp] = mPInverted[ilmp * GetMaxJM() * 2 + ilmz];
799 }
800 }
801 }
802
803 void FemtoYlmSolver::Solve(Bool_t recalc) {
804 int nbins = fCF->GetNum()->GetNbinsX();
805 for (int i = 1; i <= nbins; i++) {
806 fSlice.BuildSlice(*fCF, i);
807 if (fDebugBin == i) {
808 fDebug = kTRUE;
809 } else {
810 fDebug = kFALSE;
811 }
812 DoMath(recalc);
813 FixCF();
814 UpdateCF(i);
815 }
816 }
817
818 void FemtoYlmSolver::UpdateCF(Int_t bin) {
819 for (int j = 0; j < fMaxJM; j++) {
820 fCF->fCFReal[j]->SetBinContent(bin, fSlice.fCFReal[j]);
821 fCF->fCFImag[j]->SetBinContent(bin, fSlice.fCFImag[j]);
822 fCF->fCFReal[j]->SetBinError(bin, fSlice.fCFRealE[j]);
823 fCF->fCFImag[j]->SetBinError(bin, fSlice.fCFImagE[j]);
824 }
825 for (int a = 0; a < 2 * fMaxJM; a++) {
826 for (int b = 0; b < 2 * fMaxJM; b++) {
827 fCF->fCovCf[bin - 1][a][b] = fSlice.fCovCF[a][b];
828 // fCF->fCovNum[bin - 1][a][b] = fSlice.fCovNum[a][b]; do not update cov matrix in num
829 }
830 }
831 }
832
833 void FemtoYlmSolver::FixCF() {
834 for (int l = 0; l <= fMaxL; l++) {
835 for (int m = -l; m < 0; m++) {
836 Int_t from = fLmVals.GetIndex(l, -m);
837 Int_t to = fLmVals.GetIndex(l, m);
838 fSlice.fCFReal[to] = fSlice.fCFReal[from];
839 fSlice.fCFImag[to] = -fSlice.fCFImag[from];
840 fSlice.fCFRealE[to] = fSlice.fCFRealE[from];
841 fSlice.fCFImagE[to] = fSlice.fCFImagE[from];
842 }
843 }
844 }
845
846 void FemtoYlmSolver::SetNormalizationArea(Double_t min, Double_t max) {
847 TH1* num0 = fCF->fNumReal[0];
848 TH1* den0 = fCF->fDenReal[0];
849 fNormFactor = 1.0;
850 // TODO Fix/improvenormalization
851 double normbinmax = min; // fDenReal[0]->FindBin(fNormMax);
852 double normbinmin = max; // fDenReal[0]->FindBin(fNormMin);
853
854 if (normbinmax > 0) {
855 double sksum = 0.0;
856 double wksum = 0.0;
857
858 double sk, wk, ks;
859 if (normbinmin < 1) normbinmin = 1;
860 if (normbinmax > den0->GetNbinsX()) normbinmax = den0->GetNbinsX();
861 normbinmin = 1; // KURWA
862 normbinmax = den0->GetNbinsX(); // KURWA
863 for (int ib = normbinmin; ib <= normbinmax; ib++) {
864 ks = den0->GetXaxis()->GetBinCenter(ib);
865 sk = num0->GetBinContent(ib) / (den0->GetBinContent(ib) * (1.0 - fNormPurity / (fNormRadius * fNormBohr * ks * ks)));
866 wk = num0->GetBinContent(ib);
867 sksum += sk * wk;
868 wksum += wk;
869 }
870 fNormFactor *= sksum / wksum;
871 if (fSumw) {
872 fNormFactor /= num0->GetEffectiveEntries() / den0->GetEffectiveEntries();
873 } else
874 fNormFactor /= num0->GetEntries() / den0->GetEntries();
875 }
876 // fNumEnt = num0->GetEntries();
877 // fDenEnt = den0->GetEntries();
878 }
879
880 void FemtoYlmSolver::MakeFakeCovMatrix() {
881 double nent = fNumEnt;
882 double nentd = fDenEnt;
883 for (int ilmx = 0; ilmx < GetMaxJM(); ilmx++) {
884 const int ilmx2 = ilmx * 2;
885 for (int ilmy = 0; ilmy < GetMaxJM(); ilmy++) {
886 const int ilmy2 = ilmy * 2;
887 double t1t2rr = fSlice.fNumReal[ilmx] * fSlice.fNumReal[ilmy] / nent / nent;
888 double t1t2ri = fSlice.fNumReal[ilmx] * fSlice.fNumImag[ilmy] / nent / nent;
889 double t1t2ir = fSlice.fNumImag[ilmx] * fSlice.fNumReal[ilmy] / nent / nent;
890 double t1t2ii = fSlice.fNumImag[ilmx] * fSlice.fNumImag[ilmy] / nent / nent;
891 if (ilmx == ilmy) {
892 fSlice.fCovNum[ilmx2][ilmy2] = nent * (TMath::Power(fSlice.fNumRealE[ilmx] / nent, 2) * (nent - 1) + t1t2rr);
893 fSlice.fCovNum[ilmx2][ilmy2 + 1] = nent * t1t2ri;
894 fSlice.fCovNum[ilmx2 + 1][ilmy2] = nent * t1t2ir;
895 fSlice.fCovNum[ilmx2 + 1][ilmy2 + 1] = nent * (TMath::Power(fSlice.fNumImagE[ilmx] / nent, 2) * (nent - 1) + t1t2rr);
896 } else {
897 fSlice.fCovNum[ilmx2][ilmy2] = nent * t1t2rr;
898 fSlice.fCovNum[ilmx2][ilmy2 + 1] = nent * t1t2ri;
899 fSlice.fCovNum[ilmx2 + 1][ilmy2] = nent * t1t2ir;
900 fSlice.fCovNum[ilmx2 + 1][ilmy2 + 1] = nent * t1t2ii;
901 }
902 t1t2rr = fSlice.fDenReal[ilmx] * fSlice.fDenReal[ilmy] / nentd / nentd;
903 t1t2ri = fSlice.fDenReal[ilmx] * fSlice.fDenImag[ilmy] / nentd / nentd;
904 t1t2ir = fSlice.fDenImag[ilmx] * fSlice.fDenReal[ilmy] / nentd / nentd;
905 t1t2ii = fSlice.fDenImag[ilmx] * fSlice.fDenImag[ilmy] / nentd / nentd;
906
907
908 fSlice.fCovDen[ilmx2 + 0][ilmy2 + 0] = nentd * t1t2rr;
909 fSlice.fCovDen[ilmx2 + 0][ilmy2 + 1] = nentd * t1t2ri;
910 fSlice.fCovDen[ilmx2 + 1][ilmy2 + 0] = nentd * t1t2ir;
911 fSlice.fCovDen[ilmx2 + 1][ilmy2 + 1] = nentd * t1t2ii;
912 }
913 }
914 }
915
916} /* namespace Hal */