Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
FluctuationMath.cxx
1/*
2 * FluctuationMath.cxx
3 *
4 * Created on: 31 sie 2017
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9#include "FluctuationMath.h"
10
11#include <TMath.h>
12#include <iostream>
13
14namespace Hal {
15 FluctuationMath::FluctuationMath() :
16 fN(0),
17 fTollerance(5),
18 fNError(0),
19 fOmega(0),
20 fSSigma(0),
21 fKappaSigma(0),
22 fOmegaError(0),
23 fSigmaError(0),
24 fKappaError(0),
25 fNSample(10) {
26 fM[0] = fM[1] = 0;
27 }
28
29 void FluctuationMath::GetVals(Double_t& min, Double_t& max, Double_t val, Double_t error) {
30 min = val - TMath::Abs(error);
31 max = val + TMath::Abs(error);
32 }
33
34 void FluctuationMath::SetHistoTH1D(TH1D* h, TString opt) {
35 if (opt == "") {
36 CalcError(h);
37 } else {
38 CalcErrorMagic(h);
39 }
40 }
41
42 void FluctuationMath::Calc(TH1D* h) {
43 fN = fM[2] = fM[3] = fM[4] = 0.0;
44 Double_t entries = h->GetEntries();
45 // fHisto->Scale(1.0/entries);
46 for (int i = 1; i <= h->GetNbinsX(); i++) {
47 Double_t pN = ((Double_t) h->GetBinContent(i));
48 Double_t N = (Double_t) h->GetBinCenter(i);
49 fN += N * pN;
50 }
51 fN = fN / entries;
52 for (int i = 1; i <= h->GetNbinsX(); i++) {
53 Double_t pN = ((Double_t) h->GetBinContent(i));
54 Double_t N = (Double_t) h->GetBinCenter(i);
55 for (int j = 2; j <= 4; j++) {
56 fM[j] += TMath::Power(N - fN, j) * pN;
57 }
58 }
59 for (int j = 2; j <= 4; j++) {
60 fM[j] = fM[j] / entries;
61 }
62 fOmega = GetCentralMoment(2) / fN;
63 fSSigma = GetCentralMoment(3) / GetCentralMoment(2);
64 fKappaSigma = GetCentralMoment(4) / GetCentralMoment(2) - 3.0 * GetCentralMoment(2);
65 }
66
67 Double_t FluctuationMath::GetCentralMoment(Int_t i) const {
68 if (i < 2 || i > 4) {
69 std::cout << "wrong mu index" << std::endl;
70 return 0;
71 } else {
72 return fM[i];
73 }
74 }
75
76 void FluctuationMath::CalcError(TH1D* h) {
77 TH1D* hh = (TH1D*) h->Clone("temp_boostrap");
78 FluctuationMath* stat = new FluctuationMath();
79 Double_t range_min[4], range_max[4];
80 GetVals(range_min[0], range_max[0], GetMean(), GetMean() * 2.5);
81 GetVals(range_min[1], range_max[1], GetOmega(), GetOmega() * 2.5);
82 GetVals(range_min[2], range_max[2], GetSSigma(), GetSSigma() * 2.5);
83 GetVals(range_min[3], range_max[3], GetKappaSigma(), GetKappaSigma() * 2.5);
84 Bool_t stat_ok = kFALSE;
85 while (!stat_ok) {
86 TH1D* mean = new TH1D("mean", "mean", 1000, range_min[0], range_max[0]);
87 TH1D* omega = new TH1D("omega", "omega", 1000, range_min[1], range_max[1]);
88 TH1D* sigma = new TH1D("sigma", "sigma", 1000, range_min[2], range_max[2]);
89 TH1D* kappa = new TH1D("kapa", "kappa", 1000, range_min[3], range_max[3]);
90 for (int i = 0; i < fNSample; i++) {
91 hh->Reset();
92 for (int j = 0; j < h->GetEntries(); j++) {
93 hh->Fill(h->GetRandom());
94 }
95 stat->Calc(hh);
96 omega->Fill(stat->GetOmega());
97 sigma->Fill(stat->GetSSigma());
98 kappa->Fill(stat->GetKappaSigma());
99 mean->Fill(stat->GetMean());
100 }
101 fOmegaError = omega->GetRMS();
102 fSigmaError = sigma->GetRMS();
103 fKappaError = kappa->GetRMS();
104 fNError = mean->GetRMS();
105 stat_ok = kTRUE;
106 if (fNError > (1.0 / fTollerance) * TMath::Abs(range_min[0] - range_max[0])) {
107 stat_ok = kFALSE;
108 std::cout << " need to recalcuate because of mean" << std::endl;
109 Double_t step = TMath::Abs(range_min[0] - range_max[0]);
110 std::cout << Form(" used range %4.3f - %4.3f [%4.3f ] calculated error "
111 "was %4.3f ",
112 range_min[0],
113 range_max[0],
114 step,
115 fNError)
116 << std::endl;
117 GetVals(range_min[0], range_max[0], GetMean(), GetMeanError() * fTollerance * 2.0);
118 step = TMath::Abs(range_min[0] - range_max[0]);
119 std::cout << Form(" new range %4.3f - %4.3f [%4.3f ] ", range_min[0], range_max[0], step) << std::endl;
120 }
121 if (fOmegaError > (1.0 / fTollerance) * TMath::Abs(range_min[1] - range_max[1])) {
122 stat_ok = kFALSE;
123 std::cout << " need to recalcuate because of omega" << std::endl;
124 Double_t step = TMath::Abs(range_min[1] - range_max[1]);
125 std::cout << Form(" used range %4.3f - %4.3f [%4.3f ] calculated error "
126 "was %4.3f ",
127 range_min[1],
128 range_max[1],
129 step,
130 fOmegaError)
131 << std::endl;
132 GetVals(range_min[1], range_max[1], GetOmega(), GetOmegaError() * fTollerance * 2.0);
133 step = TMath::Abs(range_min[1] - range_max[1]);
134 std::cout << Form(" new range %4.3f - %4.3f [%4.3f ] ", range_min[1], range_max[1], step) << std::endl;
135 }
136 if (fSigmaError > (1.0 / fTollerance) * TMath::Abs(range_min[2] - range_max[2])) {
137 stat_ok = kFALSE;
138 std::cout << " need to recalcuate because of sigma" << std::endl;
139 Double_t step = TMath::Abs(range_min[2] - range_max[2]);
140 std::cout << Form(" used range %4.3f - %4.3f [%4.3f ] calculated error "
141 "was %4.3f ",
142 range_min[2],
143 range_max[2],
144 step,
145 fSigmaError)
146 << std::endl;
147 step = TMath::Abs(range_min[2] - range_max[2]);
148 GetVals(range_min[2], range_max[2], GetSSigma(), GetSSigmaError() * fTollerance * 2.0);
149 std::cout << Form(" new range %4.3f - %4.3f [%4.3f ] ", range_min[2], range_max[2], step) << std::endl;
150 }
151 if (fKappaError > (1.0 / fTollerance) * TMath::Abs(range_min[3] - range_max[3])) {
152 stat_ok = kFALSE;
153 std::cout << " need to recalcuate because of kappa" << std::endl;
154 Double_t step = TMath::Abs(range_min[3] - range_max[3]);
155 std::cout << Form(" used range %4.3f - %4.3f [%4.3f ] calculated error "
156 "was %4.3f ",
157 range_min[3],
158 range_max[3],
159 step,
160 fKappaError)
161 << std::endl;
162 GetVals(range_min[3], range_max[3], GetKappaSigma(), GetKappaSigmaError() * fTollerance * 2.0);
163 step = TMath::Abs(range_min[3] - range_max[3]);
164 std::cout << Form(" new range %4.3f - %4.3f [%4.3f ] ", range_min[3], range_max[3], step) << std::endl;
165 }
166 delete omega;
167 delete mean;
168 delete sigma;
169 delete kappa;
170 }
171 delete stat;
172 delete hh;
173 }
174
175 void FluctuationMath::CalcErrorMagic(TH1D* h) {
176 TH1D* hh = (TH1D*) h->Clone("temp_boostrap");
177 FluctuationMath* stat = new FluctuationMath();
178 Double_t range_min[4], range_max[4];
179 GetVals(range_min[0], range_max[0], GetMean(), GetMean() * 2.5);
180 GetVals(range_min[1], range_max[1], GetOmega(), GetOmega() * 2.5);
181 GetVals(range_min[2], range_max[2], GetSSigma(), GetSSigma() * 2.5);
182 GetVals(range_min[3], range_max[3], GetKappaSigma(), GetKappaSigma() * 2.5);
183 Bool_t stat_ok = kFALSE;
184 while (!stat_ok) {
185 TH1D* mean = new TH1D("mean", "mean", 1000, range_min[0], range_max[0]);
186 TH1D* omega = new TH1D("omega", "omega", 1000, range_min[1], range_max[1]);
187 TH1D* sigma = new TH1D("sigma", "sigma", 1000, range_min[2], range_max[2]);
188 TH1D* kappa = new TH1D("kapa", "kappa", 1000, range_min[3], range_max[3]);
189 for (int i = 0; i < h->GetNbinsX(); i++) {
190 hh->Reset();
191 Double_t scale = h->GetEntries();
192 Double_t bin_cont = h->GetBinContent(i);
193 scale = scale / (scale - bin_cont);
194 for (int j = 0; j <= h->GetNbinsX(); j++) {
195 hh->SetBinContent(j, h->GetBinContent(j) * scale);
196 }
197 hh->SetBinContent(i, 0);
198 stat->Calc(hh);
199 omega->Fill(stat->GetOmega(), bin_cont);
200 sigma->Fill(stat->GetSSigma(), bin_cont);
201 kappa->Fill(stat->GetKappaSigma(), bin_cont);
202 mean->Fill(stat->GetMean(), bin_cont);
203 }
204 fOmegaError = omega->GetRMS();
205 fSigmaError = sigma->GetRMS();
206 fKappaError = kappa->GetRMS();
207 fNError = mean->GetRMS();
208 stat_ok = kTRUE;
209 if (fNError > (1.0 / fTollerance) * TMath::Abs(range_min[0] - range_max[0])) {
210 stat_ok = kFALSE;
211 std::cout << " need to recalcuate because of mean" << std::endl;
212 Double_t step = TMath::Abs(range_min[0] - range_max[0]);
213 std::cout << Form(" used range %4.3f - %4.3f [%4.3f ] calculated error "
214 "was %4.3f ",
215 range_min[0],
216 range_max[0],
217 step,
218 fNError)
219 << std::endl;
220 GetVals(range_min[0], range_max[0], GetMean(), GetMeanError() * fTollerance * 2.0);
221 step = TMath::Abs(range_min[0] - range_max[0]);
222 std::cout << Form(" new range %4.3f - %4.3f [%4.3f ] ", range_min[0], range_max[0], step) << std::endl;
223 }
224 if (fOmegaError > (1.0 / fTollerance) * TMath::Abs(range_min[1] - range_max[1])) {
225 stat_ok = kFALSE;
226 std::cout << " need to recalcuate because of omega" << std::endl;
227 Double_t step = TMath::Abs(range_min[1] - range_max[1]);
228 std::cout << Form(" used range %4.3f - %4.3f [%4.3f ] calculated error "
229 "was %4.3f ",
230 range_min[1],
231 range_max[1],
232 step,
233 fOmegaError)
234 << std::endl;
235 GetVals(range_min[1], range_max[1], GetOmega(), GetOmegaError() * fTollerance * 2.0);
236 step = TMath::Abs(range_min[1] - range_max[1]);
237 std::cout << Form(" new range %4.3f - %4.3f [%4.3f ] ", range_min[1], range_max[1], step) << std::endl;
238 }
239 if (fSigmaError > (1.0 / fTollerance) * TMath::Abs(range_min[2] - range_max[2])) {
240 stat_ok = kFALSE;
241 std::cout << " need to recalcuate because of sigma" << std::endl;
242 Double_t step = TMath::Abs(range_min[2] - range_max[2]);
243 std::cout << Form(" used range %4.3f - %4.3f [%4.3f ] calculated error "
244 "was %4.3f ",
245 range_min[2],
246 range_max[2],
247 step,
248 fSigmaError)
249 << std::endl;
250 step = TMath::Abs(range_min[2] - range_max[2]);
251 GetVals(range_min[2], range_max[2], GetSSigma(), GetSSigmaError() * fTollerance * 2.0);
252 std::cout << Form(" new range %4.3f - %4.3f [%4.3f ] ", range_min[2], range_max[2], step) << std::endl;
253 }
254 if (fKappaError > (1.0 / fTollerance) * TMath::Abs(range_min[3] - range_max[3])) {
255 stat_ok = kFALSE;
256 std::cout << " need to recalcuate because of kappa" << std::endl;
257 Double_t step = TMath::Abs(range_min[3] - range_max[3]);
258 std::cout << Form(" used range %4.3f - %4.3f [%4.3f ] calculated error "
259 "was %4.3f ",
260 range_min[3],
261 range_max[3],
262 step,
263 fKappaError)
264 << std::endl;
265 GetVals(range_min[3], range_max[3], GetKappaSigma(), GetKappaSigmaError() * fTollerance * 2.0);
266 step = TMath::Abs(range_min[3] - range_max[3]);
267 std::cout << Form(" new range %4.3f - %4.3f [%4.3f ] ", range_min[3], range_max[3], step) << std::endl;
268 }
269 delete omega;
270 delete mean;
271 delete sigma;
272 delete kappa;
273 }
274 delete stat;
275 delete hh;
276 }
277
278 FluctuationMath::~FluctuationMath() {
279 // TODO Auto-generated destructor stub
280 }
281} // namespace Hal