bug fix
[u/mrichter/AliRoot.git] / FMD / flow / TestBessel.C
CommitLineData
bbcd8619 1#include <TF1.h>
2#include <TCanvas.h>
3#include <TApplication.h>
4#include <TMath.h>
5#include <TH2.h>
6#include <FMD/flow/AliFMDFlowBessel.h>
7
8//____________________________________________________________________
9double myI(double* xp, double* pp)
10{
11 double n = pp[0];
12 double x = xp[0];
13 return AliFMDFlowBessel::I(n, x);
14}
15
16//____________________________________________________________________
17double myDI(double* xp, double* pp)
18{
19 double n = pp[0];
20 double x = xp[0];
21 return AliFMDFlowBessel::DiffI(n, x);
22}
23
24//____________________________________________________________________
25double _rootI(double n, double x)
26{
27 if (fabs(x) < 0.00001) return 0;
28 if (n == 0.5)
29 return TMath::Sqrt(2 * x / TMath::Pi()) * TMath::SinH(x)/x;
30 if (n == 1.5)
31 return TMath::Sqrt(2 * x / TMath::Pi()) *
32 (TMath::CosH(x) / x - TMath::SinH(x) / (x * x));
33 if (n == 2.5) return _rootI(0.5, x) - 3 * _rootI(1.5, x) / x;
34 return TMath::BesselI(int(fabs(n)), x);
35}
36
37//____________________________________________________________________
38double rootI(double* xp, double* pp)
39{
40 return _rootI(pp[0], xp[0]);
41}
42
43//____________________________________________________________________
44const char* ntonm(double n)
45{
46 if (int(2 * fabs(n)) % 2 == 1)
47 return Form("I%d_2", int(2 * n));
48 return Form("I%d", int(n));
49}
50
51//____________________________________________________________________
52const char* ntostr(double n)
53{
54 if (int(2 * fabs(n)) % 2 == 1)
55 return Form("I_{%d/2}", int(2 * n));
56 return Form("I_{%d}", int(n));
57}
58
59//____________________________________________________________________
60TCanvas* c = 0;
61TCanvas* d = 0;
62
63//____________________________________________________________________
64void compare(double n, int col)
65{
66 Double_t min = 0.00001;
67 Double_t max = 4;
68 if (!c) {
69 c = new TCanvas("c_I", "I_{#nu}");
70 c->SetFillColor(0);
71 c->cd();
72 TH2* f = new TH2F("f_I", "I_{#nu}", 100, min, max, 100, -1, 13);
73 f->SetXTitle("x");
74 f->SetYTitle("I_{#nu}");
75 f->SetStats(0);
76 f->Draw();
77 }
78 TString str = ntostr(n);
79 TString nm = ntonm(n);
80
81 c->cd();
82 TF1* r = new TF1(Form("root_%s", nm.Data()), rootI, min, max, 1);
83 r->SetTitle(Form("%s (ROOT)", str.Data()));
84 r->SetParameter(0, n);
85 r->SetLineColor(col);
86 r->Draw("same");
87
88 c->cd();
89 TF1* m = new TF1(Form("my_%s", nm.Data()), myI, min, max, 1);
90 m->SetTitle(Form("%s (My)", str.Data()));
91 m->SetParameter(0, n);
92 m->SetLineColor(col);
93 m->SetLineStyle(2);
94 m->Draw("same");
95
96 c->Modified();
97 c->Update();
98 c->cd();
99
100}
101//____________________________________________________________________
102void dcompare(double n, int col)
103{
104 Double_t min = 0;
105 Double_t max = 4;
106 if (!d) {
107 d = new TCanvas("c_dI", "dI_{#nu}/dx");
108 d->SetFillColor(0);
109 d->cd();
110 TH2* f = new TH2F("f_dI", "dI_{#nu}/dx", 100, min, max, 100, -1, 12);
111 f->SetXTitle("x");
112 f->SetYTitle("dI_{#nu}/dx");
113 f->SetStats(0);
114 f->Draw();
115 }
116
117 TString str = ntostr(n);
118 TString nm = ntonm(n);
119
120 d->cd();
121 TF1* r = new TF1(Form("root_d%s", nm.Data()), rootI, min, max, 1);
122 r->SetTitle(Form("d%s/dx (ROOT)", str.Data()));
123 r->SetParameter(0, n);
124 r->SetLineColor(col);
125 TGraph* g = new TGraph(r, "d");
126 g->SetTitle(r->GetTitle());
127 g->Draw("same c");
128
129 d->cd();
130 TF1* m = new TF1(Form("my_d%s", nm.Data()), myDI, min, max, 1);
131 m->SetTitle(Form("d%s/dx (My)", str.Data()));
132 m->SetParameter(0, n);
133 m->SetLineColor(col);
134 m->SetLineStyle(2);
135 m->Draw("same");
136
137 d->Modified();
138 d->Update();
139 d->cd();
140}
141
142//____________________________________________________________________
143void
144TestBessel()
145{
146 gStyle->SetPalette(1);
147 // gStyle->SetOptTitle(0);
148 // double n[] = { 0, 0.5, 1, 1.5, 2, 2.5, 3 };
149 double n[] = { -1, 0, 0.5, 1, 1.5, 2, 2.5, 3 };
150 // double n[] = { 0.5, 1.5, 2.5 };
151 // double n[] = { 0.5, 1.5, 2.5 };
152 size_t m = sizeof(n)/sizeof(double);
153 Int_t nc = gStyle->GetNumberOfColors();
154 for (size_t i = 0; i < m; i++) {
155 int col = gStyle->GetColorPalette(int(float(i) / m * nc));
156 compare(n[i], col);
157 dcompare(n[i], col);
158 }
159 c->cd();
160 c->SetTopMargin(0.05);
161 c->SetRightMargin(0.05);
162 TLegend* l = c->BuildLegend(.1, .35, .4, .95);
163 l->SetFillColor(0);
164 l->SetBorderSize(1);
165 // c->Print("I.ps");
166 d->cd();
167 d->SetTopMargin(0.05);
168 d->SetRightMargin(0.05);
169 l = d->BuildLegend(.1, .35, .4, .95);
170 l->SetFillColor(0);
171 l->SetBorderSize(1);
172 // d->Print("dI.ps");
173}
174
175#ifndef __CINT__
176int
177main()
178{
179 TApplication app("app", 0, 0);
180 testBessel();
181 app.Run();
182 return 0;
183}
184#endif