Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
FemtoImaging1D.cxx
1/*
2 * FemtoImaging1D.cxx
3 *
4 * Created on: 20 lut 2016
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9#include "FemtoImaging1D.h"
10
11#include "Cout.h"
12
13#include "TMath.h"
14#ifndef DISABLE_GSL
15#include <gsl/gsl_sf_bessel.h>
16#endif
17#include <iostream>
18namespace Hal {
19 FemtoImaging1D::FemtoImaging1D() :
20 fQMax(0.5),
21 fQMin(0.0),
22 fRMax(50),
23 fRMin(0),
24 fBins(100),
25 fKinematics(Femto::EKinematics::kLCMS),
26 fCF(NULL),
27 fR(NULL),
28 fSource(NULL) {}
29
30 void FemtoImaging1D::Init() {
31 TH1* CF = fCF->GetHist(kTRUE);
32 Double_t min = CF->GetXaxis()->GetBinLowEdge(1);
33 Double_t max = CF->GetXaxis()->GetBinUpEdge(CF->GetNbinsX());
34 Double_t bins = CF->GetNbinsX();
35 fR = new TH1D("R", "R", bins, min, max);
36 for (int i = 1; i <= bins; i++) {
37 fR->SetBinContent(i, (CF->GetBinContent(i) - 1.0));
38 if (CF->GetBinContent(i) != 0) fR->SetBinError(i, CF->GetBinError(i) * fR->GetBinContent(i) / CF->GetBinContent(i));
39 }
40 delete CF;
41 // fRMax = TMath::Pi()/fR->GetBinWidth(1)/2.0/FM_TO_GEV;
42 fSource = new TH1D("Source", "Source", fBins, fRMin, fRMax);
43 if (fCF->GetFrame() == Femto::EKinematics::kPRF) {
44 fSource->GetYaxis()->SetTitle("S(r*)");
45 fSource->GetXaxis()->SetTitle("r*");
46 fKinematics = Femto::EKinematics::kPRF;
47 } else if (fCF->GetFrame() == Femto::EKinematics::kLCMS) {
48 fSource->GetYaxis()->SetTitle("S(r)");
49 fSource->GetXaxis()->SetTitle("r");
50 fKinematics = Femto::EKinematics::kLCMS;
51 }
52 }
53
54 TH1D* FemtoImaging1D::Inverse(Femto1DCF* cf) {
55 fCF = cf;
56 Init();
57 InverseTransofrm();
58 Double_t norm = fSource->Integral(1, fSource->GetNbinsX(), "width");
59 // norm = norm *fSource->GetBinWidth(1);
60 fSource->Scale(1.0 / norm);
61 TH1D* histo = (TH1D*) fSource->Clone();
62 delete fSource;
63 fSource = NULL;
64 // histo->SetDirectory(0);
65 return histo;
66 }
67
68 void FemtoImaging1D::SetRAxis(Int_t bins, Double_t min, Double_t max) {
69 fBins = bins;
70 fRMin = min;
71 fRMax = max;
72 }
73
74 TH1D* FemtoImaging1D::Inverse(TH1D* CF, Bool_t normalize) {
75 Double_t min = CF->GetXaxis()->GetBinLowEdge(1);
76 Double_t max = CF->GetXaxis()->GetBinUpEdge(CF->GetNbinsX());
77 Double_t bins = CF->GetNbinsX();
78
79 fR = new TH1D("R", "R", bins, min, max);
80 for (int i = 1; i <= bins; i++) {
81 fR->SetBinContent(i, (CF->GetBinContent(i) - 1.0));
82 if (CF->GetBinContent(i) != 0) fR->SetBinError(i, CF->GetBinError(i) * fR->GetBinContent(i) / CF->GetBinContent(i));
83 }
84 // fRMax = TMath::Pi()/fR->GetBinWidth(1)/2.0/FM_TO_GEV;
85 fSource = new TH1D("Source", "Source", fBins, fRMin, fRMax);
86 InverseTransofrm();
87 Double_t norm = fSource->Integral(1, fSource->GetNbinsX(), "width");
88 // norm = norm *fSource->GetBinWidth(1);
89 if (normalize) fSource->Scale(1.0 / norm);
90 TH1D* histo = (TH1D*) fSource->Clone();
91 delete fSource;
92 fR = NULL;
93 fSource = NULL;
94 // histo->SetDirectory(0);
95 return histo;
96 }
97
98 TH1D* FemtoImaging1D::MagicInverse(TH1D* CF, Double_t l, Bool_t /*normalize*/) {
99 Double_t min = CF->GetXaxis()->GetBinLowEdge(1);
100 Double_t max = CF->GetXaxis()->GetBinUpEdge(CF->GetNbinsX());
101 Double_t bins = CF->GetNbinsX();
102
103 fR = new TH1D("R", "R", bins, min, max);
104 for (int i = 1; i <= bins; i++) {
105 fR->SetBinContent(i, (CF->GetBinContent(i)));
106 if (CF->GetBinContent(i) != 0) fR->SetBinError(i, CF->GetBinError(i) * fR->GetBinContent(i) / CF->GetBinContent(i));
107 }
108 // fRMax = TMath::Pi()/fR->GetBinWidth(1)/2.0/FM_TO_GEV;
109 fSource = new TH1D("Source", "Source", fBins, fRMin, fRMax);
110 MagicInverseTransform(l);
111 TH1D* histo = (TH1D*) fSource->Clone();
112 delete fSource;
113 fR = NULL;
114 fSource = NULL;
115 // histo->SetDirectory(0);
116 return histo;
117 }
118
119 void FemtoImaging1D::MagicInverseTransform(Int_t l) {
120 Double_t power = TMath::Power(-1, l / 2) * TMath::Power(TMath::Pi(), 3);
121 for (double i = 1; i <= fSource->GetNbinsX(); i++) {
122 if (l % 2 == 1) continue;
123 Double_t r = fSource->GetBinCenter(i) * Femto::FmToGeV();
124 Double_t sum = 0.0;
125 Double_t errors = 0.0;
126 for (int j = fR->GetXaxis()->FindBin(fQMin); j < fR->GetXaxis()->FindBin(fQMax); j++) {
127 Double_t q = fR->GetBinCenter(j);
128 if (fKinematics == Femto::EKinematics::kPRF) q += q;
129 Double_t R = fR->GetBinContent(j);
130 // if(TMath::Abs(R)<=0.01) continue;
131 Double_t dR = fR->GetBinError(j);
132 Double_t integral = q * q * BesselJn(l, 2 * q * r) / Femto::FmToGeV() / Femto::FmToGeV();
133 sum += R * integral;
134 errors += dR * integral;
135 }
136 // Double_t end = fR->GetBinCenter(fR->FindBin(fQMax));
137 // Double_t numerical_end =TMath::Ceil( fQMax/TMath::Pi()/2.0/r);
138 if (r != 0) {
139 // Double_t correction =
140 // -1.0/r*(TMath::Cos(r*numerical_end)-TMath::Cos(r*end)); sum+=correction;
141 }
142
143 // Double_t ep = r*fR->GetBinWidth(1);
144 // empiric estimation of error
145 if (sum < 0) { sum = 0; }
146 // Double_t dFourierError = (0.05+ep*0.22+ep*ep*0.084)*sum;
147 // errors= TMath::Sqrt(errors*errors+dFourierError*dFourierError);
148 // if(r!=0){
149 fSource->SetBinContent(i, sum * power);
150 fSource->SetBinError(i, errors);
151 // }
152 }
153 delete fR;
154 fR = NULL;
155 }
156
157 Double_t FemtoImaging1D::BesselJ0(double x) const {
158#ifndef DISABLE_GSL
159 return gsl_sf_bessel_j0(x);
160#else
161 std::cout << "GSL NOT SUPPORTED" << std::endl;
162 return 0;
163#endif
164 }
165
166 Double_t FemtoImaging1D::BesselJn(int n, double x) const {
167 // gsl_sf_bessel_j
168 // return gsl_sf_bessel_Jn(n,x);
169 // gsl_sf_bessel_Inu()
170#ifndef DISABLE_GSL
171 return jn(n, x);
172#else
173 std::cout << "GSL NOT SUPPORTED" << std::endl;
174 return 0;
175#endif
176 }
177
178 double FemtoImaging1D::BesselJ1(double x) const {
179#ifndef DISABLE_GSL
180 return gsl_sf_bessel_j1(x);
181#else
182 std::cout << "GSL NOT SUPPORTED" << std::endl;
183 return 0;
184#endif
185 }
186
187 FemtoImaging1D::~FemtoImaging1D() {}
188
189 void FemtoImaging1D::InverseTransofrm() {
190 for (double i = 1; i <= fSource->GetNbinsX(); i++) {
191 Double_t r = fSource->GetBinCenter(i) * Femto::FmToGeV();
192 Double_t sum = 0.0;
193 Double_t errors = 0.0;
194 for (int j = fR->GetXaxis()->FindBin(fQMin); j < fR->GetXaxis()->FindBin(fQMax); j++) {
195 Double_t q = fR->GetBinCenter(j);
196 if (fKinematics == Femto::EKinematics::kPRF) q += q;
197 Double_t R = fR->GetBinContent(j);
198 if (R == 0) continue;
199 Double_t dR = fR->GetBinError(j);
200 Double_t integral = q * TMath::Sin(q * r * 2);
201 sum += R * integral;
202 errors += dR * integral;
203 }
204 // Double_t end = fR->GetBinCenter(fR->FindBin(fQMax));
205 // Double_t numerical_end =TMath::Ceil( fQMax/TMath::Pi()/2.0/r);
206 if (r != 0) {
207 // Double_t correction =
208 // -1.0/r*(TMath::Cos(r*numerical_end)-TMath::Cos(r*end)); sum+=correction;
209 }
210
211 Double_t ep = r * fR->GetBinWidth(1);
212 // empiric estimation of error
213 if (sum < 0) { sum = 0; }
214 Double_t dFourierError = (0.05 + ep * 0.22 + ep * ep * 0.084) * sum;
215 errors = TMath::Sqrt(errors * errors + dFourierError * dFourierError);
216 if (r != 0) {
217 fSource->SetBinContent(i, sum / r);
218 fSource->SetBinError(i, errors / r);
219 }
220 }
221 delete fR;
222 fR = NULL;
223 }
224} // namespace Hal