added UTILS dir with first macros
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / UTILS / YieldMean.C
1 /* definition of the fields in the histogram returned */
2 enum EValue_t {
3   kYield = 1,
4   kYieldStat,
5   kYieldSysHi,
6   kYieldSysLo,
7   kMean,
8   kMeanStat,
9   kMeanSysHi,
10   kMeanSysLo
11 };
12
13 TH1 *
14 YieldMean(TH1 *hstat, TH1 *hsys, TF1 *f = NULL, Double_t min = 0., Double_t max = 10., Double_t loprecision = 0.01, Double_t hiprecision = 0.1, Option_t *opt = "0q")
15 {
16   /* set many iterations when fitting the data so we don't
17      stop minimization with MAX_CALLS */
18   TVirtualFitter::SetMaxIterations(1000000);
19
20   /* create output histo */
21   Double_t integral, mean;
22   TH1 *hout = new TH1D("hout", "", 8, 0, 8);
23   TH1 *hlo, *hhi;
24   
25   /* create histo with stat+sys errors */
26   TH1 *htot = (TH1 *)hstat->Clone("htot");
27   for (Int_t ibin = 0; ibin < htot->GetNbinsX(); ibin++) {
28     htot->SetBinError(ibin + 1, TMath::Sqrt(hsys->GetBinError(ibin + 1) * hsys->GetBinError(ibin + 1) + hstat->GetBinError(ibin + 1) * hstat->GetBinError(ibin + 1)));
29   }
30
31   /*
32    *   measure the central value 
33    */
34   
35   do Int_t fitres = htot->Fit(f, opt);
36   while (fitres != 0);
37   hlo = YieldMean_LowExtrapolationHisto(htot, f, min, loprecision);
38   hhi = YieldMean_HighExtrapolationHisto(htot, f, max, hiprecision);
39   YieldMean_IntegralMean(htot, hlo, hhi, integral, mean);
40   hout->SetBinContent(kYield, integral);
41   hout->SetBinContent(kMean, mean);
42
43   /*
44    * STATISTICS
45    */
46   
47   TCanvas *cCanvasStat = new TCanvas("cCanvasStat");
48   cCanvasStat->Divide(2, 1);
49   
50   /*
51    * measure statistical error
52    */
53
54   /* fit with stat error */
55   do Int_t fitres = hstat->Fit(f, opt);
56   while (fitres != 0);
57   hlo = YieldMean_LowExtrapolationHisto(hstat, f, min, loprecision);
58   hhi = YieldMean_HighExtrapolationHisto(hstat, f, max, hiprecision);
59   
60   /* random generation with integration (coarse) */
61   TH1 *hIntegral_tmp = new TH1F("hIntegral_tmp", "", 1000, 0.75 * integral, 1.25 * integral);
62   TH1 *hMean_tmp = new TH1F("hMean_tmp", "", 1000, 0.75 * mean, 1.25 * mean);
63   for (Int_t irnd = 0; irnd < 100; irnd++) {
64     /* get random histogram */
65     TH1 *hrnd = YieldMean_ReturnRandom(hstat);
66     /* fit */
67     TH1 *hrndlo = YieldMean_ReturnCoherentRandom(hlo);
68     TH1 *hrndhi = YieldMean_ReturnCoherentRandom(hhi);
69     /* integrate */
70     YieldMean_IntegralMean(hrnd, hrndlo, hrndhi, integral, mean);
71     hIntegral_tmp->Fill(integral);
72     hMean_tmp->Fill(mean);
73     delete hrnd;
74     delete hrndlo;
75     delete hrndhi;
76   }
77   /* random generation with integration (fine) */
78   TH1 *hIntegral = new TH1F("hIntegral", "", 100, 
79                             hIntegral_tmp->GetMean() - 10. * hIntegral_tmp->GetRMS(),
80                             hIntegral_tmp->GetMean() + 10. * hIntegral_tmp->GetRMS());
81   TH1 *hMean = new TH1F("hMean", "", 100,
82                         hMean_tmp->GetMean() - 10. * hMean_tmp->GetRMS(),
83                         hMean_tmp->GetMean() + 10. * hMean_tmp->GetRMS());
84   for (Int_t irnd = 0; irnd < 1000; irnd++) {
85     /* get random histogram */
86     TH1 *hrnd = YieldMean_ReturnRandom(hstat);
87     /* fit */
88     TH1 *hrndlo = YieldMean_ReturnCoherentRandom(hlo);
89     TH1 *hrndhi = YieldMean_ReturnCoherentRandom(hhi);
90     /* integrate */
91     YieldMean_IntegralMean(hrnd, hrndlo, hrndhi, integral, mean);
92     hIntegral->Fill(integral);
93     hMean->Fill(mean);
94     delete hrnd;
95     delete hrndlo;
96     delete hrndhi;
97   }
98   TF1 *gaus = (TF1 *)gROOT->GetFunction("gaus");
99   
100   cCanvasStat->cd(1);
101   hIntegral->Fit(gaus, "q");
102   integral = hout->GetBinContent(kYield) * gaus->GetParameter(2) / gaus->GetParameter(1);
103   hout->SetBinContent(kYieldStat, integral);
104   
105   cCanvasStat->cd(2);
106   hMean->Fit(gaus, "q");
107   mean = hout->GetBinContent(kMean) * gaus->GetParameter(2) / gaus->GetParameter(1);
108   hout->SetBinContent(kMeanStat, mean);
109   
110   /*
111    * SYSTEMATICS
112    */
113
114   TCanvas *cCanvasSys = new TCanvas("cCanvasYieldSys");
115   cCanvasSys->Divide(2, 1);
116   cCanvasSys->cd(1)->DrawFrame(min, 1.e-3, max, 1.e3);
117   hsys->SetMarkerStyle(20);
118   hsys->SetMarkerColor(1);
119   hsys->SetMarkerSize(1);
120   hsys->Draw("same");
121   cCanvasSys->cd(2)->DrawFrame(min, 1.e-3, max, 1.e3);
122   hsys->Draw("same");
123   
124   /*
125    * systematic error high
126    */
127
128   TH1 *hhigh = YieldMean_ReturnExtremeHighHisto(hsys);
129   do Int_t fitres = hhigh->Fit(f, opt);
130   while (fitres != 0);
131   hlo = YieldMean_LowExtrapolationHisto(hhigh, f, min, loprecision);
132   hhi = YieldMean_HighExtrapolationHisto(hhigh, f, max, hiprecision);
133   YieldMean_IntegralMean(hhigh, hlo, hhi, integral, mean);
134   integral = TMath::Abs(integral - hout->GetBinContent(kYield));
135   hout->SetBinContent(kYieldSysHi, integral);
136
137   cCanvasSys->cd(1);
138   f->SetLineColor(2);
139   f->DrawCopy("same");
140   
141   /*
142    * systematic error hard
143    */
144
145   TH1 *hhard = YieldMean_ReturnExtremeHardHisto(hsys);
146   do Int_t fitres = hhard->Fit(f, opt);
147   while (fitres != 0);
148   hlo = YieldMean_LowExtrapolationHisto(hhard, f, min, loprecision);
149   hhi = YieldMean_HighExtrapolationHisto(hhard, f, max, hiprecision);
150   YieldMean_IntegralMean(hhard, hlo, hhi, integral, mean);
151   mean = TMath::Abs(mean - hout->GetBinContent(kMean));
152   hout->SetBinContent(kMeanSysHi, mean);
153
154   cCanvasSys->cd(2);
155   f->SetLineColor(2);
156   f->DrawCopy("same");
157   
158   /*
159    * systematic error low
160    */
161
162   TH1 *hlow = YieldMean_ReturnExtremeLowHisto(hsys);
163   do Int_t fitres = hlow->Fit(f, opt);
164   while (fitres != 0);
165   hlo = YieldMean_LowExtrapolationHisto(hlow, f, min, loprecision);
166   hhi = YieldMean_HighExtrapolationHisto(hlow, f, max, hiprecision);
167   YieldMean_IntegralMean(hlow, hlo, hhi, integral, mean);
168   integral = TMath::Abs(integral - hout->GetBinContent(kYield));
169   hout->SetBinContent(kYieldSysLo, integral);
170
171   cCanvasSys->cd(1);
172   f->SetLineColor(4);
173   f->DrawCopy("same");
174
175   /*
176    * systematic error soft
177    */
178
179   TH1 *hsoft = YieldMean_ReturnExtremeSoftHisto(hsys);
180   do Int_t fitres = hsoft->Fit(f, opt);
181   while (fitres != 0);
182   hlo = YieldMean_LowExtrapolationHisto(hsoft, f, min, loprecision);
183   hhi = YieldMean_HighExtrapolationHisto(hsoft, f, max, hiprecision);
184   YieldMean_IntegralMean(hsoft, hlo, hhi, integral, mean);
185   mean = TMath::Abs(mean - hout->GetBinContent(kMean));
186   hout->SetBinContent(kMeanSysLo, mean);
187
188   cCanvasSys->cd(2);
189   f->SetLineColor(4);
190   f->DrawCopy("same");
191
192   return hout;
193 }
194
195 TH1 *
196 YieldMean_LowExtrapolationHisto(TH1 *h, TF1 *f, Double_t min, Double_t binwidth = 0.01)
197 {
198   /* find lowest edge in histo */
199   Int_t binlo;
200   Double_t lo;
201   for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
202     if (h->GetBinContent(ibin) != 0.) {
203       binlo = ibin;
204       lo = h->GetBinLowEdge(ibin);
205       break;
206     }
207   }
208   
209   Int_t nbins = (lo - min) / binwidth;
210   TH1 *hlo = new TH1F("hlo", "", nbins, min, lo);
211   
212   /* integrate function in histogram bins */
213   Double_t cont, err, width;
214   for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) {
215     width = hlo->GetBinWidth(ibin + 1);
216     cont = f->Integral(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
217     err = f->IntegralError(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
218     hlo->SetBinContent(ibin + 1, cont / width);
219     hlo->SetBinError(ibin + 1, err / width);
220   }
221
222   return hlo;
223 }
224
225 TH1 *
226 YieldMean_HighExtrapolationHisto(TH1 *h, TF1 *f, Double_t max, Double_t binwidth = 0.1)
227 {
228   /* find highest edge in histo */
229   Int_t binhi;
230   Double_t hi;
231   for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
232     if (h->GetBinContent(ibin) != 0.) {
233       binhi = ibin + 1;
234       hi = h->GetBinLowEdge(ibin + 1);
235       break;
236     }
237   }
238   
239   Int_t nbins = (max - hi) / binwidth;
240   TH1 *hhi = new TH1F("hhi", "", nbins, hi, max);
241   
242   /* integrate function in histogram bins */
243   Double_t cont, err, width;
244   for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) {
245     width = hhi->GetBinWidth(ibin + 1);
246     cont = f->Integral(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
247     err = f->IntegralError(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
248     hhi->SetBinContent(ibin + 1, cont / width);
249     hhi->SetBinError(ibin + 1, err / width);
250   }
251
252   return hhi;
253 }
254
255 TH1 *
256 YieldMean_ReturnRandom(TH1 *hin)
257 {
258   TH1 *hout = (TH1 *)hin->Clone("hout");
259   hout->Reset();
260   Double_t cont, err;
261   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
262     if (hin->GetBinError(ibin + 1) <= 0.) continue;
263     cont = hin->GetBinContent(ibin + 1);
264     err = hin->GetBinError(ibin + 1);
265     hout->SetBinContent(ibin + 1, gRandom->Gaus(cont, err));
266     hout->SetBinError(ibin + 1, err);
267   }
268   return hout;
269 }
270
271 TH1 *
272 YieldMean_ReturnCoherentRandom(TH1 *hin)
273 {
274   TH1 *hout = (TH1 *)hin->Clone("hout");
275   hout->Reset();
276   Double_t cont, err, cohe;
277   cohe = gRandom->Gaus(0., 1.);
278   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
279     if (hin->GetBinError(ibin + 1) <= 0.) continue;
280     cont = hin->GetBinContent(ibin + 1);
281     err = hin->GetBinError(ibin + 1);
282     hout->SetBinContent(ibin + 1, cont + cohe * err);
283     hout->SetBinError(ibin + 1, err);
284   }
285   return hout;
286 }
287
288 TH1 *
289 YieldMean_ReturnExtremeHighHisto(TH1 *hin)
290 {
291   TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehigh", hin->GetName()));
292   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
293     if (hin->GetBinError(ibin + 1) <= 0.) continue;
294     Double_t val = hin->GetBinContent(ibin + 1);
295     Double_t err = hin->GetBinError(ibin + 1);
296     hout->SetBinContent(ibin + 1, val + err);
297   }
298   return hout;
299 }
300
301 TH1 *
302 YieldMean_ReturnExtremeLowHisto(TH1 *hin)
303 {
304   TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremelow", hin->GetName()));
305   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
306     if (hin->GetBinError(ibin + 1) <= 0.) continue;
307     Double_t val = hin->GetBinContent(ibin + 1);
308     Double_t err = hin->GetBinError(ibin + 1);
309     hout->SetBinContent(ibin + 1, val - err);
310   }
311   return hout;
312 }
313
314 TH1 *
315 YieldMean_ReturnExtremeSoftHisto(TH1 *hin)
316 {
317   return YieldMean_ReturnExtremeHisto(hin, -1.);
318 }
319
320 TH1 *
321 YieldMean_ReturnExtremeHardHisto(TH1 *hin)
322 {
323   return YieldMean_ReturnExtremeHisto(hin, 1.);
324 }
325
326 TH1 *
327 YieldMean_ReturnExtremeHisto(TH1 *hin, Float_t sign = 1.)
328 {
329   Double_t ptlow, pthigh;
330   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
331     if (hin->GetBinError(ibin + 1) <= 0.) continue;
332     ptlow = hin->GetBinLowEdge(ibin + 1);
333     break;
334   }
335   for (Int_t ibin = hin->GetNbinsX(); ibin >= 0; ibin--) {
336     if (hin->GetBinError(ibin + 1) <= 0.) continue;
337     pthigh = hin->GetBinLowEdge(ibin + 2);
338     break;
339   }
340
341   Double_t mean = hin->GetMean();
342   Double_t maxdiff = 0.;
343   TH1 *hmax = NULL;
344   for (Int_t inode = 0; inode < hin->GetNbinsX(); inode++) {
345
346     Double_t ptnode = hin->GetBinCenter(inode + 1);
347     TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehard", hin->GetName()));
348     
349     for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
350       if (hin->GetBinError(ibin + 1) <= 0.) continue;
351       Double_t val = hin->GetBinContent(ibin + 1);
352       Double_t err = hin->GetBinError(ibin + 1);
353       Double_t cen = hin->GetBinCenter(ibin + 1);
354       if (cen < ptnode)
355         err *= -1. + (cen - ptlow) / (ptnode - ptlow);
356       else
357         err *= (cen - ptnode) / (pthigh - ptnode);
358
359       hout->SetBinContent(ibin + 1, val + sign * err);
360     }
361
362     Double_t diff = TMath::Abs(mean - hout->GetMean());
363     if (diff > maxdiff) {
364       //      printf("found max at %f\n", ptnode);
365       if (hmax) delete hmax;
366       hmax = (TH1 *)hout->Clone("hmax");
367       maxdiff = diff;
368     }
369     delete hout;
370   }
371   return hmax;
372 }
373
374 YieldMean_IntegralMean(TH1 *hdata, TH1 *hlo, TH1 *hhi, Double_t &integral, Double_t &mean)
375 {
376   
377   /*
378    * compute integrals
379    */
380
381   Double_t cont, err, width, cent;
382   Double_t I = 0., IX = 0., Ierr = 0., IXerr = 0., Ilerr = 0., IXlerr = 0.;
383   Double_t M = 0., Merr = 0., Mlerr = 0., C;
384
385   /* integrate the data */
386   for (Int_t ibin = 0; ibin < hdata->GetNbinsX(); ibin++) {
387     cent = hdata->GetBinCenter(ibin + 1);
388     width = hdata->GetBinWidth(ibin + 1);
389     cont = width * hdata->GetBinContent(ibin + 1);
390     err = width * hdata->GetBinError(ibin + 1);
391     if (err <= 0.) continue;
392     I += cont;
393     IX += cont * cent;
394   }
395   /* integrate low */
396   for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) {
397     cent = hlo->GetBinCenter(ibin + 1);
398     width = hlo->GetBinWidth(ibin + 1);
399     cont = width * hlo->GetBinContent(ibin + 1);
400     err = width * hlo->GetBinError(ibin + 1);
401     if (err <= 0.) continue;
402     I += cont;
403     IX += cont * cent;
404   }
405   /* integrate high */
406   for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) {
407     cent = hhi->GetBinCenter(ibin + 1);
408     width = hhi->GetBinWidth(ibin + 1);
409     cont = width * hhi->GetBinContent(ibin + 1);
410     err = width * hhi->GetBinError(ibin + 1);
411     if (err <= 0.) continue;
412     I += cont;
413     IX += cont * cent;
414   }
415
416   /* set values */
417   integral = I;
418   mean = IX / I;
419 }