TOF PbPb utils macro update
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TOF / PbPb276 / macros / HistoUtils.C
1 TH1 *
2 HistoUtils_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
28 TH1 *
29 HistoUtils_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
47 TH1 *
48 HistoUtils_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
66 TH1 *
67 HistoUtils_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
83 //__________________________________________________________________
84
85 HistoUtils_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
103 HistoUtils_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
136 HistoUtils_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
183 TH1 *
184 HistoUtils_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
205 TH1 *
206 HistoUtils_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
237 HistoUtils_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
252 HistoUtils_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
271 TObjArray *
272 HistoUtils_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
330 Int_t
331 HistoUtils_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
359 void
360 HistoUtils_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 }