[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / corrs / PDFLandau.C

#include <TF1.h>
#include <TH1.h>
#include <TMath.h>
#include <cmath>

struct PDF 
{
  
  static double cheb(double v, double* p, double* q, int n) {
    double num = 0;
    double den = 0;
    for (int i = n-1; i >= 0; i--) {
      num = (p[i] + num)*v;
      den = (q[i] + den)*v;
    }
    return num/den;
  }
  static double f1(double v) { 
    static double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};

    double u = std::exp(v+1);
    if (u < 1e-10) return 0;
    
    double ue = std::exp(-1/u);
    double us = std::sqrt(u);
    return 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
  }
  static double f2(double v) {
    static double p1[] = {0.4259894875,  -0.1249762550, 
			  0.03984243700, -0.006298287635,   
			  0.001511162253};
    static double q1[] = {1.0,           -0.3388260629, 
			  0.09594393323, -0.01608042283,    
			  0.003778942063};

    double u = std::exp(-v-1);
    return std::exp(-u)*std::sqrt(u)*cheb(v,p1,q1,5);
  }
  static double f3(double v) { 
    static double p2[] = {0.1788541609,   0.1173957403, 
			  0.01488850518, -0.001394989411,   
			  0.0001283617211};
    static double q2[] = {1.0,            0.7428795082, 
			  0.3153932961,   0.06694219548,    
			  0.008790609714};
    return cheb(v,p2,q2,5);
  }
  static double f4(double v) { 
    static double p3[] = {0.1788544503,   0.09359161662,
			  0.006325387654, 0.00006611667319,
			  -0.000002031049101};
    static double q3[] = {1.0,            0.6097809921, 
			  0.2560616665,   0.04746722384,    
			  0.006957301675};
    return cheb(v,p3,q3,5);
  }
  static double f5(double v) { 
    static double p4[] = {0.9874054407, 118.6723273,  
			  849.2794360, -743.7792444,      
			  427.0262186};
    static double q4[] = {1.0,           106.8615961,  
			  337.6496214,   2016.712389,      
			  1597.063511};
    double u   = 1/v;
    return u*u*cheb(u,p4,q4,5);
  }
  static double f6(double v) { 
    static double p5[] = {1.003675074,   167.5702434,  
			   4789.711289, 21217.86767,   
			   -22324.94910};
    static double q5[] = {1.0,          156.9424537,
			   3745.310488,  9834.698876,
			   66924.28357};
    double u   = 1/v;
    return u*u*cheb(u,p5,q5,5);
  }
  static double f7(double v) { 
    static double p6[] = {1.000827619,  664.9143136,  
			  62972.92665,  475554.6998,
			  -5743609.109};
    static double q6[] = {1.0,          651.4101098,  
			  56974.73333,  165917.4725,
			  -2815759.939};
    double u   = 1/v;
    return u*u*cheb(u,p6,q6,5);
  }
  static double f8(double v) { 
    static double a2[2] = {-1.845568670,-4.284640743};
    double u = 1/(v-v*std::log(v)/(v+1));
    return u*u*(1+(a2[0]+a2[1]*u)*u);
  }
    
  static Double_t wrap(Double_t* px, Double_t* pp) { 
    Int_t    i = pp[0];
    Double_t v = px[0];
    switch (i) { 
    case 1: return f1(v);
    case 2: return f2(v);
    case 3: return f3(v);
    case 4: return f4(v);
    case 5: return f5(v);
    case 6: return f6(v);
    case 7: return f7(v);
    case 8: return f8(v);
    }
    return 0;
  }

  static TF1* getComp(Int_t i) { 
    Double_t l = -10;
    Double_t h = -5.5;
    switch (i) { 
    case 2: l =   -5.5; h =   -1; break;
    case 3: l =   -1;   h =    1; break;
    case 4: l =    1;   h =    5; break;
    case 5: l =    5;   h =   12; break;
    case 6: l =   12;   h =   50; break;
    case 7: l =   50;   h =  300; break;
    case 8: l =  300;   h = 1000; break;
    default: i = 1; break;
    }
    
    l *= (l < 0 ? 1.5 : 0.4);
    h *= (h < 0 ? 0.5 : 1.5);

    TF1* f = new TF1(Form("f%d", i), &PDF::wrap, l, h, 1);
    f->SetParameter(0, i);
    f->SetLineColor(i);
    return f;
  }
  static void draw() { 
    TH1* h = new TH1F("h","h", 1000, -10, 1000);
    h->Draw();

    TF1* l = new TF1("l", "landau", -10, 1000);
    l->SetParameters(1,0,1);
    l->SetLineStyle(2);
    l->Draw("same");

    Double_t max = 0;
    for (Int_t i = 1; i <= 8; i++) { 
      TF1* f = getComp(i);
      f->Draw("same");
      max = TMath::Max(f->GetMaximum(),max);
    }
    h->SetMaximum(1.1*max);
  }
};
Commit	Line	Data
67a218da	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 uucheb(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 uucheb(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 uucheb(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 uu(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