]>
Commit | Line | Data |
---|---|---|
4c0d7bc7 | 1 | // ---------------------------------------------------------------------- |
2 | // AliBWTools | |
3 | // | |
4 | // This class provides some tools which can be useful in the analsis | |
5 | // of spectra, to fit or transform histograms. See the comments of the | |
6 | // individual methods for details | |
7 | // | |
8 | // Author: M. Floris (CERN) | |
9 | // ---------------------------------------------------------------------- | |
10 | ||
11 | #include "AliBWTools.h" | |
12 | #include "TH1D.h" | |
13 | #include "TF1.h" | |
14 | #include "TH1.h" | |
15 | #include "TMath.h" | |
16 | #include "TGraphErrors.h" | |
17 | #include "TVirtualFitter.h" | |
18 | #include "TMinuit.h" | |
19 | #include "AliLog.h" | |
20 | #include <iostream> | |
21 | ||
22 | using namespace std; | |
23 | ||
c3c47c65 | 24 | TF1 * AliBWTools::fgFuncForNormalized = 0; |
4c0d7bc7 | 25 | |
26 | ClassImp(AliBWTools) | |
27 | ||
28 | AliBWTools::AliBWTools() { | |
29 | // ctor | |
30 | } | |
31 | ||
32 | AliBWTools::~AliBWTools(){ | |
33 | // dtor | |
34 | } | |
35 | ||
36 | TH1 * AliBWTools::GetdNdmtFromdNdpt(TH1 * hpt, Double_t mass) { | |
37 | // convert the x axis from pt to mt. Assumes you have 1/pt dNdpt in the histo you start with | |
38 | ||
39 | Int_t nbins = hpt->GetNbinsX(); | |
40 | Float_t * xbins = new Float_t[nbins+1]; | |
41 | for(Int_t ibins = 0; ibins <= nbins; ibins++){ | |
42 | xbins[ibins] = TMath::Sqrt(hpt->GetBinLowEdge(ibins+1)*hpt->GetBinLowEdge(ibins+1) + | |
43 | mass *mass) - mass; | |
44 | // xbins[ibins] = TMath::Sqrt(hpt->GetBinLowEdge(ibins+1)*hpt->GetBinLowEdge(ibins+1) + | |
45 | // mass *mass); | |
46 | // cout << ibins << " "<< xbins[ibins] << endl; | |
47 | ||
48 | } | |
49 | ||
50 | TH1D * hmt = new TH1D(TString(hpt->GetName())+"_mt", | |
51 | TString(hpt->GetName())+"_mt", | |
52 | nbins, xbins); | |
53 | for(Int_t ibins = 1; ibins <= nbins; ibins++){ | |
54 | hmt->SetBinContent(ibins, hpt->GetBinContent(ibins)); | |
55 | hmt->SetBinError(ibins, hpt->GetBinError(ibins)); | |
56 | ||
57 | } | |
58 | ||
59 | hmt->SetXTitle("m_{t} - m_{0} (GeV/c^{2})"); | |
60 | hmt->SetYTitle("1/m_{t} dN/dm_{t} (a.u.)"); | |
61 | ||
62 | return hmt; | |
63 | ||
64 | } | |
65 | ||
66 | TH1 * AliBWTools::GetdNdPtFromOneOverPt(TH1 * h1Pt) { | |
67 | ||
68 | // convert an histo from 1/pt dNdpt to dNdpt, using the pt at the center of the bin | |
69 | ||
70 | ||
71 | TH1 * hPt = (TH1 *) h1Pt->Clone((TString(h1Pt->GetName()) + "_inv").Data()); | |
72 | hPt->Reset(); | |
73 | ||
74 | Int_t nbinx = hPt->GetNbinsX(); | |
75 | ||
76 | for(Int_t ibinx = 1; ibinx <= nbinx; ibinx++){ | |
77 | ||
78 | Double_t cont = h1Pt->GetBinContent(ibinx); | |
79 | Double_t err = h1Pt->GetBinError(ibinx); | |
80 | ||
81 | Double_t pt = h1Pt->GetBinCenter(ibinx); | |
82 | ||
c3c47c65 | 83 | if(pt > 0) { |
4c0d7bc7 | 84 | cont *= pt; |
85 | err *= pt; | |
86 | } else { | |
87 | cont = 0; | |
88 | err = 0; | |
89 | } | |
90 | ||
91 | hPt->SetBinContent(ibinx, cont); | |
92 | hPt->SetBinError(ibinx, err); | |
93 | ||
94 | } | |
95 | ||
96 | hPt->SetXTitle("p_{t} (GeV)"); | |
97 | hPt->SetYTitle("dN/dp_{t} (GeV^{-2})"); | |
98 | ||
99 | return hPt; | |
100 | ||
101 | } | |
102 | ||
103 | ||
104 | ||
105 | ||
106 | TH1 * AliBWTools::GetOneOverPtdNdPt(TH1 * hPt) { | |
107 | ||
108 | // convert an histo from dNdpt to 1/pt dNdpt, using the pt at the center of the bin | |
109 | ||
110 | TH1 * h1Pt = (TH1 *) hPt->Clone((TString(hPt->GetName()) + "_inv").Data()); | |
111 | h1Pt->Reset(); | |
112 | ||
113 | Int_t nbinx = h1Pt->GetNbinsX(); | |
114 | ||
115 | for(Int_t ibinx = 1; ibinx <= nbinx; ibinx++){ | |
116 | ||
117 | Double_t cont = hPt->GetBinContent(ibinx); | |
118 | Double_t err = hPt->GetBinError(ibinx); | |
119 | ||
120 | Double_t pt = hPt->GetBinCenter(ibinx); | |
121 | ||
c3c47c65 | 122 | if(pt > 0) { |
4c0d7bc7 | 123 | cont /= pt; |
124 | err /= pt; | |
125 | } else { | |
126 | cont = 0; | |
127 | err = 0; | |
128 | } | |
129 | ||
130 | h1Pt->SetBinContent(ibinx, cont); | |
131 | h1Pt->SetBinError(ibinx, err); | |
132 | ||
133 | } | |
134 | ||
135 | h1Pt->SetXTitle("p_{t} (GeV)"); | |
136 | h1Pt->SetYTitle("1/p_{t} dN/dp_{t} (GeV^{-2})"); | |
137 | ||
138 | return h1Pt; | |
139 | ||
140 | } | |
141 | ||
142 | ||
143 | TGraphErrors * AliBWTools::GetGraphFromHisto(TH1F * h, Bool_t binWidth) { | |
144 | // Convert a histo to a graph | |
145 | // if binWidth is true ex is set to the bin width of the histos, otherwise it is set to zero | |
146 | Int_t nbin = h->GetNbinsX(); | |
147 | ||
148 | TGraphErrors * g = new TGraphErrors(); | |
149 | Int_t ipoint = 0; | |
150 | for(Int_t ibin = 1; ibin <= nbin; ibin++){ | |
151 | Double_t xerr = binWidth ? h->GetBinWidth(ibin)/2 : 0; | |
152 | if (h->GetBinContent(ibin)!=0) { | |
153 | g->SetPoint (ipoint, h->GetBinCenter(ibin), h->GetBinContent(ibin)); | |
154 | g->SetPointError(ipoint, xerr, h->GetBinError(ibin)); | |
155 | ipoint++; | |
156 | } | |
157 | } | |
158 | ||
159 | g->SetMarkerStyle(h->GetMarkerStyle()); | |
160 | g->SetMarkerColor(h->GetMarkerColor()); | |
161 | g->SetLineColor(h->GetLineColor()); | |
162 | ||
163 | g->SetTitle(h->GetTitle()); | |
164 | g->SetName(TString("g_")+h->GetName()); | |
165 | ||
166 | return g; | |
167 | ||
168 | } | |
169 | ||
170 | TH1F * AliBWTools::GetHistoFromGraph(TGraphErrors * g, TH1F* hTemplate) { | |
171 | ||
172 | // convert a graph to histo with the binning of hTemplate. | |
173 | // Warning: the template should be chosen | |
174 | // properly: if you have two graph points in the same histo bin this | |
175 | // won't work! | |
176 | ||
177 | TH1F * h = (TH1F*) hTemplate->Clone(TString("h_")+g->GetName()); | |
178 | h->Reset(); | |
179 | Int_t npoint = g->GetN(); | |
180 | // g->Print(); | |
181 | for(Int_t ipoint = 0; ipoint < npoint; ipoint++){ | |
182 | Float_t x = g->GetX() [ipoint]; | |
183 | Float_t y = g->GetY() [ipoint]; | |
184 | Float_t ey = g->GetEY()[ipoint]; | |
185 | Int_t bin = h->FindBin(x); | |
186 | // cout << "bin: "<< bin << " -> " << x << ", "<< y <<", " << ey << endl; | |
187 | ||
188 | h->SetBinContent(bin,y); | |
189 | h->SetBinError (bin,ey); | |
190 | } | |
191 | ||
192 | h->SetMarkerStyle(g->GetMarkerStyle()); | |
193 | h->SetMarkerColor(g->GetMarkerColor()); | |
194 | h->SetLineColor (g->GetLineColor()); | |
195 | ||
196 | ||
197 | return h; | |
198 | } | |
199 | ||
200 | TGraphErrors * AliBWTools::ConcatenateGraphs(TGraphErrors * g1,TGraphErrors * g2){ | |
201 | ||
202 | // concatenates two graphs | |
203 | ||
204 | Int_t npoint1=g1->GetN(); | |
205 | Int_t npoint2=g2->GetN(); | |
206 | ||
207 | TGraphErrors * gClone = (TGraphErrors*) g1->Clone(); | |
208 | ||
209 | // for(Int_t ipoint = 0; ipoint < npoint1; ipoint++){ | |
210 | // gClone->SetPointError(ipoint,0,g1->GetEY()[ipoint]); | |
211 | ||
212 | // } | |
213 | for(Int_t ipoint = 0; ipoint < npoint2; ipoint++){ | |
214 | gClone->SetPoint(ipoint+npoint1,g2->GetX()[ipoint],g2->GetY()[ipoint]); | |
215 | gClone->SetPointError(ipoint+npoint1,g2->GetEX()[ipoint],g2->GetEY()[ipoint]); | |
216 | // gClone->SetPointError(ipoint+npoint1,0,g2->GetEY()[ipoint]); | |
217 | } | |
218 | ||
219 | gClone->GetHistogram()->GetXaxis()->SetTimeDisplay(1); | |
220 | gClone->SetTitle(TString(gClone->GetTitle())+" + "+g2->GetTitle()); | |
221 | gClone->SetName(TString(gClone->GetName())+"_"+g2->GetName()); | |
222 | ||
223 | return gClone; | |
224 | } | |
225 | ||
226 | ||
ee08b77d | 227 | TH1F * AliBWTools::Combine3HistosWithErrors(TH1 * h1, TH1 * h2, TH1* h3, |
228 | TH1 * he1, TH1 * he2, TH1 * he3, | |
229 | TH1* htemplate, Int_t statFrom, | |
230 | Float_t renorm1, Float_t renorm2, Float_t renorm3) { | |
231 | ||
232 | // Combines 3 histos (h1,h2,h3), weighting by the errors provided in | |
233 | // he1,he2,he3, supposed to be the independent systematic errors. | |
234 | // he1,he2,he3 are also assumed to have the same binning as h1,h2,h3 | |
235 | // The combined histo must fit the template provided (no check is performed on this) | |
236 | // The histogram are supposed to come from the same (nearly) sample | |
237 | // of tracks. statFrom tells the stat error of which of the 3 is | |
238 | // going to be assigned to the combined | |
239 | // Optionally, it is possible to rescale any of the histograms. | |
240 | ||
241 | TH1F * hcomb = (TH1F*) htemplate->Clone(TString("hComb_")+h1->GetName()+"_"+h2->GetName()+h3->GetName()); | |
242 | ||
243 | // TODO: I should have used an array for h*local... | |
244 | ||
245 | // Clone histos locally to rescale them | |
246 | TH1F * h1local = (TH1F*) h1->Clone(); | |
247 | TH1F * h2local = (TH1F*) h2->Clone(); | |
248 | TH1F * h3local = (TH1F*) h3->Clone(); | |
249 | h1local->Scale(renorm1); | |
250 | h2local->Scale(renorm2); | |
251 | h3local->Scale(renorm3); | |
252 | ||
253 | TH1 * hStatError = 0; | |
254 | if (statFrom == 0) hStatError = h1; | |
255 | else if (statFrom == 1) hStatError = h2; | |
256 | else if (statFrom == 2) hStatError = h3; | |
257 | else Printf("AliBWTools::Combine3HistosWithErrors: wrong value of the statFrom parameter"); | |
258 | Printf("AliBWTools::Combine3HistosWithErrors: improve error on combined"); | |
259 | // Loop over all bins and take weighted mean of all points | |
260 | Int_t nBinComb = hcomb->GetNbinsX(); | |
261 | for(Int_t ibin = 1; ibin <= nBinComb; ibin++){ | |
262 | Int_t ibin1 = h1local->FindBin(hcomb->GetBinCenter(ibin)); | |
263 | Int_t ibin2 = h2local->FindBin(hcomb->GetBinCenter(ibin)); | |
264 | Int_t ibin3 = h3local->FindBin(hcomb->GetBinCenter(ibin)); | |
265 | Int_t ibinError = -1; // bin used to get stat error | |
266 | ||
267 | if (statFrom == 0) ibinError = ibin1; | |
268 | else if (statFrom == 1) ibinError = ibin2; | |
269 | else if (statFrom == 2) ibinError = ibin3; | |
270 | else Printf("AliBWTools::Combine3HistosWithErrors: wrong value of the statFrom parameter"); | |
271 | ||
272 | ||
273 | Double_t y[3]; | |
274 | Double_t ye[3]; | |
275 | y[0] = h1local->GetBinContent(ibin1); | |
276 | y[1] = h2local->GetBinContent(ibin2); | |
277 | y[2] = h3local->GetBinContent(ibin3); | |
278 | ye[0] = he1->GetBinError(ibin1); | |
279 | ye[1] = he2->GetBinError(ibin2); | |
280 | ye[2] = he3->GetBinError(ibin3); | |
281 | ||
282 | // Set error to 0 if content is 0 (means it was not filled) | |
283 | if(h1local->GetBinContent(ibin1)==0) ye[0] = 0; | |
284 | if(h2local->GetBinContent(ibin2)==0) ye[1] = 0; | |
285 | if(h3local->GetBinContent(ibin3)==0) ye[2] = 0; | |
286 | ||
287 | // Compute weighted mean | |
288 | Double_t mean, err; | |
289 | WeightedMean(3,y,ye,mean,err); | |
290 | ||
291 | ||
292 | // Fill combined | |
293 | // TODO: return error from weighted mean somehow... | |
294 | hcomb->SetBinContent(ibin,mean); | |
295 | Double_t statError = 0; | |
296 | if (hStatError->GetBinContent(ibinError) != 0) { | |
297 | // cout << "def" << endl; | |
298 | statError = hStatError->GetBinError(ibinError)/hStatError->GetBinContent(ibinError)*hcomb->GetBinContent(ibin); | |
299 | } | |
300 | else if (h1local->GetBinContent(ibin1) != 0) { | |
301 | // cout << "1" << endl; | |
302 | statError = h1local->GetBinError(ibin1)/h1local->GetBinContent(ibin1)*hcomb->GetBinContent(ibin); | |
303 | } | |
304 | else if (h2local->GetBinContent(ibin2) != 0) { | |
305 | // cout << "2" << endl; | |
306 | statError = h2local->GetBinError(ibin2)/h2local->GetBinContent(ibin2)*hcomb->GetBinContent(ibin); | |
307 | } | |
308 | else if (h3local->GetBinContent(ibin3) != 0) { | |
309 | // cout << "3" << endl; | |
310 | statError = h3local->GetBinError(ibin3)/h3local->GetBinContent(ibin3)*hcomb->GetBinContent(ibin); | |
311 | } | |
312 | hcomb->SetBinError (ibin,statError); | |
313 | ||
314 | // cout << "BIN " << ibin << " " << mean << " " << statError << endl; | |
315 | ||
316 | } | |
317 | ||
318 | hcomb->SetMarkerStyle(hStatError->GetMarkerStyle()); | |
319 | hcomb->SetMarkerColor(hStatError->GetMarkerColor()); | |
320 | hcomb->SetLineColor (hStatError->GetLineColor()); | |
321 | ||
322 | hcomb->SetXTitle(hStatError->GetXaxis()->GetTitle()); | |
323 | hcomb->SetYTitle(hStatError->GetYaxis()->GetTitle()); | |
324 | ||
325 | delete h1local; | |
326 | delete h2local; | |
327 | delete h3local; | |
328 | ||
329 | return hcomb; | |
330 | } | |
331 | ||
332 | ||
4c0d7bc7 | 333 | TH1F * AliBWTools::CombineHistos(TH1 * h1, TH1 * h2, TH1* htemplate, Float_t renorm1){ |
334 | ||
335 | // Combine two histos. This assumes the histos have the same binning | |
336 | // in the overlapping region. It computes the arithmetic mean in the | |
337 | // overlapping region and assigns as an error the relative error | |
338 | // h2. TO BE IMPROVED | |
339 | ||
340 | TH1F * hcomb = (TH1F*) htemplate->Clone(TString(h1->GetName())+"_"+h2->GetName()); | |
341 | ||
342 | TH1F * h1local = (TH1F*) h1->Clone(); | |
343 | h1local->Scale(renorm1); | |
344 | ||
345 | Int_t nBinComb = hcomb->GetNbinsX(); | |
346 | for(Int_t ibin = 1; ibin <= nBinComb; ibin++){ | |
347 | Int_t ibin1 = h1local->FindBin(hcomb->GetBinCenter(ibin)); | |
348 | Int_t ibin2 = h2->FindBin(hcomb->GetBinCenter(ibin)); | |
349 | ||
350 | if (h1local->GetBinContent(ibin1) == 0 && h2->GetBinContent(ibin2) == 0) { | |
351 | // None has data: go to next bin | |
352 | hcomb->SetBinContent(ibin,0); | |
353 | hcomb->SetBinError (ibin,0); | |
354 | } else if(h1local->GetBinContent(ibin1) != 0 && h2->GetBinContent(ibin2) == 0) { | |
355 | // take data from h1local: | |
356 | hcomb->SetBinContent(ibin,h1local->GetBinContent(ibin1)); | |
357 | hcomb->SetBinError (ibin,h1local->GetBinError(ibin1)); | |
358 | } else if(h1local->GetBinContent(ibin1) == 0 && h2->GetBinContent(ibin2) != 0) { | |
359 | // take data from h2: | |
360 | hcomb->SetBinContent(ibin,h2->GetBinContent(ibin2)); | |
361 | hcomb->SetBinError (ibin,h2->GetBinError(ibin2)); | |
362 | } else { | |
363 | hcomb->SetBinContent(ibin,(h1local->GetBinContent(ibin1) +h2->GetBinContent(ibin2))/2); | |
364 | // hcomb->SetBinError (ibin,h1local->GetBinError(ibin1)/h1local->GetBinContent(ibin1)*hcomb->GetBinContent(ibin)); | |
365 | hcomb->SetBinError (ibin,h2->GetBinError(ibin2)/h2->GetBinContent(ibin2)*hcomb->GetBinContent(ibin)); | |
366 | } | |
367 | ||
368 | ||
369 | } | |
370 | ||
371 | ||
372 | hcomb->SetMarkerStyle(h1local->GetMarkerStyle()); | |
373 | hcomb->SetMarkerColor(h1local->GetMarkerColor()); | |
374 | hcomb->SetLineColor (h1local->GetLineColor()); | |
375 | ||
376 | hcomb->SetXTitle(h1local->GetXaxis()->GetTitle()); | |
377 | hcomb->SetYTitle(h1local->GetYaxis()->GetTitle()); | |
378 | delete h1local; | |
379 | return hcomb; | |
380 | ||
381 | } | |
382 | ||
383 | ||
384 | void AliBWTools::GetFromHistoGraphDifferentX(TH1F * h, TF1 * f, TGraphErrors ** gBarycentre, TGraphErrors ** gXlw) { | |
385 | ||
386 | // Computes the baycentre in each bin and the xlw as defined in NIMA | |
387 | // 355 - 541 using f. Returs 2 graphs with the same y content of h | |
388 | // but with a different x (barycentre and xlw) | |
389 | ||
390 | Int_t nbin = h->GetNbinsX(); | |
391 | ||
392 | (*gBarycentre) = new TGraphErrors(); | |
393 | (*gXlw) = new TGraphErrors(); | |
394 | ||
395 | Int_t ipoint = 0; | |
396 | for(Int_t ibin = 1; ibin <= nbin; ibin++){ | |
397 | Float_t min = h->GetBinLowEdge(ibin); | |
398 | Float_t max = h->GetBinLowEdge(ibin+1); | |
399 | Double_t xerr = 0; | |
400 | Double_t xbar = f->Mean(min,max); | |
c3c47c65 | 401 | // compute xLW |
4c0d7bc7 | 402 | Double_t temp = 1./(max-min) * f->Integral(min,max); |
403 | Double_t epsilon = 0.000000001; | |
404 | Double_t increment = 0.0000000001; | |
c3c47c65 | 405 | Double_t xLW = min; |
4c0d7bc7 | 406 | |
c3c47c65 | 407 | while ((f->Eval(xLW)- temp) > epsilon) { |
408 | xLW += increment; | |
409 | if(xLW > max) { | |
410 | Printf("Cannot find xLW"); | |
4c0d7bc7 | 411 | break; |
412 | } | |
413 | } | |
414 | ||
415 | if (h->GetBinContent(ibin)!=0) { | |
416 | (*gBarycentre)->SetPoint (ipoint, xbar, h->GetBinContent(ibin)); | |
417 | (*gBarycentre)->SetPointError(ipoint, xerr, h->GetBinError(ibin)); | |
c3c47c65 | 418 | (*gXlw) ->SetPoint (ipoint, xLW, h->GetBinContent(ibin)); |
4c0d7bc7 | 419 | (*gXlw) ->SetPointError(ipoint, xerr, h->GetBinError(ibin)); |
420 | ipoint++; | |
421 | } | |
422 | } | |
423 | ||
424 | (*gBarycentre)->SetMarkerStyle(h->GetMarkerStyle()); | |
425 | (*gBarycentre)->SetMarkerColor(h->GetMarkerColor()); | |
426 | (*gBarycentre)->SetLineColor(h->GetLineColor()); | |
427 | ||
428 | (*gBarycentre)->SetTitle(h->GetTitle()); | |
429 | (*gBarycentre)->SetName(TString("g_")+h->GetName()); | |
430 | ||
431 | (*gXlw)->SetMarkerStyle(h->GetMarkerStyle()); | |
432 | (*gXlw)->SetMarkerColor(h->GetMarkerColor()); | |
433 | (*gXlw)->SetLineColor(h->GetLineColor()); | |
434 | (*gXlw)->SetTitle(h->GetTitle()); | |
435 | (*gXlw)->SetName(TString("g_")+h->GetName()); | |
436 | ||
437 | ||
438 | } | |
439 | ||
440 | ||
441 | Float_t AliBWTools::GetMean(TH1F * h, Float_t min, Float_t max) { | |
442 | ||
443 | // Get the mean of histo in a range; root is not reliable in sub ranges with variable binning. | |
444 | ||
445 | Int_t minBin = h->FindBin(min); | |
446 | Int_t maxBin = h->FindBin(max); | |
447 | ||
448 | Float_t mean = 0 ; | |
449 | Float_t integral = 0; | |
450 | for(Int_t ibin = minBin; ibin < maxBin; ibin++){ | |
451 | mean = mean + (h->GetBinCenter(ibin) * h->GetBinWidth(ibin)* h->GetBinContent(ibin)); | |
452 | integral = integral + h->GetBinContent(ibin) * h->GetBinWidth(ibin); | |
453 | } | |
454 | ||
455 | return mean/integral; | |
456 | ||
457 | ||
458 | } | |
459 | ||
460 | void AliBWTools::GetMean(TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) { | |
461 | ||
462 | // Get the mean of function in a range; | |
463 | ||
464 | return GetMoment("fMean", "x*", func, mean, error, min,max); | |
465 | ||
466 | } | |
467 | ||
468 | void AliBWTools::GetMeanSquare(TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) { | |
469 | ||
470 | // Get the mean2 of function in a range; | |
471 | ||
472 | return GetMoment("fMean2", "x*x*", func, mean, error, min,max); | |
473 | ||
474 | ||
475 | } | |
476 | ||
477 | void AliBWTools::GetMoment(TString name, TString var, TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) { | |
478 | ||
479 | // returns the integral of a function derived from func by prefixing some operation. | |
480 | // Used as a base method for mean and mean2 | |
481 | Printf("AliBWTools::GetMoment: Error on <pt> is not correct!!! It is overestimated, fix required"); | |
482 | Int_t npar = func->GetNpar(); | |
e50c2fcb | 483 | Double_t integr = func->Integral(min,max); |
4c0d7bc7 | 484 | |
485 | TF1 * f = new TF1(name, var+func->GetName()); // FIXME | |
c3c47c65 | 486 | // fgFuncForNormalized = func;// TMP: computing mean pt |
4c0d7bc7 | 487 | // TF1 * f = new TF1(name,GetNormalizedFunc,0,10,npar);// FIXME |
488 | // for(Int_t ipar = 0; ipar < npar; ipar++){ // FIXME | |
489 | // f->SetParameter(ipar,func->GetParameter(ipar)); // FIXME | |
490 | // } | |
491 | ||
492 | ||
493 | // The parameter of the function used to compute the mean should be | |
494 | // the same as the parent function: fixed if needed and they should | |
495 | // also have the same errors. | |
496 | ||
497 | // cout << "npar :" << npar << endl; | |
498 | ||
499 | for(Int_t ipar = 0; ipar < npar; ipar++){ | |
500 | Double_t parmin, parmax; | |
501 | Double_t value = func->GetParameter(ipar); | |
502 | func->GetParLimits(ipar, parmin, parmax); | |
503 | if ( parmin == parmax ) { | |
504 | if ( parmin != 0 || (parmin == 1 && value == 0) ) { | |
505 | f->FixParameter(ipar,func->GetParameter(ipar)); | |
506 | // cout << "Fixing " << ipar << "("<<value<<","<<parmin<<","<<parmax<<")"<<endl; | |
507 | } | |
508 | else { | |
509 | f->SetParError (ipar,func->GetParError(ipar) ); | |
510 | // cout << "Setting Err" << ipar << "("<<func->GetParError(ipar)<<")"<<endl; | |
511 | } | |
512 | } | |
513 | else { | |
514 | f->SetParError (ipar,func->GetParError(ipar) ); | |
515 | // cout << "Setting Err" << ipar << "("<<func->GetParError(ipar)<<")"<<endl; | |
516 | } | |
517 | } | |
518 | // f->Print(); | |
519 | // mean = f->Integral (min,max)/func->Integral(min,max); | |
520 | // error = f->IntegralError(min,max)/func->Integral(min,max); | |
e50c2fcb | 521 | mean = f->Integral (min,max)/integr; |
522 | error = f->IntegralError(min,max)/integr; | |
4c0d7bc7 | 523 | // cout << "Mean " << mean <<"+-"<< error<< endl; |
524 | // cout << "Integral Error " << error << endl; | |
525 | ||
526 | } | |
527 | ||
528 | Double_t AliBWTools::GetNormalizedFunc(double * x, double* p){ | |
529 | ||
530 | // Static function used to provide normalized pointer to a function | |
531 | ||
c3c47c65 | 532 | Int_t npar = fgFuncForNormalized->GetNpar(); |
4c0d7bc7 | 533 | for(Int_t ipar = 0; ipar < npar; ipar++){ // FIXME |
c3c47c65 | 534 | fgFuncForNormalized->SetParameter(ipar,p[ipar]); // FIXME |
4c0d7bc7 | 535 | } |
536 | ||
c3c47c65 | 537 | return x[0]*fgFuncForNormalized->Eval(x[0])/fgFuncForNormalized->Integral(0,100); |
4c0d7bc7 | 538 | |
539 | } | |
540 | ||
541 | ||
542 | Bool_t AliBWTools::Fit (TH1 * h1, TF1* func, Float_t min, Float_t max) { | |
543 | ||
544 | // Fits h1 with func, preforms several checks on the quality of the | |
545 | // fit and tries to improve it. If the fit is not good enough, it | |
546 | // returs false. | |
547 | ||
548 | Double_t amin; Double_t edm; Double_t errdef; Int_t nvpar; Int_t nparx; | |
549 | TVirtualFitter *fitter; | |
550 | cout << "--- Fitting : " << h1->GetName() << " ["<< h1->GetTitle() <<"] ---"<< endl; | |
551 | ||
552 | h1-> Fit(func,"IME0","",min,max); | |
553 | Int_t fitResult = h1-> Fit(func,"IME0","",min,max); | |
554 | // h1-> Fit(func,"0","",min,max); | |
555 | // Int_t fitResult = h1-> Fit(func,"0IE","",min,max); | |
556 | ||
557 | ||
558 | // From TH1: | |
559 | // The fitStatus is 0 if the fit is OK (i.e no error occurred). The | |
560 | // value of the fit status code is negative in case of an error not | |
561 | // connected with the minimization procedure, for example when a wrong | |
562 | // function is used. Otherwise the return value is the one returned | |
563 | // from the minimization procedure. When TMinuit (default case) or | |
564 | // Minuit2 are used as minimizer the status returned is : fitStatus = | |
565 | // migradResult + 10*minosResult + 100*hesseResult + | |
566 | // 1000*improveResult. TMinuit will return 0 (for migrad, minos, | |
567 | // hesse or improve) in case of success and 4 in case of error (see | |
568 | // the documentation of TMinuit::mnexcm). So for example, for an error | |
569 | // only in Minos but not in Migrad a fitStatus of 40 will be returned. | |
570 | // Minuit2 will return also 0 in case of success and different values | |
571 | // in migrad minos or hesse depending on the error. See in this case | |
572 | // the documentation of Minuit2Minimizer::Minimize for the | |
573 | // migradResult, Minuit2Minimizer::GetMinosError for the minosResult | |
574 | // and Minuit2Minimizer::Hesse for the hesseResult. If other | |
575 | // minimizers are used see their specific documentation for the status | |
576 | // code returned. For example in the case of Fumili, for the status | |
577 | // returned see TFumili::Minimize. | |
578 | ||
579 | ||
580 | if( gMinuit->fLimset ) { | |
581 | Printf("ERROR: AliBWTools: Parameters at limits"); | |
582 | return kFALSE; | |
583 | } | |
584 | ||
585 | ||
586 | /// | |
587 | fitter = TVirtualFitter::GetFitter(); | |
588 | Int_t fitStat = fitter->GetStats(amin, edm, errdef, nvpar, nparx); | |
589 | ||
590 | if( ( (fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ")|| (edm > 1e6) || (fitResult !=0 && fitResult < 4000) ) && | |
591 | TString(gMinuit->fCstatu) != "SUCCESSFUL" && | |
592 | TString(gMinuit->fCstatu) != "CONVERGED " ) { | |
593 | if(fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ") { | |
594 | Printf("WARNING: AliBWTools: Cannot properly compute errors"); | |
595 | if (fitStat == 0) Printf(" not calculated at all"); | |
596 | if (fitStat == 1) Printf(" approximation only, not accurate"); | |
597 | if (fitStat == 2) Printf(" full matrix, but forced positive-definite"); | |
598 | } | |
599 | ||
600 | if (edm > 1e6) { | |
601 | ||
602 | Printf("WARNING: AliBWTools: Huge EDM (%f): Fit probably failed!", edm); | |
603 | } | |
604 | if (fitResult != 0) { | |
605 | Printf("WARNING: AliBWTools: Fit Result (%d)", fitResult); | |
606 | } | |
607 | ||
608 | Printf ("AliBWTools: Trying Again with Strategy = 2"); | |
609 | ||
610 | gMinuit->Command("SET STRATEGY 2"); // more effort | |
611 | fitResult = 0; | |
612 | fitResult = h1-> Fit(func,"0","",min,max); | |
613 | fitResult = h1-> Fit(func,"IME0","",min,max); | |
614 | fitResult = h1-> Fit(func,"IME0","",min,max); | |
615 | ||
616 | fitter = TVirtualFitter::GetFitter(); | |
617 | ||
618 | fitStat = fitter->GetStats(amin, edm, errdef, nvpar, nparx); | |
619 | ||
620 | if(fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ") { | |
621 | Printf("ERROR: AliBWTools: Cannot properly compute errors"); | |
622 | if (fitStat == 0) Printf(" not calculated at all"); | |
623 | if (fitStat == 1) Printf(" approximation only, not accurate"); | |
624 | if (fitStat == 2) Printf(" full matrix, but forced positive-definite"); | |
625 | cout << "[" <<gMinuit->fCstatu<<"]" << endl; | |
626 | return kFALSE; | |
627 | } | |
628 | ||
629 | if (edm > 1e6) { | |
630 | Printf("ERROR: AliBWTools: Huge EDM (%f): Fit probably failed!", edm); | |
631 | ||
632 | return kFALSE; | |
633 | } | |
634 | if (fitResult != 0) { | |
635 | Printf("ERROR: AliBWTools: Fit Result (%d)", fitResult); | |
636 | ||
637 | return kFALSE; | |
638 | } | |
639 | ||
640 | gMinuit->Command("SET STRATEGY 1"); // back to normal value | |
641 | ||
642 | } | |
643 | ||
644 | cout << "---- FIT OK ----" << endl; | |
645 | ||
646 | return kTRUE; | |
647 | ||
648 | } | |
649 | ||
650 | Int_t AliBWTools::GetLowestNotEmptyBin(TH1*h) { | |
651 | ||
652 | // Return the index of the lowest non empty bin in the histo h | |
653 | ||
654 | Int_t nbin = h->GetNbinsX(); | |
655 | for(Int_t ibin = 1; ibin <= nbin; ibin++){ | |
656 | if(h->GetBinContent(ibin)>0) return ibin; | |
657 | } | |
658 | ||
659 | return -1; | |
660 | ||
661 | } | |
662 | ||
663 | Int_t AliBWTools::GetHighestNotEmptyBin(TH1*h) { | |
664 | ||
665 | // Return the index of the highest non empty bin in the histo h | |
666 | ||
667 | Int_t nbin = h->GetNbinsX(); | |
668 | for(Int_t ibin = nbin; ibin > 0; ibin--){ | |
669 | if(h->GetBinContent(ibin)>0) return ibin; | |
670 | } | |
671 | ||
672 | return -1; | |
673 | ||
674 | } | |
675 | ||
676 | void AliBWTools::GetResiduals(TGraphErrors * gdata, TF1 * func, TH1F ** hres, TGraphErrors ** gres) { | |
677 | ||
678 | // Returns a graph of residuals vs point and the res/err distribution | |
679 | ||
680 | Int_t npoint = gdata->GetN(); | |
681 | ||
682 | (*gres) =new TGraphErrors(); | |
683 | (*hres) = new TH1F(TString("hres_")+gdata->GetName()+"-"+func->GetName(), | |
684 | TString("hres_")+gdata->GetName()+"-"+func->GetName(), | |
685 | 20,-5,5); | |
686 | ||
687 | ||
688 | for(Int_t ipoint = 0; ipoint < npoint; ipoint++){ | |
689 | Float_t x = gdata->GetX()[ipoint]; | |
690 | Float_t res = (gdata->GetY()[ipoint] - func->Eval(x))/func->Eval(x); | |
691 | Float_t err = gdata->GetEY()[ipoint]/func->Eval(x); | |
692 | (*hres)->Fill(res/err); | |
693 | (*gres)->SetPoint(ipoint, x, res/err); | |
694 | // (*gres)->SetPointError(ipoint, 0, err); | |
695 | ||
696 | } | |
697 | ||
698 | (*gres)->SetMarkerStyle(gdata->GetMarkerStyle()); | |
699 | (*gres)->SetMarkerColor(gdata->GetMarkerColor()); | |
700 | (*gres)->SetLineColor (gdata->GetLineColor()); | |
701 | (*gres)->GetHistogram()->GetYaxis()->SetTitle("(data-function)/function"); | |
702 | (*hres)->SetMarkerStyle(gdata->GetMarkerStyle()); | |
703 | (*hres)->SetMarkerColor(gdata->GetMarkerColor()); | |
704 | (*hres)->SetLineColor (gdata->GetLineColor()); | |
705 | ||
706 | ||
707 | ||
708 | } | |
709 | ||
710 | void AliBWTools::GetResiduals(TH1F* hdata, TF1 * func, TH1F ** hres, TH1F ** hresVsBin) { | |
711 | ||
712 | // Returns an histo of residuals bin by bin and the res/err distribution | |
713 | ||
714 | if (!func) { | |
715 | Printf("AliBWTools::GetResiduals: No function provided"); | |
716 | return; | |
717 | } | |
718 | if (!hdata) { | |
719 | Printf("AliBWTools::GetResiduals: No data provided"); | |
720 | return; | |
721 | } | |
722 | ||
723 | (*hresVsBin) = (TH1F*) hdata->Clone(TString("hres_")+hdata->GetName()); | |
724 | (*hresVsBin)->Reset(); | |
725 | (*hres) = new TH1F(TString("hres_")+hdata->GetName()+"-"+func->GetName(), | |
726 | TString("hres_")+hdata->GetName()+"-"+func->GetName(), | |
727 | 20,-5,5); | |
728 | ||
729 | Int_t nbin = hdata->GetNbinsX(); | |
730 | for(Int_t ibin = 1; ibin <= nbin; ibin++){ | |
731 | if(hdata->GetBinContent(ibin)==0) continue; | |
732 | Float_t res = (hdata->GetBinContent(ibin) - func->Eval(hdata->GetBinCenter(ibin)) ) / | |
733 | func->Eval(hdata->GetBinCenter(ibin)); | |
734 | Float_t err = hdata->GetBinError (ibin) / func->Eval(hdata->GetBinCenter(ibin)); | |
735 | (*hresVsBin)->SetBinContent(ibin,res); | |
736 | (*hresVsBin)->SetBinError (ibin,err); | |
737 | (*hres)->Fill(res/err); | |
738 | ||
739 | } | |
740 | ||
741 | (*hresVsBin)->SetMarkerStyle(hdata->GetMarkerStyle()); | |
742 | (*hresVsBin)->SetMarkerColor(hdata->GetMarkerColor()); | |
743 | (*hresVsBin)->SetLineColor (hdata->GetLineColor() ); | |
744 | (*hresVsBin)->GetYaxis()->SetTitle("(data-function)/function"); | |
745 | (*hres)->SetMarkerStyle(hdata->GetMarkerStyle()); | |
746 | (*hres)->SetMarkerColor(hdata->GetMarkerColor()); | |
747 | (*hres)->SetLineColor (hdata->GetLineColor() ); | |
748 | ||
749 | } | |
750 | ||
751 | void AliBWTools::GetYield(TH1* h, TF1 * f, Double_t &yield, Double_t &yieldError, Float_t min, Float_t max, | |
752 | Double_t *partialYields, Double_t *partialYieldsErrors){ | |
753 | ||
754 | // Returns the yield extracted from the data in the histo where | |
755 | // there are points and from the fit to extrapolate, in the given | |
756 | // range. | |
757 | ||
758 | // Partial yields are also returned if the corresponding pointers are non null | |
759 | ||
760 | Int_t bin1 = h->FindBin(min); | |
761 | Int_t bin2 = h->FindBin(max); | |
762 | Float_t bin1Edge = GetLowestNotEmptyBinEdge (h); | |
763 | Float_t bin2Edge = GetHighestNotEmptyBinEdge(h); | |
764 | ||
765 | Double_t integralFromHistoError ; | |
766 | Double_t integralFromHisto = h->IntegralAndError(bin1,bin2,integralFromHistoError,"width"); | |
767 | ||
768 | Double_t integralBelow = min < bin1Edge ? f->Integral(min,bin1Edge) : 0; | |
769 | Double_t integralBelowError = min < bin1Edge ? f->IntegralError(min,bin1Edge) : 0; | |
770 | Double_t integralAbove = max > bin2Edge ? f->Integral(bin2Edge,max) : 0; | |
771 | Double_t integralAboveError = max > bin2Edge ? f->Integral(bin2Edge,max) : 0; | |
772 | ||
773 | cout << "GetYield INFO" << endl; | |
774 | cout << " " << bin1Edge << " " << bin2Edge << endl; | |
775 | cout << " " << integralFromHisto << " " << integralBelow << " " << integralAbove << endl; | |
776 | cout << " " << integralFromHistoError << " " << integralBelowError << " " << integralAboveError << endl; | |
777 | ||
778 | if(partialYields) { | |
779 | partialYields[0] = integralFromHisto; | |
780 | partialYields[1] = integralBelow; | |
781 | partialYields[2] = integralAbove; | |
782 | } | |
783 | if(partialYieldsErrors) { | |
784 | partialYieldsErrors[0] = integralFromHistoError; | |
785 | partialYieldsErrors[1] = integralBelowError; | |
786 | partialYieldsErrors[2] = integralAboveError; | |
787 | } | |
788 | yield = integralFromHisto+integralBelow+integralAbove; | |
789 | yieldError = TMath::Sqrt(integralFromHistoError*integralFromHistoError+ | |
790 | integralBelowError*integralBelowError+ | |
791 | integralAboveError*integralAboveError); | |
792 | ||
793 | } | |
794 | ||
795 | TGraphErrors * AliBWTools::DivideGraphByFunc(TGraphErrors * g, TF1 * f, Bool_t invert){ | |
796 | ||
797 | // Divides g/f. If invert == true => f/g | |
798 | ||
799 | TGraphErrors * gRatio = new TGraphErrors(); | |
800 | Int_t npoint = g->GetN(); | |
801 | for(Int_t ipoint = 0; ipoint < npoint; ipoint++){ | |
802 | Double_t x = g->GetX()[ipoint]; | |
803 | Double_t ratio = invert ? f->Eval(x)/g->GetY()[ipoint] :g->GetY()[ipoint]/f->Eval(x); | |
804 | gRatio->SetPoint (ipoint, x, ratio); | |
805 | gRatio->SetPointError(ipoint, 0, g->GetEY()[ipoint]/f->Eval(x)); | |
806 | // cout << x << " " << g->GetY()[ipoint] << " " << f->Eval(x) << endl; | |
807 | ||
808 | } | |
809 | gRatio->SetMarkerStyle(20); | |
810 | //gRatio->Print(); | |
811 | return gRatio; | |
812 | ||
813 | } | |
814 | ||
815 | TGraphErrors * AliBWTools::DivideGraphByHisto(TGraphErrors * g, TH1 * h, Bool_t invert){ | |
816 | ||
817 | // Divides g/h. If invert == true => h/g | |
818 | ||
819 | ||
820 | TGraphErrors * gRatio = new TGraphErrors(); | |
821 | Int_t npoint = g->GetN(); | |
822 | for(Int_t ipoint = 0; ipoint < npoint; ipoint++){ | |
823 | Double_t xj = g->GetX()[ipoint]; | |
824 | Double_t yj = g->GetY()[ipoint]; | |
825 | Double_t yje = g->GetEY()[ipoint]; | |
826 | ||
827 | Int_t binData = h->FindBin(xj); | |
828 | Double_t yd = h->GetBinContent(binData); | |
829 | Double_t yde = h->GetBinError(binData); | |
830 | Double_t xd = h->GetBinCenter(binData); | |
831 | ||
832 | // cout << TMath::Abs((xd-xj)/xd) << endl; | |
833 | ||
834 | ||
835 | ||
836 | if (yd == 0) continue; | |
837 | // if (binData == 28 ) cout << TMath::Abs(xd-xj)/TMath::Abs(xd) << endl; | |
838 | ||
839 | if (TMath::Abs(xd-xj)/TMath::Abs(xd) > 0.01){ | |
840 | Printf( "WARNING: bin center (%f) and x graph (%f) are more than 1 %% away, skipping",xd,xj ); | |
841 | continue; | |
842 | ||
843 | } | |
844 | ||
845 | Double_t ratio = invert ? yd/yj : yj/yd; | |
846 | ||
847 | gRatio->SetPoint(ipoint, xj, ratio); | |
848 | gRatio->SetPointError(ipoint, 0, TMath::Sqrt(yde*yde/yd/yd + yje*yje/yj/yj)*ratio); | |
849 | // gRatio->SetPointError(ipoint, 0, yje/yj * ratio); | |
850 | } | |
851 | ||
852 | return gRatio; | |
853 | ||
854 | ||
855 | } | |
856 | ||
857 | TH1F * AliBWTools::DivideHistoByFunc(TH1F * h, TF1 * f, Bool_t invert){ | |
858 | ||
859 | // Divides h/f. If invert == true => f/g | |
860 | // Performs the integral of f on the bin range to perform the ratio | |
861 | // Returns a histo with the same binnig as h | |
862 | ||
863 | // Prepare histo for ratio | |
864 | TH1F * hRatio = (TH1F*) h->Clone(TString("hRatio_")+h->GetName()+"_"+f->GetName()); | |
865 | hRatio->Reset(); | |
866 | // Set y title | |
867 | if(!invert) hRatio->SetYTitle(TString(h->GetName())+"/"+f->GetName()); | |
868 | else hRatio->SetYTitle(TString(f->GetName())+"/"+h->GetName()); | |
869 | ||
870 | // Loop over all bins | |
871 | Int_t nbin = hRatio->GetNbinsX(); | |
872 | ||
873 | for(Int_t ibin = 1; ibin <= nbin; ibin++){ | |
874 | Double_t yhisto = h->GetBinContent(ibin); | |
875 | Double_t yerror = h->GetBinError(ibin); | |
876 | Double_t xmin = h->GetBinLowEdge(ibin); | |
877 | Double_t xmax = h->GetBinLowEdge(ibin+1); | |
878 | Double_t yfunc = f->Integral(xmin,xmax)/(xmax-xmin); | |
879 | Double_t ratio = invert ? yfunc/yhisto : yhisto/yfunc ; | |
880 | Double_t error = yerror/yfunc ; | |
881 | hRatio->SetBinContent(ibin,ratio); | |
882 | hRatio->SetBinError (ibin,error); | |
883 | } | |
884 | ||
885 | return hRatio; | |
886 | ||
887 | } | |
ee08b77d | 888 | |
889 | void AliBWTools::WeightedMean(Int_t npoints, const Double_t *x, const Double_t *xerr, Double_t &mean, Double_t &meanerr){ | |
890 | ||
c3c47c65 | 891 | // Performs the weighted mean of npoints number in x with errors in xerr |
ee08b77d | 892 | |
893 | mean = 0; | |
894 | meanerr = 0; | |
895 | ||
896 | Double_t sumweight = 0; | |
897 | ||
898 | for (Int_t ipoint = 0; ipoint < npoints; ipoint++){ | |
899 | ||
900 | Double_t xerr2 = xerr[ipoint]*xerr[ipoint]; | |
901 | if(xerr2>0){ | |
902 | // cout << "xe2 " << xerr2 << endl; | |
903 | Double_t weight = 1. / xerr2; | |
904 | sumweight += weight; | |
905 | mean += weight * x[ipoint]; | |
906 | } | |
907 | ||
908 | } | |
909 | ||
910 | ||
911 | if(sumweight){ | |
912 | mean /= sumweight; | |
913 | meanerr = TMath::Sqrt(1./ sumweight); | |
914 | } | |
915 | else { | |
916 | mean = 0; | |
917 | meanerr = 0; | |
918 | } | |
919 | ||
920 | ||
921 | } |