Added test script of Bessel functions

diff --git a/FMD/flow/AliFMDFlowBessel.cxx b/FMD/flow/AliFMDFlowBessel.cxx
index 89488ba30e6a3ee3652d82ab9bf7cacc53ff9784..2ebe5e9f76a78fe1d2461fd884dfc5d4ca94fef9 100644 (file)
--- a/FMD/flow/AliFMDFlowBessel.cxx
+++ b/FMD/flow/AliFMDFlowBessel.cxx
@@ -285,6 +285,25 @@ AliFMDFlowBessel::Iwhole(UInt_t in, Double_t x, Double_t* bi, Double_t* di)
   return mn;
 }
 
+//____________________________________________________________________
+namespace 
+{
+  Double_t* EnsureSize(Double_t*& a, UInt_t& old, const UInt_t size)
+  {
+    if (a && old < size) { 
+      // Will delete, and reallocate. 
+      if (old != 1) delete [] a;
+      a = 0;
+    }
+    if (!a && size > 0) {
+      // const UInt_t n = (size <= 1 ? 2 : size);
+      a              = new Double_t[size];
+      old            = size; // n;
+    }
+    return a;
+  }
+}
+      
 //____________________________________________________________________
 UInt_t 
 AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x, 
@@ -296,6 +315,10 @@ AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x,
   UInt_t  in1 = UInt_t(fabs(n1));
   UInt_t  in2 = UInt_t(fabs(n2));
   UInt_t  nt  = UInt_t(n2 - n1 + 1);
