AliPhysics  608b256 (608b256)
PlotUtils.C
Go to the documentation of this file.
1 
12 #if !defined(__CINT__) || defined(__MAKECINT__)
13 
14 #include <TH1D.h>
15 #include <TF1.h>
16 #include <TMath.h>
17 #include <TGraphErrors.h>
18 
19 #endif
20 
24 //-----------------------------------------------------------------------------
26 {
27  Double_t N = par[0];
28  Double_t width = par[1];
29 
30  Double_t mean = par[2];
31  Double_t sigma = par[3];
32  Double_t alpha = par[4];
33  Double_t n = par[5];
34 
35  Double_t A = pow(n/fabs(alpha),n)*exp(-pow(alpha,2)/2);
36  Double_t B = n/fabs(alpha) - fabs(alpha);
37 
38  if ((x[0]-mean)/sigma>-alpha)
39  return N*width*TMath::Gaus(x[0],mean,sigma,1);
40  else
41  return N/(sqrt(TMath::TwoPi())*sigma)*width*A*pow(B-(x[0]-mean)/sigma,-n);
42 }
43 
47 //-----------------------------------------------------------------------------
49 {
50  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
51  (2*par[2]*par[2]) );
52  Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0] + par[6]*x[0]*x[0]*x[0];
53  return gaus+back;
54 }
55 
59 //-----------------------------------------------------------------------------
61 {
62  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
63  (2*par[2]*par[2]) );
64  Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
65  return gaus+back;
66 }
67 
71 //-----------------------------------------------------------------------------
73 {
74  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
75  (2*par[2]*par[2]) );
76  Double_t back = par[3] + par[4]*x[0];
77  return gaus+back;
78 }
79 
83 //-----------------------------------------------------------------------------
85 {
86  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
87  (2*par[2]*par[2]) );
88  Double_t back = 0;//par[3];
89  return gaus+back;
90 }
91 
95 //-----------------------------------------------------------------------------
97 {
98  if ( x[0] > 0.425 && x[0] < 0.65 )
99  {
100  TF1::RejectPoint();
101  return 0;
102  }
103 
104  return par[0] + par[1]*x[0];
105  // + x[0]*x[0]*par[2] + x[0]*x[0]*x[0]*par[3];//+ x[0]*x[0]*x[0]*x[0]*par[4];
106 }
107 
111 //-----------------------------------------------------------------------------
113 {
114  //if ( x[0] > 0.07 && x[0] < 0.2 )
115  if ( x[0] < 0.25 )
116  {
117  TF1::RejectPoint();
118  return 0;
119  }
120 
121  return par[0] + par[1]*x[0];
122  // + x[0]*x[0]*par[2] + x[0]*x[0]*x[0]*par[3] ;//+ x[0]*x[0]*x[0]*x[0]*par[4];
123 }
124 
125 //-----------------------------------------------------------------------------
130 //-----------------------------------------------------------------------------
132 (Double_t num , Double_t den,
133  Double_t numErr, Double_t denErr)
134 {
135  if ( num == 0 || den == 0 ) return 0.;
136 
137  //printf("\t num %e den %e numErr %e denErr %e\n",num,den,numErr,denErr);
138 
139  return num/den * TMath::Sqrt( ( numErr * numErr ) / ( num * num ) +
140  ( denErr * denErr ) / ( den * den ) );
141 }
142 
143 //-----------------------------------------------------------------------------
149 //-----------------------------------------------------------------------------
150 static TGraphErrors * DivideGraphs
152 {
153  if ( !gDen || !gNum )
154  {
155  printf("Graph num %p or graph den %p not available\n",gNum,gDen);
156  return 0x0;
157  }
158  const Int_t nBins = gNum->GetN();
159  if ( nBins != gDen->GetN() )
160  {
161  printf("Cannot divide %s with %d bins and %s with %d bins!\n",
162  gNum->GetName(),nBins,gDen->GetName(),gDen->GetN());
163  return 0x0;
164  }
165 
166  Double_t ratio [nBins];
167  Double_t ratioErr[nBins];
168  Double_t x [nBins];
169  Double_t xErr [nBins];
170  for (Int_t ibin = 0; ibin < nBins; ibin++)
171  {
172  Double_t num = gNum->GetY ()[ibin];
173  Double_t den = gDen->GetY ()[ibin];
174  Double_t numErr = gNum->GetEY()[ibin];
175  Double_t denErr = gDen->GetEY()[ibin];
176 
177  x [ibin] = gNum->GetX ()[ibin];
178  xErr[ibin] = gNum->GetEX()[ibin];
179 
180  if ( num == 0 || den == 0 )
181  {
182  ratio [ibin] = 0;
183  ratioErr[ibin] = 0;
184  continue;
185  }
186 
187  ratio [ibin] = num / den ;
188  ratioErr[ibin] = GetFractionError(num,den,numErr,denErr);
189 
190  // printf("bin %d, x %f (%f) num %f (%f), den %f (%f), ratio %f (%f) \n",
191  // ibin,x[ibin],xErr[ibin],num,numErr,den,denErr,ratio[ibin],ratioErr[ibin]);
192  } // do the ratio to sum
193 
194  return new TGraphErrors(nBins,x,ratio,xErr,ratioErr);
195 }
196 
197 
198 //-----------------------------------------------------------------------------
206 //-----------------------------------------------------------------------------
207 static void GetRangeIntegralAndError
208 (TH1D* h, Int_t binMin, Int_t binMax,
209  Double_t & integral, Double_t & integralErr )
210 {
211  integral = 0;
212  integralErr = 0;
213  for(Int_t ibin = binMin; ibin <= binMax; ibin++)
214  {
215  //if ( h->GetBinContent(ibin) == 0 ) continue ;
216 
217  integral += h->GetBinContent(ibin);
218  integralErr += h->GetBinError (ibin) * h->GetBinError(ibin);
219  }
220 // printf("\t bin range [%d,%d], n %d, sum %2.2e err %2.2e\n",
221 // binMin,binMax,n,integral,integralErr);
222  if ( integralErr > 0 ) integralErr = TMath::Sqrt(integralErr);
223 }
224 
225 //-----------------------------------------------------------------------------
233 //-----------------------------------------------------------------------------
234 static void GetRangeIntegralAndError
235 (TH1D* h, Float_t minX, Float_t maxX,
236  Double_t & integral, Double_t & integralErr )
237 {
238  Int_t binMin = h->FindBin(minX);
239  Int_t binMax = h->FindBin(maxX)-1;
240  return GetRangeIntegralAndError(h, binMin, binMax,
241  integral, integralErr );
242 }
243 
244 //-----------------------------------------------------------------------------
255 //-----------------------------------------------------------------------------
257 (TGraphErrors* graph,
258  Double_t min, Double_t max,
259  Float_t errFraction = 0.5,
260  Float_t value = -1, Float_t valueErr = 0 )
261 {
262  for(Int_t ibin = 0; ibin < graph->GetN(); ibin++)
263  {
264  Float_t x = graph->GetX()[ibin];
265  Float_t y = graph->GetY()[ibin];
266  Float_t yE = graph->GetEY()[ibin];
267  if(x < min || x > max )
268  {
269  graph->SetPoint (ibin,x, value );
270  graph->SetPointError(ibin,x, valueErr);
271  }
272  else
273  if(yE > errFraction*y) // Remove points with more than 50% error
274  {
275  //printf("ibin %d, y %f, yE %f\n",ibin,y,yE);
276  graph->SetPoint (ibin,x,value );
277  graph->SetPointError(ibin,x,valueErr);
278  }
279  }
280 }
281 
282 //-----------------------------------------------------------------------------
284 //-----------------------------------------------------------------------------
285 static void ScaleBinBySize(TH1D* h)
286 {
287  for(Int_t ibin = 1; ibin < h->GetNbinsX();ibin++)
288  {
289  Double_t width = h->GetBinWidth(ibin);
290  Double_t content = h->GetBinContent(ibin);
291  //printf("bin %d, width %f, content %e\n",ibin,width,content);
292  h->SetBinContent(ibin,content/width);
293  }
294 }
295 
296 //-----------------------------------------------------------------------------
303 //-----------------------------------------------------------------------------
304 void GetCanvasColRowNumber(Int_t npad, Int_t & ncol, Int_t & nrow,
305  Bool_t square = kFALSE)
306 {
307  if(npad <= 0)
308  {
309  ncol = 0;
310  nrow = 0;
311  return;
312  }
313 
314  ncol = TMath::Sqrt(npad);
315  nrow = ncol;
316 
317  if ( square )
318  return;
319 
320  if ( npad == 1 )
321  {
322  ncol = 1;
323  nrow = 1;
324  }
325  else if ( npad == 2 )
326  {
327  ncol = 2;
328  nrow = 1;
329  }
330  else if ( npad == 3 )
331  {
332  ncol = 3;
333  nrow = 1;
334  }
335  else if ( npad == 4 )
336  {
337  ncol = 2;
338  nrow = 2;
339  }
340  else if ( npad < 7 )
341  {
342  ncol = 3;
343  nrow = 2;
344  }
345  else if ( npad < 10 )
346  {
347  ncol = 3;
348  nrow = 3;
349  }
350  else if ( npad < 13 )
351  {
352  ncol = 4;
353  nrow = 3;
354  }
355  else if ( npad < 17 )
356  {
357  ncol = 4;
358  nrow = 4;
359  }
360  else if ( npad < 21 )
361  {
362  ncol = 5;
363  nrow = 4;
364  }
365  else if ( npad < 26 )
366  {
367  ncol = 5;
368  nrow = 5;
369  }
370  else if ( npad < 31 )
371  {
372  ncol = 6;
373  nrow = 5;
374  }
375  else if ( npad < 37 )
376  {
377  ncol = 6;
378  nrow = 6;
379  }
380 
381 }
382 
static Double_t FunctionCrystalBall(Double_t *x, Double_t *par)
Definition: PlotUtils.C:25
static Double_t FunctionTruncatedPolPi0(Double_t *x, Double_t *par)
Definition: PlotUtils.C:112
double Double_t
Definition: External.C:58
static void ScaleBinBySize(TH1D *h)
Scale histogram bins by its size.
Definition: PlotUtils.C:285
static Double_t FunctionGaussPol0(Double_t *x, Double_t *par)
Definition: PlotUtils.C:84
static Double_t FunctionGaussPol1(Double_t *x, Double_t *par)
Definition: PlotUtils.C:72
static TGraphErrors * DivideGraphs(TGraphErrors *gNum, TGraphErrors *gDen)
Definition: PlotUtils.C:151
static Double_t GetFractionError(Double_t num, Double_t den, Double_t numErr, Double_t denErr)
Definition: PlotUtils.C:132
void RemovePointsOutOfRangeOrLargeErrorFromGraph(TGraphErrors *graph, Double_t min, Double_t max, Float_t errFraction=0.5, Float_t value=-1, Float_t valueErr=0)
Definition: PlotUtils.C:257
Double_t * sigma
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Definition: External.C:212
static Double_t FunctionGaussPol3(Double_t *x, Double_t *par)
Definition: PlotUtils.C:48
static void GetRangeIntegralAndError(TH1D *h, Int_t binMin, Int_t binMax, Double_t &integral, Double_t &integralErr)
Definition: PlotUtils.C:208
static Double_t FunctionGaussPol2(Double_t *x, Double_t *par)
Definition: PlotUtils.C:60
static Double_t FunctionTruncatedPolEta(Double_t *x, Double_t *par)
Definition: PlotUtils.C:96
void GetCanvasColRowNumber(Int_t npad, Int_t &ncol, Int_t &nrow, Bool_t square=kFALSE)
Definition: PlotUtils.C:304
bool Bool_t
Definition: External.C:53