TOF PbPb utils macro update
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TOF / PbPb276 / macros / HistoUtils.C
CommitLineData
aab90846 1TH1 *
2HistoUtils_weightedmean(TH1 *h1, TH1 *h2)
3{
4
5 TH1 *ho = (TH1 *)h1->Clone("ho");
6 ho->Reset();
7 Double_t val1, val2, w1, w2, mean, meane;
8 for (Int_t i = 0; i < ho->GetNbinsX(); i++) {
9 val1 = h1->GetBinContent(i + 1);
10 w1 = 1. / (h1->GetBinError(i + 1) * h1->GetBinError(i + 1));
11 val2 = h2->GetBinContent(i + 1);
12 w2 = 1. / (h2->GetBinError(i + 1) * h2->GetBinError(i + 1));
13
14 if (val1 == 0 && val2 == 0) continue;
15
16 mean = (w1 * val1 + w2 * val2) / (w1 + w2);
17 meane = TMath::Sqrt(1. / (w1 + w2));
18
19 ho->SetBinContent(i + 1, mean);
20 ho->SetBinError(i + 1, meane);
21 }
22
23 return ho;
24}
25
26//__________________________________________________________________
27
28TH1 *
29HistoUtils_smartdifference(TH1 *hnum, TH1 *hden)
30{
31
32 TH1 *hr = (TH1 *)hnum->Clone("hr");
33 hr->Reset();
34 Double_t ref;
35 for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
36 if (hnum->GetBinError(i + 1) <= 0.) continue;
37 ref = hden->Interpolate(hr->GetBinCenter(i + 1));
38 if (ref <= 0.) continue;
39 hr->SetBinContent(i + 1, (hnum->GetBinContent(i + 1) - ref) / ref);
40 hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
41 }
42 return hr;
43}
44
45//__________________________________________________________________
46
47TH1 *
48HistoUtils_smartratio(TH1 *hnum, TH1 *hden)
49{
50
51 TH1 *hr = (TH1 *)hnum->Clone("hr");
52 hr->Reset();
53 Double_t ref;
54 for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
55 if (hnum->GetBinError(i + 1) <= 0.) continue;
56 ref = hden->Interpolate(hr->GetBinCenter(i + 1));
57 if (ref <= 0.) continue;
58 hr->SetBinContent(i + 1, hnum->GetBinContent(i + 1) / ref);
59 hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
60 }
61 return hr;
62}
63
64//__________________________________________________________________
65
66TH1 *
67HistoUtils_smartratio(TH1 *hnum, TGraph *hden)
68{
69
70 TH1 *hr = (TH1 *)hnum->Clone("hr");
71 hr->Reset();
72 Double_t ref;
73 for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
74 if (hnum->GetBinError(i + 1) <= 0.) continue;
75 ref = hden->Eval(hr->GetBinCenter(i + 1));
76 if (ref <= 0.) continue;
77 hr->SetBinContent(i + 1, hnum->GetBinContent(i + 1) / ref);
78 hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
79 }
80 return hr;
81}
82
48a73d39 83//__________________________________________________________________
84
85HistoUtils_drawthemall(const Char_t *filename, Int_t sleepms = 100)
86{
87 TFile *f1 = TFile::Open(filename);
88 TList *l1 = f1->GetListOfKeys();
89 TObject *o;
90 Char_t *name;
91 for (Int_t i = 0; i < l1->GetEntries(); i++) {
92 name = l1->At(i)->GetName();
93 o = f1->Get(name);
94 o->Draw();
95 gPad->Update();
96 gSystem->Sleep(sleepms);
97 }
98 f1->Close();
99}
100
101//__________________________________________________________________
102
103HistoUtils_autoratio(const Char_t *f1name, const Char_t *f2name, const Char_t *outname, const Char_t *title = NULL)
104{
105
106 TFile *f1 = TFile::Open(f1name);
107 TFile *f2 = TFile::Open(f2name);
108 TFile *fo = TFile::Open(outname, "RECREATE");
109 TList *l1 = f1->GetListOfKeys();
110
111
112 TH1D *h1, *h2, *hr;
113 Char_t *name;
114 for (Int_t i = 0; i < l1->GetEntries(); i++) {
115 name = l1->At(i)->GetName();
116 h1 = (TH1D *)f1->Get(name);
117 h2 = (TH1D *)f2->Get(name);
118 if (!h1 || !h2) continue;
119 hr = new TH1D(*h1);
120 hr->Divide(h2);
121 if (title)
122 hr->SetTitle(title);
123 fo->cd();
124 hr->Write();
125 }
126
127 delete hr;
128 f1->Close();
129 f2->Close();
130 fo->Close();
131
132}
133
134//__________________________________________________________________
135
136HistoUtils_autosystematics(const Char_t *f1name, const Char_t *f2name, const Char_t *outname, const Char_t *title = NULL)
137{
138
139 TFile *f1 = TFile::Open(f1name);
140 TFile *f2 = TFile::Open(f2name);
141 if (!f1 || !f1->IsOpen() || !f2 || !f2->IsOpen()) return;
142 TFile *fo = TFile::Open(outname, "RECREATE");
143 TList *l1 = f1->GetListOfKeys();
144
145
146 TH1D *h1, *h2, *hd;
147 Char_t *name;
148 for (Int_t i = 0; i < l1->GetEntries(); i++) {
149 name = l1->At(i)->GetName();
150 h1 = (TH1D *)f1->Get(name);
151 h2 = (TH1D *)f2->Get(name);
152 if (!h1 || !h2) continue;
153
154 hd = new TH1D(*h1);
155 Double_t val1, val2, vald, vale1, vale2, valde;
156 for (Int_t ii = 0; ii < hd->GetNbinsX(); ii++) {
157 val1 = h1->GetBinContent(ii + 1);
158 vale1 = h1->GetBinError(ii + 1);
159 val2 = h2->GetBinContent(ii + 1);
160 vale2 = h2->GetBinError(ii + 1);
161 if (val2 == 0.) continue;
162 vald = (val1 - val2) / val2;
163 valde = TMath::Sqrt(TMath::Abs((vale1 * vale1 - vale2 * vale2))) / val2;
164 hd->SetBinContent(ii + 1, vald);
165 hd->SetBinError(ii + 1, valde);
166 }
167
168 if (title)
169 hd->SetTitle(title);
170 fo->cd();
171 hd->Write();
172 delete hd;
173 }
174
175 f1->Close();
176 f2->Close();
177 fo->Close();
178
179}
180
181//__________________________________________________________________
182
183TH1 *
184HistoUtils_ratio(const Char_t *f1name, const Char_t *f2name, const Char_t *h1name, const Char_t *h2name)
185{
186
187 if (!f2name) f2name = f1name;
188 TFile *f1 = TFile::Open(f1name);
189 TFile *f2 = TFile::Open(f2name);
190
191 if (!h2name) h2name = h1name;
192 TH1 *h1 = (TH1 *)f1->Get(h1name);
193 TH1 *h2 = (TH1 *)f2->Get(h2name);
194
195 TH1 *hr = h1->Clone("hr");
196 hr->Sumw2();
197 hr->Divide(h2);
198
199 return hr;
200
201}
202
203//__________________________________________________________________
204
205TH1 *
206HistoUtils_systematics(const Char_t *f1name, const Char_t *f2name, const Char_t *h1name, const Char_t *h2name)
207{
208
209 if (!f2name) f2name = f1name;
210 TFile *f1 = TFile::Open(f1name);
211 TFile *f2 = TFile::Open(f2name);
212
213 if (!h2name) h2name = h1name;
214 TH1 *h1 = (TH1 *)f1->Get(h1name);
215 TH1 *h2 = (TH1 *)f2->Get(h2name);
216
217 TH1 *hd = h1->Clone("hd");
218 Double_t val1, val2, vald, vale1, vale2, valde;
219 for (Int_t i = 0; i < h1->GetNbinsX(); i++) {
220 val1 = h1->GetBinContent(i + 1);
221 vale1 = h1->GetBinError(i + 1);
222 val2 = h2->GetBinContent(i + 1);
223 vale2 = h2->GetBinError(i + 1);
224 if (val2 == 0.) continue;
225 vald = (val1 - val2) / val2;
226 valde = TMath::Sqrt(TMath::Abs((vale1 * vale1 - vale2 * vale2))) / val2;
227 hd->SetBinContent(i + 1, vald);
228 hd->SetBinError(i + 1, valde);
229 }
230
231 return hd;
232
233}
234
235//__________________________________________________________________
236
237HistoUtils_BinLogX(TH1 *h)
238{
239 TAxis *axis = h->GetXaxis();
240 Int_t bins = axis->GetNbins();
241 Axis_t from = axis->GetXmin();
242 Axis_t to = axis->GetXmax();
243 Axis_t width = (to - from) / bins;
244 Axis_t *new_bins = new Axis_t[bins + 1];
245 for (int i = 0; i <= bins; i++) new_bins[i] = TMath::Power(10, from + i * width);
246 axis->Set(bins, new_bins);
247 delete new_bins;
248}
249
250//__________________________________________________________________
251
252HistoUtils_BinNormX(TH1 *h)
253{
254 TAxis *axis = h->GetXaxis();
255 Int_t bins = axis->GetNbins();
256 Double_t c, ec;
257 Double_t w;
258 for (Int_t i = 0; i < bins; i++) {
259 c = h->GetBinContent(i + 1);
260 ec = h->GetBinError(i + 1);
261 w = axis->GetBinWidth(i + 1);
262 c /= w;
263 ec /= w;
264 h->SetBinContent(i + 1, c);
265 h->SetBinError(i + 1, ec);
266 }
267}
268
269//__________________________________________________________________
270
271TObjArray *
272HistoUtils_FitPeak(TF1 *fitFunc, TH2 *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax, Int_t minIntegral = 100., const Char_t *basename = "hParam", Bool_t monitor = kFALSE)
273{
274
275 /* gaus function if not specified */
276 if (!fitFunc)
277 fitFunc = (TF1*)gROOT->GetFunction("gaus");
278
279 /* prepare output histos */
280 TObjArray *outArray = new TObjArray();
281 Int_t npars = fitFunc->GetNpar();
282 TH1D *hpx = h->ProjectionX("hpx");
283 TH1D *hParam;
284 for (Int_t ipar = 0; ipar < npars; ipar++) {
285 hParam = new TH1D(*hpx);
286 hParam->SetName(Form("%s_%d", basename, ipar));
287 hParam->Reset();
288 outArray->Add(hParam);
289 }
290
291 /* loop over x-bins */
292 for (Int_t ibin = 0; ibin < hpx->GetNbinsX(); ibin++) {
293
294 /* check integral */
295 if (hpx->GetBinContent(ibin + 1) < minIntegral) continue;
296 /* projection y */
297 TH1D *hpy = h->ProjectionY("hpy", ibin + 1, ibin + 1);
298 /* fit peak */
299 if (HistoUtils_FitPeak(fitFunc, hpy, startSigma, nSigmaMin, nSigmaMax) != 0) {
300 delete hpy;
301 continue;
302 }
303 /* setup output histos */
304 for (Int_t ipar = 0; ipar < npars; ipar++) {
305 hParam = (TH1D *)outArray->At(ipar);
306 hParam->SetBinContent(ibin + 1, fitFunc->GetParameter(ipar));
307 hParam->SetBinError(ibin + 1, fitFunc->GetParError(ipar));
308 }
309 /* monitor */
310 if (monitor) {
311 hpy->SetMarkerStyle(20);
312 hpy->SetMarkerColor(4);
313 hpy->Draw("E1");
314 fitFunc->Draw("same");
315 gPad->Update();
316 getchar();
317 }
318 /* delete */
319 delete hpy;
320
321 }
322 /* delete */
323 delete hpx;
324 /* return output array */
325 return outArray;
326}
327
328//__________________________________________________________________
329
330Int_t
331HistoUtils_FitPeak(TF1 *fitFunc, TH1 *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax)
332{
333
334 Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
335 Double_t fitMin = fitCent - nSigmaMin * startSigma;
336 Double_t fitMax = fitCent + nSigmaMax * startSigma;
337 if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
338 if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
339 fitFunc->SetParameter(1, fitCent);
340 fitFunc->SetParameter(2, startSigma);
341 Int_t fitres = h->Fit(fitFunc, "WWq0", "", fitMin, fitMax);
342 if (fitres != 0) return fitres;
343 /* refit with better range */
344 for (Int_t i = 0; i < 3; i++) {
345 fitCent = fitFunc->GetParameter(1);
346 fitMin = fitCent - nSigmaMin * fitFunc->GetParameter(2);
347 fitMax = fitCent + nSigmaMax * fitFunc->GetParameter(2);
348 if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
349 if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
350 fitres = h->Fit(fitFunc, "q0", "", fitMin, fitMax);
351 if (fitres != 0) return fitres;
352 }
353 return fitres;
354
355}
356
357//__________________________________________________________________
358
359void
360HistoUtils_Function2Profile(TF1 *fin, TProfile *p, Int_t ntry = 10000)
361{
362
363 TF1 *f = new TF1(*fin);
364 Int_t npars = f->GetNpar();
365 Double_t par[1000];
366 Double_t pare[1000];
367 for (Int_t ipar = 0; ipar < npars; ipar++) {
368 par[ipar] = f->GetParameter(ipar);
369 pare[ipar] = f->GetParError(ipar);
370 }
371
372 Double_t x;
373 for (Int_t ibin = 0; ibin < p->GetNbinsX(); ibin++) {
374 for (Int_t itry = 0; itry < ntry; itry++) {
375 for (Int_t ipar = 0; ipar < npars; ipar++)
376 f->SetParameter(ipar, gRandom->Gaus(par[ipar], pare[ipar]));
377 x = gRandom->Uniform(p->GetXaxis()->GetBinLowEdge(ibin + 1), p->GetXaxis()->GetBinUpEdge(ibin + 1));
378 p->Fill(x, f->Eval(x));
379 }
380 }
381}