AliPhysics  fde8a9f (fde8a9f)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliPWGHistoTools.cxx
Go to the documentation of this file.
1 // ----------------------------------------------------------------------
2 // AliPWGHistoTools
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 "AliPWGHistoTools.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 "TSpline.h"
21 #include <iostream>
22 #include "TGraphAsymmErrors.h"
23 
24 using namespace std;
25 
26 
28 
29 TF1 * AliPWGHistoTools::fdNdptForETCalc = 0;
30 
32  // ctor
33 }
34 
36  // dtor
37 }
38 
40  // convert the x axis from pt to mt. Assumes you have 1/pt dNdpt in the histo you start with
41 
42  Int_t nbins = hpt->GetNbinsX();
43  Float_t * xbins = new Float_t[nbins+1];
44  for(Int_t ibins = 0; ibins <= nbins; ibins++){
45  xbins[ibins] = TMath::Sqrt(hpt->GetBinLowEdge(ibins+1)*hpt->GetBinLowEdge(ibins+1) +
46  mass *mass) - mass;
47  // // xbins[ibins] = TMath::Sqrt(hpt->GetBinLowEdge(ibins+1)*hpt->GetBinLowEdge(ibins+1) +
48  // // mass *mass);
49  // cout << ibins << " "<< xbins[ibins] << endl;
50 
51  }
52 
53  TH1D * hmt = new TH1D(TString(hpt->GetName())+"_mt",
54  TString(hpt->GetName())+"_mt",
55  nbins, xbins);
56  for(Int_t ibins = 1; ibins <= nbins; ibins++){
57  hmt->SetBinContent(ibins, hpt->GetBinContent(ibins));
58  hmt->SetBinError(ibins, hpt->GetBinError(ibins));
59 
60  }
61 
62  hmt->SetXTitle("m_{t} - m_{0} (GeV/c^{2})");
63  hmt->SetYTitle("1/m_{t} dN/dm_{t} (a.u.)");
64  hmt->SetMarkerStyle(hpt->GetMarkerStyle());
65  hmt->SetMarkerColor(hpt->GetMarkerColor());
66  hmt->SetLineColor(hpt->GetLineColor());
67 
68  return hmt;
69 
70 }
71 
73  // convert the x axis from mt to pt. Assumes you have 1/mt dNdmt in the histo you start with
74 
75  Int_t nbins = hmt->GetNbinsX();
76  Float_t * xbins = new Float_t[nbins+1];
77  for(Int_t ibins = 0; ibins <= nbins; ibins++){
78  xbins[ibins] = TMath::Sqrt((hmt->GetBinLowEdge(ibins+1)+mass)*(hmt->GetBinLowEdge(ibins+1)+mass) -
79  mass *mass);
80  xbins[ibins] = Float_t(TMath::Nint(xbins[ibins]*100))/100;
81  // // xbins[ibins] = TMath::Sqrt(hmt->GetBinLowEdge(ibins+1)*hmt->GetBinLowEdge(ibins+1) +
82  // // mass *mass);
83  cout << ibins << " "<< xbins[ibins] << endl;
84 
85  }
86 
87  TH1D * hptL = new TH1D(TString(hmt->GetName())+"_pt",
88  TString(hmt->GetName())+"_pt",
89  nbins, xbins);
90  for(Int_t ibins = 1; ibins <= nbins; ibins++){
91  hptL->SetBinContent(ibins, hmt->GetBinContent(ibins));
92  hptL->SetBinError(ibins, hmt->GetBinError(ibins));
93 
94  }
95 
96  hptL->SetXTitle("p_{t} (GeV/c)");
97  hptL->SetYTitle("1/p_{t} dN/dp_{t} (a.u.)");
98  hptL->SetMarkerStyle(hmt->GetMarkerStyle());
99  hptL->SetMarkerColor(hmt->GetMarkerColor());
100  hptL->SetLineColor(hmt->GetLineColor());
101 
102  return hptL;
103 
104 }
105 
106 
108 
109  // convert an histo from 1/pt dNdpt to dNdpt, using the pt at the center of the bin
110 
111 
112  TH1 * hPt = (TH1 *) h1Pt->Clone((TString(h1Pt->GetName()) + "_inv").Data());
113  hPt->Reset();
114 
115  Int_t nbinx = hPt->GetNbinsX();
116 
117  for(Int_t ibinx = 1; ibinx <= nbinx; ibinx++){
118 
119  Double_t cont = h1Pt->GetBinContent(ibinx);
120  Double_t err = h1Pt->GetBinError(ibinx);
121 
122  Double_t pt = h1Pt->GetBinCenter(ibinx);
123 
124  if(pt > 0) {
125  cont *= pt;
126  err *= pt;
127  } else {
128  cont = 0;
129  err = 0;
130  }
131 
132  hPt->SetBinContent(ibinx, cont);
133  hPt->SetBinError(ibinx, err);
134 
135  }
136 
137  hPt->SetXTitle("p_{t} (GeV)");
138  hPt->SetYTitle("dN/dp_{t} (GeV^{-2})");
139 
140  return hPt;
141 
142 }
143 
144 
145 
146 
148 
149  // convert an histo from dNdpt to 1/pt dNdpt, using the pt at the center of the bin
150 
151  TH1 * h1Pt = (TH1 *) hPt->Clone((TString(hPt->GetName()) + "_inv").Data());
152  h1Pt->Reset();
153 
154  Int_t nbinx = h1Pt->GetNbinsX();
155 
156  for(Int_t ibinx = 1; ibinx <= nbinx; ibinx++){
157 
158  Double_t cont = hPt->GetBinContent(ibinx);
159  Double_t err = hPt->GetBinError(ibinx);
160 
161  Double_t pt = hPt->GetBinCenter(ibinx);
162 
163  if(pt > 0) {
164  cont /= pt;
165  err /= pt;
166  } else {
167  cont = 0;
168  err = 0;
169  }
170 
171  h1Pt->SetBinContent(ibinx, cont);
172  h1Pt->SetBinError(ibinx, err);
173 
174  }
175 
176  h1Pt->SetXTitle("p_{t} (GeV)");
177  h1Pt->SetYTitle("1/p_{t} dN/dp_{t} (GeV^{-2})");
178 
179  return h1Pt;
180 
181 }
182 
183 
185  // Convert a histo to a graph
186  // if binWidth is true ex is set to the bin width of the histos, otherwise it is set to zero
187  Int_t nbin = h->GetNbinsX();
188 
189  TGraphErrors * g = new TGraphErrors();
190  Int_t ipoint = 0;
191  for(Int_t ibin = 1; ibin <= nbin; ibin++){
192  Double_t xerr = binWidth ? h->GetBinWidth(ibin)/2 : 0;
193  if (h->GetBinContent(ibin)) {
194  g->SetPoint (ipoint, h->GetBinCenter(ibin), h->GetBinContent(ibin));
195  g->SetPointError(ipoint, xerr, h->GetBinError(ibin));
196  ipoint++;
197  }
198  }
199 
200  g->SetMarkerStyle(h->GetMarkerStyle());
201  g->SetMarkerColor(h->GetMarkerColor());
202  g->SetLineColor(h->GetLineColor());
203  g->SetLineStyle(h->GetLineStyle());
204  g->SetLineWidth(h->GetLineWidth());
205 
206  g->SetTitle(h->GetTitle());
207  g->SetName(TString("g_")+h->GetName());
208 
209  return g;
210 
211 }
212 
213 TH1F * AliPWGHistoTools::GetHistoFromGraph(const TGraphErrors * g, const TH1F* hTemplate) {
214 
215  // convert a graph to histo with the binning of hTemplate.
216  // Warning: the template should be chosen
217  // properly: if you have two graph points in the same histo bin this
218  // won't work!
219 
220  TH1F * h = (TH1F*) hTemplate->Clone(TString("h_")+g->GetName());
221  h->Reset();
222  Int_t npoint = g->GetN();
223  // g->Print();
224  for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
225  Float_t x = g->GetX() [ipoint];
226  Float_t y = g->GetY() [ipoint];
227  Float_t ey = 0;
228  if(g->InheritsFrom("TGraphErrors"))
229  ey = g->GetEY()[ipoint];
230  Int_t bin = h->FindBin(x);
231  // cout << "bin: "<< bin << " -> " << x << ", "<< y <<", " << ey << endl;
232 
233  h->SetBinContent(bin,y);
234  h->SetBinError (bin,ey);
235  }
236 
237  h->SetMarkerStyle(g->GetMarkerStyle());
238  h->SetMarkerColor(g->GetMarkerColor());
239  h->SetLineColor (g->GetLineColor());
240 
241 
242  return h;
243 }
244 
246 
247  // concatenates two graphs
248 
249  Int_t npoint1=g1->GetN();
250  Int_t npoint2=g2->GetN();
251 
252  TGraphErrors * gClone = (TGraphErrors*) g1->Clone();
253 
254 // for(Int_t ipoint = 0; ipoint < npoint1; ipoint++){
255 // gClone->SetPointError(ipoint,0,g1->GetEY()[ipoint]);
256 
257 // }
258  for(Int_t ipoint = 0; ipoint < npoint2; ipoint++){
259  gClone->SetPoint(ipoint+npoint1,g2->GetX()[ipoint],g2->GetY()[ipoint]);
260  gClone->SetPointError(ipoint+npoint1,g2->GetEX()[ipoint],g2->GetEY()[ipoint]);
261  // gClone->SetPointError(ipoint+npoint1,0,g2->GetEY()[ipoint]);
262  }
263 
264  gClone->GetHistogram()->GetXaxis()->SetTimeDisplay(1);
265  gClone->SetTitle(TString(gClone->GetTitle())+" + "+g2->GetTitle());
266  gClone->SetName(TString(gClone->GetName())+"_"+g2->GetName());
267 
268  return gClone;
269 }
270 
271 
272 TH1F * AliPWGHistoTools::Combine3HistosWithErrors(const TH1 * h1, const TH1 * h2, const TH1* h3,
273  TH1 * he1, TH1 * he2, TH1 * he3,
274  const TH1* htemplate, Int_t statFrom,
275  Float_t renorm1, Float_t renorm2, Float_t renorm3,
276  TH1 ** hSyst, Bool_t errorFromBinContent) {
277 
278  // Combines 3 histos (h1,h2,h3), weighting by the errors provided in
279  // he1,he2,he3, supposed to be the independent systematic errors.
280  // he1,he2,he3 are also assumed to have the same binning as h1,h2,h3
281  // The combined histo must fit the template provided (no check is performed on this)
282  // The histogram are supposed to come from the same (nearly) sample
283  // of tracks. statFrom tells the stat error of which of the 3 is
284  // going to be assigned to the combined
285  // Optionally, it is possible to rescale any of the histograms.
286  // if hSyst is give, the histo is filled with combined syst error vs pt
287  // if errorFromBinContent is true, the weights are taken from the he* content rather than errors
288  TH1F * hcomb = (TH1F*) htemplate->Clone(TString("hComb_")+h1->GetName()+"_"+h2->GetName()+h3->GetName());
289 
290  // TODO: I should have used an array for h*local...
291 
292  // Clone histos locally to rescale them
293  TH1F * h1local = (TH1F*) h1->Clone();
294  TH1F * h2local = (TH1F*) h2->Clone();
295  TH1F * h3local = (TH1F*) h3->Clone();
296  h1local->Scale(renorm1);
297  h2local->Scale(renorm2);
298  h3local->Scale(renorm3);
299 
300  const TH1 * hStatError = 0;
301  if (statFrom == 0) hStatError = h1;
302  else if (statFrom == 1) hStatError = h2;
303  else if (statFrom == 2) hStatError = h3;
304  else {
305  Printf("AliPWGHistoTools::Combine3HistosWithErrors: wrong value of the statFrom parameter");
306  return NULL;
307  }
308  Printf("AliPWGHistoTools::Combine3HistosWithErrors: improve error on combined");
309  // Loop over all bins and take weighted mean of all points
310  Int_t nBinComb = hcomb->GetNbinsX();
311  for(Int_t ibin = 1; ibin <= nBinComb; ibin++){
312  Int_t ibin1 = h1local->FindBin(hcomb->GetBinCenter(ibin));
313  Int_t ibin2 = h2local->FindBin(hcomb->GetBinCenter(ibin));
314  Int_t ibin3 = h3local->FindBin(hcomb->GetBinCenter(ibin));
315  Int_t ibinError = -1; // bin used to get stat error
316 
317  if (statFrom == 0) ibinError = ibin1;
318  else if (statFrom == 1) ibinError = ibin2;
319  else if (statFrom == 2) ibinError = ibin3;
320  else Printf("AliPWGHistoTools::Combine3HistosWithErrors: wrong value of the statFrom parameter");
321 
322 
323  Double_t y[3];
324  Double_t ye[3];
325  y[0] = h1local->GetBinContent(ibin1);
326  y[1] = h2local->GetBinContent(ibin2);
327  y[2] = h3local->GetBinContent(ibin3);
328  if (errorFromBinContent) {
329  ye[0] = he1->GetBinContent(he1->FindBin(hcomb->GetBinCenter(ibin)));
330  ye[1] = he2->GetBinContent(he2->FindBin(hcomb->GetBinCenter(ibin)));
331  ye[2] = he3->GetBinContent(he3->FindBin(hcomb->GetBinCenter(ibin)));
332  } else {
333  ye[0] = he1->GetBinError(ibin1);
334  ye[1] = he2->GetBinError(ibin2);
335  ye[2] = he3->GetBinError(ibin3);
336  }
337  // Set error to 0 if content is 0 (means it was not filled)
338  if(!h1local->GetBinContent(ibin1)) ye[0] = 0;
339  if(!h2local->GetBinContent(ibin2)) ye[1] = 0;
340  if(!h3local->GetBinContent(ibin3)) ye[2] = 0;
341 
342  // Compute weighted mean
343  // cout << "Bins: "<< ibin1 << " " << ibin2 << " " << ibin3 << endl;
344  Double_t mean, err;
345  WeightedMean(3,y,ye,mean,err);
346 
347 
348  // Fill combined
349  hcomb->SetBinContent(ibin,mean);
350  Double_t statError = 0;
351  if (hStatError->GetBinContent(ibinError)) {
352  // cout << "def" << endl;
353  statError = hStatError->GetBinError(ibinError)/hStatError->GetBinContent(ibinError)*hcomb->GetBinContent(ibin);
354  }
355  else if (h1local->GetBinContent(ibin1)) {
356  // cout << "1" << endl;
357  statError = h1local->GetBinError(ibin1)/h1local->GetBinContent(ibin1)*hcomb->GetBinContent(ibin);
358  }
359  else if (h2local->GetBinContent(ibin2)) {
360  // cout << "2" << endl;
361  statError = h2local->GetBinError(ibin2)/h2local->GetBinContent(ibin2)*hcomb->GetBinContent(ibin);
362  }
363  else if (h3local->GetBinContent(ibin3)) {
364  // cout << "3" << endl;
365  statError = h3local->GetBinError(ibin3)/h3local->GetBinContent(ibin3)*hcomb->GetBinContent(ibin);
366  }
367  hcomb->SetBinError (ibin,statError);
368  if(hSyst) (*hSyst)->SetBinContent(ibin,err);
369  // cout << "BIN " << ibin << " " << mean << " " << statError << endl;
370 
371  }
372 
373  hcomb->SetMarkerStyle(hStatError->GetMarkerStyle());
374  hcomb->SetMarkerColor(hStatError->GetMarkerColor());
375  hcomb->SetLineColor (hStatError->GetLineColor());
376 
377  hcomb->SetXTitle(hStatError->GetXaxis()->GetTitle());
378  hcomb->SetYTitle(hStatError->GetYaxis()->GetTitle());
379 
380  delete h1local;
381  delete h2local;
382  delete h3local;
383 
384  return hcomb;
385 }
386 
387 
388 void AliPWGHistoTools::GetMeanDataAndExtrapolation(const TH1 * hData, TF1 * fExtrapolation, Double_t &mean, Double_t &error, Float_t min, Float_t max){
389  // Computes the mean of the combined data and extrapolation in a
390  // given range, use data where they are available and the function
391  // where data are not available
392  // ERROR from DATA ONLY is returned in this version!
393  //
394  Printf("AliPWGHistoTools::GetMeanDataAndExtrapolation: WARNING from data only");
395  Float_t minData = GetLowestNotEmptyBinEdge (hData);
396  Int_t minDataBin = GetLowestNotEmptyBin (hData);
397  Float_t maxData = GetHighestNotEmptyBinEdge(hData);
398  Int_t maxDataBin = GetHighestNotEmptyBin (hData);
399  Double_t integral = 0;
400  mean = 0;
401  error = 0;
402  if (min < minData) {
403  // Compute integral exploiting root function to calculate moments, "unnormalizing" them
404  mean += fExtrapolation->Mean(min,minData)*fExtrapolation->Integral(min,minData);
405  integral += fExtrapolation->Integral(min,minData);
406  cout << " Low "<< mean << " " << integral << endl;
407 
408  }
409 
410  if (max > maxData) {
411  // Compute integral exploiting root function to calculate moments, "unnormalizing" them
412  mean += fExtrapolation->Mean(maxData,max)*fExtrapolation->Integral(maxData,max);
413  integral += fExtrapolation->Integral(maxData,max);
414  cout << " Hi "<< mean << " " << integral << endl;
415  }
416  Float_t err2 = 0;
417 
418  for(Int_t ibin = minDataBin; ibin <= maxDataBin; ibin++){
419  if(hData->GetBinCenter(ibin) < min) continue;
420  if(hData->GetBinCenter(ibin) > max) continue;
421  if(hData->GetBinContent(ibin)) {
422  mean = mean + (hData->GetBinCenter(ibin) * hData->GetBinWidth(ibin)* hData->GetBinContent(ibin));
423  err2 = err2 + TMath::Power(hData->GetBinError(ibin) * hData->GetBinCenter(ibin) * hData->GetBinWidth(ibin),2);
424  integral = integral + hData->GetBinContent(ibin) * hData->GetBinWidth(ibin);
425  }
426  else {
427  Double_t locMin = hData->GetBinLowEdge(ibin);
428  Double_t locMax = hData->GetBinLowEdge(ibin) + hData->GetBinWidth(ibin);
429  mean += fExtrapolation->Mean(locMin,locMax)*fExtrapolation->Integral(locMin,locMax);
430  Double_t deltaIntegral = fExtrapolation->Integral(locMin,locMax);
431  integral += deltaIntegral;
432  cout << "WARNING: bin " << ibin << " is empty, patching with function (+" << deltaIntegral<<")" << endl;
433  }
434  }
435  cout << " Data "<< mean << " " << integral << endl;
436 
437  mean = mean/integral;
438  error = TMath::Sqrt(err2)/integral;
439 
440 
441 }
442 
443 TH1F * AliPWGHistoTools::CombineHistos(const TH1 * h1, TH1 * h2, const TH1* htemplate, Float_t renorm1){
444  // Combine two histos. This assumes the histos have the same binning
445  // in the overlapping region. It computes the arithmetic mean in the
446  // overlapping region and assigns as an error the relative error
447  // h2. TO BE IMPROVED
448 
449  TH1F * hcomb = (TH1F*) htemplate->Clone(TString(h1->GetName())+"_"+h2->GetName());
450 
451  TH1F * h1local = (TH1F*) h1->Clone();
452  h1local->Scale(renorm1);
453 
454  Int_t nBinComb = hcomb->GetNbinsX();
455  for(Int_t ibin = 1; ibin <= nBinComb; ibin++){
456  Int_t ibin1 = h1local->FindBin(hcomb->GetBinCenter(ibin));
457  Int_t ibin2 = h2->FindBin(hcomb->GetBinCenter(ibin));
458 
459  if (!h1local->GetBinContent(ibin1) && !h2->GetBinContent(ibin2) ) {
460  // None has data: go to next bin
461  hcomb->SetBinContent(ibin,0);
462  hcomb->SetBinError (ibin,0);
463  } else if(h1local->GetBinContent(ibin1) && !h2->GetBinContent(ibin2)) {
464  // take data from h1local:
465  hcomb->SetBinContent(ibin,h1local->GetBinContent(ibin1));
466  hcomb->SetBinError (ibin,h1local->GetBinError(ibin1));
467  } else if(!h1local->GetBinContent(ibin1) && h2->GetBinContent(ibin2)) {
468  // take data from h2:
469  hcomb->SetBinContent(ibin,h2->GetBinContent(ibin2));
470  hcomb->SetBinError (ibin,h2->GetBinError(ibin2));
471  } else {
472  hcomb->SetBinContent(ibin,(h1local->GetBinContent(ibin1) +h2->GetBinContent(ibin2))/2);
473  // hcomb->SetBinError (ibin,h1local->GetBinError(ibin1)/h1local->GetBinContent(ibin1)*hcomb->GetBinContent(ibin));
474  hcomb->SetBinError (ibin,h2->GetBinError(ibin2)/h2->GetBinContent(ibin2)*hcomb->GetBinContent(ibin));
475  }
476 
477 
478  }
479 
480 
481  hcomb->SetMarkerStyle(h1local->GetMarkerStyle());
482  hcomb->SetMarkerColor(h1local->GetMarkerColor());
483  hcomb->SetLineColor (h1local->GetLineColor());
484 
485  hcomb->SetXTitle(h1local->GetXaxis()->GetTitle());
486  hcomb->SetYTitle(h1local->GetYaxis()->GetTitle());
487  delete h1local;
488  return hcomb;
489 
490 }
491 
492 
493 void AliPWGHistoTools::GetFromHistoGraphDifferentX(const TH1F * h, TF1 * f, TGraphErrors ** gBarycentre, TGraphErrors ** gXlw) {
494 
495  // Computes the baycentre in each bin and the xlw as defined in NIMA
496  // 355 - 541 using f. Returs 2 graphs with the same y content of h
497  // but with a different x (barycentre and xlw)
498 
499  Int_t nbin = h->GetNbinsX();
500 
501  (*gBarycentre) = new TGraphErrors();
502  (*gXlw) = new TGraphErrors();
503 
504  Int_t ipoint = 0;
505  for(Int_t ibin = 1; ibin <= nbin; ibin++){
506  Float_t min = h->GetBinLowEdge(ibin);
507  Float_t max = h->GetBinLowEdge(ibin+1);
508  Double_t xerr = 0;
509  Double_t xbar = f->Mean(min,max);
510  // compute xLW
511  Double_t temp = 1./(max-min) * f->Integral(min,max);
512  Double_t epsilon = 0.000000001;
513  Double_t increment = 0.0000000001;
514  Double_t xLW = min;
515 
516  while ((f->Eval(xLW)- temp) > epsilon) {
517  xLW += increment;
518  if(xLW > max) {
519  Printf("Cannot find xLW");
520  break;
521  }
522  }
523 
524  if (h->GetBinContent(ibin)!=0) {
525  (*gBarycentre)->SetPoint (ipoint, xbar, h->GetBinContent(ibin));
526  (*gBarycentre)->SetPointError(ipoint, xerr, h->GetBinError(ibin));
527  (*gXlw) ->SetPoint (ipoint, xLW, h->GetBinContent(ibin));
528  (*gXlw) ->SetPointError(ipoint, xerr, h->GetBinError(ibin));
529  ipoint++;
530  }
531  }
532 
533  (*gBarycentre)->SetMarkerStyle(h->GetMarkerStyle());
534  (*gBarycentre)->SetMarkerColor(h->GetMarkerColor());
535  (*gBarycentre)->SetLineColor(h->GetLineColor());
536 
537  (*gBarycentre)->SetTitle(h->GetTitle());
538  (*gBarycentre)->SetName(TString("g_")+h->GetName());
539 
540  (*gXlw)->SetMarkerStyle(h->GetMarkerStyle());
541  (*gXlw)->SetMarkerColor(h->GetMarkerColor());
542  (*gXlw)->SetLineColor(h->GetLineColor());
543  (*gXlw)->SetTitle(h->GetTitle());
544  (*gXlw)->SetName(TString("g_")+h->GetName());
545 
546 
547 }
548 
549 
551 
552  // Get the mean of histo in a range; root is not reliable in sub
553  // ranges with variable binning.
554  Int_t minBin = h->FindBin(min);
555  Int_t maxBin = h->FindBin(max-0.00001);
556 
557  Float_t mean = 0 ;
558  Float_t integral = 0;
559  Float_t err2 = 0;
560  for(Int_t ibin = minBin; ibin <= maxBin; ibin++){
561  mean = mean + (h->GetBinCenter(ibin) * h->GetBinWidth(ibin)* h->GetBinContent(ibin));
562  err2 = err2 + TMath::Power(h->GetBinError(ibin) * h->GetBinCenter(ibin) * h->GetBinWidth(ibin),2);
563  integral = integral + h->GetBinContent(ibin) * h->GetBinWidth(ibin);
564  }
565 
566  Float_t value = mean/integral;
567  if (error) (*error) = TMath::Sqrt(err2);
568  return value;
569 
570 
571 }
572 
573 void AliPWGHistoTools::GetMean(TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max, Int_t normPar) {
574 
575  // Get the mean of function in a range; If normPar is >= 0, it means
576  // that the function is defined such that par[normPar] is its
577  // integral. In this case the error on meanpt can be calculated
578  // correctly. Otherwise, the function is normalized in get moment,
579  // but the error is not computed correctly.
580 
581  return GetMoment("fMean", TString("x*")+func->GetExpFormula(), func, mean, error, min, max, normPar);
582 
583 }
584 
585 void AliPWGHistoTools::GetMeanSquare(TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max, Int_t normPar) {
586 
587  // Get the mean2 of function in a range; If normPar is >= 0, it means
588  // that the function is defined such that par[normPar] is its
589  // integral. In this case the error on meanpt can be calculated
590  // correctly. Otherwise, the function is normalized in get moment,
591  // but the error is not computed correctly.
592 
593  return GetMoment("fMean2", TString("x*x*")+func->GetExpFormula(), func, mean, error, min, max, normPar);
594 
595 
596 }
597 
598 void AliPWGHistoTools::GetMoment(TString name, TString var, TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max, Int_t normPar) {
599 
600  // returns the integral of a function derived from func by prefixing some operation.
601  // the derived function MUST have the same parameter in the same order
602  // Used as a base method for mean and mean2
603  // If normPar is >= 0, it means that the function is defined such
604  // that par[normPar] is its integral. In this case the error on
605  // meanpt can be calculated correctly. Otherwise, the function is
606  // normalized using its numerical integral, but the error is not computed
607  // correctly.
608 
609  // TODO:
610  // - improve to propagate error also in the case you need the
611  // numerical integrals (will be rather slow)
612  // - this version assumes that func is defined using a
613  // TFormula. Generalize to the case of a C++ function
614 
615  if (normPar<0) Printf("AliPWGHistoTools::GetMoment: Warning: If normalization is required, the error may bot be correct");
616  if (!strcmp(func->GetExpFormula(),"")) Printf("AliPWGHistoTools::GetMoment: Warning: Empty formula in the base function");
617  Int_t npar = func->GetNpar();
618 
619  // Definition changes according to the value of normPar
620  TF1 * f = normPar < 0 ?
621  new TF1(name, var) : // not normalized
622  new TF1(name, var+Form("/[%d]",normPar)); // normalized with par normPar
623 
624  // integr is used to normalize if no parameter is provided
625  Double_t integr = normPar < 0 ? func->Integral(min,max) : 1;
626 
627  // The parameter of the function used to compute the mean should be
628  // the same as the parent function: fixed if needed and they should
629  // also have the same errors.
630 
631  // cout << "npar :" << npar << endl;
632 
633  for(Int_t ipar = 0; ipar < npar; ipar++){
634  Double_t parmin, parmax;
635  Double_t value = func->GetParameter(ipar);
636  f->SetParameter(ipar,value);
637  func->GetParLimits(ipar, parmin, parmax);
638  if ( parmin == parmax ) {
639  // if ( parmin || (parmin == 1 && !value) ) { // not sure we I check parmin == 1 here.
640  if ( parmin || (TMath::Abs(parmin-1)<0.000001 && !value) ) { // not sure we I check parmin == 1 here. Changed like this because of coding conventions. Does it still work? FIXME
641  f->FixParameter(ipar,func->GetParameter(ipar));
642  // cout << "Fixing " << ipar << "("<<value<<","<<parmin<<","<<parmax<<")"<<endl;
643  }
644  else {
645  f->SetParError (ipar,func->GetParError(ipar) );
646  // cout << "Setting Err" << ipar << "("<<func->GetParError(ipar)<<")"<<endl;
647  }
648  }
649  else {
650  f->SetParError (ipar,func->GetParError(ipar) );
651  // cout << "Setting Err" << ipar << "("<<func->GetParError(ipar)<<")"<<endl;
652  }
653  }
654 
655 // f->Print();
656 // cout << "----" << endl;
657 // func->Print();
658 
659  mean = normPar < 0 ? f->Integral (min,max)/integr : f->Integral (min,max);
660  error = normPar < 0 ? f->IntegralError(min,max)/integr : f->IntegralError(min,max);
661 // cout << "Mean " << mean <<"+-"<< error<< endl;
662 // cout << "Integral Error " << error << endl;
663 
664 }
665 
666 
667 
668 Bool_t AliPWGHistoTools::Fit (TH1 * h1, TF1* func, Float_t min, Float_t max) {
669 
670  // Fits h1 with func, preforms several checks on the quality of the
671  // fit and tries to improve it. If the fit is not good enough, it
672  // returs false.
673 
674  Double_t amin; Double_t edm; Double_t errdef; Int_t nvpar; Int_t nparx;
675  TVirtualFitter *fitter;
676  cout << "--- Fitting : " << h1->GetName() << " ["<< h1->GetTitle() <<"] ---"<< endl;
677 
678  h1-> Fit(func,"IME0","",min,max);
679  Int_t fitResult = h1-> Fit(func,"IME0","",min,max);
680 // h1-> Fit(func,"0","",min,max);
681 // Int_t fitResult = h1-> Fit(func,"0IE","",min,max);
682 
683 
684 // From TH1:
685 // The fitStatus is 0 if the fit is OK (i.e no error occurred). The
686 // value of the fit status code is negative in case of an error not
687 // connected with the minimization procedure, for example when a wrong
688 // function is used. Otherwise the return value is the one returned
689 // from the minimization procedure. When TMinuit (default case) or
690 // Minuit2 are used as minimizer the status returned is : fitStatus =
691 // migradResult + 10*minosResult + 100*hesseResult +
692 // 1000*improveResult. TMinuit will return 0 (for migrad, minos,
693 // hesse or improve) in case of success and 4 in case of error (see
694 // the documentation of TMinuit::mnexcm). So for example, for an error
695 // only in Minos but not in Migrad a fitStatus of 40 will be returned.
696 // Minuit2 will return also 0 in case of success and different values
697 // in migrad minos or hesse depending on the error. See in this case
698 // the documentation of Minuit2Minimizer::Minimize for the
699 // migradResult, Minuit2Minimizer::GetMinosError for the minosResult
700 // and Minuit2Minimizer::Hesse for the hesseResult. If other
701 // minimizers are used see their specific documentation for the status
702 // code returned. For example in the case of Fumili, for the status
703 // returned see TFumili::Minimize.
704 
705 
706  if( gMinuit->fLimset ) {
707  Printf("ERROR: AliPWGHistoTools: Parameters at limits");
708  return kFALSE;
709  }
710 
711 
713  fitter = TVirtualFitter::GetFitter();
714  Int_t fitStat = fitter->GetStats(amin, edm, errdef, nvpar, nparx);
715 
716  if( ( (fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ")|| (edm > 1e6) || (fitResult !=0 && fitResult < 4000) ) &&
717  TString(gMinuit->fCstatu) != "SUCCESSFUL" &&
718  TString(gMinuit->fCstatu) != "CONVERGED " ) {
719  if(fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ") {
720  Printf("WARNING: AliPWGHistoTools: Cannot properly compute errors");
721  if (fitStat == 0) Printf(" not calculated at all");
722  if (fitStat == 1) Printf(" approximation only, not accurate");
723  if (fitStat == 2) Printf(" full matrix, but forced positive-definite");
724  }
725 
726  if (edm > 1e6) {
727 
728  Printf("WARNING: AliPWGHistoTools: Huge EDM (%f): Fit probably failed!", edm);
729  }
730  if (fitResult != 0) {
731  Printf("WARNING: AliPWGHistoTools: Fit Result (%d)", fitResult);
732  }
733 
734  Printf ("AliPWGHistoTools: Trying Again with Strategy = 2");
735 
736  gMinuit->Command("SET STRATEGY 2"); // more effort
737  fitResult = 0;
738  fitResult = h1-> Fit(func,"0","",min,max);
739  fitResult = h1-> Fit(func,"IME0","",min,max);
740  fitResult = h1-> Fit(func,"IME0","",min,max);
741 
742  fitter = TVirtualFitter::GetFitter();
743 
744  fitStat = fitter->GetStats(amin, edm, errdef, nvpar, nparx);
745 
746  if(fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ") {
747  Printf("ERROR: AliPWGHistoTools: Cannot properly compute errors");
748  if (fitStat == 0) Printf(" not calculated at all");
749  if (fitStat == 1) Printf(" approximation only, not accurate");
750  if (fitStat == 2) Printf(" full matrix, but forced positive-definite");
751  cout << "[" <<gMinuit->fCstatu<<"]" << endl;
752  return kFALSE;
753  }
754 
755  if (edm > 1e6) {
756  Printf("ERROR: AliPWGHistoTools: Huge EDM (%f): Fit probably failed!", edm);
757 
758  return kFALSE;
759  }
760  if (fitResult != 0) {
761  Printf("ERROR: AliPWGHistoTools: Fit Result (%d)", fitResult);
762 
763  return kFALSE;
764  }
765 
766  gMinuit->Command("SET STRATEGY 1"); // back to normal value
767 
768  }
769 
770  cout << "---- FIT OK ----" << endl;
771 
772  return kTRUE;
773 
774 }
775 
777 
778  // Return the index of the lowest non empty bin in the histo h
779 
780  Int_t nbin = h->GetNbinsX();
781  for(Int_t ibin = 1; ibin <= nbin; ibin++){
782  if(h->GetBinContent(ibin)>0) return ibin;
783  }
784 
785  return -1;
786 
787 }
788 
790 
791  // Return the index of the highest non empty bin in the histo h
792 
793  Int_t nbin = h->GetNbinsX();
794  for(Int_t ibin = nbin; ibin > 0; ibin--){
795  if(h->GetBinContent(ibin)>0) return ibin;
796  }
797 
798  return -1;
799 
800 }
801 
802 void AliPWGHistoTools::GetResiduals(const TGraphErrors * gdata, const TF1 * func, TH1F ** hres, TGraphErrors ** gres) {
803 
804  // Returns a graph of residuals vs point and the res/err distribution
805 
806  Int_t npoint = gdata->GetN();
807 
808  (*gres) =new TGraphErrors();
809  (*hres) = new TH1F(TString("hres_")+gdata->GetName()+"-"+func->GetName(),
810  TString("hres_")+gdata->GetName()+"-"+func->GetName(),
811  20,-5,5);
812 
813 
814  for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
815  Float_t x = gdata->GetX()[ipoint];
816  Float_t res = (gdata->GetY()[ipoint] - func->Eval(x))/func->Eval(x);
817  Float_t err = gdata->GetEY()[ipoint]/func->Eval(x);
818  (*hres)->Fill(res/err);
819  (*gres)->SetPoint(ipoint, x, res/err);
820  // (*gres)->SetPointError(ipoint, 0, err);
821 
822  }
823 
824  (*gres)->SetMarkerStyle(gdata->GetMarkerStyle());
825  (*gres)->SetMarkerColor(gdata->GetMarkerColor());
826  (*gres)->SetLineColor (gdata->GetLineColor());
827  (*gres)->GetHistogram()->GetYaxis()->SetTitle("(data-function)/function");
828  (*hres)->SetMarkerStyle(gdata->GetMarkerStyle());
829  (*hres)->SetMarkerColor(gdata->GetMarkerColor());
830  (*hres)->SetLineColor (gdata->GetLineColor());
831 
832 
833 
834 }
835 
836 void AliPWGHistoTools::GetResiduals(const TH1F* hdata, const TF1 * func, TH1F ** hres, TH1F ** hresVsBin) {
837 
838  // Returns an histo of residuals bin by bin and the res/err distribution
839 
840  if (!func) {
841  Printf("AliPWGHistoTools::GetResiduals: No function provided");
842  return;
843  }
844  if (!hdata) {
845  Printf("AliPWGHistoTools::GetResiduals: No data provided");
846  return;
847  }
848 
849  (*hresVsBin) = (TH1F*) hdata->Clone(TString("hres_")+hdata->GetName());
850  (*hresVsBin)->Reset();
851  (*hres) = new TH1F(TString("hres_")+hdata->GetName()+"-"+func->GetName(),
852  TString("hres_")+hdata->GetName()+"-"+func->GetName(),
853  20,-5,5);
854 
855  Int_t nbin = hdata->GetNbinsX();
856  for(Int_t ibin = 1; ibin <= nbin; ibin++){
857  if(!hdata->GetBinContent(ibin)) continue;
858  Float_t res = (hdata->GetBinContent(ibin) - func->Eval(hdata->GetBinCenter(ibin)) ) /
859  func->Eval(hdata->GetBinCenter(ibin));
860  Float_t err = hdata->GetBinError (ibin) / func->Eval(hdata->GetBinCenter(ibin));
861  (*hresVsBin)->SetBinContent(ibin,res);
862  (*hresVsBin)->SetBinError (ibin,err);
863  (*hres)->Fill(res/err);
864 
865  }
866 
867  (*hresVsBin)->SetMarkerStyle(hdata->GetMarkerStyle());
868  (*hresVsBin)->SetMarkerColor(hdata->GetMarkerColor());
869  (*hresVsBin)->SetLineColor (hdata->GetLineColor() );
870  (*hresVsBin)->GetYaxis()->SetTitle("(data-function)/function");
871  (*hres)->SetMarkerStyle(hdata->GetMarkerStyle());
872  (*hres)->SetMarkerColor(hdata->GetMarkerColor());
873  (*hres)->SetLineColor (hdata->GetLineColor() );
874 
875 }
876 
877 void AliPWGHistoTools::GetYield(TH1* h, TF1 * f, Double_t &yield, Double_t &yieldError, Float_t min, Float_t max,
878  Double_t *partialYields, Double_t *partialYieldsErrors){
879 
880  // Returns the yield extracted from the data in the histo where
881  // there are points and from the fit to extrapolate, in the given
882  // range.
883 
884  // Partial yields are also returned if the corresponding pointers are non null
885 
886  Int_t bin1 = h->FindBin(min);
887  Int_t bin2 = h->FindBin(max);
888  Float_t bin1Edge = GetLowestNotEmptyBinEdge (h);
889  Float_t bin2Edge = GetHighestNotEmptyBinEdge(h);
890  Int_t lowestNE = GetLowestNotEmptyBin (h);
891  Int_t highestNE = GetHighestNotEmptyBin(h);
892 
893  Double_t integralFromHistoError = 0;
894  Double_t integralFromHisto = 0;
895  // Double_t integralFromHisto = DoIntegral(h,bin1,bin2,-1,-1,-1,-1,integralFromHistoError,"width",1);
896 
897  for(Int_t ibin = TMath::Max(bin1,lowestNE); ibin <= TMath::Min(bin2,highestNE); ibin++){
898  if(h->GetBinContent(ibin)) {
899  integralFromHisto = integralFromHisto + h->GetBinContent(ibin) * h->GetBinWidth(ibin);
900  integralFromHistoError = integralFromHistoError + h->GetBinError(ibin)*h->GetBinError(ibin) * h->GetBinWidth(ibin);
901  }
902  else {
903  Double_t locMin = h->GetBinLowEdge(ibin);
904  Double_t locMax = h->GetBinLowEdge(ibin) + h->GetBinWidth(ibin);
905  Double_t deltaIntegral = f->Integral(locMin,locMax);
906  integralFromHisto += deltaIntegral;
907  integralFromHistoError = integralFromHistoError + f->IntegralError(locMin,locMax)*f->IntegralError(locMin,locMax);
908  cout << "WARNING: bin " << ibin << " is empty, patching with function (+" << deltaIntegral<<")" << endl;
909  }
910  }
911 
912  integralFromHistoError = TMath::Sqrt(integralFromHistoError);// ok, this is a bit dumb. Just in case the code above is changed again.
913 
914 
915  Double_t integralBelow = min < bin1Edge ? f->Integral(min,bin1Edge) : 0;
916  Double_t integralBelowError = min < bin1Edge ? f->IntegralError(min,bin1Edge) : 0;
917  Double_t integralAbove = max > bin2Edge ? f->Integral(bin2Edge,max) : 0;
918  Double_t integralAboveError = max > bin2Edge ? f->IntegralError(bin2Edge,max) : 0;
919 
920 // cout << "GetYield INFO" << endl;
921 // cout << " " << bin1Edge << " " << bin2Edge << endl;
922 // cout << " " << integralFromHisto << " " << integralBelow << " " << integralAbove << endl;
923 // cout << " " << integralFromHistoError << " " << integralBelowError << " " << integralAboveError << endl;
924 
925  if(partialYields) {
926  partialYields[0] = integralFromHisto;
927  partialYields[1] = integralBelow;
928  partialYields[2] = integralAbove;
929  }
930  if(partialYieldsErrors) {
931  partialYieldsErrors[0] = integralFromHistoError;
932  partialYieldsErrors[1] = integralBelowError;
933  partialYieldsErrors[2] = integralAboveError;
934  }
935  yield = integralFromHisto+integralBelow+integralAbove;
936  yieldError = TMath::Sqrt(integralFromHistoError*integralFromHistoError+
937  integralBelowError*integralBelowError+
938  integralAboveError*integralAboveError);
939 
940 }
941 
943 
944  // Divides g/f. If invert == true => f/g
945 
946  TGraphErrors * gRatio = new TGraphErrors();
947  Int_t npoint = g->GetN();
948  for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
949  Double_t x = g->GetX()[ipoint];
950  Double_t ratio = invert ? f->Eval(x)/g->GetY()[ipoint] :g->GetY()[ipoint]/f->Eval(x);
951  gRatio->SetPoint (ipoint, x, ratio);
952  if(f->Eval(x) && strcmp(g->ClassName(),"TGraphAsymmErrors")) gRatio->SetPointError(ipoint, 0, g->GetEY()[ipoint]/f->Eval(x));
953  // cout << x << " " << g->GetY()[ipoint] << " " << f->Eval(x) << endl;
954 
955  }
956  gRatio->SetMarkerStyle(20);
957  //gRatio->Print();
958  return gRatio;
959 
960 }
961 
962 TGraphErrors * AliPWGHistoTools::Divide2Graphs(const TGraph * g1, const TGraph * g2, Int_t strategy){
963 
964  // Divides g1/g2,
965  // 2 strategies are possible:
966  // - strategy = 0 looks for point with very close centers (also propagates error in this case)
967  // - strategy = 1 Loop over all points from g1 (xi,yi), interpolates the TGraphs with a spline and computes the ratio as yi / SplineG2(xi)
968 
969  TGraphErrors * gRatio = new TGraphErrors();
970 
971  if(strategy == 0) {
972  Int_t ipoint=0;
973  Int_t npoint1 = g1->GetN();
974  Int_t npoint2 = g2->GetN();
975  for(Int_t ipoint1 = 0; ipoint1 < npoint1; ipoint1++){
976  Double_t x1 = g1->GetX()[ipoint1];
977  for(Int_t ipoint2 = 0; ipoint2 < npoint2; ipoint2++){
978  Double_t x2 = g2->GetX()[ipoint2];
979  if((TMath::Abs(x1-x2)/(x1+x2)*2)<0.01) {
980  Double_t ratio = g2->GetY()[ipoint2] ? g1->GetY()[ipoint1]/g2->GetY()[ipoint2] : 0;
981  Double_t eratio = 0;
982  if(g1->InheritsFrom("TGraphErrors") && g2->InheritsFrom("TGraphErrors")) {
983  eratio = g2->GetY()[ipoint2] ?
984  TMath::Sqrt(g1->GetEY()[ipoint1]*g1->GetEY()[ipoint1]/g1->GetY()[ipoint1]/g1->GetY()[ipoint1] +
985  g2->GetEY()[ipoint2]/g2->GetY()[ipoint2]/g2->GetY()[ipoint2] ) * ratio
986  : 0;
987  }
988  gRatio->SetPoint (ipoint, x1, ratio);
989  gRatio->SetPointError(ipoint, 0, eratio);
990  ipoint++;
991  cout << ipoint << " [" << x1 << "] " << g1->GetY()[ipoint1] << "/" << g2->GetY()[ipoint2] << " = " << ratio <<"+-"<<eratio<< endl;
992 
993  // cout << x << " " << g->GetY()[ipoint] << " " << f->Eval(x) << endl;
994  }
995 
996  }
997  }
998  }
999  else if (strategy == 1) {
1000  TSpline3 * splG2 = new TSpline3("fGraphSplG2",g2);
1001  Int_t npoint1 = g1->GetN();
1002  for(Int_t ipoint1 = 0; ipoint1 < npoint1; ipoint1++){
1003  Double_t x1 = g1->GetX()[ipoint1];
1004  Double_t y1 = g1->GetY()[ipoint1];
1005  Double_t ratio = y1/splG2->Eval(x1);
1006  gRatio->SetPoint(ipoint1, x1, ratio);
1007 
1008  }
1009  delete splG2;
1010  }
1011  else Printf("AliPWGHistoTools::Divide2Graphs(): Invalid strategy [%d]", strategy);
1012 
1013  gRatio->SetMarkerStyle(20);
1014  //gRatio->Print();
1015  return gRatio;
1016 
1017 }
1019 
1020  // Sums g1/g2, looks for point with very close centers
1021  Int_t ipoint=0;
1022  TGraphErrors * gSum = new TGraphErrors();
1023  Int_t npoint1 = g1->GetN();
1024  Int_t npoint2 = g2->GetN();
1025  for(Int_t ipoint1 = 0; ipoint1 < npoint1; ipoint1++){
1026  Double_t x1 = g1->GetX()[ipoint1];
1027  for(Int_t ipoint2 = 0; ipoint2 < npoint2; ipoint2++){
1028  Double_t x2 = g2->GetX()[ipoint2];
1029  if((TMath::Abs(x1-x2)/(x1+x2)*2)<0.01) {
1030  Double_t sum = g1->GetY()[ipoint1]+g2->GetY()[ipoint2];
1031  Double_t esum = 0;
1032  if(g1->InheritsFrom("TGraphErrors") && g2->InheritsFrom("TGraphErrors")) {
1033  esum =
1034  TMath::Sqrt(g1->GetEY()[ipoint1]*g1->GetEY()[ipoint1] +
1035  g2->GetEY()[ipoint2]*g2->GetEY()[ipoint2] ) ;
1036 
1037  }
1038  gSum->SetPoint (ipoint, x1, sum);
1039  gSum->SetPointError(ipoint, 0, esum);
1040  ipoint++;
1041  cout << ipoint << " [" << x1 << "] " << g1->GetY()[ipoint1] << "+" << g2->GetY()[ipoint2] << " = " << sum <<"+-"<<esum<< endl;
1042 
1043  // cout << x << " " << g->GetY()[ipoint] << " " << f->Eval(x) << endl;
1044  }
1045 
1046  }
1047  }
1048  gSum->SetMarkerStyle(20);
1049  //gSum->Print();
1050  return gSum;
1051 
1052 }
1053 
1055 
1056  // Divides g/h. If invert == true => h/g
1057 
1058  Bool_t skipError = kFALSE;
1059  if(!strcmp(g->ClassName(),"TGraph")) skipError = kTRUE;
1060  if(!strcmp(g->ClassName(),"TGraphAsymmErrors")) skipError = kTRUE;
1061  if(skipError) {
1062  Printf("WARNING: Skipping graph errors");
1063  }
1064  TGraphErrors * gRatio = new TGraphErrors();
1065  Int_t npoint = g->GetN();
1066  for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
1067  Double_t xj = g->GetX()[ipoint];
1068  Double_t yj = g->GetY()[ipoint];
1069  Double_t yje = skipError ? 0 : g->GetEY()[ipoint];
1070 
1071  Int_t binData = h->FindBin(xj);
1072  Double_t yd = h->GetBinContent(binData);
1073  Double_t yde = h->GetBinError(binData);
1074  Double_t xd = h->GetBinCenter(binData);
1075 
1076 
1077 
1078  if (!yd) continue;
1079 
1080  if (TMath::Abs(xd-xj)/TMath::Abs(xd) > 0.01){
1081  Printf( "WARNING: bin center (%f) and x graph (%f) are more than 1 %% away, skipping",xd,xj );
1082  continue;
1083 
1084  }
1085 
1086  Double_t ratio = invert ? yd/yj : yj/yd;
1087 
1088  gRatio->SetPoint(ipoint, xj, ratio);
1089  gRatio->SetPointError(ipoint, 0, TMath::Sqrt(yde*yde/yd/yd + yje*yje/yj/yj)*ratio);
1090  // gRatio->SetPointError(ipoint, 0, yje/yj * ratio);
1091  }
1092 
1093  return gRatio;
1094 
1095 
1096 }
1097 
1098 TH1F * AliPWGHistoTools::DivideHistoByFunc(TH1F * h, TF1 * f, Bool_t invert){
1099 
1100  // Divides h/f. If invert == true => f/g
1101  // Performs the integral of f on the bin range to perform the ratio
1102  // Returns a histo with the same binnig as h
1103 
1104  // Prepare histo for ratio
1105  TH1F * hRatio = (TH1F*) h->Clone(TString("hRatio_")+h->GetName()+"_"+f->GetName());
1106  hRatio->Reset();
1107  // Set y title
1108  if(!invert) hRatio->SetYTitle(TString(h->GetName())+"/"+f->GetName());
1109  else hRatio->SetYTitle(TString(f->GetName())+"/"+h->GetName());
1110 
1111  // Loop over all bins
1112  Int_t nbin = hRatio->GetNbinsX();
1113 
1114  for(Int_t ibin = 1; ibin <= nbin; ibin++){
1115  Double_t yhisto = h->GetBinContent(ibin);
1116  Double_t yerror = h->GetBinError(ibin);
1117  Double_t xmin = h->GetBinLowEdge(ibin);
1118  Double_t xmax = h->GetBinLowEdge(ibin+1);
1119  Double_t yfunc = f->Integral(xmin,xmax)/(xmax-xmin);
1120  Double_t ratio = yfunc&&yhisto ? (invert ? yfunc/yhisto : yhisto/yfunc) : 0 ;
1121  Double_t error = yfunc ? yerror/yfunc : 0 ;
1122  hRatio->SetBinContent(ibin,ratio);
1123  hRatio->SetBinError (ibin,error);
1124  }
1125 
1126  return hRatio;
1127 
1128 }
1129 
1130 void AliPWGHistoTools::WeightedMean(Int_t npoints, const Double_t *x, const Double_t *xerr, Double_t &mean, Double_t &meanerr){
1131 
1132  // Performs the weighted mean of npoints numbers in x with errors in xerr
1133 
1134  mean = 0;
1135  meanerr = 0;
1136 
1137  Double_t sumweight = 0;
1138 
1139  for (Int_t ipoint = 0; ipoint < npoints; ipoint++){
1140 
1141  Double_t xerr2 = xerr[ipoint]*xerr[ipoint];
1142  if(xerr2>0){
1143  // cout << "xe2 " << xerr2 << endl;
1144  Double_t weight = 1. / xerr2;
1145  sumweight += weight;
1146  mean += weight * x[ipoint];
1147  }// else cout << " Skipping " << ipoint << endl;
1148 
1149  }
1150 
1151 
1152  if(sumweight){
1153  mean /= sumweight;
1154  meanerr = TMath::Sqrt(1./ sumweight);
1155  }
1156  else {
1157  // cout << " No sumweight" << endl;
1158  mean = 0;
1159  meanerr = 0;
1160  }
1161 
1162 
1163 }
1164 
1166  // Returns an histogram with the same binning as h, filled with the relative error bin by bin
1167  TH1 * hrel = (TH1*) h->Clone(TString(h->GetName())+"_rel");
1168  hrel->Reset();
1169  Int_t nbin = hrel->GetNbinsX();
1170  for(Int_t ibin = 1; ibin <= nbin; ibin++){
1171  if(h->GetBinError(ibin)){
1172  hrel->SetBinContent(ibin,h->GetBinError(ibin)/h->GetBinContent(ibin));
1173  } else {
1174  hrel->SetBinContent(ibin,0);
1175  }
1176  hrel->SetBinError(ibin,0);
1177  }
1178 
1179  return hrel;
1180 }
1181 
1182 
1183 void AliPWGHistoTools::GetValueAndError(TH1 * hdest, const TH1 * hvalue, const TH1 * herror, Bool_t isPercentError) {
1184 
1185  // Put into source, bin-by-bin, the values from hvalue and the
1186  // errors from content from herror.
1187  // Used mainly to combine histos of systemati errors with their spectra
1188  // Set isPercentError to kTRUE if the error is given in %
1189 
1190  if(hdest == NULL){
1191  Printf("AliPWGHistoTools::GetValueAndError Errror: hdest is null");
1192  return;
1193  }
1194 
1195 
1196  Int_t nbin = hdest->GetNbinsX();
1197  Int_t nBinSourceVal = hvalue->GetNbinsX();
1198  Int_t nBinSourceErr = herror->GetNbinsX();
1199 
1200  for(Int_t iBinDest = 1; iBinDest <= nbin; iBinDest++){
1201  Float_t lowPtDest=hdest->GetBinLowEdge(iBinDest);
1202  Float_t binWidDest=hdest->GetBinWidth(iBinDest);
1203  // Loop over Source bins and find overlapping bins to Dest
1204  // First value then error
1205  // Value
1206  Bool_t foundValue = kFALSE;
1207  for(Int_t iBinSourceVal=1; iBinSourceVal<=nBinSourceVal; iBinSourceVal++){
1208  Float_t lowPtSource= hvalue->GetBinLowEdge(iBinSourceVal) ;
1209  Float_t binWidSource= hvalue->GetBinWidth(iBinSourceVal);
1210  if(TMath::Abs(lowPtDest-lowPtSource)<0.001 && TMath::Abs(binWidSource-binWidDest)<0.001){
1211  Double_t content = hvalue->GetBinContent(iBinSourceVal);
1212  hdest->SetBinContent(iBinDest, content);
1213  foundValue = kTRUE;
1214  break;
1215  }
1216  }
1217  // if (!foundValue){
1218  // Printf("AliPWGHistoTools::GetValueAndError: Error: cannot find matching value source bin for destination %d",iBinDest);
1219  // }
1220 
1221  // Error
1222  Bool_t foundError = kFALSE;
1223  for(Int_t iBinSourceErr=1; iBinSourceErr<=nBinSourceErr; iBinSourceErr++){
1224  Float_t lowPtSource= herror->GetBinLowEdge(iBinSourceErr) ;
1225  Float_t binWidSource= herror->GetBinWidth(iBinSourceErr);
1226  if(TMath::Abs(lowPtDest-lowPtSource)<0.001 && TMath::Abs(binWidSource-binWidDest)<0.001){
1227  Double_t error = herror->GetBinContent(iBinSourceErr);
1228  // cout << "-> " << iBinDest << " " << error << " " << hdest->GetBinContent(iBinDest) << endl;
1229 
1230  hdest->SetBinError(iBinDest, isPercentError ? error * hdest->GetBinContent(iBinDest) : error);
1231  foundError=kTRUE;
1232  break;
1233  }
1234  }
1235  // if (!foundError ){
1236  // Printf("AliPWGHistoTools::GetValueAndError: Error: cannot find matching error source bin for destination %d",iBinDest);
1237  // }
1238  }
1239 
1240 
1241 }
1242 
1243 void AliPWGHistoTools::AddHisto(TH1 * hdest, const TH1* hsource, Bool_t getMirrorBins) {
1244 
1245  // Adds hsource to hdest bin by bin, even if they have a different
1246  // binning If getMirrorBins is true, it takes the negative bins
1247  // (Needed because sometimes the TPC uses the positive axis for
1248  // negative particles and the possitive axis for positive
1249  // particles).
1250 
1251 
1252  if (hdest == NULL) {
1253  Printf("Error: hdest is NULL\n");
1254  return;
1255  }
1256  if (hsource == NULL) {
1257  Printf("Error: hsource is NULL\n");
1258  return;
1259  }
1260 
1261  Int_t nBinSource = hsource->GetNbinsX();
1262  Int_t nBinDest = hdest->GetNbinsX();
1263 
1264  // Loop over destination bins,
1265  for(Int_t iBinDest=1; iBinDest<=nBinDest; iBinDest++){
1266  Float_t lowPtDest=hdest->GetBinLowEdge(iBinDest);
1267  Float_t binWidDest=hdest->GetBinWidth(iBinDest);
1268  // Loop over Source bins and find overlapping bins to Dest
1269  Bool_t found = kFALSE;
1270  for(Int_t iBinSource=1; iBinSource<=nBinSource; iBinSource++){
1271  Float_t lowPtSource= getMirrorBins ? -hsource->GetBinLowEdge(iBinSource)+hsource->GetBinWidth(iBinSource) : hsource->GetBinLowEdge(iBinSource) ;
1272  Float_t binWidSource= hsource->GetBinWidth(iBinSource) ;
1273  if(TMath::Abs(lowPtDest-lowPtSource)<0.001 && TMath::Abs(binWidSource-binWidDest)<0.001){
1274  Float_t dest=hdest->GetBinContent(iBinDest);
1275  Float_t source=hsource->GetBinContent(iBinSource);
1276  Float_t edest=hdest->GetBinError(iBinDest);
1277  Float_t esource=hsource->GetBinError(iBinSource);
1278  Double_t cont=dest+source;
1279  Double_t econt=TMath::Sqrt(edest*edest+esource*esource);
1280  hdest->SetBinContent(iBinDest,cont);
1281  hdest->SetBinError (iBinDest,econt);
1282  found = kTRUE;
1283 
1284  break;
1285  }
1286  }
1287  // if (!found){
1288  // Printf("Error: cannot find matching source bin for destination %d",iBinDest);
1289  // }
1290  }
1291 
1292 
1293 }
1294 
1296 
1297  // Combine the errors of hdest with the errors of h1, summing in
1298  // quadrature. Results are put in hdest. Histograms are assumed to
1299  // have the same binning
1300 
1301  Int_t nbin = hdest->GetNbinsX();
1302  for(Int_t ibin = 1; ibin <= nbin; ibin++){
1303  Double_t e1 = hdest->GetBinError(ibin);
1304  Double_t e2 = h1->GetBinError(ibin);
1305  hdest->SetBinError(ibin, TMath::Sqrt(e1*e1+e2*e2));
1306  }
1307 
1308 
1309 }
1310 
1311 TH1F * AliPWGHistoTools::DivideHistosDifferentBins(const TH1F* h1, const TH1F* h2) {
1312  // Divides 2 histos even if they have a different binning. Finds
1313  // overlapping bins and divides them
1314 
1315  // 1. clone histo
1316  TH1F * hRatio = new TH1F(*h1);
1317  Int_t nBinsH1=h1->GetNbinsX();
1318  Int_t nBinsH2=h2->GetNbinsX();
1319  // Loop over H1 bins,
1320  for(Int_t iBin=1; iBin<=nBinsH1; iBin++){
1321  hRatio->SetBinContent(iBin,0.);
1322  hRatio->SetBinContent(iBin,0.);
1323  Float_t lowPtH1=h1->GetBinLowEdge(iBin);
1324  Float_t binWidH1=h1->GetBinWidth(iBin);
1325  // Loop over H2 bins and find overlapping bins to H1
1326  for(Int_t jBin=1; jBin<=nBinsH2; jBin++){
1327  Float_t lowPtH2=h2->GetBinLowEdge(jBin);
1328  Float_t binWidH2=h2->GetBinWidth(jBin);
1329  if(TMath::Abs(lowPtH1-lowPtH2)<0.001 && TMath::Abs(binWidH2-binWidH1)<0.001){
1330  Float_t numer=h1->GetBinContent(iBin);
1331  Float_t denom=h2->GetBinContent(jBin);
1332  Float_t enumer=h1->GetBinError(iBin);
1333  Float_t edenom=h2->GetBinError(jBin);
1334  Double_t ratio=0.;
1335  Double_t eratio=0.;
1336  if(numer>0. && denom>0.){
1337  ratio=numer/denom;
1338  eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1339  }
1340  hRatio->SetBinContent(iBin,ratio);
1341  hRatio->SetBinError(iBin,eratio);
1342  break;
1343  }
1344  }
1345  }
1346  return hRatio;
1347 }
1348 
1349 Double_t AliPWGHistoTools::DoIntegral(TH1* h, Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Int_t binz1, Int_t binz2, Double_t & error ,
1350  Option_t *option, Bool_t doError)
1351 {
1352  // function to compute integral and optionally the error between the limits
1353  // specified by the bin number values working for all histograms (1D, 2D and 3D)
1354  // copied from TH! to fix a bug still present in 5-27-06b
1355  Int_t nbinsx = h->GetNbinsX();
1356  if (binx1 < 0) binx1 = 0;
1357  if (binx2 > nbinsx+1 || binx2 < binx1) binx2 = nbinsx+1;
1358  if (h->GetDimension() > 1) {
1359  Int_t nbinsy = h->GetNbinsY();
1360  if (biny1 < 0) biny1 = 0;
1361  if (biny2 > nbinsy+1 || biny2 < biny1) biny2 = nbinsy+1;
1362  } else {
1363  biny1 = 0; biny2 = 0;
1364  }
1365  if (h->GetDimension() > 2) {
1366  Int_t nbinsz = h->GetNbinsZ();
1367  if (binz1 < 0) binz1 = 0;
1368  if (binz2 > nbinsz+1 || binz2 < binz1) binz2 = nbinsz+1;
1369  } else {
1370  binz1 = 0; binz2 = 0;
1371  }
1372 
1373  // - Loop on bins in specified range
1374  TString opt = option;
1375  opt.ToLower();
1376  Bool_t width = kFALSE;
1377  if (opt.Contains("width")) width = kTRUE;
1378 
1379 
1380  Double_t dx = 1.;
1381  Double_t dy = 1.;
1382  Double_t dz = 1.;
1383  Double_t integral = 0;
1384  Double_t igerr2 = 0;
1385  for (Int_t binx = binx1; binx <= binx2; ++binx) {
1386  if (width) dx = h->GetXaxis()->GetBinWidth(binx);
1387  for (Int_t biny = biny1; biny <= biny2; ++biny) {
1388  if (width) dy = h->GetYaxis()->GetBinWidth(biny);
1389  for (Int_t binz = binz1; binz <= binz2; ++binz) {
1390  if (width) dz = h->GetZaxis()->GetBinWidth(binz);
1391  Int_t bin = h->GetBin(binx, biny, binz);
1392  if (width) integral += h->GetBinContent(bin)*dx*dy*dz;
1393  else integral += h->GetBinContent(bin);
1394  if (doError) {
1395  if (width) igerr2 += h->GetBinError(bin)*h->GetBinError(bin)*dx*dy*dz*dx*dy*dz;
1396  else igerr2 += h->GetBinError(bin)*h->GetBinError(bin);
1397  }
1398  // cout << h->GetBinContent(bin) << " " << h->GetBinError(bin) << " " << dx*dy*dz << " " << integral << " +- " << igerr2 << endl;
1399 
1400  }
1401  }
1402  }
1403 
1404  if (doError) error = TMath::Sqrt(igerr2);
1405  return integral;
1406 }
1407 
1409 
1410  // Computes the dmt/dptdeta function using the dN/dpt function
1411  // This is a protected function used internally by GetdMtdy to integrate dN/dpt function using mt as a weight
1412  // The mass of the particle is given as p[0]
1413  Double_t pt = x[0];
1414  Double_t mass = p[0];
1415  Double_t mt = TMath::Sqrt(pt*pt + mass*mass);
1416  Double_t jacobian = pt/mt;
1417  if(!fdNdptForETCalc){
1418  Printf("AliPWGHistoTools::dMtdptFunction: ERROR: fdNdptForETCalc not defined");
1419  return 0;
1420  }
1421  Double_t dNdpt = fdNdptForETCalc->Eval(pt);
1422  return dNdpt*mt*jacobian; // FIXME: do I have to normalize somehow?
1423 
1424 }
1425 
1426 Double_t AliPWGHistoTools::GetdMtdEta(TH1 *hData, TF1 * fExtrapolation, Double_t mass) {
1427  // Computes dMtdEta integrating dN/dptdy with the proper weights and jacobian.
1428  Printf("WARNING ALIBWTOOLS::GetdMtdEta: ONLY USING FUNCTION FOR THE TIME BEING, hData");
1429  if(!hData) {
1430  Printf("hData not set");
1431  }
1432 
1433  // Assign the fiunction used internally by dMtdptFunction
1434  fdNdptForETCalc = fExtrapolation;
1435  // Create the function to be integrated
1436  TF1 * funcdMtdPt = new TF1 ("funcdMtdPt", dMtdptFunction, 0.0, 20, 1);
1437  funcdMtdPt->SetParameter(0,mass);
1438  // Integrate it
1439  Double_t dMtdEta = funcdMtdPt->Integral(0,100);
1440  // Clean up
1441  fdNdptForETCalc=0;
1442  delete funcdMtdPt;
1443  //return
1444  return dMtdEta;
1445 
1446 }
1447 
1448 
1449 
1450 
1452 
1453  // Scale a grh by the factor scale. If g1 is a TGraph(Asym)Error, errors are also scaled
1454 
1455  if(!g1) return;
1456  Bool_t isErrors = TString(g1->ClassName()) == "TGraphErrors";
1457  Bool_t isAErrors = TString(g1->ClassName()) == "TGraphAsymmErrors";
1458 
1459  Int_t npoint = g1->GetN();
1460  for (Int_t ipoint = 0; ipoint<npoint; ipoint++) {
1461  Double_t x = g1->GetX()[ipoint];
1462  Double_t y = g1->GetY()[ipoint];
1463  Double_t val = y*scale;
1464  g1->SetPoint(ipoint, x, val);
1465  // if it is aclass with errors, also scale those
1466  if (isAErrors) {
1467  Double_t errxlow = g1->GetEXlow() [ipoint];
1468  Double_t errxhigh = g1->GetEXhigh()[ipoint];
1469  Double_t errylow = g1->GetEYlow() [ipoint]*scale;
1470  Double_t erryhigh = g1->GetEYhigh()[ipoint]*scale;
1471  ((TGraphAsymmErrors*)g1)->SetPointError(ipoint, errxlow, errxhigh, errylow, erryhigh);
1472  }
1473  else if (isErrors) {
1474  Double_t erry = g1->GetEY()[ipoint]*scale;
1475  ((TGraphErrors*)g1)->SetPointError(ipoint, g1->GetEX()[ipoint], erry);
1476  }
1477  }
1478 
1479 
1480 }
static void GetValueAndError(TH1 *hdest, const TH1 *hvalue, const TH1 *herror, Bool_t isPercentError)
static TH1 * GetdNdPtFromOneOverPt(const TH1 *h1Pt)
static TGraphErrors * ConcatenateGraphs(const TGraphErrors *g1, const TGraphErrors *g2)
static void GetMeanDataAndExtrapolation(const TH1 *hData, TF1 *fExtrapolation, Double_t &mean, Double_t &error, Float_t min=0, Float_t max=100)
double Double_t
Definition: External.C:58
static Int_t GetLowestNotEmptyBin(const TH1 *h)
static TH1 * GetRelativeError(TH1 *h)
Int_t nbinsy
static Bool_t Fit(TH1 *h, TF1 *f, Float_t min, Float_t max)
Double_t mass
static TGraphErrors * Divide2Graphs(const TGraph *g1, const TGraph *g2, Int_t strategy)
static Int_t GetHighestNotEmptyBin(const TH1 *h)
static Double_t DoIntegral(TH1 *h, Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Int_t binz1, Int_t binz2, Double_t &error, Option_t *option, Bool_t doError)
static TH1 * GetdNdmtFromdNdpt(const TH1 *hpt, Double_t mass)
static void GetMoment(TString name, TString var, TF1 *func, Float_t &mean, Float_t &error, Float_t min, Float_t max, Int_t normPar=-1)
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Definition: External.C:212
static void GetMeanSquare(TF1 *func, Float_t &mean, Float_t &error, Float_t min=0, Float_t max=100, Int_t normPar=-1)
static void AddHisto(TH1 *hdest, const TH1 *hsource, Bool_t getMirrorBins=kFALSE)
static TH1F * DivideHistosDifferentBins(const TH1F *h1, const TH1F *h2)
static TH1F * CombineHistos(const TH1 *h1, TH1 *h2, const TH1 *htemplate, Float_t renorm1=1.)
TList * fitter
Definition: DrawAnaELoss.C:26
static void GetResiduals(const TGraphErrors *gdata, const TF1 *func, TH1F **hres, TGraphErrors **gres)
static TH1 * GetOneOverPtdNdPt(const TH1 *hPt)
static TGraphErrors * GetGraphFromHisto(const TH1F *h, Bool_t binWidth=kTRUE)
static void GetFromHistoGraphDifferentX(const TH1F *h, TF1 *f, TGraphErrors **gBarycentre, TGraphErrors **gXlw)
static TH1F * GetHistoFromGraph(const TGraphErrors *g, const TH1F *hTemplate)
static Double_t dMtdptFunction(Double_t *x, Double_t *p)
Int_t nbinsx
static TGraphErrors * DivideGraphByFunc(const TGraphErrors *g, const TF1 *f, Bool_t invert=kFALSE)
static void WeightedMean(Int_t npoints, const Double_t *x, const Double_t *xerr, Double_t &mean, Double_t &meanerr)
static TH1F * DivideHistoByFunc(TH1F *h, TF1 *f, Bool_t invert=kFALSE)
static void GetYield(TH1 *h, TF1 *f, Double_t &yield, Double_t &yieldError, Float_t min=0, Float_t max=100, Double_t *partialYields=0, Double_t *partialYieldsErrors=0)
static TGraphErrors * Add2Graphs(const TGraphErrors *g1, const TGraphErrors *g2)
static TH1F * Combine3HistosWithErrors(const TH1 *h1, const TH1 *h2, const TH1 *h3, TH1 *he1, TH1 *he2, TH1 *he3, const TH1 *htemplate, Int_t statFrom=0, Float_t renorm1=1., Float_t renorm2=1., Float_t renorm3=1., TH1 **hSyst=0, Bool_t errorFromBinContent=kFALSE)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
const char Option_t
Definition: External.C:48
static void ScaleGraph(TGraph *g1, Double_t scale)
const Int_t nbins
bool Bool_t
Definition: External.C:53
static TGraphErrors * DivideGraphByHisto(const TGraphErrors *g, TH1 *h, Bool_t invert=kFALSE)
static void GetHistoCombinedErrors(TH1 *hdest, const TH1 *h1)
static TH1 * GetdNdptFromdNdmt(const TH1 *hmt, Double_t mass)
Definition: External.C:196
static Double_t GetdMtdEta(TH1 *hData, TF1 *fExtrapolation, Double_t mass)
static Float_t GetMean(TH1F *h, Float_t min, Float_t max, Float_t *error=NULL)