coverity 15108 fixed
[u/mrichter/AliRoot.git] / FMD / flow / TestBessel.C
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 //____________________________________________________________________
9 double 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 //____________________________________________________________________
17 double 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 //____________________________________________________________________
25 double _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 //____________________________________________________________________
38 double rootI(double* xp, double* pp)
39 {
40   return _rootI(pp[0], xp[0]);
41 }
42
43 //____________________________________________________________________
44 const 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 //____________________________________________________________________
52 const 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 //____________________________________________________________________
60 TCanvas* c = 0;
61 TCanvas* d = 0;
62
63 //____________________________________________________________________
64 void 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 //____________________________________________________________________
102 void 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 //____________________________________________________________________
143 void
144 TestBessel()
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__
176 int 
177 main()
178 {
179   TApplication app("app", 0, 0);
180   testBessel();
181   app.Run();
182   return 0;
183 }
184 #endif