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