16#include <TCollection.h>
22#include <TObjString.h>
27#include <gsl/gsl_blas.h>
28#include <gsl/gsl_cblas.h>
29#include <gsl/gsl_matrix_double.h>
30#include <gsl/gsl_vector_double.h>
34#include "CorrFitSHCF.h"
39#include "FemtoSHCFPainter.h"
40#include "FemtoSHSlice.h"
41#include "FemtoSerializationInterfaceSH.h"
42#include "FemtoYlmSolver.h"
43#include "HistoStyle.h"
46#include "HtmlObject.h"
49#include "MarkerStyle.h"
62 fFrame(Femto::EKinematics::kLCMS),
67 gSystem->Load(
"libgsl.so");
68 gSystem->Load(
"libgslcblas.so");
71 FemtoSHCF::FemtoSHCF(TString name, Int_t maxL, Int_t bins, Double_t min, Double_t max, Femto::EKinematics kinematics) :
73 fMaxJM((TMath::Min(maxL, 5) + 1) * (TMath::Min(maxL, 5) + 1)),
82 fNumReal =
new TH1D*[fMaxJM];
83 fNumImag =
new TH1D*[fMaxJM];
84 fDenReal =
new TH1D*[fMaxJM];
85 fDenImag =
new TH1D*[fMaxJM];
87 TString x_name = Hal::Femto::KinematicsToAxisLabel(fFrame, 0, 1);
88 TString y_name = Hal::Femto::KinematicsToAxisLabel(fFrame, 1, 1);
90 for (
int ihist = 0; ihist < fMaxJM; ihist++) {
91 TString hname = Form(
"NumReYlm%i%i", fLmVals.
GetElsi(ihist), fLmVals.
GetEmsi(ihist));
92 fNumReal[ihist] =
new TH1D(hname, hname, bins, min, max);
93 hname = Form(
"NumImYlm%i%i", fLmVals.
GetElsi(ihist), fLmVals.
GetEmsi(ihist));
94 fNumImag[ihist] =
new TH1D(hname, hname, bins, min, max);
95 hname = Form(
"DenReYlm%i%i", fLmVals.
GetElsi(ihist), fLmVals.
GetEmsi(ihist));
96 fDenReal[ihist] =
new TH1D(hname, hname, bins, min, max);
97 hname = Form(
"DenImYlm%i%i", (
int) fLmVals.
GetElsi(ihist), fLmVals.
GetEmsi(ihist));
98 fDenImag[ihist] =
new TH1D(hname, hname, bins, min, max);
99 fNumReal[ihist]->GetXaxis()->SetTitle(x_name);
100 fNumReal[ihist]->GetYaxis()->SetTitle(y_name);
101 fNumImag[ihist]->GetXaxis()->SetTitle(x_name);
102 fNumImag[ihist]->GetYaxis()->SetTitle(y_name);
103 fDenReal[ihist]->GetXaxis()->SetTitle(x_name);
104 fDenReal[ihist]->GetYaxis()->SetTitle(y_name);
105 fDenImag[ihist]->GetXaxis()->SetTitle(x_name);
106 fDenImag[ihist]->GetYaxis()->SetTitle(y_name);
107 fNumReal[ihist]->Sumw2();
108 fNumImag[ihist]->Sumw2();
109 fDenReal[ihist]->Sumw2();
110 fDenImag[ihist]->Sumw2();
113 fCovNum.
MakeBigger(bins, fMaxJM * 2, fMaxJM * 2);
114 fCovDen.
MakeBigger(bins, fMaxJM * 2, fMaxJM * 2);
115 fCovCf.
MakeBigger(bins, fMaxJM * 2, fMaxJM * 2);
120 AddLabel(Hal::Femto::KinematicsToLabel(fFrame));
126 fMaxJM(other.fMaxJM),
127 fFrame(other.fFrame),
128 fNormPurity(other.fNormPurity),
129 fNormRadius(other.fNormRadius),
130 fNormBohr(other.fNormBohr),
131 fLmVals(other.fLmVals) {
132 if (other.fNumReal) {
133 fNumReal = new TH1D*[fMaxJM];
134 fNumImag = new TH1D*[fMaxJM];
135 fDenReal = new TH1D*[fMaxJM];
136 fDenImag = new TH1D*[fMaxJM];
137 for (int i = 0; i < fMaxJM; i++) {
138 fNumReal[i] = (TH1D*) other.fNumReal[i]->Clone();
139 fNumImag[i] = (TH1D*) other.fNumImag[i]->Clone();
140 fDenReal[i] = (TH1D*) other.fDenReal[i]->Clone();
141 fDenImag[i] = (TH1D*) other.fDenImag[i]->Clone();
146 fCFReal = new TH1D*[fMaxJM];
147 fCFImag = new TH1D*[fMaxJM];
148 for (int i = 0; i < fMaxJM; i++) {
149 fCFReal[i] = (TH1D*) other.fCFReal[i]->Clone();
150 fCFImag[i] = (TH1D*) other.fCFImag[i]->Clone();
153 fCovNum = other.fCovNum;
154 fCovDen = other.fCovDen;
155 fCovCf = other.fCovCf;
156 if (other.fCfcov) { fCfcov = (TH3D*) other.fCfcov->Clone(); }
161 TH1D* FemtoSHCF::GetCFRe(
int el,
int em)
const {
162 if (fLmVals.GetIndex(el, em) >= 0 && fCFReal !=
nullptr) {
163 return fCFReal[fLmVals.GetIndex(el, em)];
169 TH1D* FemtoSHCF::GetCFIm(
int el,
int em)
const {
170 if (fLmVals.GetIndex(el, em) >= 0 && fCFImag !=
nullptr) {
171 return fCFImag[fLmVals.GetIndex(el, em)];
177 TH1D* FemtoSHCF::GetNumRe(
int el,
int em)
const {
178 if (fLmVals.GetIndex(el, em) >= 0)
179 return fNumReal[fLmVals.GetIndex(el, em)];
184 TH1D* FemtoSHCF::GetNumIm(
int el,
int em)
const {
185 if (fLmVals.GetIndex(el, em) >= 0)
186 return fNumImag[fLmVals.GetIndex(el, em)];
191 TH1D* FemtoSHCF::GetDenRe(
int el,
int em)
const {
192 if (fLmVals.GetIndex(el, em) >= 0)
193 return fDenReal[fLmVals.GetIndex(el, em)];
198 TH1D* FemtoSHCF::GetDenIm(
int el,
int em)
const {
199 if (fLmVals.GetIndex(el, em) >= 0)
200 return fDenImag[fLmVals.GetIndex(el, em)];
205 void FemtoSHCF::SetNumRe(TH1D** histograms, Bool_t clone) {
207 if (fNumReal ==
nullptr) fNumReal =
new TH1D*[fMaxJM];
208 for (
int i = 0; i < fMaxJM; i++) {
209 fNumReal[i] = (TH1D*) histograms[i]->Clone();
212 fNumReal = histograms;
216 void FemtoSHCF::SetNumIm(TH1D** histograms, Bool_t clone) {
218 if (fNumImag ==
nullptr) fNumImag =
new TH1D*[fMaxJM];
219 for (
int i = 0; i < fMaxJM; i++) {
220 fNumImag[i] = (TH1D*) histograms[i]->Clone();
223 fNumImag = histograms;
227 void FemtoSHCF::SetDenRe(TH1D** histograms, Bool_t clone) {
229 if (fDenReal ==
nullptr) fDenReal =
new TH1D*[fMaxJM];
230 for (
int i = 0; i < fMaxJM; i++) {
231 fDenReal[i] = (TH1D*) histograms[i]->Clone();
234 fDenReal = histograms;
238 void FemtoSHCF::SetDenIm(TH1D** histograms, Bool_t clone) {
240 if (fDenImag ==
nullptr) fDenImag =
new TH1D*[fMaxJM];
241 for (
int i = 0; i < fMaxJM; i++) {
242 fDenImag[i] = (TH1D*) histograms[i]->Clone();
245 fDenImag = histograms;
249 void FemtoSHCF::Draw(Option_t* opt) {
250 TString option = opt;
252 fPainter->SetOption(option);
256 fPainter->SetOption(option);
261 void FemtoSHCF::PackCfcCovariance() {
263 if (fCfcov)
delete fCfcov;
264 bufName = Form(
"CovCfc%s", fNumReal[0]->GetName());
265 fCfcov =
new TH3D(bufName,
267 fCFReal[0]->GetNbinsX(),
268 fCFReal[0]->GetXaxis()->GetXmin(),
269 fCFReal[0]->GetXaxis()->GetXmax(),
272 GetMaxJM() * 2 - 0.5,
275 GetMaxJM() * 2 - 0.5);
278 for (
int ibin = 1; ibin <= fCfcov->GetNbinsX(); ibin++) {
279 const int binCov = ibin - 1;
280 for (
int ilmz = 0; ilmz < GetMaxJM() * 2; ilmz++)
281 for (
int ilmp = 0; ilmp < GetMaxJM() * 2; ilmp++) {
282 tK = fCovCf[binCov][ilmz][ilmp];
283 fCfcov->SetBinContent(ibin, ilmz + 1, ilmp + 1, tK);
286 fCfcov->SetBinContent(0, 0, 0, 1.0);
289 void FemtoSHCF::Browse(TBrowser* b) {
290 if (fCFReal ==
nullptr) {
291 Cout::PrintInfo(
"No correlation functions, call RecalculateCF", EInfo::kError);
294 for (
int i = 0; i < fMaxJM; i++) {
300 for (
int i = 0; i < fMaxJM; i++) {
311 void FemtoSHCF::NormalizeBy(
int el,
int em, Option_t* opt) {
312 TString option = opt;
315 if (option.EqualTo(
"re")) {
316 num = GetNumRe(el, em);
317 den = GetDenRe(el, em);
319 num = GetNumIm(el, em);
320 den = GetDenIm(el, em);
321 if (num) num->SetLineColor(kRed);
325 Double_t n_integral, d_integral;
326 if (fNormMax == fNormMin) {
327 n_integral = num->Integral(1, num->GetNbinsX());
328 d_integral = den->Integral(1, den->GetNbinsX());
330 n_integral = num->Integral(num->GetXaxis()->FindBin(fNormMin[0]), num->GetXaxis()->FindBin(fNormMax[0]));
331 d_integral = den->Integral(den->GetXaxis()->FindBin(fNormMin[0]), den->GetXaxis()->FindBin(fNormMax[0]));
333 fScale = d_integral / n_integral;
336 void FemtoSHCF::AddLabel(TString label) {
337 if (fLabels ==
nullptr) {
338 fLabels =
new TList();
339 fLabels->SetName(
"Labels");
341 fLabels->AddLast(
new TObjString(label));
344 void FemtoSHCF::Add(
const Object* pack) {
346 if (cf->fMaxJM != fMaxJM) {
347 Cout::PrintInfo(
"Not compatible maxL numbers during merging FemtoSHCF", EInfo::kError);
350 if (cf->fFrame != this->fFrame) {
351 Cout::PrintInfo(
"Incompatible frames !", EInfo::kError);
354 for (
int i = 0; i < fMaxJM; i++) {
355 fNumReal[i]->Add(cf->fNumReal[i]);
356 fNumImag[i]->Add(cf->fNumImag[i]);
357 fDenReal[i]->Add(cf->fDenReal[i]);
358 fDenImag[i]->Add(cf->fDenImag[i]);
362 for (
int i = 0; i < fCovNum.GetSize(); i++) {
363 for (
int j = 0; j < fCovNum[0].GetSize(); j++) {
364 for (
int k = 0; k < fCovNum[0][0].GetSize(); k++) {
365 fCovNum[i][j][k] = fCovNum[i][j][k] + cf->fCovNum[i][j][k];
366 fCovDen[i][j][k] = fCovDen[i][j][k] + cf->fCovDen[i][j][k];
373 TH1D* FemtoSHCF::GetHisto(
int el,
int em, Bool_t norm, Option_t* opt)
const {
374 TString option = opt;
376 if (option.EqualTo(
"re")) {
377 cf = GetCFRe(el, em);
379 cf = GetCFIm(el, em);
381 cf = (TH1D*) cf->Clone();
382 cf->GetYaxis()->SetTitle(Form(
"C^{%i}_{%i}", el, em));
383 cf->SetTitle(Form(
"CF^{%i}_{%i}", el, em));
385 cf->SetMaximum(cf->GetBinContent(cf->GetMaximumBin()) * 1.1);
386 if (norm && fScale != 0) cf->Scale(fScale);
390 TH1D* FemtoSHCF::GetCubic(Option_t* opt, Bool_t R)
const {
391 if (fCFReal ==
nullptr) {
392 Cout::PrintInfo(
"No correlation functions, call RecalculateCF", EInfo::kError);
395 TString option = opt;
396 if (option.Length() > GetLMax())
return nullptr;
402 if (option.Length() == 1) {
403 Double_t scale1 = 2.0 * Sqr(1, 3);
404 Double_t scale2 = Sqr(2, 3);
407 cf = Histo(1, 0,
"re", scale1);
408 }
else if (option ==
"y") {
409 cf = Histo(1, 1,
"im", scale2);
410 cf2 = Histo(1, -1,
"im", scale2);
411 }
else if (option ==
"x") {
412 cf = Histo(1, 1,
"re", -scale2);
413 cf2 = Histo(1, -1,
"re", scale2);
414 }
else if (option ==
"0") {
415 cf = Histo(0, 0,
"re", 1.0);
418 if (option.Length() == 2) {
419 Double_t scale = Sqr(2, 15);
420 if (option ==
"zz") {
421 cf = Histo(2, 0,
"re", 4.0 / 3.0 * Sqr(1, 5));
422 }
else if (option ==
"yz") {
423 cf = Histo(2, 1,
"im", scale);
424 cf2 = Histo(2, -1,
"im", scale);
425 }
else if (option ==
"yy") {
426 cf = Histo(2, 2,
"re", -scale);
427 cf2 = Histo(2, 0,
"re", -2.0 / 3.0 * Sqr(1, 5));
428 cf3 = Histo(2, -2,
"re", -scale);
429 }
else if (option ==
"xz") {
430 cf = Histo(2, 1,
"re", -scale);
431 cf2 = Histo(2, -1,
"re", scale);
432 }
else if (option ==
"xy") {
433 cf = Histo(2, 2,
"im", -scale);
434 cf2 = Histo(2, -2,
"im", scale);
435 }
else if (option ==
"xx") {
436 cf = Histo(2, 2,
"re", scale);
437 cf2 = Histo(2, 0,
"re", -2.0 / 3.0 * Sqr(1, 5));
438 cf3 = Histo(2, -2,
"re", scale);
441 if (option.Length() == 3) {
442 Double_t scale7 = 1.0 / 5.0 * Sqr(1, 7);
443 Double_t scale21 = 1.0 / 5.0 * Sqr(1, 21);
444 Double_t scale105 = Sqr(2, 105);
445 Double_t scale35 = Sqr(1, 35);
446 if (option ==
"zzz") {
447 cf = Histo(3, 0,
"re", 4.0 * scale7);
448 }
else if (option ==
"yzz") {
449 cf = Histo(3, 1,
"im", 4.0 * scale21);
450 cf2 = Histo(3, -1,
"im", 4.0 * scale21);
451 }
else if (option ==
"yyz") {
452 cf = Histo(3, 2,
"re", -scale105);
453 cf2 = Histo(3, 0,
"re", -2.0 * scale7);
454 cf3 = Histo(3, -2,
"re", -scale105);
455 }
else if (option ==
"yyy") {
456 cf = Histo(3, 3,
"im", -scale35);
457 cf2 = Histo(3, 1,
"im", -Sqr(3, 7) / 5.0);
458 cf3 = Histo(3, -1,
"im", -Sqr(3, 7) / 5.0);
459 cf4 = Histo(3, -3,
"im", -scale35);
460 }
else if (option ==
"xzz") {
461 cf = Histo(3, 1,
"re", -4.0 * scale21);
462 cf2 = Histo(3, -1,
"re", 4.0 * scale21);
463 }
else if (option ==
"xyz") {
464 cf = Histo(3, 2,
"im", -scale105);
465 cf2 = Histo(3, -2,
"im", scale105);
466 }
else if (option ==
"xyy") {
467 cf = Histo(3, 3,
"re", scale35);
468 cf2 = Histo(3, 1,
"re", scale21);
469 cf3 = Histo(3, -1,
"re", -scale21);
470 cf4 = Histo(3, -3,
"re", -scale35);
471 }
else if (option ==
"xxz") {
472 cf = Histo(3, 2,
"re", scale105);
473 cf2 = Histo(3, 0,
"re", -2.0 * scale7);
474 cf3 = Histo(3, -2,
"re", scale105);
475 }
else if (option ==
"xxy") {
476 cf = Histo(3, 3,
"im", scale35);
477 cf2 = Histo(3, 1,
"im", -scale21);
478 cf3 = Histo(3, -1,
"im", -scale21);
479 cf4 = Histo(3, -3,
"im", scale35);
480 }
else if (option ==
"xxx") {
481 cf = Histo(3, 3,
"re", -scale35);
482 cf2 = Histo(3, 1,
"re", scale7);
483 cf3 = Histo(3, -1,
"re", -scale7);
484 cf4 = Histo(3, -3,
"re", scale35);
487 if (option.Length() == 4) {
488 Double_t scale5 = Sqr(1, 5) / 21.0;
489 Double_t scale25 = Sqr(2, 5) / 21.0;
490 Double_t scale35 = Sqr(1, 35) / 3.0;
491 Double_t scale235 = Sqr(2, 35) / 3.0;
492 Double_t scale57 = Sqr(1, 5) / 7.0;
493 Double_t pis = TMath::Sqrt(TMath::Pi());
494 if (option ==
"zzzz") {
495 cf = Histo(4, 0,
"re", TMath::Sqrt(TMath::Pi()) * 16.0 / 105.0);
496 }
else if (option ==
"yzzz") {
497 cf = Histo(4, 1,
"im", scale5 * 4);
498 cf2 = Histo(4, -1,
"im", scale5 * 4.0);
499 }
else if (option ==
"yyzz") {
500 cf = Histo(4, 2,
"re", -2.0 * scale25);
501 cf2 = Histo(4, 0,
"re", -8.0 / 105.0 * pis);
502 cf3 = Histo(4, -2,
"re", -2.0 * scale25);
503 }
else if (option ==
"yyyz") {
504 cf = Histo(4, 3,
"im", -scale35);
505 cf2 = Histo(4, 1,
"im", -scale57);
506 cf3 = Histo(4, -1,
"im", -scale57);
507 cf4 = Histo(4, -3,
"im", -scale35);
508 }
else if (option ==
"yyyy") {
509 cf = Histo(4, 4,
"re", scale235);
510 cf2 = Histo(4, 2,
"re", scale25 * 2);
511 cf3 = Histo(4, 0,
"re", 2.0 / 35.0 * pis);
512 cf4 = Histo(4, 2,
"re", scale25 * 2);
513 cf5 = Histo(4, -4,
"re", scale235);
514 }
else if (option ==
"xzzz") {
515 cf = Histo(4, 1,
"re", -4.0 * scale5);
516 cf2 = Histo(4, -1,
"re", 4.0 * scale5);
517 }
else if (option ==
"xyzz") {
518 cf = Histo(4, 2,
"im", -2.0 * scale25);
519 cf2 = Histo(4, -2,
"im", 2.0 * scale25);
520 }
else if (option ==
"xyzz") {
521 cf = Histo(4, 3,
"re", scale35);
522 cf2 = Histo(4, 1,
"re", scale5);
523 cf3 = Histo(4, -1,
"re", -scale5);
524 cf4 = Histo(4, -3,
"re", -scale35);
525 }
else if (option ==
"xyyz") {
526 cf = Histo(4, 3,
"re", scale35);
527 cf2 = Histo(4, 1,
"re", scale5);
528 cf3 = Histo(4, -1,
"re", -scale5);
529 cf4 = Histo(4, -3,
"re", -scale35);
530 }
else if (option ==
"xyyy") {
531 cf = Histo(4, 4,
"im", scale235);
532 cf2 = Histo(4, 2,
"im", scale25);
533 cf3 = Histo(4, -2,
"im", -scale25);
534 cf4 = Histo(4, -4,
"im", -scale235);
535 }
else if (option ==
"xxzz") {
536 cf = Histo(4, 2,
"re", 2.0 * scale25);
537 cf2 = Histo(4, 0,
"re", -8.0 / 105.0 * pis);
538 cf3 = Histo(4, -2,
"re", 2.0 * scale25);
539 }
else if (option ==
"xxyz") {
540 cf = Histo(4, 3,
"im", scale35);
541 cf2 = Histo(4, 1,
"im", -scale5);
542 cf3 = Histo(4, -1,
"im", -scale5);
543 cf4 = Histo(4, -3,
"im", scale35);
544 }
else if (option ==
"xxyy") {
545 cf = Histo(4, 4,
"re", -scale235);
546 cf2 = Histo(4, 0,
"re", 2.0 / 105. * pis);
547 cf3 = Histo(4, 4,
"re", -scale235);
548 }
else if (option ==
"xxxz") {
549 cf = Histo(4, 3,
"re", -scale35);
550 cf2 = Histo(4, 1,
"re", scale57);
551 cf3 = Histo(4, -1,
"re", -scale57);
552 cf4 = Histo(4, -3,
"re", scale35);
553 }
else if (option ==
"xxxy") {
554 cf = Histo(4, 4,
"im", -scale235);
555 cf2 = Histo(4, 2,
"im", scale25);
556 cf3 = Histo(4, -2,
"im", -scale25);
557 cf4 = Histo(4, -4,
"im", scale235);
558 }
else if (option ==
"xxxx") {
559 cf = Histo(4, 4,
"re", scale235);
560 cf2 = Histo(4, 2,
"re", -2.0 * scale25);
561 cf3 = Histo(4, 0,
"re", 2.0 / 35.0 * pis);
562 cf4 = Histo(4, -2,
"re", -2.0 * scale25);
563 cf5 = Histo(4, -4,
"re", scale235);
566 if (cf ==
nullptr)
return nullptr;
583 cf->SetTitle(option);
586 cf->SetTitle(option +
" R(q)");
587 for (
int i = 1; i <= cf->GetNbinsX(); i++) {
588 cf->SetBinContent(i, cf->GetBinContent(i) - fScale);
594 Double_t FemtoSHCF::Sqr(Double_t val1, Double_t val2)
const {
return TMath::Sqrt(TMath::Pi() * val1 / val2); }
596 TH1D* FemtoSHCF::Histo(
int lm,
int em, Option_t* opt, Double_t scale)
const {
598 TString option = opt;
599 if (option.EqualTo(
"im")) {
600 cf = (TH1D*) GetCFIm(lm, TMath::Abs(em))->Clone();
602 cf = (TH1D*) GetCFRe(lm, TMath::Abs(em))->Clone();
604 if (option ==
"im") {
612 TH3D* FemtoSHCF::GetCovCF()
const {
return fCfcov; }
614 void FemtoSHCF::RecalculateCF(Int_t debug, Bool_t suwm) {
616 if (fDenImag ==
nullptr) {
617 Cout::PrintInfo(
"No imaginary denominators!", EInfo::kError);
620 if (fDenReal ==
nullptr) {
621 Cout::PrintInfo(
"No real denominators", EInfo::kError);
624 if (fNumReal ==
nullptr) {
625 Cout::PrintInfo(
"No real numeraotrs", EInfo::kError);
628 if (fNumImag ==
nullptr) {
629 Cout::PrintInfo(
"No imaginary numerators", EInfo::kError);
633 Cout::PrintInfo(
"Computing SF correlation functions", EInfo::kDebugInfo);
635 for (
int i = 0; i < fMaxJM; i++)
640 for (
int i = 0; i < fMaxJM; i++)
644 fCFReal =
new TH1D*[fMaxJM];
645 fCFImag =
new TH1D*[fMaxJM];
647 for (
int i = 0; i < fMaxJM; i++) {
648 TString name = fNumReal[i]->GetName();
649 name.ReplaceAll(
"Num",
"CF");
652 Std::GetAxisPar(*fNumReal[i], nbins, min, max,
"x");
653 fCFReal[i] =
new TH1D(name, name, nbins, min, max);
654 name = fNumImag[i]->GetName();
655 name.ReplaceAll(
"Num",
"CF");
656 fCFImag[i] =
new TH1D(name, name, nbins, min, max);
657 fCFImag[i]->SetLineColor(kRed);
658 fCFImag[i]->SetMarkerColor(kRed);
660 el = fLmVals.GetElsi(i);
661 em = fLmVals.GetEmsi(i);
662 fCFReal[i]->GetYaxis()->SetTitle(Form(
"C^{%i}_{%i}", el, em));
663 fCFReal[i]->GetXaxis()->SetTitle(Hal::Femto::KinematicsToAxisLabel(fFrame, 0, 1));
664 fCFReal[i]->SetTitle(Form(
"CF^{%i}_{%i} Re", el, em));
665 fCFImag[i]->GetYaxis()->SetTitle(Form(
"C^{%i}_{%i}", el, em));
666 fCFImag[i]->SetTitle(Form(
"CF^{%i}_{%i} Im", el, em));
667 fCFImag[i]->GetXaxis()->SetTitle(Hal::Femto::KinematicsToAxisLabel(fFrame, 0, 1));
669 if ((fCovNum[0][0][0] > 0.0)) {
670 std::cout <<
"Detected calculated covariance matrix. Do not recalculate !!!" << std::endl;
673 if (fNum->GetEntries() == 0) {
674 std::cout <<
"Emtpy numerator, stop calculating cF" << std::endl;
680 solver.SetDebugBin(debug);
681 solver.SetNormalizationArea(GetNormMin(0), GetNormMax(0));
685 std::cout <<
"GLS NOT ENABLED!" << std::endl;
689 Long64_t FemtoSHCF::Merge(TCollection* collection) {
692 TIter iterator(collection);
693 while ((pack = (
FemtoSHCF*) iterator())) {
703 if (cf->fMaxJM != fMaxJM) {
704 Cout::PrintInfo(
"Not compatible maxL numbers during merging FemtoSHCF", EInfo::kError);
707 if (cf->fFrame != this->fFrame) {
708 Cout::PrintInfo(
"Incompatible frames !", EInfo::kError);
711 for (
int i = 0; i < fMaxJM; i++) {
712 fNumReal[i]->Add(cf->fNumReal[i]);
713 fNumImag[i]->Add(cf->fNumImag[i]);
714 fDenReal[i]->Add(cf->fDenReal[i]);
715 fDenImag[i]->Add(cf->fDenImag[i]);
717 fCovCf += obj->fCovCf;
718 fCovNum += obj->fCovNum;
719 fCovDen += obj->fCovDen;
724 FemtoSHCF::~FemtoSHCF() {
725 auto cleanArray = [&](TH1D** x) {
726 if (x ==
nullptr)
return;
727 for (
int i = 0; i < fMaxJM; i++)
731 cleanArray(fNumImag);
732 cleanArray(fNumReal);
733 cleanArray(fDenReal);
734 cleanArray(fDenImag);
737 if (fCfcov)
delete fCfcov;
738 if (fPainter)
delete fPainter;
741 TString FemtoSHCF::HTMLExtract(Int_t counter, TString dir)
const {
742 gSystem->MakeDirectory(Form(
"%s/divided_%i", dir.Data(), counter));
743 TString filename = Form(
"%s/divided_%i/divided.html", dir.Data(), counter);
746 HtmlRow row1(
"",
"dark_blue",
"");
750 HtmlRow row2(
"",
"dark_blue",
"");
754 HtmlRow row3(
"",
"dark_blue",
"");
759 HtmlRow row4(
"",
"dark_blue",
"");
764 HtmlRow row4(
"",
"dark_blue",
"");
770 file.AddContent(table);
772 HtmlRow row5(
"",
"dark_blue",
"");
778 for (
int i = 0; i < fLabels->GetEntries(); i++) {
779 TObjString* obj = (TObjString*) fLabels->At(i);
780 HtmlRow row6(
"",
"light_blue",
"");
785 file.AddContent(table2);
789 Bool_t batch = gROOT->IsBatch();
790 gROOT->SetBatch(kTRUE);
791 Int_t jmax = GetLMax() * 2 + 1;
792 Int_t width = 966 / jmax;
793 Int_t height = 966 / GetLMax();
794 if (height > width) height = width;
795 for (
int i = 0; i <= GetLMax(); i++) {
796 HtmlRow row7(
"",
"dark_blue",
"");
799 for (
int j = 0; j <= GetLMax(); j++) {
802 TH1D* cfr = GetHisto(I, J, kFALSE,
"re");
803 TH1D* cfi = GetHisto(I, J, kFALSE,
"im");
804 if (cfi) cfi->SetLineColor(kRed);
807 TCanvas* c1 =
new TCanvas();
810 c1->SaveAs(Form(
"%s/divided_%i/C_%i%i.png", dir.Data(), counter, I, J));
811 TString cell_content = Form(
"<a href=\"C_%i%i.png\">", I, J);
812 cell_content = cell_content
813 + Form(
"<img border=\"0\" src=\"C_%i%i.png\" "
814 "width=\"%i\" height=\"%i\"></a>",
820 cell.SetStringContent(cell_content);
825 file.AddContent(row7);
828 file.AddContent(table3);
830 gROOT->SetBatch(batch);
831 return HtmlCore::GetUrl(Form(
"divided_%i/divided.html", counter), this->ClassName());
834 void FemtoSHCF::FillNumObj(TObject* obj) {
837 std::complex<double>* YlmBuffer = fLmMath.YlmUpToL(fLmVals.GetMaxL(), pair->
GetX(), pair->
GetY(), pair->
GetZ());
838 for (
int ilm = 0; ilm < GetMaxJM(); ilm++) {
839 fNumReal[ilm]->Fill(kv, real(YlmBuffer[ilm]) * pair->
GetWeight());
840 fNumImag[ilm]->Fill(kv, -imag(YlmBuffer[ilm]) * pair->
GetWeight());
845 int nqbin = fNum->GetXaxis()->FindFixBin(kv) - 1;
846 if (nqbin < fNum->GetNbinsX()) {
848 for (
int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++) {
849 const int twoilmzero = ilmzero * 2;
850 for (
int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
851 const int twoilmprim = ilmprim * 2;
852 fCovNum[nqbin][twoilmzero][twoilmprim] += real(YlmBuffer[ilmzero]) * real(YlmBuffer[ilmprim]) * weight2;
853 fCovNum[nqbin][twoilmzero][twoilmprim + 1] += real(YlmBuffer[ilmzero]) * -imag(YlmBuffer[ilmprim]) * weight2;
854 fCovNum[nqbin][twoilmzero + 1][twoilmprim] -= imag(YlmBuffer[ilmzero]) * real(YlmBuffer[ilmprim]) * weight2;
855 fCovNum[nqbin][twoilmzero + 1][twoilmprim + 1] -= imag(YlmBuffer[ilmzero]) * -imag(YlmBuffer[ilmprim]) * weight2;
861 void FemtoSHCF::FillDenObj(TObject* obj) {
864 std::complex<double>* YlmBuffer = fLmMath.YlmUpToL(fLmVals.GetMaxL(), pair->
GetX(), pair->
GetY(), pair->
GetZ());
865 for (
int ilm = 0; ilm < GetMaxJM(); ilm++) {
866 fDenReal[ilm]->Fill(kv, real(YlmBuffer[ilm]) * pair->
GetWeight());
867 fDenImag[ilm]->Fill(kv, -imag(YlmBuffer[ilm]) * pair->
GetWeight());
872 int nqbin = fDen->GetXaxis()->FindFixBin(kv) - 1;
873 if (nqbin < fDen->GetNbinsX()) {
875 for (
int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++) {
876 const int twoilmzero = ilmzero * 2;
877 for (
int ilmprim = 0; ilmprim < GetMaxJM(); ilmprim++) {
878 const int twoilmprim = ilmprim * 2;
879 fCovDen[nqbin][twoilmzero][twoilmprim] += real(YlmBuffer[ilmzero]) * real(YlmBuffer[ilmprim]) * weight2;
880 fCovDen[nqbin][twoilmzero][twoilmprim + 1] += real(YlmBuffer[ilmzero]) * -imag(YlmBuffer[ilmprim]) * weight2;
881 fCovDen[nqbin][twoilmzero + 1][twoilmprim] -= imag(YlmBuffer[ilmzero]) * real(YlmBuffer[ilmprim]) * weight2;
882 fCovDen[nqbin][twoilmzero + 1][twoilmprim + 1] -= imag(YlmBuffer[ilmzero]) * -imag(YlmBuffer[ilmprim]) * weight2;
890 ExportIntoToFlatNum(data);
894 void FemtoSHCF::ExportIntoToFlatNum(Array_1<Float_t>* output)
const {
896 Int_t nbins = GetNum()->GetNbinsX();
897 output->
MakeBigger(fNum->GetNbinsX() * fMaxJM * 2);
898 for (
int ibin = 1; ibin <= nbins; ibin++) {
899 for (
int l = 0; l <= GetLMax(); l++) {
900 for (
int m = -l; m <= l; m++) {
901 TH1D* h = GetCFRe(l, m);
902 (*output)[bin++] = h->GetBinContent(ibin);
904 (*output)[bin++] = h->GetBinContent(ibin);
910 void FemtoSHCF::ExportIntoToFlatNumValkyria(Array_1<Float_t>* output)
const {
912 Int_t nbins = GetNum()->GetNbinsX();
913 output->MakeBigger(fNum->GetNbinsX() * fMaxJM * 2);
914 for (
int ibin = 1; ibin <= nbins; ibin++) {
915 for (
int l = 0; l <= GetLMax(); l++) {
916 for (
int m = -l; m <= l; m++) {
917 TH1D* h = GetNumRe(l, m);
918 (*output)[bin++] = h->GetBinContent(ibin);
920 (*output)[bin++] = h->GetBinContent(ibin);
926 void FemtoSHCF::MakeDummyCov() {
927 Cout::Text(
"FemtoSHCF::MakeDummyCov",
"M", kRed);
928 double nent = fNum->GetEntries();
929 double nentd = fDen->GetEntries();
930 for (
int iBin = 1; iBin <= fNum->GetNbinsX(); iBin++) {
931 const int bin = iBin - 1;
932 for (
int ilmx = 0; ilmx < GetMaxJM(); ilmx++) {
933 const int ilmx2 = ilmx * 2;
934 for (
int ilmy = 0; ilmy < GetMaxJM(); ilmy++) {
935 const int ilmy2 = ilmy * 2;
936 double t1t2rr = fNumReal[ilmx]->GetBinContent(iBin) * fNumReal[ilmy]->GetBinContent(iBin) / nent / nent;
937 double t1t2ri = fNumReal[ilmx]->GetBinContent(iBin) * fNumImag[ilmy]->GetBinContent(iBin) / nent / nent;
938 double t1t2ir = fNumImag[ilmx]->GetBinContent(iBin) * fNumReal[ilmy]->GetBinContent(iBin) / nent / nent;
939 double t1t2ii = fNumImag[ilmx]->GetBinContent(iBin) * fNumImag[ilmy]->GetBinContent(iBin) / nent / nent;
941 fCovNum[bin][ilmx2][ilmy2] = nent * (TMath::Power(fNumReal[ilmx]->GetBinError(iBin) / nent, 2) * (nent - 1) + t1t2rr);
942 fCovNum[bin][ilmx2][ilmy2 + 1] = nent * t1t2ri;
943 fCovNum[bin][ilmx2 + 1][ilmy2] = nent * t1t2ir;
944 fCovNum[bin][ilmx2 + 1][ilmy2 + 1] =
945 nent * (TMath::Power(fNumImag[ilmx]->GetBinError(iBin) / nent, 2) * (nent - 1) + t1t2rr);
947 fCovNum[bin][ilmx2][ilmy2] = nent * t1t2rr;
948 fCovNum[bin][ilmx2][ilmy2 + 1] = nent * t1t2ri;
949 fCovNum[bin][ilmx2 + 1][ilmy2] = nent * t1t2ir;
950 fCovNum[bin][ilmx2 + 1][ilmy2 + 1] = nent * t1t2ii;
952 t1t2rr = fDenReal[ilmx]->GetBinContent(iBin) * fDenReal[ilmy]->GetBinContent(iBin) / nentd / nentd;
953 t1t2ri = fDenReal[ilmx]->GetBinContent(iBin) * fDenImag[ilmy]->GetBinContent(iBin) / nentd / nentd;
954 t1t2ir = fDenImag[ilmx]->GetBinContent(iBin) * fDenReal[ilmy]->GetBinContent(iBin) / nentd / nentd;
955 t1t2ii = fDenImag[ilmx]->GetBinContent(iBin) * fDenImag[ilmy]->GetBinContent(iBin) / nentd / nentd;
958 fCovDen[bin][ilmx2 + 0][ilmy2 + 0] = nentd * t1t2rr;
959 fCovDen[bin][ilmx2 + 0][ilmy2 + 1] = nentd * t1t2ri;
960 fCovDen[bin][ilmx2 + 1][ilmy2 + 0] = nentd * t1t2ir;
961 fCovDen[bin][ilmx2 + 1][ilmy2 + 1] = nentd * t1t2ii;
968 for (
int i = 0; i < fNum->GetNbinsX(); i++) {
969 for (
int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++) {
970 for (
int ilmzero2 = 0; ilmzero2 < GetMaxJM(); ilmzero2++) {
975 for (
int ilmzero = 0; ilmzero < GetMaxJM(); ilmzero++) {
982 void FemtoSHCF::Rebin(Int_t ngroup, Option_t* opt) {
984 int oldBins = fNumReal[0]->GetNbinsX();
985 for (
int i = 0; i < fMaxJM; i++) {
986 fNumReal[i]->Rebin(ngroup);
987 fNumImag[i]->Rebin(ngroup);
988 fDenReal[i]->Rebin(ngroup);
989 fDenImag[i]->Rebin(ngroup);
991 Int_t newBins = oldBins / ngroup;
993 CovNum.
MakeBigger(oldBins, fMaxJM * 2, fMaxJM * 2);
994 CovDen.
MakeBigger(oldBins, fMaxJM * 2, fMaxJM * 2);
995 CovCf.
MakeBigger(oldBins, fMaxJM * 2, fMaxJM * 2);
996 for (
int i = 0; i < newBins; i++) {
997 for (
int a = 0; a < fMaxJM * 2; a++) {
998 for (
int b = 0; b < fMaxJM * 2; b++) {
999 for (
int j = 0; j < ngroup; j++) {
1000 CovNum[i][a][b] = CovNum[i][a][b] + fCovNum[i * ngroup + j][a][b];
1001 CovDen[i][a][b] = CovDen[i][a][b] + fCovDen[i * ngroup + j][a][b];
1002 CovCf[i][a][b] = CovCf[i][a][b] + fCovCf[i * ngroup + j][a][b];
1010 if (fCFReal[0]) { RecalculateCF(); }
1013 void FemtoSHCF::Fit(CorrFitSHCF* fit) { fit->Fit(
this); }
1015 void FemtoSHCF::FitDummy(CorrFitSHCF* fit) { fit->FitDummy(
this); }
1017 CorrFitMaskSH FemtoSHCF::MakeEmptyMask()
const {
1024 void FemtoSHCF::ImportSlice(Array_1<Float_t>* array, Int_t bin) {
1029 fNum->SetBinContent(bin, array->Get(count++));
1030 fDen->SetBinContent(bin, array->Get(count++));
1031 for (
int i = 0; i < fMaxJM; i++) {
1032 fNumReal[i]->SetBinContent(bin, array->Get(count++));
1033 fDenReal[i]->SetBinContent(bin, array->Get(count++));
1034 fNumImag[i]->SetBinContent(bin, array->Get(count++));
1035 fDenImag[i]->SetBinContent(bin, array->Get(count++));
1037 for (
int i = 0; i < fMaxJM * 2; i++) {
1038 for (
int j = 0; j < fMaxJM * 2; j++) {
1039 fCovNum[bin - 1][i][j] = array->Get(count++);
1050 for (
int i = 0; i < fMaxJM; i++) {
1051 h.Apply(*fNumReal[i]);
1052 h.Apply(*fDenReal[i]);
1053 anti.Apply(*fNumImag[i]);
1054 anti.Apply(*fDenImag[i]);
1055 if (fCFReal[i]) h.Apply(*fCFReal[i]);
1056 if (fCFImag[i]) anti.Apply(*fCFImag[i]);
1058 && (copy.GetMarkerStyle().
Find(Hal::MarkerStyle::kColor)
1059 || copy.GetLineStyle().
Find(Hal::LineStyle::kColor))) {
1061 fColRe = fCFReal[i]->GetMarkerColor();
1062 fColIm = fCFImag[i]->GetMarkerColor();
1067 TObject* FemtoSHCF::GetSpecial(TString opt)
const {
1069 if (opt ==
"painter")
return fPainter;
void MakeBigger(Int_t new_dim)
void MakeBigger(Int_t sizeA, Int_t sizeB, Int_t sizeC)
static void PrintInfo(TString text, Hal::EInfo status)
virtual void ApplyStyle(const HistoStyle &h)
void SetNorm(Double_t min, Double_t max, Int_t axis=0)
Double_t GetWeight() const
void RecalculateCF(Int_t debugBin=-1, Bool_t suwm=kFALSE)
void AddLabel(TString label)
Int_t GetEmsi(Int_t i) const
Int_t GetElsi(Int_t i) const
void SetAntiColor(Bool_t safe=kFALSE)
virtual void AddContent(const HtmlObject &obj)
virtual void AddContent(const HtmlObject &obj)
virtual void Add(const Object *pack)
Bool_t Find(Int_t bit) const