AliPhysics  2c6b7ad (2c6b7ad)
PDFLandau.C
Go to the documentation of this file.
1 #include <TF1.h>
2 #include <TH1.h>
3 #include <TMath.h>
4 #include <cmath>
5 
6 struct PDF
7 {
8 
9  static double cheb(double v, double* p, double* q, int n) {
10  double num = 0;
11  double den = 0;
12  for (int i = n-1; i >= 0; i--) {
13  num = (p[i] + num)*v;
14  den = (q[i] + den)*v;
15  }
16  return num/den;
17  }
18  static double f1(double v) {
19  static double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
20 
21  double u = std::exp(v+1);
22  if (u < 1e-10) return 0;
23 
24  double ue = std::exp(-1/u);
25  double us = std::sqrt(u);
26  return 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
27  }
28  static double f2(double v) {
29  static double p1[] = {0.4259894875, -0.1249762550,
30  0.03984243700, -0.006298287635,
31  0.001511162253};
32  static double q1[] = {1.0, -0.3388260629,
33  0.09594393323, -0.01608042283,
34  0.003778942063};
35 
36  double u = std::exp(-v-1);
37  return std::exp(-u)*std::sqrt(u)*cheb(v,p1,q1,5);
38  }
39  static double f3(double v) {
40  static double p2[] = {0.1788541609, 0.1173957403,
41  0.01488850518, -0.001394989411,
42  0.0001283617211};
43  static double q2[] = {1.0, 0.7428795082,
44  0.3153932961, 0.06694219548,
45  0.008790609714};
46  return cheb(v,p2,q2,5);
47  }
48  static double f4(double v) {
49  static double p3[] = {0.1788544503, 0.09359161662,
50  0.006325387654, 0.00006611667319,
51  -0.000002031049101};
52  static double q3[] = {1.0, 0.6097809921,
53  0.2560616665, 0.04746722384,
54  0.006957301675};
55  return cheb(v,p3,q3,5);
56  }
57  static double f5(double v) {
58  static double p4[] = {0.9874054407, 118.6723273,
59  849.2794360, -743.7792444,
60  427.0262186};
61  static double q4[] = {1.0, 106.8615961,
62  337.6496214, 2016.712389,
63  1597.063511};
64  double u = 1/v;
65  return u*u*cheb(u,p4,q4,5);
66  }
67  static double f6(double v) {
68  static double p5[] = {1.003675074, 167.5702434,
69  4789.711289, 21217.86767,
70  -22324.94910};
71  static double q5[] = {1.0, 156.9424537,
72  3745.310488, 9834.698876,
73  66924.28357};
74  double u = 1/v;
75  return u*u*cheb(u,p5,q5,5);
76  }
77  static double f7(double v) {
78  static double p6[] = {1.000827619, 664.9143136,
79  62972.92665, 475554.6998,
80  -5743609.109};
81  static double q6[] = {1.0, 651.4101098,
82  56974.73333, 165917.4725,
83  -2815759.939};
84  double u = 1/v;
85  return u*u*cheb(u,p6,q6,5);
86  }
87  static double f8(double v) {
88  static double a2[2] = {-1.845568670,-4.284640743};
89  double u = 1/(v-v*std::log(v)/(v+1));
90  return u*u*(1+(a2[0]+a2[1]*u)*u);
91  }
92 
93  static Double_t wrap(Double_t* px, Double_t* pp) {
94  Int_t i = pp[0];
95  Double_t v = px[0];
96  switch (i) {
97  case 1: return f1(v);
98  case 2: return f2(v);
99  case 3: return f3(v);
100  case 4: return f4(v);
101  case 5: return f5(v);
102  case 6: return f6(v);
103  case 7: return f7(v);
104  case 8: return f8(v);
105  }
106  return 0;
107  }
108 
109  static TF1* getComp(Int_t i) {
110  Double_t l = -10;
111  Double_t h = -5.5;
112  switch (i) {
113  case 2: l = -5.5; h = -1; break;
114  case 3: l = -1; h = 1; break;
115  case 4: l = 1; h = 5; break;
116  case 5: l = 5; h = 12; break;
117  case 6: l = 12; h = 50; break;
118  case 7: l = 50; h = 300; break;
119  case 8: l = 300; h = 1000; break;
120  default: i = 1; break;
121  }
122 
123  l *= (l < 0 ? 1.5 : 0.4);
124  h *= (h < 0 ? 0.5 : 1.5);
125 
126  TF1* f = new TF1(Form("f%d", i), &PDF::wrap, l, h, 1);
127  f->SetParameter(0, i);
128  f->SetLineColor(i);
129  return f;
130  }
131  static void draw() {
132  TH1* h = new TH1F("h","h", 1000, -10, 1000);
133  h->Draw();
134 
135  TF1* l = new TF1("l", "landau", -10, 1000);
136  l->SetParameters(1,0,1);
137  l->SetLineStyle(2);
138  l->Draw("same");
139 
140  Double_t max = 0;
141  for (Int_t i = 1; i <= 8; i++) {
142  TF1* f = getComp(i);
143  f->Draw("same");
144  max = TMath::Max(f->GetMaximum(),max);
145  }
146  h->SetMaximum(1.1*max);
147  }
148 };
149 
150 
151 
double Double_t
Definition: External.C:58
static double f7(double v)
Definition: PDFLandau.C:77
static double f3(double v)
Definition: PDFLandau.C:39
Definition: PDFLandau.C:6
static Double_t wrap(Double_t *px, Double_t *pp)
Definition: PDFLandau.C:93
int Int_t
Definition: External.C:63
static void draw()
Definition: PDFLandau.C:131
static double f8(double v)
Definition: PDFLandau.C:87
static double f1(double v)
Definition: PDFLandau.C:18
static double f2(double v)
Definition: PDFLandau.C:28
static double f5(double v)
Definition: PDFLandau.C:57
static double f6(double v)
Definition: PDFLandau.C:67
static double cheb(double v, double *p, double *q, int n)
Definition: PDFLandau.C:9
static double f4(double v)
Definition: PDFLandau.C:48
static TF1 * getComp(Int_t i)
Definition: PDFLandau.C:109
Definition: External.C:196