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