+  static Double_t* tbi = 0;
+  static Double_t* tdi = 0;
+  static UInt_t    tn  = 0;
+
   if (Int_t(2 * fabs(n1)) % 2 == 1) { // Half-integer argument 
     if (Int_t(2 * fabs(n2)) % 2 != 1) { 
       std::cout << "If n1 is half-integer (" << n1 << ") then n2 "
@@ -310,16 +333,8 @@ AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x,
     Int_t l2    = (s1 > 0 ? in1 : 0);
 
     // Temporary buffers - only grows in size. 
-    static Double_t* tbi = 0;
-    static Double_t* tdi = 0;
-    static UInt_t    tn  = 0;
-    if (!tbi || tn < nt+1) { 
-      if (tbi) delete [] tbi;
-      if (tdi) delete [] tdi;
-      tn  = nt+1;
-      tbi = new Double_t[tn];
-      tdi = new Double_t[tn];
-    }
+    tbi = EnsureSize(tbi, tn, nt+1);
+    tdi = EnsureSize(tdi, tn, nt+1);
     // Double_t tbi[nt+1], tdi[nt+1];
 
     if (s1 < 0) { 
@@ -348,18 +363,12 @@ AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x,
   }
 
   UInt_t  n   = UInt_t(in1 > in2 ? in1 : in2);
-  static Double_t* tbi = 0;
-  static Double_t* tdi = 0;
-  static UInt_t    tn  = 0;
-  if (!tbi || tn < n+1) { 
-    if (tbi) delete [] tbi;
-    if (tdi) delete [] tdi;
-    tn  = n+1;
-    tbi = new Double_t[tn];
-    tdi = new Double_t[tn];
-  }
+  // Temporary buffers - only grows in size. 
+  tbi = EnsureSize(tbi, tn, n+1);
+  tdi = EnsureSize(tdi, tn, n+1);
   // Double_t  tbi[n+1];
   // Double_t  tdi[n+1];
+
   UInt_t  r   = Iwhole(n, x, tbi, tdi);
   if (r < n) 
     std::cerr << "Only got " << r << "/" << n << " values" 

diff --git a/FMD/flow/TestBessel.C b/FMD/flow/TestBessel.C
new file mode 100644 (file)
index 0000000..5b95061
--- /dev/null
+++ b/FMD/flow/TestBessel.C
@@ -0,0 +1,184 @@
+#include <TF1.h>
+#include <TCanvas.h>
+#include <TApplication.h>
+#include <TMath.h>
+#include <TH2.h>
+#include <FMD/flow/AliFMDFlowBessel.h>
+
+//____________________________________________________________________
+double myI(double* xp, double* pp)
+{
+  double n = pp[0];
+  double x = xp[0];
+  return AliFMDFlowBessel::I(n, x);
+}
+
+//____________________________________________________________________
+double myDI(double* xp, double* pp)
+{
+  double n = pp[0];
+  double x = xp[0];
+  return AliFMDFlowBessel::DiffI(n, x);
+}
+
+//____________________________________________________________________
+double _rootI(double n, double x)
+{
+  if (fabs(x) < 0.00001) return 0;
+  if (n == 0.5) 
+    return TMath::Sqrt(2 * x / TMath::Pi()) * TMath::SinH(x)/x;
+  if (n == 1.5) 
+    return TMath::Sqrt(2 * x / TMath::Pi()) * 
+      (TMath::CosH(x) / x - TMath::SinH(x) / (x * x));
+  if (n == 2.5) return _rootI(0.5, x) - 3 * _rootI(1.5, x) / x;
+  return TMath::BesselI(int(fabs(n)), x);
+}
+
+//____________________________________________________________________
+double rootI(double* xp, double* pp)
+{
+  return _rootI(pp[0], xp[0]);
+}
+
+//____________________________________________________________________
+const char* ntonm(double n)
+{
+  if (int(2 * fabs(n)) % 2 == 1) 
+    return Form("I%d_2", int(2 * n));
+  return Form("I%d", int(n));
+}
+
+//____________________________________________________________________
+const char* ntostr(double n)
+{
+  if (int(2 * fabs(n)) % 2 == 1) 
+    return Form("I_{%d/2}", int(2 * n));
+  return Form("I_{%d}", int(n));
+}
+
+//____________________________________________________________________
+TCanvas* c = 0;
+TCanvas* d = 0;
+
+//____________________________________________________________________
+void compare(double n, int col)
+{
+  Double_t min = 0.00001;
+  Double_t max = 4;
+  if (!c) {
+    c = new TCanvas("c_I", "I_{#nu}");
+    c->SetFillColor(0);
+    c->cd();
+    TH2* f = new TH2F("f_I", "I_{#nu}", 100, min, max, 100, -1, 13);
+    f->SetXTitle("x");
+    f->SetYTitle("I_{#nu}"); 
+    f->SetStats(0);
+    f->Draw();
+  }
+  TString  str = ntostr(n);
+  TString  nm  = ntonm(n);
+
+  c->cd();
+  TF1* r = new TF1(Form("root_%s", nm.Data()), rootI, min, max, 1);
+  r->SetTitle(Form("%s (ROOT)", str.Data()));
+  r->SetParameter(0, n);
+  r->SetLineColor(col);
+  r->Draw("same");
+
+  c->cd();
+  TF1* m = new TF1(Form("my_%s", nm.Data()), myI, min, max, 1);
+  m->SetTitle(Form("%s (My)", str.Data()));
+  m->SetParameter(0, n);
+  m->SetLineColor(col);
+  m->SetLineStyle(2);
+  m->Draw("same");
+
+  c->Modified();
+  c->Update();
+  c->cd();
+
+}
+//____________________________________________________________________
+void dcompare(double n, int col)
+{
+  Double_t min = 0;
+  Double_t max = 4;
+  if (!d) { 
+    d = new TCanvas("c_dI", "dI_{#nu}/dx");
+    d->SetFillColor(0);
+    d->cd();
+    TH2* f = new TH2F("f_dI", "dI_{#nu}/dx", 100, min, max, 100, -1, 12);
+    f->SetXTitle("x");
+    f->SetYTitle("dI_{#nu}/dx");
+    f->SetStats(0);
+    f->Draw();
+  }
+
+  TString  str = ntostr(n);
+  TString  nm  = ntonm(n);
+
+  d->cd();
+  TF1* r = new TF1(Form("root_d%s", nm.Data()), rootI, min, max, 1);
+  r->SetTitle(Form("d%s/dx (ROOT)", str.Data()));
+  r->SetParameter(0, n);
+  r->SetLineColor(col);
+  TGraph* g = new TGraph(r, "d");
+  g->SetTitle(r->GetTitle());
+  g->Draw("same c");
+
+  d->cd();
+  TF1* m = new TF1(Form("my_d%s", nm.Data()), myDI, min, max, 1);
+  m->SetTitle(Form("d%s/dx (My)", str.Data()));
+  m->SetParameter(0, n);
+  m->SetLineColor(col);
+  m->SetLineStyle(2);
+  m->Draw("same");
+
+  d->Modified();
+  d->Update();
+  d->cd();
+}
+  
+//____________________________________________________________________
+void
+TestBessel()
+{
+  gStyle->SetPalette(1);
+  // gStyle->SetOptTitle(0);
+  // double n[] = { 0, 0.5, 1, 1.5, 2, 2.5, 3 };
+  double n[] = { -1, 0, 0.5, 1, 1.5, 2, 2.5, 3 };
+  // double n[] = { 0.5, 1.5, 2.5 };
+  // double n[] = { 0.5, 1.5, 2.5 };
+  size_t  m  = sizeof(n)/sizeof(double);
+  Int_t   nc = gStyle->GetNumberOfColors();
+  for (size_t i = 0; i < m; i++) { 
+    int col = gStyle->GetColorPalette(int(float(i) / m * nc));
+    compare(n[i], col);
+    dcompare(n[i], col);
+  }
+  c->cd();
+  c->SetTopMargin(0.05);
+  c->SetRightMargin(0.05);
+  TLegend* l = c->BuildLegend(.1, .35, .4, .95);
+  l->SetFillColor(0);
+  l->SetBorderSize(1);
+  // c->Print("I.ps");
+  d->cd();
+  d->SetTopMargin(0.05);
+  d->SetRightMargin(0.05);
+  l = d->BuildLegend(.1, .35, .4, .95);
+  l->SetFillColor(0);
+  l->SetBorderSize(1);
+  // d->Print("dI.ps");
+}
+
+#ifndef __CINT__
+int 
+main()
+{
+  TApplication app("app", 0, 0);
+  testBessel();
+  app.Run();
+  return 0;
+}
+#endif