Added code to do flow analysis.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 18 Sep 2007 14:51:23 +0000 (14:51 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 18 Sep 2007 14:51:23 +0000 (14:51 +0000)
The core of the code is in the class AliFMDFlowBin.   It contains an
AliFMDFlowHarmonic object to calculate the harmonic, 3
AliFMDFlowEventPlane objects to calculate the full event event-plane,
and two sub-event event-planes, as well as three AliFMDFlowResolution
objects for calculating the event-plane resolution in 3 different
ways.

The classes AliFMDFlowBinned1D and AliFMDFlowBinned2D are "histograms"
of AliFMDFlowBin objects.

AliFMDFlowAxis is an axis in an AliFMDFlowBinned1D or AliFMDFlowBinned2D
"histogram".

AliFMDFlowTrueBin and AliFMDFlowTrue1D are for calculating the flow
from a known event plane angle - e.g., when reading data from an
event generator.

The class AliFMDFlowStat is single variable statistics that does not
suffer from rounding errors when adding small numbers to a large sum.

AliFMDFlowBessel is  namespace for a number of Bessel functions.  These
allow one to efficiently calculate the value of several Modified Bessel
functions of the 1st kind in one go.

AliFMDFlowUtil contains some utility functions used through out the code.

Example scripts and programs will follow later.

24 files changed:
FMD/flow/AliFMDFlowAxis.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowAxis.h [new file with mode: 0644]
FMD/flow/AliFMDFlowBessel.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowBessel.h [new file with mode: 0644]
FMD/flow/AliFMDFlowBin.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowBin.h [new file with mode: 0644]
FMD/flow/AliFMDFlowBinned1D.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowBinned1D.h [new file with mode: 0644]
FMD/flow/AliFMDFlowBinned2D.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowBinned2D.h [new file with mode: 0644]
FMD/flow/AliFMDFlowEfficiency.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowEfficiency.h [new file with mode: 0644]
FMD/flow/AliFMDFlowEventPlane.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowEventPlane.h [new file with mode: 0644]
FMD/flow/AliFMDFlowHarmonic.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowHarmonic.h [new file with mode: 0644]
FMD/flow/AliFMDFlowResolution.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowResolution.h [new file with mode: 0644]
FMD/flow/AliFMDFlowStat.h [new file with mode: 0644]
FMD/flow/AliFMDFlowTrue.cxx [new file with mode: 0644]
FMD/flow/AliFMDFlowTrue.h [new file with mode: 0644]
FMD/flow/AliFMDFlowUtil.h [new file with mode: 0644]
FMD/flow/FMDflowLinkDef.h [new file with mode: 0644]
FMD/libFMDflow.pkg [new file with mode: 0644]

diff --git a/FMD/flow/AliFMDFlowAxis.cxx b/FMD/flow/AliFMDFlowAxis.cxx
new file mode 100644 (file)
index 0000000..51e1a4f
--- /dev/null
@@ -0,0 +1,96 @@
+/** @file 
+    @brief Implementation of an Axis in a Flow "histogram" */
+#include "flow/AliFMDFlowAxis.h"
+#include <cmath>
+
+//====================================================================
+AliFMDFlowAxis::AliFMDFlowAxis(UShort_t n, Double_t* bins) 
+  : fN(n), 
+    fBins(0)
+{
+  fBins = new Double_t[fN+1];
+  for (UInt_t i = 0; i <= fN; i++) fBins[i] = bins[i];
+}
+//____________________________________________________________________
+AliFMDFlowAxis::AliFMDFlowAxis(UShort_t n, Double_t min, Double_t max)
+  : fN(n), 
+    fBins(0)
+{
+  fBins     = new Double_t[fN+1];
+  Double_t dx = (max-min)/fN;
+  for (UInt_t i = 0; i <= fN; i++) fBins[i] = min + i * dx;
+}
+//____________________________________________________________________
+AliFMDFlowAxis::AliFMDFlowAxis(const AliFMDFlowAxis& other)
+  : fN(other.fN), 
+    fBins(0)
+{
+  fBins = new Double_t[fN+1];
+  for (UInt_t i = 0; i <= fN; i++) fBins[i] = other.fBins[i];
+}
+//____________________________________________________________________
+AliFMDFlowAxis& 
+AliFMDFlowAxis::operator=(const AliFMDFlowAxis& other)
+{
+  if (fBins) delete [] fBins;
+  fN    = other.fN;
+  fBins = new Double_t[fN+1];
+  for (UInt_t i = 0; i <= fN; i++) fBins[i] = other.fBins[i];
+  return *this;
+}
+//____________________________________________________________________
+AliFMDFlowAxis::~AliFMDFlowAxis()
+{
+  if (fBins) delete [] fBins;
+  fBins = 0;
+}
+//____________________________________________________________________
+Int_t 
+AliFMDFlowAxis::FindBin(Double_t x) const 
+{
+  if (x < fBins[0])  return -1;
+  if (x > fBins[fN]) return -1;
+  UInt_t above, below, middle;
+  above = fN+2;
+  below = 0;
+  while((above - below) > 1) {
+    middle = (above + below) / 2;
+    if (x == fBins[middle-1]) return middle-1;
+    if (x  < fBins[middle-1]) above = middle;
+    else                      below = middle;
+  }
+  return below-1;
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowAxis::BinWidth(UShort_t i) const
+{
+  if (i >= fN) return 0;
+  return (fBins[i+1] - fBins[i]);
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowAxis::BinCenter(UShort_t i) const
+{
+  if (i >= fN) return 0;
+  return fBins[i] + BinWidth(i) / 2;
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowAxis::BinLower(UShort_t i) const
+{
+  if (i >= fN) return 0;
+  return fBins[i];
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowAxis::BinUpper(UShort_t i) const
+{
+  if (i >= fN) return 0;
+  return fBins[i+1];
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowAxis.h b/FMD/flow/AliFMDFlowAxis.h
new file mode 100644 (file)
index 0000000..81a5f6f
--- /dev/null
@@ -0,0 +1,73 @@
+// -*- mode: C++ -*-
+/** @file 
+    @brief Declaration of an Axis in a Flow "histogram" */
+#ifndef FLOW_AXIS_H
+#define FLOW_AXIS_H
+#ifndef ROOT_TObject
+# include <TObject.h>
+#endif
+
+//______________________________________________________
+/** @class AliFMDFlowAxis flow/AliFMDFlowAxis.h <flow/AliFMDFlowAxis.h>
+    @brief Axis object for the AliFMDFlowBinned1D and
+    AliFMDFlowBinned2D "histograms" of  objects of class AliFMDFlowBin. 
+    @ingroup c_binned 
+*/
+class AliFMDFlowAxis : public TObject
+{
+public:
+  /** Constructor 
+      @param n    Number of bins 
+      @param bins Bin limits (@a n+1 entries) */
+  AliFMDFlowAxis(UShort_t n, Double_t* bins);
+  /** Constructor 
+      @param n    Number of bins 
+      @param min  Minimum 
+      @param max  Maximum */
+  AliFMDFlowAxis(UShort_t n, Double_t min, Double_t max);
+  /** Copy constructor 
+      @param other Object to copy from */
+  AliFMDFlowAxis(const AliFMDFlowAxis& other);
+  /** Assignement operator 
+      @param other Object to assign from */ 
+  AliFMDFlowAxis& operator=(const AliFMDFlowAxis& other);
+  /** Destructor */ 
+  virtual ~AliFMDFlowAxis();
+  /** Find a bin 
+      @param v Value to look for 
+      @return bin number of @a v */
+  Int_t FindBin(Double_t v) const;
+  /** Get the width of the @a i th bin 
+      @param i  Bin to get width of */
+  Double_t BinWidth(UShort_t i) const;
+  /** Get the center of a bin 
+      @param i Bin to get center of 
+      @return Center of the @a i th bin */
+  Double_t BinCenter(UShort_t i) const;
+  /** Get the lower limit of a bin 
+      @param i Bin to get the lower limit of 
+      @return lower limit of bin @a i */
+  Double_t BinLower(UShort_t i) const;
+  /** Get the upper limit of a bin 
+      @param i Bin to get the upper limit of 
+      @return upper limit of bin @a i */
+  Double_t BinUpper(UShort_t i) const;
+  /** Get a pointer to the bins 
+      @return pointer to the bins */ 
+  Double_t* Bins() const { return fBins; }
+  /** Get the number of bins */ 
+  UShort_t N() const { return fN; }
+protected:
+  /** Number of bins */ 
+  UShort_t fN;
+  /** Borders of the bins */ 
+  Double_t* fBins; //[fN]
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowAxis,1);
+};  
+
+
+#endif
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowBessel.cxx b/FMD/flow/AliFMDFlowBessel.cxx
new file mode 100644 (file)
index 0000000..193e7bd
--- /dev/null
@@ -0,0 +1,319 @@
+/** @file 
+    @brief Implementation of Bessel functions */
+#include "flow/AliFMDFlowBessel.h"
+#include <cmath>
+#include <iostream>
+
+/** Namespace for utility functions */
+namespace 
+{ 
+  //__________________________________________________________________
+  /** Utility function to compute 
+      @f[ 
+      envj_n(x) = \frac12 \log_{10}(6.28 n) - n \log_{10}(1.36 x / n)
+      @f]
+      @param n Order 
+      @param x Argument 
+      @return @f$ envj_n(x)@f$ */
+  Double_t Envj(UInt_t n, Double_t x) { 
+    return .5 * log10(6.28 * n) - n * log10(1.36 * x / n);
+  }
+  //__________________________________________________________________
+  /** Utility function to do loop in Msta1 and Msta2 */ 
+  UInt_t MstaLoop(Int_t n0, Double_t a0, Float_t mp) {
+    Double_t f0 = Envj(n0, a0) - mp;
+    Int_t    n1 = n0 + 5;
+    Double_t f1 = Envj(n1, a0) - mp;
+    Int_t    nn = 0;
+    for (UInt_t i = 0; i < 20; i++) { 
+      nn = n1 - (n1 - n0) / (1. - f0 / f1);
+      if (fabs(nn - n1) < 1) break;
+      n0 = n1;
+      f0 = f1;
+      n1 = nn;
+      f1 = Envj(nn, a0) - mp;
+    }
+    return nn;
+  }
+  //__________________________________________________________________
+  /** Determine the starting point for backward recurrence such that
+      the magnitude of @f$J_n(x)@f$ at that point is about @f$
+      10^{-mp}@f$  
+      @param x   Argument to @f$ J_n(x)@f$ 
+      @param mp  Value of magnitude @f$ mp@f$ 
+      @return The starting point */
+  UInt_t Msta1(Double_t x, UInt_t mp) { 
+    Double_t a0 = fabs(x);
+    Int_t    n0 = Int_t(1.1 * a0) + 1;
+    return MstaLoop(n0, a0, mp);
+  }
+  //__________________________________________________________________
+  /** Determine the starting point for backward recurrence such that
+      all @f$ J_n(x)@f$ has @a mp significant digits.
+      @param x   Argument to @f$ J_n(x)@f$ 
+      @param n   Order of @f$ J_n@f$
+      @param mp  Number of significant digits 
+      @reutnr starting point */
+  UInt_t Msta2(Double_t x, Int_t n, Int_t mp) { 
+    Double_t a0  = fabs(x);
+    Double_t hmp = 0.5 * mp;
+    Double_t ejn = Envj(n, a0);
+    Double_t obj = 0;
+    Int_t    n0  = 0;
+    if (ejn <= hmp) {
+      obj = mp;
+      n0  = Int_t(1.1 * a0) + 1; //Bug for x<0.1-vl, 2-8.2002
+    }
+    else { 
+      obj = hmp + ejn;
+      n0  = n;
+    }
+    return MstaLoop(n0, a0, obj) + 10;
+  }
+}
+//____________________________________________________________________
+void 
+AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0, Double_t& bi1, Double_t& di1)
+{
+  Double_t       x2 =  x * x;
+  if (x == 0) { 
+    bi0 = 1;
+    bi1 = 0;
+    di0 = 0;
+    di1 = .5;
+    return;
+  }
+    
+  if (x <= 18) { 
+    bi0      = 1;
+    Double_t r = 1;
+    for (UInt_t k = 1; k <= 50; k++) { 
+      r   =  .25 * r * x2 / (k * k);
+      bi0 += r;
+      if (fabs(r / bi0) < 1e-15) break;
+    }
+    bi1 = 1;
+    r   = 1;
+    for (UInt_t k = 1; k <= 50; k++) { 
+      r   =  .25 * r * x2 / (k * (k + 1));
+      bi1 += r;
+    }
+    bi1 *= .5 * x;
+  }
+  else { 
+    const Double_t a[] = { 0.125,            0.0703125, 
+                        0.0732421875,     0.11215209960938, 
+                        0.22710800170898, 0.57250142097473, 
+                        1.7277275025845,  6.0740420012735, 
+                        24.380529699556,  110.01714026925,
+                        551.33589612202,  3038.0905109224 };
+    const Double_t b[] = { -0.375,          -0.1171875, 
+                        -0.1025390625,    -0.14419555664063,
+                        -0.2775764465332, -0.67659258842468,
+                        -1.9935317337513, -6.8839142681099,
+                        -27.248827311269, -121.59789187654,
+                        -603.84407670507, -3302.2722944809 };
+    UInt_t k0 = 12;
+    if (x >= 35) k0 = 9;
+    if (x >= 50) k0 = 7;
+      
+    Double_t ca = exp(x) / sqrt(2 * M_PI * x);
+    Double_t xr = 1. / x;
+    bi0       = 1;
+      
+    for (UInt_t k = 1; k <= k0; k++) {
+      bi0 += a[k-1] * pow(xr, Int_t(k));
+      bi1 += b[k-1] * pow(xr, Int_t(k));
+    }
+      
+    bi0 *= ca;
+    bi1 *= ca;
+  }
+  di0 = bi1;
+  di1 = bi0 - bi1 / x;
+}
+//____________________________________________________________________
+UInt_t 
+AliFMDFlowBessel::Ihalf(Int_t n, Double_t x, Double_t* bi, Double_t* di) 
+{
+  typedef Double_t (*fun_t)(Double_t);
+  Int_t   p = (n > 0 ? 1     : -1);
+  fun_t s = (n > 0 ? &sinh : &cosh);
+  fun_t c = (n > 0 ? &cosh : &sinh);
+  
+  // Temporary buffer 
+  Double_t tbi[p*n/2+2];
+  Double_t tdi[p*n/2+2];
+  
+  // Calculate I(-1/2) for n>0 or I(1/2) for n<0 
+  Double_t bim = sqrt(2) * (*c)(x) / sqrt(M_PI * x);
+
+  // We calculate one more than n, that is we calculate 
+  //   I(1/2), I(3/2),...,I(n/2+1) 
+  // because we need that for the negative orders 
+  for (Int_t i = 1; i <= p * n + 2; i += 2) { 
+    switch (i) { 
+    case 1: 
+      // Explicit formula for 1/2
+      tbi[0] = sqrt(2 * x / M_PI) * (*s)(x) / x;
+      break;
+    case 3:
+      // Explicit formula for 3/2
+      tbi[1] = sqrt(2 * x / M_PI) * ((*c)(x) / x - (*s)(x) / (x * x));
+      break;
+    default:
+      // Recursive formula for n/2 in terms of (n/2-2) and (n/2-1) 
+      tbi[i/2] = tbi[i/2-2] - (i-2) * tbi[i/2-1] / x;
+      break;
+    }
+  }
+  // We calculate the first one dI(1/2) here, so that we can use it in
+  // the loop. 
+  tdi[0] = bim - tbi[0] / (2 * x);
+
+  // Now, calculate dI(3/2), ..., dI(n/2)
+  for (Int_t i = 3; i <= p * n; i += 2) {
+    tdi[i/2] = tbi[i/2-p*1] - p * i * tbi[i/2] / (2 * x);
+  }
+
+  // Finally, store the output 
+  for (Int_t i = 0; i < p * n / 2 + 1; i++) { 
+    UInt_t j = (p < 0 ? p * n / 2 - i : i);
+    bi[i]    = tbi[j];
+    di[i]    = tdi[j];
+  }
+  return n / 2;
+}
+               
+//____________________________________________________________________
+UInt_t 
+AliFMDFlowBessel::Iwhole(UInt_t in, Double_t x, Double_t* bi, Double_t* di) 
+{
+  UInt_t mn = in;
+  if (x < 1e-100) { 
+    for (UInt_t k = 0; k <= in; k++) bi[k] = di[k]  = 0;
+    bi[0] = 1;
+    di[1] = .5;
+    return 1;
+  }
+  Double_t bi0, di0, bi1, di1;
+  I01(x, bi0, di0, bi1, di1);
+  bi[0] = bi0; di[0] = di0; bi[1] = bi1; di[1] = di1; 
+
+  if (in <= 1) return 2;
+    
+  if (x > 40 and Int_t(in) < Int_t(.25 * x)) { 
+    Double_t h0 = bi0;
+    Double_t h1 = bi1;
+    for (UInt_t k = 2; k <= in; k++) { 
+      bi[k] = -2 * (k - 1) / x * h1 + h0;
+      h0    = h1;
+      h1    = bi[k];
+    }
+  }
+  else { 
+    UInt_t      m  = Msta1(x, 100);
+    if (m < in) mn = m;
+    else        m  = Msta2(x, in, 15);
+      
+    Double_t f0 = 0;
+    Double_t f1 = 1e-100;
+    Double_t f  = 0;
+    for (Int_t k = m; k >= 0; k--) { 
+      f  = 2 * (k + 1) * f1 / x + f0;
+      if (k <= Int_t(mn)) bi[k] = f;
+      f0 = f1;
+      f1 = f;
+    }
+    Double_t s0 = bi0 / f;
+    for (UInt_t k = 0; k <= mn; k++) 
+      bi[k] *= s0;
+  }
+  for (UInt_t k = 2; k <= mn; k++) 
+    di[k] =  bi[k - 1] - k / x * bi[k];
+  return mn;
+}
+
+//____________________________________________________________________
+UInt_t 
+AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x, Double_t* bi, Double_t* di)
+{
+  UInt_t  in1 = UInt_t(fabs(n1));
+  UInt_t  in2 = UInt_t(fabs(n2));
+  UInt_t  nt  = UInt_t(n2 - n1 + 1);
+  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 "
+               << "must also be half integer" << std::endl;
+      return 0;
+    }
+    // std::cout << "Half-integers " << n1 << "," << n2 << std::endl;
+
+    Int_t s1    = (n1 < 0 ? -1 : 1);
+    Int_t s2    = (n2 < 0 ? -1 : 1);
+    Int_t l1    = (s1 < 0 ? (2*in1+1)/2 : 0);
+    Int_t l2    = (s1 > 0 ? in1 : 0);
+    Double_t tbi[nt+1], tdi[nt+1];
+    if (s1 < 0) { 
+      Ihalf(2 * n1, x, tbi, tdi);
+      for (Int_t i = 0; i < l1 && i < Int_t(nt); i++) { 
+       bi[i] = tbi[i];
+       di[i] = tdi[i];
+      }
+    }
+    if (s2 > 0) { 
+      // std::cout << "Evaluating Ihalf(" << 2 * n2 << "," << x << ",...)";
+      Ihalf(2 * n2, x, tbi, tdi);
+      for (Int_t i = l1; i <= 2 * Int_t(in2) && i < Int_t(nt); i++) { 
+       UInt_t j = i + l2 - l1;
+       bi[i] = tbi[j];
+       di[i] = tdi[j];
+      }
+    }
+    return nt;
+  }
+  if (Int_t(n1) != n1 || Int_t(n2) != n2) { 
+    std::cerr << "n1 (" << n1 << "/" << in1 << ") and "
+             << "n2 (" << n2 << "/" << in2 << ") must be "
+             << "half-integer or integer" << std::endl;
+    return 0; // Not integer!
+  }
+
+  UInt_t  n   = UInt_t(in1 > in2 ? in1 : in2);
+  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" 
+             << std::endl;
+  for (UInt_t i = 0; i < nt; i++) { 
+    UInt_t j = (i + n1 < 0 ? -n1-i : n1+i);
+    bi[i]    = tbi[j];
+    di[i]    = tdi[j];
+  }
+  return nt;
+}
+
+    
+//____________________________________________________________________
+Double_t 
+AliFMDFlowBessel::I(Double_t n, Double_t x) 
+{ 
+  Double_t i[2], di[2];
+  Inu(n, n, x, i, di);
+  return i[0];
+}
+      
+//____________________________________________________________________
+Double_t 
+AliFMDFlowBessel::DiffI(Double_t n, Double_t x) 
+{ 
+  Double_t i[2], di[2];
+  Inu(n, n, x, i, di);
+  return di[0];
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowBessel.h b/FMD/flow/AliFMDFlowBessel.h
new file mode 100644 (file)
index 0000000..a1011b2
--- /dev/null
@@ -0,0 +1,110 @@
+// -*- mode: C++ -*-
+/** @file 
+    @brief Declaration of Bessel functions */
+#ifndef Bessel_h
+#define Bessel_h
+#include <cctype>
+#include <Rtypes.h>
+
+/** @defgroup z_bessel Bessel utility functions.  
+    @brief This group contains a number of functions to calculate
+    the modified Bessel function of the first kind @f$
+    I_{\nu}(x)@f$ for integer and half-integer values of @f$
+    \nu@f$, of any real number @f$ x@f$.
+    
+    The main entry of these functions is the function Inu, which
+    in turn call I01, Iwhole, or Ihalf.  Inu always returns a list
+    of the function values for the specified values of @f$ \nu@f$.
+    
+    The convinience functions I and DiffI returns a single value
+    of @f$ I_{\nu}(x)@f$ or @f$ dI_{\nu}(x)/dx@f$ for a specified
+    value of @f$ \nu@f$. 
+    
+    @example test_bessel.cxx */
+
+/** Namespace for Bessel functions 
+    @ingroup z_bessel */
+namespace AliFMDFlowBessel
+{
+  /** @{ 
+      @ingroup z_bessel */
+  /** Compute the modified Bessel functions 
+      @f[
+      I_0(x) = \sum_{k=1}^\infty\frac{\left(\frac{x^2}{4}\right)^k}{(k!)^2}
+      @f]
+      and @f$ I_1(x)@f$, and their first derivatives  
+      @param x   Argument @f$ x@f$ 
+      @param bi0 On return, @f$ I_0(x)@f$ 
+      @param di0 On return, @f$ I_0'(x)@f$ 
+      @param bi1 On return, @f$ I_1(x)@f$ 
+      @param di1 On return, @f$ I_1'(x)@f$ 
+  */
+  void I01(Double_t x, 
+          Double_t& bi0, Double_t& di0, Double_t& bi1, Double_t& di1);
+  
+  /** Compute the modified Bessel functions @f$ I_{n/2}(x)@f$ and their
+      derivatives.  
+      @param x  Argument. 
+      @param n  Order
+      @param bi On output, @f$ I_{1/2}(x), \ldots, I_{n/2}(x)@f$ for
+      @f$ n > 0@f$ and @f$ I_{-n/2}(x), \ldots, I_{-1/2}(x)@f$ for
+      @f$ n < 0@f$ 
+      @param di On output, @f$ I_{1/2}'(x), \ldots, I_{n/2}'(x)@f$ for
+      @f$ n > 0@f$ and @f$ I_{-n/2}'(x), \ldots, I_{-1/2}'(x)@f$ for
+      @f$ n < 0@f$ 
+      @return number of valid entries in @a bi and @a di */
+  UInt_t Ihalf(Int_t n, Double_t x, Double_t* bi, Double_t* di);
+  /** Compute the modified Bessel functions @f$ I_n(x)@f$ and their
+      derivatives.  Note, that @f$ I_{-n}(x) = I_n(x)@f$ and 
+      @f$ dI_{-n}(x)/dx = dI_{n}(x)/dx@f$ so this function can be used
+      to evaluate @f$ I_n \forall n \in \mathcal{Z}@f$ 
+      @param x  Argument. 
+      @param n  Order
+      @param bi On output, @f$ I_0(x), \ldots, I_{mn}(x)@f$ 
+      @param di On output, @f$ I_0'(x), \ldots, I_{mn}'(x)@f$ 
+      @return number of valid entries in @a bi and @a di */
+  UInt_t Iwhole(UInt_t n, Double_t x, Double_t* bi, Double_t* di);
+
+  /** Compute the modified Bessel functions @f$ I_{\nu}(x)@f$ and
+      their derivatives for @f$ \nu = n_1, n_1+1, \ldots, n_2@f$ for
+      and integer @f$ n_1, n_2@f$ or any half-integer @f$ n_1,
+      n_2@f$. 
+      @param x  Argument. 
+      @param n1 Lower order
+      @param n2 Upper order
+      @param bi On output, @f$ I_{n_1}(x), \ldots, I_{n_2}(x)@f$ 
+      @param di On output, @f$ I_{n_1}'(x), \ldots, I_{n_2}'(x)@f$ 
+      @return number of valid entries in @a bi and @a di */
+  UInt_t Inu(Double_t n1, Double_t n2, Double_t x, Double_t* bi, Double_t* di);
+
+  /** Compute the modified Bessel function of the first kind 
+      @f[
+      I_n(x) = \left(\frac{x}{2}\right)^n
+      \sum_{k=0}^\infty\frac{\left(\frac{x^2}{4}\right)^k}
+      {k!\Gamma(n+k+1)}
+      @f]
+      for arbirary integer order @f$ n@f$ 
+      @param n  Order
+      @param x  Argument 
+      @return @f$ I_n(x)@f$ */ 
+  Double_t I(Double_t n, Double_t x);
+      
+  /** Compute the derivative of the modified Bessel function of the
+      first kind 
+      @f[ 
+      \frac{dI_n(x)}{dx} = I_{n-1}(x) - \frac{n}{x} I_{n}(x)
+      @f]
+      for arbirary integer order @f$ n@f$
+      @param n Order 
+      @param x Argument 
+      @return @f$ \frac{dI_n(x)}{dx}@f$ */ 
+  Double_t DiffI(Double_t n, Double_t x);
+  /** @} */
+}
+
+#endif
+//
+// EOF
+//
+
+
diff --git a/FMD/flow/AliFMDFlowBin.cxx b/FMD/flow/AliFMDFlowBin.cxx
new file mode 100644 (file)
index 0000000..eeb7842
--- /dev/null
@@ -0,0 +1,150 @@
+/** @file 
+    @brief implementation of a Bin in a Flow histogram */
+#include "flow/AliFMDFlowBin.h"
+#include <cmath>
+#include <iostream>
+#include <iomanip>
+
+//====================================================================
+void 
+AliFMDFlowBin::Begin() 
+{
+  // Clear event plane calculators 
+  fPsi.Clear();
+  fPsiA.Clear();
+  fPsiB.Clear();
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBin::AddToEventPlane(Double_t phi, Double_t w, Bool_t a) 
+{
+  fPsi.Add(phi, w);
+  if (a) fPsiA.Add(phi, w);
+  else   fPsiB.Add(phi, w);
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBin::AddToHarmonic(Double_t phi, Double_t w)
+{
+  // Disregard the obervation of phi from the event plane angle. 
+  Double_t psi   = fPsi.Psi(phi, w);
+  fHarmonic.Add(phi, psi);
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBin::End()
+{
+  Double_t psi_A = fPsiA.Psi();
+  Double_t psi_B = fPsiB.Psi();
+
+  // Update the resolutions 
+  fRes.Add(psi_A, psi_B);
+  fResStar.Add(psi_A, psi_B);
+  fResTdr.Add(psi_A, psi_B);
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBin::Event(Double_t* phis, Double_t* ws, UInt_t n) 
+{ 
+  Begin();
+  
+  // Calculate split. 
+  UInt_t split = n / 2;
+  // First sub-event. 
+  for (UInt_t i = 0; i < split; i++) 
+    AddToEventPlane(phis[i], (ws ? ws[i] : 1), kTRUE);
+  // Second sub-event. 
+  for (UInt_t i = split; i < n; i++) 
+    AddToEventPlane(phis[i], (ws ? ws[i] : 1), kFALSE);
+  // Add contributions to the harmonic. 
+  for (UInt_t i = 0; i < n; i++)     
+    AddToHarmonic(phis[i], (ws ? ws[i] : 1));
+
+  End();
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowBin::Value(CorType t) const
+{ 
+  Double_t e;
+  return Value(e, t);
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowBin::EValue(CorType t) const 
+{ 
+  Double_t e2;
+  Value(e2, t);
+  return sqrt(e2);
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowBin::Value(Double_t& e2, CorType t) const
+{ 
+  Double_t r, er2;
+  r = Correction(er2, t);
+  return fHarmonic.Value(r, er2, e2);
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowBin::Correction(Double_t& er2, CorType t) const
+{
+  Double_t r = 1;
+  UShort_t k = fHarmonic.Order()/fRes.Order();
+  switch (t) { 
+  case naive: r = fRes.Correction(k, er2);     break;
+  case star:  r = fResStar.Correction(k, er2); break;
+  case tdr:   r = fResTdr.Correction(k, er2);  break;
+  default:    r = 1; er2 = 0;                  break;
+  }
+  return r;
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBin::Finish() 
+{}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBin::Print(Option_t*) const
+{
+  Double_t e2v[4], v[4], r[4], e2r[4];
+  const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
+  v[0] = 100 * Value(e2v[0], AliFMDFlowBin::none);
+  v[1] = 100 * Value(e2v[1], AliFMDFlowBin::naive);
+  v[2] = 100 * Value(e2v[2], AliFMDFlowBin::star);
+  v[3] = 100 * Value(e2v[3], AliFMDFlowBin::tdr);
+  r[0] = 100 * Correction(e2r[0], AliFMDFlowBin::none);
+  r[1] = 100 * Correction(e2r[1], AliFMDFlowBin::naive);
+  r[2] = 100 * Correction(e2r[2], AliFMDFlowBin::star);
+  r[3] = 100 * Correction(e2r[3], AliFMDFlowBin::tdr);
+  
+  std::streamsize         old_prec  = std::cout.precision(3);
+  std::ios_base::fmtflags old_flags = std::cout.setf(std::ios_base::fixed, 
+                                                    std::ios_base::floatfield);
+  std::cout << "  v" << std::setw(1) << fHarmonic.Order() << ": ";
+  for (UInt_t i = 0; i < 4; i++) 
+    std::cout << std::setw(6+(i == 0 ? 0 : 6)) << names[i] << ": " 
+             << std::setw(6) << v[i] << " +/- " 
+             << std::setw(6) << 100*sqrt(e2v[i]) << " ["
+             << std::setw(7) << r[i] << " +/- " 
+             << std::setw(7) << 100*sqrt(e2r[i]) << "]\n";
+  std::cout << std::flush;
+  std::cout.precision(old_prec);
+  std::cout.setf(old_flags, std::ios_base::floatfield);
+}
+
+
+//____________________________________________________________________
+//
+// EOF
+// 
diff --git a/FMD/flow/AliFMDFlowBin.h b/FMD/flow/AliFMDFlowBin.h
new file mode 100644 (file)
index 0000000..67205ea
--- /dev/null
@@ -0,0 +1,135 @@
+// -*- mode: C++ -*-
+/** @file 
+    @brief Declaration of a Bin in a Flow "histogram" */
+#ifndef FLOW_BIN_H
+#define FLOW_BIN_H
+#include <flow/AliFMDFlowEventPlane.h>
+#include <flow/AliFMDFlowHarmonic.h>
+#include <flow/AliFMDFlowResolution.h>
+#include <TObject.h>
+
+/** @defgroup c_binned Binned flow 
+    @brief This group contains code for binned flow analysis.  Two
+    kinds of "histograms" are defined - a 1 dimensional and a 2
+    dimensional set of binned objects of class AliFMDFlowBin.   
+
+    Objects of class AliFMDFlowBin contains all the code needed to compute
+    flow in a given bin.   
+
+    The class AliFMDFlowAxis encodes look-up of a object of class
+    AliFMDFlowBin in a flow "Histogram"
+*/
+//______________________________________________________
+/** @class AliFMDFlowBin flow/AliFMDFlowBin.h <flow/AliFMDFlowBin.h>
+    @brief A bin of flow.   
+
+    This contains an of class AliFMDFlowHarmonic and an object of
+    class AliFMDFlowEventPlane to calculate @f$ v_n@f$ and
+    @f$\Psi_k@f$.  It contain two objects of class
+    AliFMDFlowEventPlane to calculate the sub-event event planes
+    @f$\Psi_A, \Psi_B@f$.  It also contain 3 objects of class
+    AliFMDFlowResolution to calculate the event plane angle
+    resolution.
+
+    @ingroup c_binned 
+*/
+class AliFMDFlowBin : public TObject
+{
+public:
+  /** Correction type */
+  enum CorType {
+    /** No correction */
+    none, 
+    /** Naive, using the formulas in Voloshins paper */
+    naive,
+    /** STARs way */
+    star, 
+    /** The way used in the TDR */
+    tdr
+  };
+  /** Constructor */
+  AliFMDFlowBin(UShort_t order, UShort_t k=1) 
+    : fPsi(order / k), 
+      fPsiA(order / k), 
+      fPsiB(order / k), 
+      fRes(order / k), 
+      fResStar(order / k), 
+      fResTdr(order / k),
+      fHarmonic(order) 
+  {}
+  /** Destructor */
+  virtual ~AliFMDFlowBin() {} 
+  /** Should be called at the start of an event */ 
+  virtual void Begin();
+  /** Called to add a contribution to the event plane 
+      @param phi The angle @f$ \varphi \in[0,2\pi]@f$ 
+      @param w   Weight
+      @param a   If true, add to sub-event A, otherwise to sub-event
+      B. */
+  virtual void AddToEventPlane(Double_t phi, Double_t w=1, Bool_t a=kTRUE);
+  /** Called to add a contribution to the harmonic. 
+      @param phi The angle @f$ \varphi \in[0,2\pi]@f$
+      @param w   Weight of @a phi (only used in the calculation of
+      the event plane). */
+  virtual void AddToHarmonic(Double_t phi, Double_t w=1);
+  /** Should be called at the end of an event */ 
+  virtual void End();
+  /** Analyse events 
+      @param phis @f$ (\varphi_i, \ldots, \varphi_n)@f$ 
+      @param ws   Weights (optional)
+      @param n    Size of @a phis and possibly @a ws */
+  virtual void Event(Double_t* phis, Double_t* ws, UInt_t n);
+  /** Finish run */
+  virtual void Finish();
+  /** Get the value in this bin 
+      @param t  Which type of correction
+      @return the value of the harmonic */
+  virtual Double_t Value(CorType t=naive) const;
+  /** Get the value in this bin 
+      @param t  Which type of correction 
+      @return the error on the value of the harmonic */
+  virtual Double_t EValue(CorType t=naive) const;
+  /** Get the value in this bin 
+      @param e2 On return, the square error. 
+      @param t  Which type of correction
+      @return the value of the harmonic */
+  virtual Double_t Value(Double_t& e2, CorType t=naive) const;
+  /** Get the value in this bin 
+      @param e2 On return, the square error. 
+      @param t  Which type  of correction
+      @return the value of the harmonic */
+  virtual Double_t Correction(Double_t& e2, CorType t=naive) const;
+  /** Print summary to standard output */ 
+  virtual void Print(Option_t* option="") const;
+    
+  /** Get the event plane angle */
+  virtual Double_t Psi() const { return fPsi.Psi(); } 
+  /** Get the sub-event A plane angle */
+  virtual Double_t PsiA() const { return fPsiA.Psi(); } 
+  /** Get the sub-event B plane angle */
+  virtual Double_t PsiB() const { return fPsiB.Psi(); } 
+
+protected:
+  /** Major event plane */
+  AliFMDFlowEventPlane fPsi;
+  /** Sub-event A event plane */
+  AliFMDFlowEventPlane fPsiA;
+  /** Sub-event B event plane */
+  AliFMDFlowEventPlane fPsiB;
+  /** Resolution */
+  AliFMDFlowResolution fRes;
+  /** Resolution */
+  AliFMDFlowResolutionStar fResStar;
+  /** Resolution */
+  AliFMDFlowResolutionTDR fResTdr;
+  /** The harmonic */
+  AliFMDFlowHarmonic fHarmonic;
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowBin,1);
+};
+
+
+#endif
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowBinned1D.cxx b/FMD/flow/AliFMDFlowBinned1D.cxx
new file mode 100644 (file)
index 0000000..113fa4e
--- /dev/null
@@ -0,0 +1,201 @@
+/** @file 
+    @brief Implementation of a 1-dimensional Flow "histogram" */
+#include "flow/AliFMDFlowBinned1D.h"
+#include "flow/AliFMDFlowBin.h"
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
+#include <TBrowser.h>
+#include <TString.h>
+
+//====================================================================
+AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order, 
+                                      UShort_t nxbins, 
+                                      Double_t* xbins) 
+  : fXAxis(nxbins, xbins),
+    fBins(0)
+{
+  UShort_t n = fXAxis.N();
+  fBins   = new AliFMDFlowBin*[n];
+  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
+}
+//____________________________________________________________________
+AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order, 
+                                      const AliFMDFlowAxis& xaxis)
+  : fXAxis(xaxis)
+{
+  UShort_t n = fXAxis.N();
+  fBins   = new AliFMDFlowBin*[n];
+  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
+}
+//____________________________________________________________________
+AliFMDFlowBinned1D::AliFMDFlowBinned1D(const AliFMDFlowBinned1D& o)
+  : fXAxis(o.fXAxis)
+{
+  UShort_t n = fXAxis.N();
+  fBins   = new AliFMDFlowBin*[n];
+  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
+}
+//____________________________________________________________________
+AliFMDFlowBinned1D&
+AliFMDFlowBinned1D::operator=(const AliFMDFlowBinned1D& o)
+{
+  if (fBins) { 
+    for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
+    delete [] fBins;
+  }
+  fXAxis     = o.fXAxis;
+  UShort_t n = fXAxis.N();
+  fBins   = new AliFMDFlowBin*[n];
+  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
+  return *this;
+}
+  
+//____________________________________________________________________
+AliFMDFlowBinned1D::~AliFMDFlowBinned1D()
+{
+  if (fBins) { 
+    for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
+    delete [] fBins;
+  }
+}
+
+//____________________________________________________________________
+AliFMDFlowBin* 
+AliFMDFlowBinned1D::GetBin(UShort_t i) const
+{
+  if (i >= fXAxis.N()) return 0;
+  return fBins[i];
+}
+//____________________________________________________________________
+AliFMDFlowBin* 
+AliFMDFlowBinned1D::GetBin(Double_t x) const
+{
+  Int_t i = fXAxis.FindBin(x);
+  if (i < 0) return 0;
+  UShort_t j = i;
+  return GetBin(j);
+}
+  
+//____________________________________________________________________
+void 
+AliFMDFlowBinned1D::Begin()
+{
+  for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Begin();
+}
+//____________________________________________________________________
+void 
+AliFMDFlowBinned1D::End()
+{
+  for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->End();
+}
+//____________________________________________________________________
+void 
+AliFMDFlowBinned1D::Finish()
+{
+  for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Finish();
+}
+//____________________________________________________________________
+Bool_t 
+AliFMDFlowBinned1D::AddToEventPlane(Double_t x, Double_t phi, Double_t w, Bool_t a)
+{
+  AliFMDFlowBin* bin = GetBin(x);
+  if (!bin) return kFALSE;
+  bin->AddToEventPlane(phi, w, a);
+  return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t 
+AliFMDFlowBinned1D::AddToHarmonic(Double_t x, Double_t phi)
+{
+  AliFMDFlowBin* bin = GetBin(x);
+  if (!bin) return kFALSE;
+  bin->AddToHarmonic(phi);
+  return kTRUE;
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBinned1D::Event(Double_t* phis, Double_t* xs, Double_t* ws, ULong_t n)
+{
+  Begin();
+  for (UInt_t i = 0; i < n; i++) 
+    AddToEventPlane(xs[i], phis[i], (ws ? ws[i] : 1), 
+                   Float_t(rand()) / RAND_MAX > 0.5);
+  for (UInt_t i = 0; i < n; i++) 
+    AddToHarmonic(xs[i], phis[i]);
+  End();
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBinned1D::Browse(TBrowser* b)
+{
+  b->Add(&fXAxis, "xaxis");
+  for (UInt_t i = 0; i < fXAxis.N(); i++) 
+    b->Add(fBins[i], Form("bin_%03d", i));
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBinned1D::Print(Option_t* option) const
+{
+  TString opt(option);
+  opt.ToLower();
+  Bool_t det = opt.Contains("d");
+  Bool_t sum = opt.Contains("s");
+  if (det) { 
+    for (UShort_t i = 0; i < fXAxis.N(); i++) { 
+      Double_t x = fXAxis.BinCenter(i);
+      std::streamsize         old_p = std::cout.precision(3);
+      std::ios_base::fmtflags old_f = std::cout.setf(std::ios_base::fixed, 
+                                                    std::ios_base::floatfield);
+      std::cout << "x=" << std::setw(5) << x << std::endl;
+      fBins[i]->Print();
+      std::cout.precision(old_p);
+      std::cout.setf(old_f, std::ios_base::floatfield);
+    }
+  }
+  
+  if (sum) { 
+    UInt_t       nType = 4;
+    const char*  names[] = { "Bare",    "Naive",    "STAR",    "TDR" };
+    AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::none, 
+                                      AliFMDFlowBin::naive, 
+                                      AliFMDFlowBin::star, 
+                                      AliFMDFlowBin::tdr };
+    std::cout << "    x";
+    for (UInt_t i = 0; i < nType; i++) 
+      std::cout << " | " << std::setw(6+6+5) << names[i];
+    std::cout << "\n-----" << std::setfill('-');
+    for (UInt_t i = 0; i < nType; i++) 
+      std::cout << "-+-" <<  std::setw(6+6+5) << "-";
+    std::cout << std::setfill(' ') << std::endl;
+    
+    std::streamsize         old_p = std::cout.precision(2);
+    std::ios_base::fmtflags old_f = std::cout.setf(std::ios_base::fixed, 
+                                                  std::ios_base::floatfield);
+    for (UShort_t i = 0; i < fXAxis.N(); i++) { 
+      Double_t x = fXAxis.BinCenter(i);
+      std::cout << std::setprecision(2) << std::setw(5) << x << std::flush;
+      for (UShort_t j = 0; j < nType; j++) { 
+       Double_t e2v;
+       Double_t v  = fBins[i]->Value(e2v, types[j]);
+       std::cout << std::setprecision(3)    << " | "
+                 << std::setw(6) << 100 * v << " +/- " 
+                 << std::setw(6) << 100 * sqrt(e2v);
+      }
+      std::cout << std::endl;
+    }
+    std::cout.precision(old_p);
+    std::cout.setf(old_f, std::ios_base::floatfield);
+  }
+}
+
+    
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowBinned1D.h b/FMD/flow/AliFMDFlowBinned1D.h
new file mode 100644 (file)
index 0000000..a145f61
--- /dev/null
@@ -0,0 +1,95 @@
+// -*- mode: C++ -*-
+/** @file 
+    @brief Declaration of a 1-dimensional Flow "histogram" */
+#ifndef FLOW_BINNED1D_H
+#define FLOW_BINNED1D_H
+#include <flow/AliFMDFlowAxis.h>
+
+// Forward declaration 
+class AliFMDFlowBin;
+class TBrowser;
+
+//______________________________________________________
+/** @class AliFMDFlowBinned1D flow/AliFMDFlowBinned1D.h <flow/AliFMDFlowBinned1D.h>
+    @brief A 1 dimensional "histogram" of objects of class AliFMDFlowBin. 
+    @ingroup c_binned 
+    @example test_flow.cxx 
+    @example ana_flow.cxx 
+*/
+class AliFMDFlowBinned1D : public TObject
+{
+public:
+  /** Flags for print */
+  enum {
+    /** Show details */ 
+    details = 0x1, 
+    /** Show summary */ 
+    summary = 0x2
+  };
+  /** Constructor 
+      @param order    Order of the harmonic
+      @param nxbins   Number of X bins.
+      @param xbins    Borders of X bins (@a nxbins+1 entries) */
+  AliFMDFlowBinned1D(UShort_t order, UShort_t nxbins, Double_t* xbins);
+  /** Constructor 
+      @param order    Order of the harmonic
+      @param xaxis Axis object  */
+  AliFMDFlowBinned1D(UShort_t order, const AliFMDFlowAxis& xaxis);
+  /** Copy constructor */
+  AliFMDFlowBinned1D(const AliFMDFlowBinned1D& other);
+  /** Copy constructor */
+  AliFMDFlowBinned1D& operator=(const AliFMDFlowBinned1D& other);
+  /** Destructor */
+  virtual ~AliFMDFlowBinned1D();
+  
+  /** Called at the beginning of an event */
+  virtual void Begin();
+  /** Called at the end of an event */ 
+  virtual void End();
+  /** Called to add a contribution to the event plane 
+      @param x   Bin to fill into 
+      @param w   Weight
+      @param phi The angle @f$ \varphi@f$ in radians 
+      @param a   If true, add to sub-event A, otherwise sub-event B */ 
+  virtual Bool_t AddToEventPlane(Double_t x, Double_t phi, 
+                                Double_t w, Bool_t a);
+  /** Called to add a contribution to the harmonic
+      @param x   Bin to fill into 
+      @param phi The angle @f$ \varphi@f$ in radians */
+  virtual Bool_t AddToHarmonic(Double_t x, Double_t phi);
+  /** Process a full event. 
+      @param phis List of @a n @f$ \varphi=[0,2\pi]@f$ angles 
+      @param xs   List of @a n @f$ x@f$ values. 
+      @param ws   Weights
+      @param n    Size of @a phis and @a xs */
+  virtual void Event(Double_t* phis, Double_t* xs, Double_t* ws, 
+                    ULong_t n);
+  /** Called at the end of a run */
+  virtual void Finish();
+  /** Get the bin at @a x 
+      @param x The bin value to find a flow object for. 
+      @return The flow object at @a x or 0 if out of range */ 
+  virtual AliFMDFlowBin* GetBin(Double_t x) const;
+  /** Get the bin @a i 
+      @param i The bin number to get
+      @return The flow object in bin @a i or 0 if out of range */ 
+  virtual AliFMDFlowBin* GetBin(UShort_t i) const;
+  /** Print to standard out */ 
+  virtual void Print(Option_t* option="s") const;
+  /** Whether this is to be considered a folder */
+  Bool_t IsFolder() const { return kTRUE; }
+  /** Browse this object */ 
+  void Browse(TBrowser* b);
+protected:
+  /** X axis */ 
+  AliFMDFlowAxis fXAxis;
+  /** Array of the flow objects */ 
+  AliFMDFlowBin** fBins;
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowBinned1D,1);
+};
+
+#endif
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowBinned2D.cxx b/FMD/flow/AliFMDFlowBinned2D.cxx
new file mode 100644 (file)
index 0000000..fa655c0
--- /dev/null
@@ -0,0 +1,159 @@
+/** @file 
+    @brief Implementation of a 2-dimensional Flow "histogram" */
+#include "flow/AliFMDFlowBinned2D.h"
+#include "flow/AliFMDFlowBin.h"
+#include <cmath>
+#include <cstdlib>
+#include <TString.h>
+#include <TBrowser.h>
+
+//====================================================================
+AliFMDFlowBinned2D::AliFMDFlowBinned2D(UShort_t order, 
+                        UShort_t nxbins, Double_t* xbins,
+                        UShort_t nybins, Double_t* ybins) 
+  : fXAxis(nxbins, xbins),
+    fYAxis(nybins, ybins),
+    fBins(0)
+{
+  UInt_t n = fXAxis.N() * fYAxis.N();
+  fBins   = new AliFMDFlowBin*[n];
+  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
+}
+//____________________________________________________________________
+AliFMDFlowBinned2D::AliFMDFlowBinned2D(UShort_t order, 
+                        const AliFMDFlowAxis&    xaxis, 
+                        const AliFMDFlowAxis&    yaxis)
+  : fXAxis(xaxis), 
+    fYAxis(yaxis)
+{
+  UShort_t n = fXAxis.N() * fYAxis.N();
+  fBins   = new AliFMDFlowBin*[n];
+  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
+}
+//____________________________________________________________________
+AliFMDFlowBinned2D::AliFMDFlowBinned2D(const AliFMDFlowBinned2D& o)
+  : TObject(o), 
+    fXAxis(o.fXAxis), 
+    fYAxis(o.fYAxis)
+{
+  UShort_t n = fXAxis.N() * fYAxis.N();
+  fBins   = new AliFMDFlowBin*[n];
+  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
+}
+//____________________________________________________________________
+AliFMDFlowBinned2D&
+AliFMDFlowBinned2D::operator=(const AliFMDFlowBinned2D& o)
+{
+  if (fBins) { 
+    UInt_t n = fXAxis.N() * fYAxis.N();
+    for (UInt_t i = 0; i < n; i++) delete fBins[i];
+    delete [] fBins;
+  }
+  fXAxis     = o.fXAxis;
+  UShort_t n = fXAxis.N() * fYAxis.N();
+  fBins   = new AliFMDFlowBin*[n];
+  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
+  return *this;
+}
+
+//____________________________________________________________________
+AliFMDFlowBinned2D::~AliFMDFlowBinned2D()
+{
+  if (fBins) { 
+    UInt_t n = fXAxis.N() * fYAxis.N();
+    for (UInt_t i = 0; i < n; i++) delete fBins[i];
+    delete [] fBins;
+  }
+}
+//____________________________________________________________________
+AliFMDFlowBin* 
+AliFMDFlowBinned2D::GetBin(UShort_t i, UShort_t j) const
+{
+  if (i >= fXAxis.N() || j >= fYAxis.N()) return 0;
+  return fBins[i * fYAxis.N() + j];
+}
+//____________________________________________________________________
+AliFMDFlowBin* 
+AliFMDFlowBinned2D::GetBin(Double_t x, Double_t y) const
+{
+  Int_t i = fXAxis.FindBin(x);
+  if (i < 0) return 0;
+  Int_t j = fYAxis.FindBin(y);
+  if (j < 0) return 0;
+  UShort_t k = i;
+  UShort_t l = j;
+  return GetBin(k, l);
+}
+//____________________________________________________________________
+void 
+AliFMDFlowBinned2D::Begin()
+{
+  UInt_t n = fXAxis.N() * fYAxis.N();
+  for (UInt_t i = 0; i < n; i++) fBins[i]->Begin();
+}
+//____________________________________________________________________
+void 
+AliFMDFlowBinned2D::End()
+{
+  UInt_t n = fXAxis.N() * fYAxis.N();
+  for (UInt_t i = 0; i < n; i++) fBins[i]->End();
+}
+//____________________________________________________________________
+void 
+AliFMDFlowBinned2D::Finish()
+{
+  UInt_t n = fXAxis.N() * fYAxis.N();
+  for (UInt_t i = 0; i < n; i++) fBins[i]->Finish();
+}
+//____________________________________________________________________
+Bool_t 
+AliFMDFlowBinned2D::AddToEventPlane(Double_t x, Double_t y, Double_t phi, 
+                               Double_t w, Bool_t a)
+{
+  AliFMDFlowBin* bin = GetBin(x, y);
+  if (!bin) return kFALSE;
+  bin->AddToEventPlane(phi, w, a);
+  return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t 
+AliFMDFlowBinned2D::AddToHarmonic(Double_t x, Double_t y, Double_t phi)
+{
+  AliFMDFlowBin* bin = GetBin(x, y);
+  if (!bin) return kFALSE;
+  bin->AddToHarmonic(phi);
+  return kTRUE;
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBinned2D::Event(Double_t* phis, Double_t* xs, Double_t* ys, 
+                         Double_t* ws, ULong_t n)
+{
+  Begin();
+  for (UInt_t i = 0; i < n; i++) 
+    AddToEventPlane(xs[i], ys[i], phis[i], (ws ? ws[i] : 1), 
+                   Float_t(rand()) / RAND_MAX > 0.5);
+  for (UInt_t i = 0; i < n; i++) 
+    AddToHarmonic(xs[i], ys[i], phis[i]);
+  End();
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowBinned2D::Browse(TBrowser* b)
+{
+  b->Add(&fXAxis, "xaxis");
+  b->Add(&fYAxis, "yaxis");
+  for (UInt_t i = 0; i < fXAxis.N(); i++) { 
+    for (UInt_t j = 0; i < fYAxis.N(); j++) { 
+      b->Add(fBins[i*fXAxis.N()+j], Form("bin_%03d_%03d", i, j));
+    }
+  }
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowBinned2D.h b/FMD/flow/AliFMDFlowBinned2D.h
new file mode 100644 (file)
index 0000000..08a47f1
--- /dev/null
@@ -0,0 +1,100 @@
+qq// -*- mode: C++ -*-
+/** @file 
+    @brief Declaration of a 2-dimensional Flow "histogram" */
+#ifndef FLOW_BINNED2D_H
+#define FLOW_BINNED2D_H
+#include <flow/AliFMDFlowAxis.h>
+#include <TObject.h>
+
+// Forward declaration 
+class AliFMDFlowBin;
+class TBrowser;
+
+//______________________________________________________
+/** @class AliFMDFlowBinned2D flow/AliFMDFlowBinned2D.h <flow/AliFMDFlowBinned2D.h>
+    @brief A 2 dimensional "histogram" of objects of class AliFMDFlowBin. 
+    @ingroup c_binned */
+class AliFMDFlowBinned2D : public TObject
+{
+public:
+  /** Constructor 
+      @param order    Order of the harmonic 
+      @param nxbins   Number of X bins.
+      @param xbins    Borders of X bins (@a nxbins+1 entries) 
+      @param nybins   Number of Y bins.
+      @param ybins    Borders of Y bins (@a nybins+1 entries) */
+  AliFMDFlowBinned2D(UShort_t order, 
+                    UShort_t nxbins, Double_t* xbins,
+                    UShort_t nybins, Double_t* ybins);
+  /** Constructor 
+      @param order Order of the harmonic 
+      @param xaxis Axis object 
+      @param yaxis Axis object */
+  AliFMDFlowBinned2D(UShort_t order, const AliFMDFlowAxis& xaxis, 
+                    const AliFMDFlowAxis& yaxis);
+      
+  /** Copy constructor */
+  AliFMDFlowBinned2D(const AliFMDFlowBinned2D& other);
+  /** Copy constructor */
+  AliFMDFlowBinned2D& operator=(const AliFMDFlowBinned2D& other);
+  /** Destructor */
+  virtual ~AliFMDFlowBinned2D();
+  
+  /** Called at the beginning of an event */
+  virtual void Begin();
+  /** Called at the end of an event */ 
+  virtual void End();
+  /** Called to add a contribution to the event plane 
+      @param x   Bin to fill into 
+      @param y   Bin to fill into 
+      @param phi The angle @f$ \varphi@f$ in radians 
+      @param w   Weight
+      @param a   If true, add to sub-event A, otherwise sub-event B */ 
+  virtual Bool_t AddToEventPlane(Double_t x, Double_t y, Double_t phi, 
+                                Double_t w, Bool_t a);
+  /** Called to add a contribution to the harmonic
+      @param x   Bin to fill into 
+      @param y   Bin to fill into 
+      @param phi The angle @f$ \varphi@f$ in radians */
+  virtual Bool_t AddToHarmonic(Double_t x, Double_t y, Double_t phi);
+  /** Process a full event. 
+      @param phis List of @a n @f$ \varphi=[0,2\pi]@f$ angles 
+      @param xs   List of @a n @f$ x@f$ values. 
+      @param ys   List of @a n @f$ y@f$ values. 
+      @param ws   Weights
+      @param n    Size of @a phis and @a xs */
+  virtual void Event(Double_t* phis, Double_t* xs, Double_t* ys, 
+                    Double_t* ws, ULong_t n);
+  /** Called at the end of a run */
+  virtual void Finish();
+  /** Get the bin at @a x 
+      @param x The bin value to find a flow object for. 
+      @param y The bin value to find a flow object for. 
+      @return The flow object at @a x,y or 0 if out of range */ 
+  virtual AliFMDFlowBin* GetBin(Double_t x, Double_t y) const;
+  /** Get the bin at @a x 
+      @param i The bin index to find a flow object for. 
+      @param j The bin index to find a flow object for. 
+      @return The flow object at @a i,j or 0 if out of range */ 
+  virtual AliFMDFlowBin* GetBin(UShort_t i, UShort_t j) const;
+  /** Whether this is to be considered a folder */
+  Bool_t IsFolder() const { return kTRUE; }
+  /** Browse this object */ 
+  void Browse(TBrowser* b);
+protected:
+  /** X axis */ 
+  AliFMDFlowAxis fXAxis;
+  /** Y axis */ 
+  AliFMDFlowAxis fYAxis;
+  /** Array of the flow objects */ 
+  AliFMDFlowBin** fBins;
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowBinned2D,1);
+};
+
+
+#endif
+//
+// EOF
+//
+
diff --git a/FMD/flow/AliFMDFlowEfficiency.cxx b/FMD/flow/AliFMDFlowEfficiency.cxx
new file mode 100644 (file)
index 0000000..150e772
--- /dev/null
@@ -0,0 +1,284 @@
+#include <flow/AliFMDFlowEfficiency.h>
+#include <cmath>
+#include <iostream>
+template <typename T>
+T sign(const T& x, const T& y) 
+{
+  return (y >= 0 ? fabs(x) : -fabs(x));
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDFlowEfficiency::LnGamma(Double_t z)
+{
+  if (z <= 0) return 0;
+  Double_t c[] = { 2.5066282746310005, 76.18009172947146, 
+                  -86.50532032941677, 24.01409824083091,  
+                  -1.231739572450155, 0.1208650973866179e-2, 
+                  -0.5395239384953e-5};
+  Double_t x = z;
+  Double_t y = x;
+  Double_t t = x + 5.5;
+  t        = (x+.5 * log(t) - t);
+  Double_t s = 1.000000000190015;
+  for (UInt_t i = 1; i < 7; i++) { 
+    y  += 1;
+    s  += c[i] / y;
+  }
+  Double_t v = t + log(c[0]*s/x);
+  return v;
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDFlowEfficiency::BetaCf(Double_t x, Double_t a, Double_t b)
+{
+  Int_t    itmax = 500;   // Maximum # of iterations
+  Double_t eps   = 3e-14; // Precision. 
+  Double_t fpmin = 1e-30;
+  
+  Double_t qab = a + b;
+  Double_t qap = a + 1;
+  Double_t qam = a - 1;
+  Double_t c   = 1;
+  Double_t d   = 1 - qab * x  / qap;
+  if (fabs(d) < fpmin) d = fpmin;
+  d          = 1/d;
+  Double_t h   = d;
+  Int_t m;
+  for (m = 1; m <= itmax; m++) { 
+    Int_t    m2 =  m * 2;
+    Double_t aa =  m * (b - m) * x / ((qam + m2) * (a + m2));
+    d         =  1 + aa * d;
+    if (fabs(d) < fpmin) d = fpmin;
+    c         =  1 + aa / c;
+    if (fabs(c) < fpmin) d = fpmin;
+    d         =  1 / d;
+    h         *= d * c;
+    aa        =  -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
+    d         = 1 + aa * d;
+    if (fabs(d) < fpmin) d = fpmin;
+    c         =  1 + aa / c;
+    if (fabs(c) < fpmin) d = fpmin;
+    d         =  1 / d;
+    Double_t dd = d * c;
+    h         *= dd;
+    if (fabs(dd - 1) <= eps) break;
+  }
+  if (m > itmax) { 
+    std::cerr << "a or b too big, or max # iterations too small"
+             << std::endl;
+  }
+  return h;
+}
+
+  
+
+//____________________________________________________________________
+Double_t
+AliFMDFlowEfficiency::IBetaI(Double_t a, Double_t b, Double_t x)
+{
+  Double_t bt = 0;
+  if (x < 0 || x > 1) { 
+    std::cerr << "x not in [0,1]" << std::endl;
+    return 0;
+  }
+  if (x == 0 || x == 1) bt = 0;
+  else 
+    bt = (exp(LnGamma(a + b) - LnGamma(a) - LnGamma(b) 
+             + a * log(x) + b * log(1 - x)));
+  if (x < (a + 1) / (a + b + 2)) 
+    return bt * BetaCf(x, a, b) / a;
+  return 1 - bt * BetaCf(1-x, b, a)/ b;
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDFlowEfficiency::BetaAB(Double_t a, Double_t b, Int_t k, Int_t n)
+{
+  if (a == b) return 0;
+  Int_t c1 = k + 1;
+  Int_t c2 = n - k + 1;
+  return IBetaI(c1, c2, b) - IBetaI(c1, c2, a);
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDFlowEfficiency::SearchUpper(Double_t low, Int_t k, Int_t n, Double_t c)
+{
+  Double_t integral = BetaAB(low, 1, k, n);
+  if (integral == c) return 1; // lucky -- this is the solution
+  if (integral <  c) return -1;// no solution exists
+  Double_t too_high = 1;         // upper edge estimate
+  Double_t too_low  = low;
+  Double_t test     = 0;
+
+  // Use a bracket-and-bisect search.  Loop 20 times, to end up with a
+  // root guaranteed accurate to better than 0.1%.  
+  for (UInt_t i = 0; i < 20; i++) { 
+    test     = 0.5 * (too_low + too_high);
+    integral = BetaAB(low, test, k, n);
+    if (integral > c) too_high = test;
+    else              too_low  = test;
+  }
+  return test;
+}
+//____________________________________________________________________
+Double_t
+AliFMDFlowEfficiency::SearchLower(Double_t high, Int_t k, Int_t n, Double_t c)
+{
+  Double_t integral = BetaAB(0, high, k, n);
+  if (integral == c) return 0; // lucky -- this is the solution
+  if (integral <  c) return -1;// no solution exists
+  Double_t too_low  = 0;         // lower edge estimate
+  Double_t too_high = high;
+  Double_t test     = 0;
+
+  // Use a bracket-and-bisect search.  Loop 20 times, to end up with a
+  // root guaranteed accurate to better than 0.1%.  
+  for (UInt_t i = 0; i < 20; i++) { 
+    test     = 0.5 * (too_high + too_low);
+    integral = BetaAB(test, high, k, n);
+    if (integral > c) too_low  = test;
+    else              too_high = test;
+  }
+  return test;
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowEfficiency::Brent(Double_t ax, Double_t bx, Double_t cx, 
+                           Double_t& xmin, Int_t k, Int_t n, Double_t c, 
+                           Double_t tol)
+{
+  const Int_t    itmax = 100;
+  const Double_t gold  = 0.3819660;
+  const Double_t eps   = 1e-10;
+  
+  Int_t iter;
+  Double_t a  = (ax < cx ? ax : cx);
+  Double_t b  = (ax > cx ? ax : cx);
+  Double_t d  = 0;
+  Double_t e  = 0;
+  Double_t x  = bx;
+  Double_t w  = bx;
+  Double_t v  = bx;
+  Double_t fw = Length(x, k, n, c);
+  Double_t fv = fw;
+  Double_t fx = fv;
+  
+  for (iter = 1; iter <= itmax; iter++) { 
+    Double_t xm   = .5 * (a + b);
+    Double_t tol1 = tol * fabs(x) + eps;
+    Double_t tol2 = 2 * tol1;
+    if (fabs(x - xm) < (tol2 - .5 * (b-a))) {
+      xmin = x;
+      return fx;
+    }
+    if (fabs(e) > tol1) { 
+      Double_t r = (x - w) * (fx - fv);
+      Double_t q = (x - v) * (fx - fw);
+      Double_t p = (x - v) * q - (x - w) * r;
+      q        = 2 * (q - r);
+      if (q > 0) p *= -1;
+      q        = fabs(q);
+      Double_t t = e;
+      e        = d;
+      if (fabs(p) >= fabs(.5 * q * t) || 
+         p <= q * (a - x)            || 
+         p >= q * (b - x)) {
+       e = (x > xm ? a - x : b - x);
+       d = gold * e;
+      }
+      else {
+       d        = p / q;
+       Double_t u = x + d;
+       if (u - a < tol2 || b - u < tol2) 
+         d = sign(tol1, xm - x);
+      }
+    }
+    else { 
+      e = (x >= xm ? a - x : b - x);
+      d = gold * e;
+    }
+    Double_t u  = (fabs(d) >= tol1 ? x + d : x + sign(tol1, d));
+    Double_t fu = Length(u, k, n, c);
+    if (fu <= fx) { 
+      if (u >= x) a = x;
+      else        b = x;
+      v  = w;
+      w  = x;
+      x  = u;
+      fv = fw;
+      fw = fx;
+      fx = fu;
+    }
+    else { 
+      if (u < x) a = u;
+      else       b = u;
+      if (fu <= fw || w == x) { 
+       v  = w;
+       w  = u;
+       fv = fw;
+       fw = fu;
+      }
+      else if (fu <= fv || v == x || v == w) { 
+       v  = u;
+       fv = fu;
+      }
+    }
+  }
+  std::cerr << "Brett: too many iterations" << std::endl;
+  xmin = x;
+  return fx;
+}
+  
+//____________________________________________________________________
+Double_t 
+AliFMDFlowEfficiency::Length(Double_t l, Int_t k, Int_t n, Double_t c)
+{
+  Double_t h = SearchUpper(l, k, n, c);
+  if (h == -1) return 2;
+  return (h-l);
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowEfficiency::Interval(Int_t k, Int_t n, 
+                              Double_t& low, Double_t& high, 
+                              Double_t c)
+{
+  if (n == 0) { 
+    low  = 0;
+    high = 1;
+    return .5;
+  }
+  
+  Double_t e = Double_t(k) / n;
+  Double_t l, h;
+  
+  if (k == 0) { 
+    l = 0;
+    h = SearchUpper(l, k, n, c);
+  }
+  else if (k == n) { 
+    h = 1;
+    l = SearchLower(h, k, n, c);
+  }
+  else { 
+    Brent(0, .5, 1, l, n, k, 1e-9, c);
+    h = l + Length(l, n, k, c);
+  }
+  low  = l;
+  high = h;
+  return e;
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+    
+  
+  
diff --git a/FMD/flow/AliFMDFlowEfficiency.h b/FMD/flow/AliFMDFlowEfficiency.h
new file mode 100644 (file)
index 0000000..59bd9f6
--- /dev/null
@@ -0,0 +1,128 @@
+// -*- mode: C++ -*-
+#ifndef Efficiency_h
+#define Efficiency_h
+#include <Rtypes.h>
+
+/** @defgroup z_eff Efficiency calculations 
+    @brief Functions to do efficiency calculations based on a
+    Baysian analysis. 
+*/
+/** Namespace for efficency calculations 
+    @ingroup z_eff */
+namespace AliFMDFlowEfficiency 
+{
+  /** @{ 
+      @ingroup z_eff */
+  /** Calculate @f$ \log(\Gamma(z))\ \forall z>0@f$ 
+      @param z Argument. 
+      @return  @f$ \log(\Gamma(z))@f$ */ 
+  Double_t LnGamma(Double_t z);
+  
+  /** Continued fraction evaluation by modified Lentz's method 
+      used in calculation of incomplete  Beta function. 
+      @param x argument. 
+      @param a lower limit
+      @param b upper limit
+      @return incomplete Beta function evaluated at x */ 
+  Double_t BetaCf(Double_t x, Double_t a, Double_t b);
+  
+  /** Calculates the incomplete Beta function @f$ I_x(a,b)@f$;
+      this is the incomplete Beta function divided by the
+      complete Beta function. 
+      @param a Lower bound 
+      @param b Upper bound 
+      @param x Order 
+      @return  @f$ I_x(a,b)@f$ */
+  Double_t IBetaI(Double_t a, Double_t b, Double_t x);
+
+  /** Calculates the fraction of the area under the curve 
+      @f$ x^k (1-x)^{n-k}@f$ between @f$ x=a@f$ and @f$ x=b@f$ 
+      @param a lower limit 
+      @param b upper limit 
+      @param k Parameter @f$ k@f$ 
+      @param n Parameter @f$ n@f$ 
+      @return The fraction under the curve */ 
+  Double_t BetaAB(Double_t a, Double_t b, Int_t k, Int_t n);
+  
+  /** Integrates the Binomial distribution with parameters @a k and @a
+      n, and determines the upper edge of the integration region, 
+      starting at @a low, which contains probability content @a c.  If
+      an upper limit is found, the value is returned. If no solution
+      is found, -1 is returned. Check to see if there is any solution
+      by verifying that the integral up to the maximum upper limit (1)
+      is greater than c 
+      @param low Where to start the integration from. 
+      @param k   @a k parameter of the Binomial distribution. 
+      @param n   @a N parameter of the Binomial distribution. 
+      @param c   Wanted confidence limit (defaults to 68% - similar to
+      @f$ 1\sigma@f$ of a Gaussian distribution)
+      @return the upper limit of the confidence interval, or -1 in
+      case of failure */
+  Double_t SearchUpper(Double_t low, Int_t k, Int_t n, Double_t c=0.683);
+
+  /** Integrates the Binomial distribution with parameters @a k and @a
+      n, and determines the lower edge of the integration region, 
+      ending at @a high, which contains probability content @a c.  If
+      a lower limit is found, the value is returned. If no solution
+      is found, -1 is returned. Check to see if there is any solution
+      by verifying that the integral up to the maximum upper limit (1)
+      is greater than c 
+      @param high Where to end the integration at. 
+      @param k    @a k parameter of the Binomial distribution. 
+      @param n    @a N parameter of the Binomial distribution. 
+      @param c    Wanted confidence limit (defaults to 68% - similar to
+      @f$ 1\sigma@f$ of a Gaussian distribution) 
+      @return the upper limit of the confidence interval, or -1 in
+      case of failure */
+  Double_t SearchLower(Double_t high, Int_t k, Int_t n, Double_t c=0.683);
+
+  /** Numerical equation solver.  This includes root finding and
+      minimum finding algorithms. Adapted from Numerical Recipes in
+      C, 2nd edition. Translated to C++ by Marc Paterno
+      @param ax   Left side of interval 
+      @param bx   Middle of interval 
+      @param cx   Right side of interval 
+      @param tol  Tolerance 
+      @param xmin On return, the value of @f$ x@f$ such that @f$
+      f(x)@f$ i minimal. 
+      @param k    Parameter of the Binomial distribution 
+      @param n    Parameter of the Binomial distribution 
+      @param c    Confidence level. 
+      @return Minimum of @f$ f = f(x)@f$. */
+  Double_t Brent(Double_t ax, Double_t bx, Double_t cx, Double_t& xmin, 
+              Int_t k, Int_t n, Double_t c=.683, Double_t tol=1e-9);
+  
+  /** Return the length of the interval starting at @a l that contains
+      @a c of the @f$ x^k (1-x)^{n-k}@f$ distribution. If there is no
+      sufficient interval starting at @a l, we return 2.0 
+      @param l Lower bound 
+      @param k Binomial parameter k 
+      @param n Binomial parameter n 
+      @param c Condifience level 
+      @return Legnth of interval */
+  Double_t Length(Double_t l, Int_t k, Int_t n, Double_t c=0.683);
+    
+
+  /** Calculate the shortest central confidence interval containing
+      the required probability content @a c. Interval(low) returns the
+      length of the interval starting at low that contains @a c
+      probability. We use Brent's method, except in two special cases:
+      when @a k=0, or when @a k=n 
+      @author Marc Paterno
+      @param k    Binomial parameter @a k 
+      @param n    Binomial parameter @a n 
+      @param low  On return, the lower limit
+      @param high On return, the upper limit 
+      @param c    Required confidence level (defaults to 68% - similar
+      to @f$ 1\sigma@f$ of a Gaussian distribution)
+      @return The mode */ 
+  Double_t Interval(Int_t k, Int_t n, Double_t& low, Double_t& high, Double_t c=0.683);
+  /** @} */
+}
+
+
+#endif
+//
+// EOF
+//
+
diff --git a/FMD/flow/AliFMDFlowEventPlane.cxx b/FMD/flow/AliFMDFlowEventPlane.cxx
new file mode 100644 (file)
index 0000000..bfdecb6
--- /dev/null
@@ -0,0 +1,77 @@
+/** @file 
+    @brief Implementation of an EventPlane class */
+#include "flow/AliFMDFlowEventPlane.h"
+#include "flow/AliFMDFlowUtil.h"
+// #include <cmath>
+#ifndef _GNU_SOURCE
+extern "C" 
+{
+  /** Function to caculate @f$ \sin(a), \cos(a)@f$ in one go.  Note,
+      that with GCC, this is a built-in, and we only have this to
+      serve as a replacement function 
+      @ingroup utils */
+  inline void sincos(Double_t a, Double_t* sina, Double_t* cosa) {
+    *sina = sin(a);
+    *cosa = cos(a);
+  }
+}
+#endif
+
+//====================================================================
+void 
+AliFMDFlowEventPlane::Clear(Option_t*) 
+{ 
+  fSumSinMPhi = 0;
+  fSumCosMPhi = 0;
+  fCache      = -1;
+}
+//____________________________________________________________________
+void 
+AliFMDFlowEventPlane::Add(Double_t phi, Double_t weight) 
+{ 
+  Double_t a = NormalizeAngle(fOrder * phi);
+  Double_t s, c;
+  sincos(a, &s, &c);
+  if (isnan(s) || isinf(s) || isnan(c) || isinf(s)) return;
+  fSumSinMPhi += weight * s;
+  fSumCosMPhi += weight * c;
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowEventPlane::Psi() const  
+{ 
+  if (fCache < 0) fCache = DoPsi(fSumSinMPhi, fSumCosMPhi);
+  return fCache;
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowEventPlane::Psi(Double_t phi, Double_t w) const  
+{ 
+  Double_t a = NormalizeAngle(fOrder * phi);
+  Double_t s, c;
+  sincos(a, &s, &c);
+  if (isnan(s) || isinf(s) || isnan(c) || isinf(s)) return Psi();
+  Double_t psi = DoPsi(fSumSinMPhi - w * s, fSumCosMPhi - w * c);
+  return psi;
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowEventPlane::DoPsi(Double_t sumsin, Double_t sumcos) const
+{
+  Double_t psi = 0;
+  // Make sure we get an angle everywhere 
+  if      (sumcos != 0) psi =  atan2(sumsin, sumcos);
+  else if (sumsin == 0) psi =  0;
+  else if (sumsin >  0) psi =  M_PI / 2;
+  else                  psi = -M_PI / 2;
+  psi =  NormalizeAngle(psi);
+  psi /= fOrder;
+  return psi;
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
diff --git a/FMD/flow/AliFMDFlowEventPlane.h b/FMD/flow/AliFMDFlowEventPlane.h
new file mode 100644 (file)
index 0000000..8258eab
--- /dev/null
@@ -0,0 +1,79 @@
+// -*- C++ -*- 
+/** @file 
+    @brief Declaration of an EventPlane class */
+#ifndef FLOW_EVENTPLANE_H
+#define FLOW_EVENTPLANE_H
+#include <TObject.h>
+
+//______________________________________________________
+/** @class AliFMDFlowEventPlane flow/AliFMDFlowEventPlane.h <flow/AliFMDFlowEventPlane.h>
+    @brief Class to determine the event plane 
+
+    The event plane is calculated as 
+    @f[ 
+    \Psi_n = \frac1n\tan^{-1}\left[\frac{\sum_i(w_i\sin(n\varphi_i))}
+    {\sum_i(w_i\cos(n\varphi_i))}\right]
+    @f]
+    where @f$ i @f$ runs over all observations of @f$\varphi@f$ in an
+    event, and @f$ w_i@f$ is the weight of the @f$ i@f$ observation of
+    @f$ \varphi@f$
+
+    @ingroup a_basic
+*/
+class AliFMDFlowEventPlane : public TObject
+{
+public:
+  /** Constructor 
+      @param m Harmonic number */
+  AliFMDFlowEventPlane(UShort_t m) 
+    : fSumSinMPhi(0), 
+      fSumCosMPhi(0),
+      fOrder(m),
+      fCache(0)
+  { Clear(); }
+  /** Destructor */
+  ~AliFMDFlowEventPlane() {} 
+  /** Clear it */
+  void Clear(Option_t* option="");
+  /** Add a data point 
+      @param phi The angle @f$\varphi\in[0,2\pi]@f$
+      @param weight The weight */
+  void Add(Double_t phi, Double_t weight=1);
+  /** Get the event plane 
+      @return @f$ \Psi_k@f$ */
+  Double_t Psi() const;
+  /** Get the event plane angle @f$ \Psi_k@f$ @e disregarding the
+      contribution from the observation @f$ \varphi_i@f$ with weight
+      @f$ w_i@f$.  This is to avoid auto-correlations 
+      @param phi The observation  @f$ \varphi_i@f$
+      @param w   The weight @f$ w_i@f$ of the obervation. 
+      @return The event plane angle @f$ \Psi_k@f$ with out the
+      contribution from @f$ \varphi_i@f$ */
+  Double_t Psi(Double_t phi, Double_t w=1) const;
+  /** Get the harmnic order 
+      @return @f$ k@f$  */
+  UShort_t Order() const { return fOrder; }
+protected:
+  /** Utility function to calculate @f$ \Psi@f$ from the sum of
+      sines and cosines. 
+      @param sumsin Sum of sines 
+      @param sumcos Sum of cosine. 
+      @return @f$ \Psi@f$ */
+  Double_t DoPsi(Double_t sumsin, Double_t sumcos) const;
+  /** @f$ \sum_i w_i \sin(k \varphi_i)@f$ */
+  Double_t fSumSinMPhi;
+  /** @f$ \sum_i w_i \cos(k \varphi_i)@f$ */
+  Double_t fSumCosMPhi;
+  /** Order */
+  UShort_t fOrder;
+  /** Cache of Psi */
+  mutable Double_t fCache;
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowEventPlane,1); 
+};
+
+
+#endif
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowHarmonic.cxx b/FMD/flow/AliFMDFlowHarmonic.cxx
new file mode 100644 (file)
index 0000000..68c3c0d
--- /dev/null
@@ -0,0 +1,52 @@
+/** @file 
+    @brief Implementation of a Harmonic class */
+#include "flow/AliFMDFlowHarmonic.h"
+#include "flow/AliFMDFlowUtil.h"
+// #include <cmath>
+
+//====================================================================
+void 
+AliFMDFlowHarmonic::Add(Double_t phi, Double_t psi) 
+{ 
+  Double_t a       = NormalizeAngle(fOrder * (phi - psi));
+  Double_t contrib = cos(a);
+  AliFMDFlowStat::Add(contrib);
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowHarmonic::Value(Double_t r, Double_t er2, Double_t& e2) const 
+{
+  // The corrected value is given by 
+  //
+  //          v_n^obs
+  //    v_n = -------
+  //             R
+  // 
+  // where 
+  // 
+  //              1
+  //    v_n^obs = - \sum_i(cos(n(\phi_i - \Psi)))
+  //              N 
+  // 
+  // and R is the resolution 
+  // 
+  // The error on the corrected value is given by 
+  //
+  //                dv_n                    dv_n
+  //    d^2v_n = (--------)^2 d^2v_n^obs + (----)^2 d^2R
+  //              dv_n^obs                   dR
+  //
+  //             d^2v_n^obs R^2 + d^2R v_n^obs^2
+  //           = -------------------------------
+  //                         R^4 
+  // 
+  Double_t a = fAverage;
+  Double_t v = a / r;
+  Double_t s = fSqVar / fN;
+  e2       = (s * r * r + er2 * a * a) / pow(r, 4);
+  return v;
+}
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowHarmonic.h b/FMD/flow/AliFMDFlowHarmonic.h
new file mode 100644 (file)
index 0000000..967651c
--- /dev/null
@@ -0,0 +1,65 @@
+// -*- mode: C++ -*-
+/** @file 
+    @brief Declaration of a Harmonic class */
+#ifndef FLOW_HARMONIC_H
+#define FLOW_HARMONIC_H
+#include <flow/AliFMDFlowStat.h>
+
+/** @defgroup a_basic Basic classes for doing Flow analysis. 
+    @brief This group of class handles the low-level stuff to do
+    flow analysis. */
+//______________________________________________________
+/** @class AliFMDFlowHarmonic flow/AliFMDFlowHarmonic.h <flow/AliFMDFlowHarmonic.h>
+    @brief Calculate the @f$ n^{th}@f$ order harmonic
+    @ingroup a_basic 
+
+    Calculate the @f$ n^{th}@f$ order harmonic, given by 
+    @f{eqnarray*}
+    v_n       &=& \frac{v_n^{obs}}{R}\\
+    v_n^{obs} &=& \frac1M\sum_i^M\cos(n (\varphi_i - \Psi))
+    @f}
+    where @f$ R@f$ is the resolution, @f$ i@f$ runs over all @f$
+    M@f$ observations of @f$ \varphi_i@f$ in all events, and
+    @f$\Psi@f$ is the estimated event plane. 
+
+    The error on the corrected value is given by 
+    @f{eqnarray*} 
+    \delta^2v_n & = & \left(\frac{dv_n}{dv_n^{obs}}\right)^2\delta^2
+    v_n^{obs} + \left(\frac{dv_n}{dR}\right)^2\delta^2R \\ 
+    & = & \frac{\delta^2v_n^{obs} R^2 + \delta^2R (v_n^{obs})^2}
+    {R^4}
+    @f}
+*/    
+class AliFMDFlowHarmonic : public AliFMDFlowStat
+{
+public:
+  /** Constructor 
+      @param n Order of the harmonic */
+  AliFMDFlowHarmonic(UShort_t n) : fOrder(n) {} 
+  /** Destructor */ 
+  virtual ~AliFMDFlowHarmonic() {}
+  /** Add a data point 
+      @param phi The absolute angle @f$ \varphi \in[0,2\pi]@f$ 
+      @param psi The event plane angle @f$ \Psi \in[0,2\pi] @f$ */
+  void Add(Double_t phi, Double_t psi);
+  /** Get the harmonic. 
+      @param r   Event plane resolution 
+      @param er2 Square error on event plane resolution 
+      @param e2  On return the square error 
+      @return The harmonic value */
+  Double_t Value(Double_t r, Double_t er2, Double_t& e2) const;
+  /** Get the order of the harmonic */
+  UShort_t Order() const { return fOrder; }
+protected:
+  /** the order */ 
+  UShort_t fOrder;
+  /** Define for ROOT I/O */ 
+  ClassDef(AliFMDFlowHarmonic,1);
+};
+
+
+#endif
+//
+// EOF
+//
+
diff --git a/FMD/flow/AliFMDFlowResolution.cxx b/FMD/flow/AliFMDFlowResolution.cxx
new file mode 100644 (file)
index 0000000..57fdc0e
--- /dev/null
@@ -0,0 +1,266 @@
+/** @file 
+    @brief Implementation of an Resolution class */
+#include "flow/AliFMDFlowResolution.h"
+#include "flow/AliFMDFlowUtil.h"
+#include "flow/AliFMDFlowBessel.h"
+#include <iostream>
+//#include <cmath>
+
+//====================================================================
+void 
+AliFMDFlowResolution::Add(Double_t psiA, Double_t psiB) 
+{ 
+  Double_t diff    = NormalizeAngle(fOrder * (psiA - psiB));
+  Double_t contrib = cos(diff);
+  AliFMDFlowStat::Add(contrib);
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowResolution::Correction(UShort_t k) const 
+{ 
+  Double_t e;
+  return Correction(k, e); 
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowResolution::Correction(UShort_t, Double_t& e2) const 
+{ 
+  e2 = fSqVar / fN;
+  return sqrt(2) * sqrt(fabs(fAverage));
+}
+
+//====================================================================
+Double_t 
+AliFMDFlowResolutionStar::Correction(UShort_t k, Double_t& e2) const 
+{ 
+  if (k > 4) return 0;
+  Double_t delta = 0;
+  Double_t chi   = Chi(fAverage, k, delta);
+  Double_t dr    = 0;
+  Double_t res   = Res(sqrt(2) * chi, k, dr);
+  e2           = pow(dr * delta,2);
+  return res;
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowResolutionStar::Correction(UShort_t k) const 
+{ 
+  Double_t e;
+  return Correction(k, e); 
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k, 
+                             Double_t& delta) const 
+{
+  delta        = 1;
+  Double_t chi   = 2;
+  Double_t dr    = 0;
+  for (UInt_t i = 0; i < 15; i++) { 
+    if (Res(chi, k, dr) < res) chi += delta;
+    else                       chi -= delta;
+    delta /= 2;
+  }
+  return chi;
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, 
+                                   Double_t& dr) const 
+{ 
+  // The resolution function is 
+  // 
+  //          sqrt(pi) x exp(-x/4) (f1(x^2/4) + f2(x^2/4))
+  //   r(x) = --------------------------------------------
+  //                          2 sqrt(2) 
+  // 
+  //        
+  //        = c x (f1(y) - f2(y))
+  //
+  // where f1 is the modified Bessel function first kind I_{(k-1)/2}, 
+  // and f2 is the modified Bessel function of the first kind
+  // I_{(k+1)/2}, and 
+  // 
+  //          sqrt(pi) exp(-x^2/4) 
+  //      c = --------------------,   y = x^2/4
+  //              2 sqrt(2)
+  // 
+  // The derivative of the resolution function is 
+  //
+  //            c 
+  //    r'(y) = - (4 sqrt(y) (f1'(y) - f2'(y)) - (4 y - 2)(f1(y) - f2(y)))
+  //            2
+  // 
+  //            c                                    r(y)   
+  //          = - (4 sqrt(y) (f1'(y) - f2'(y))) + --------- - sqrt(y) r(y)
+  //           2                                 2 sqrt(y)       
+  // 
+  // Since dI_n(x)/dx = I_(n-1)(x) - n / x I_n(x), and substituting 
+  // f3(y) = I_((k-3)/2)(y) 
+  // 
+  //            c  
+  //    r'(y) = - ((4 - 2 k) f1(y) - (4 y + 2 k) f2(y) + 4 y f3(y))
+  //            2   
+  // 
+  Double_t y  = chi * chi / 4;
+  Double_t c  = sqrt(M_PI) * exp(-y) / (2 * sqrt(2));   
+  
+  // i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
+  Double_t i[3], di[3];
+  AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
+  
+  Double_t r  = c * chi * (i[2] + i[1]);
+  dr        = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]);
+
+  return r;  
+}
+
+//====================================================================
+void 
+AliFMDFlowResolutionTDR::Clear() 
+{
+  fN = 0;
+  fLarge = 0;
+}
+//____________________________________________________________________
+void 
+AliFMDFlowResolutionTDR::Add(Double_t psi_a, Double_t psi_b)
+{ 
+  Double_t a = fabs(psi_a - psi_b);
+  if (a >= .5 * M_PI) fLarge++;
+  fN++;
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowResolutionTDR::Correction(UShort_t k, Double_t& e2) const 
+{ 
+  // From nucl-ex/9711003v2 
+  //
+  //
+  //   <cos n Delta phi> = 
+  //
+  //         sqrt(pi)
+  //         -------- chi exp(-z) (I_((n-1)/2)(z) + I_((n+1)/2)(z))
+  //            2
+  // 
+  // where z = chi^2 / 2
+  //
+  Double_t echi2 = 0;
+  Double_t y     = Chi2Over2(echi2);
+  Double_t chi   = sqrt(2 * y);
+  Double_t c     = sqrt(M_PI) * exp(-y) / 2;
+  
+  // i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
+  Double_t i[3], di[3];
+  AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
+  
+  Double_t r  = c * chi * (i[2] + i[1]);
+  Double_t dr = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]);
+  e2        = dr * dr * echi2;
+  return r;
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowResolutionTDR::Correction(UShort_t k) const 
+{ 
+  Double_t e;
+  return Correction(k, e); 
+}
+//____________________________________________________________________
+Double_t 
+AliFMDFlowResolutionTDR::Chi2Over2(Double_t& e2) const 
+{
+  // From nucl-ex/9711003v2 
+  //
+  // N equal to the total number of events, and k equal to the
+  // number of events where |Psi_A - Psi_B| > pi/2, we get 
+  //
+  //   k   exp(-chi^2/2)                               k
+  //   - = -------------     =>    chi^2 / 2 = - log 2 -
+  //   N          2                                   N
+  //
+  //                                               N
+  //                         =>    chi^2 / 2 = log --   (*)
+  //                                               k
+  //            
+  //                                                   2 k
+  //                         =>    chi^2     = - 2 log ---
+  //                                                    N
+  //                                                            
+  //                         =>    chi       = -/+ sgrt (- 2 log (2k/N))
+  //
+  // (*) this is what is returned.  
+  //
+  //  The error on chi is given by 
+  //
+  //              dchi                           1
+  //    d^2chi = (------)^2 d^2(k/N) = - -------------------- d^2(k/N)
+  //              d(k/N)                 4 (k/N)^2 log(2 k/N)
+  // 
+  // where 
+  //                      k         k^2
+  //               (1 - 2 -) dk^2 + --- dN^2
+  //                      N        N^2
+  //   d^2(k/N) = --------------------------
+  //                           N^2 
+  // 
+  // 
+  // with dk = sqrt(k) and dN = sqrt(N), this reduces to 
+  //
+  //                     k      k^2        N k     k^2   k^2
+  //              (1 - 2 -) k + ----- N    --- - 2 --- + ---
+  //                     N      N^2         N       N     N
+  //   d^2(k/N) = --------------------- = --------------------------
+  //                           N^2                 N^2 
+  //
+  //              k N - k^2   k (N - k)   k (1 - k / N)
+  //            = --------- = - ------- = - -----------
+  //                 N^3      N   N^2     N      N
+  //
+  // Identifying r = k/N, we get 
+  //
+  //            
+  //   d(k/N)   = sqrt(r (1 - r) / N) 
+  //
+  // which is the well-known result for Binomial errors.
+  // Alternatively, one could compute these error using an Baysian
+  // confidence level of say 68.3% (see for example 
+  // http://home.fnal.gov/~paterno/images/effic.pdf). 
+  // 
+  // Our final expression for the error on chi then becomes 
+  //
+  //                      1         
+  //   d^2chi = - -------------------- r (1 - r) / N
+  //             4 r^2 log(2 r)
+  //
+  //                 1 - r            r - 1
+  //          = - -------------- = ------------
+  //              4 r N log(2 r)   4 k log(2 r)  
+
+  if (fLarge == 0) { 
+    std::cerr << "TDR: Large = 0" << std::endl;
+    return -1;
+  }
+  Double_t ratio = Double_t(fN) / (2 * fLarge);
+  if (ratio <= 0) {
+    std::cerr << "TDR: Large/All = " << ratio << " <= 0!" << std::endl;
+    return -1;
+  }
+  Double_t chi2over2  = log(ratio);
+  if (chi2over2 < 0) { 
+    std::cerr << "TDR: log(" << ratio << ") = " << chi2over2 
+             << " < 0" << std::endl; 
+    return -1;
+  }
+  Double_t r = Double_t(fLarge) / fN;
+  e2 = (r - 1) / (4 * fLarge * log(2 * r));
+  return chi2over2;
+}
+
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowResolution.h b/FMD/flow/AliFMDFlowResolution.h
new file mode 100644 (file)
index 0000000..e3954ff
--- /dev/null
@@ -0,0 +1,238 @@
+// -*- mode: C++ -*-
+/** @file 
+    @brief Declaration of an Resolution class */
+#ifndef FLOW_RESOLUTION_H
+#define FLOW_RESOLUTION_H
+#include <flow/AliFMDFlowStat.h>
+
+//______________________________________________________
+/** @class AliFMDFlowResolution flow/AliFMDFlowResolution.h <flow/AliFMDFlowResolution.h>
+    @brief Class to calculate the event plane resolution based on two
+    sub-events 
+    @ingroup a_basic
+    
+    This class calculates the event plane resolution based on the
+    basic formulas given in Phys. Rev. @b C58, 1671.   That is, the
+    resolution is given by 
+    @f[ 
+    R_{k} = \langle\cos(km(\Psi_m-\Psi_R))\rangle 
+    @f] 
+    where @f$ \Psi_R@f$ is the unknown @e true event plane angle, and
+    @f$ n = km@f$
+
+    Using two random sub-events, @f$ A, B@f$ we get that 
+    @f[
+    \langle\cos(n(\Psi^A_m-\Psi^B_m))\rangle = 
+    \langle\cos(n(\Psi^A_m-\Psi_R))\rangle 
+    \langle\cos(n(\Psi^B_m-\Psi_R))\rangle 
+    @f]
+    If the sub-events are of equal size, and randomly chosen, then we
+    get that 
+    @f[ 
+    \langle\cos(n(\Psi^A_m-\Psi_R))\rangle = 
+    \sqrt{\langle\cos(n(\Psi^A_m-\Psi^B_m))}
+    @f]
+    and it follows that 
+    @f{eqnarray*}
+    \langle\cos(km(\Psi_m-\Psi_R))\rangle & = & 
+    \sqrt{2}\langle\cos(n(\Psi^A_m-\Psi_R))\rangle\\
+    \sqrt{2}\sqrt{\langle\cos(n(\Psi^A_m-\Psi^B_m))}
+    @f}
+
+    Hence, the event-plane resolution is simply the square root of the
+    scaled average distance between the two sub-events, multiplied by
+    @f$ \sqrt{s}@f$ 
+
+    The error is therefor @f$ \sqrt{s}@f$ times the variance of the
+    cosine of the distance between the two sub-events. 
+*/
+class AliFMDFlowResolution : public AliFMDFlowStat
+{
+public:
+  /** Constructor
+      @param n Harmonic order */
+  AliFMDFlowResolution(UShort_t n) : fOrder(n) {}
+  /** Destructor */
+  virtual ~AliFMDFlowResolution() {}
+  /** add data point 
+      @param psiA A sub-event plane angle @f$ \Psi_A \in[0,2\pi]@f$
+      @param psiB B sub-event plane angle @f$ \Psi_B \in[0,2\pi]@f$ */
+  virtual void Add(Double_t psiA, Double_t psiB);
+  /** Get the correction for harmonic strength of order @a k 
+      @param k  The harminic strenght order to get the correction for
+      @param e2 The square error on the correction 
+      @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
+  virtual Double_t Correction(UShort_t k, Double_t& e2) const;
+  /** Get the correction for harmonic strength of order @a k 
+      @param k The harminic strenght order to get the correction for
+      @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
+  virtual Double_t Correction(UShort_t k) const;
+  /** Get the harmnic order */
+  UShort_t Order() const { return fOrder; }
+protected:
+  /** Order */
+  UShort_t fOrder;
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowResolution,1);
+};
+
+//______________________________________________________
+/** 
+    @ingroup a_basic
+    The event plane angle resolution function is 
+    @f[
+    R_k(\chi) = \frac{\pi} \chi e^{-\chi^2/4}\left( 
+    I_{\frac{k-1}2}(\chi^2/4) + I_{\frac{k+1}2}(\chi^2/4)\right)
+    @f]
+    Where @f$ I_n(x)@f$ is the modified Bessel function of the first
+    kind.   Identifying 
+    @f[
+    y = \chi^2/4\quad C=\frac{\sqrt{\pi}e^{-y}}{2\sqrt{2}}
+    @f]
+    and the short hands @f$ f2(y) = I_{\frac{k-1}2}@f$ and @f$ f3(y) =
+    I_{\frac{k+1}2}(y)@f$ we can write this more compact as 
+    @f[ 
+    R_k(y) = C y (f2(y) + f3(y))
+    @f] 
+    The derivative of the resolution function is 
+    @f[ 
+    R_k'(y) = \frac{C}{2}\left[4\sqrt{y}\left(f2'(y)-f3'(y)\right) 
+    - (4 y - 2)\left(f2(y) + f3(y)\right)\right]
+    @f]
+    Since 
+    @f[
+    I_\nu'(x) = I_{\nu-1}(x) - \frac{\nu}{x} I_\nu(x)\quad,
+    @f]
+    and setting @f$ f1(y) = I_{\frac{k-3}2}(y)@f$, we get 
+    @f[ 
+    R_k'(y) = \frac{C}{2}\left[4y f1(y) + (4-2k) f2(y) - (4y+2k)
+    f3(y)\right]
+    @f]
+
+    In this class, the argument @f$ \chi@f$ is estimated by finding
+    the minima of @f$ R_k(\chi)@f$ near the average of
+    @f$\cos(n(\Psi_A-\Psi_B))@f$.  The error @f$ \delta\chi@f$ is
+    estimated as the largest step size in the minimisation. 
+
+    The total error on the correction is then 
+    @f[ 
+    \delta R_k = R_k'(\chi) \delta\chi
+    @f] 
+*/
+
+class AliFMDFlowResolutionStar : public AliFMDFlowResolution
+{
+public:
+  /** Constructor
+      @param n Harmonic order */
+  AliFMDFlowResolutionStar(UShort_t n) 
+    : AliFMDFlowResolution(n) {}
+  /** Destructor */
+  ~AliFMDFlowResolutionStar() {}
+  /** Get the correction for harmonic strength of order @a k 
+      @param k The harminic strenght order to get the correction for
+      @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
+  virtual Double_t Correction(UShort_t k) const;
+
+  /** Get the correction for harmonic strength of order @a k 
+      @param k  The harminic strenght order to get the correction for
+      @param e2 The square error on the correction 
+      @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
+  virtual Double_t Correction(UShort_t k, Double_t& e2) const;
+  /** Get @f$ \chi@f$ 
+      @param res  First shot at the resolution. 
+      @param k    Order 
+      @param delta On return, the last step size in @f$ \chi@f$ -
+      which is taken to be @f$ \delta\chi@f$  
+      @return @f$\chi@f$ */
+  virtual Double_t Chi(Double_t res, UShort_t k, Double_t& delta) const;
+protected:
+  /** Calculate resolution 
+      @param chi @f$ \chi@f$
+      @param k   Order factor 
+      @param dr  On return, the derivative of @f$ R(\chi)@f$
+      @return 
+      @f[ 
+      \frac{\sqrt{\pi/2}}{2}\chi e^{-\chi^2/4}
+      (I_{\frac{(k-1)}{2}}(\chi^2/4)+ I_{\frac{(k+1)}{2}}(\chi^2/4))
+      @f] */
+  Double_t Res(Double_t chi, UShort_t k, Double_t& dr) const;
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowResolutionStar,1);
+};
+
+//______________________________________________________
+/** 
+    @ingroup a_basic
+
+    For more on the event plane angle resolution function, please
+    refer to the class description of ResolutionStar. 
+
+    In this class @f$ \chi@f$ is calculated from the ratio @f$ k/N@f$
+    of events with @f$ |\Psi_A - \Psi_B| > \pi/2@f$ to the total
+    number of events. 
+
+    The pre-print @c nucl-ex/9711003v2 gives the formula 
+    @f[ 
+    \frac{k}{N} = \frac{e^{-\chi^2/2}}{2}
+    @f] 
+    for @f$ \chi@f$.  Note, that this differs from the @f$ \chi@f$
+    used in ResolutionStar by a factor of @f$ 1/sqrt{2}@f$. 
+    
+    We can isolate @f$ \chi = \mp\sqrt{-2\log(2k/N)}@f$ from the
+    equation above.  
+
+    Since @f$ r=k/N@f$ is obviously a efficiency-like ratio, we get
+    that error @f$ \delta r@f$ is given by Binomial errors 
+    @f[ 
+    \delta r =  \sqrt{r\frac{1 - r}{N}}\quad.
+    @f] 
+    The total error @f$ \delta\chi@f$ then becomes 
+    @f[ 
+    \delta^2\chi = \left(\frac{d\chi}{dr}\right)^2 \delta^2r 
+    = \frac{r - 1}{4 k log(2 r)}
+    @f]
+
+    The total error on the correction is, as in
+    ResolutionStar, then given by 
+    @f[ 
+    \delta R_k = R_k'(\chi) \delta\chi
+    @f]
+*/
+class AliFMDFlowResolutionTDR : public AliFMDFlowResolution
+{
+public:
+  /** Constructor
+      @param n Harmonic order */
+  AliFMDFlowResolutionTDR(UShort_t n) 
+    : AliFMDFlowResolution(n), fLarge(0) {}
+  /** DTOR */
+  ~AliFMDFlowResolutionTDR() {}
+  virtual void Clear();
+  /** add a data  point */
+  virtual void Add(Double_t psi_a, Double_t psi_b);
+  /** Get the correction for harmonic strength of order @a k 
+      @param k The harminic strenght order to get the correction for
+      @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
+  virtual Double_t Correction(UShort_t k) const;
+  /** Get the correction for harmonic strength of order @a k 
+      @param e2 The square error on the correction 
+      @param k  The harminic strenght order to get the correction for
+      @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
+  virtual Double_t Correction(UShort_t k, Double_t& e2) const;
+  /** Get @f$ \chi^2/2@f$ 
+      @param e2 The square error on the correction 
+      @return @f$ \chi^2/2@f$ */
+  virtual Double_t Chi2Over2(Double_t& e2) const;
+protected:
+  /** Number of events with large diviation */
+  ULong_t fLarge;
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowResolutionTDR,1);
+};
+
+
+#endif
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowStat.h b/FMD/flow/AliFMDFlowStat.h
new file mode 100644 (file)
index 0000000..60a0071
--- /dev/null
@@ -0,0 +1,82 @@
+// -*- mode: C++ -*- 
+#ifndef FLOW_STAT_H
+#define FLOW_STAT_H
+#include <cmath>
+#include <TObject.h>
+
+//______________________________________________________
+/** @class AliFMDFlowStat flow/AliFMDFlowStat.h <flow/AliFMDFlowStat.h>
+    @brief Single variable statistics. 
+    @ingroup u_utils 
+
+    This class calculates the single variable statistics of a random
+    variable @f$ x@f$. 
+    
+    The average @f$\bar{x}@f$ after @f$n@f$ observations is calculated as 
+    @f[ 
+    \bar{x}^{(n)} = \left\{\begin{array}{rl} 
+    x^{(n)} & n = 1\\
+    \left(1-\frac1n\right)\bar{x}^{(n-1)} + \frac{x^{(n)}}{n} & n > 1
+    \end{array}\right. 
+    @f]
+    and the square variance @f$s^2@f$ after @f$ n@f$ observations is
+    given by  
+    @f[ 
+    {s^{2}}^{(n)} = \left\{\begin{array}{rl}
+    0 & n = 1\\ 
+    \left(1-\frac1n\right){s^{2}}^{(n-1)} + 
+    \frac{\left(x^{(n)}-\bar{x}^{(n)}\right)^2}{n-1} & n > 1
+    \end{array}\right.
+    @f] 
+*/
+class AliFMDFlowStat
+{
+public:
+  /** Constructor */
+  AliFMDFlowStat() : fAverage(0), fSqVar(0), fN(0) {}
+  /** Destructor */
+  virtual ~AliFMDFlowStat() {}
+  /** Reset */
+  void Clear() { fAverage = fSqVar = 0; fN = 0; } 
+  /** Add another data point */
+  void Add(Double_t x);
+  /** Get the average */
+  Double_t Average() const { return fAverage; }
+  /** Get the square variance */
+  Double_t SqVar() const { return fSqVar; } 
+  /** Get the number of data points */
+  ULong_t N() const { return fN; }
+protected:
+  /** Average */
+  Double_t        fAverage;
+  /** Square variance */
+  Double_t        fSqVar;
+  /** Number of data points */
+  ULong_t fN;
+  /** Define for ROOT I/O */
+  ClassDef(AliFMDFlowStat,1);
+};
+
+//__________________________________________________________________
+inline void 
+AliFMDFlowStat::Add(Double_t x) 
+{ 
+  if (isnan(x) || isinf(x)) return;
+  fN++;
+  if (fN == 1) { 
+    fAverage = x;
+    fSqVar   = 0;
+    return;
+  }
+  Double_t c = (1. - 1. / fN);
+  fAverage  *= c;
+  fAverage  += x / fN;
+  fSqVar    *= c;
+  fSqVar    += pow(x - fAverage, 2) / (fN - 1);
+}
+
+
+#endif
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowTrue.cxx b/FMD/flow/AliFMDFlowTrue.cxx
new file mode 100644 (file)
index 0000000..d85072b
--- /dev/null
@@ -0,0 +1,112 @@
+#include "flow/AliFMDFlowTrue.h"
+#include "flow/AliFMDFlowUtil.h"
+#include <iostream>
+#include <iomanip>
+#include <TString.h>
+
+//====================================================================
+void
+AliFMDFlowTrueBin::End()
+{
+  Double_t psi  = fPsi.Psi();
+  Double_t dpsi = NormalizeAngle(fPsi.Order() * (psi-fPsiR));
+  fResReal.Add(cos(dpsi));
+}
+
+//____________________________________________________________________
+Double_t 
+AliFMDFlowTrueBin::Value(CorType t) const
+{ 
+  Double_t e;
+  return Value(e, t);
+}
+//____________________________________________________________________
+Double_t
+AliFMDFlowTrueBin::Value(Double_t& e2, CorType) const 
+{
+  return fHarmonic.Value(1, 0, e2);
+}
+//____________________________________________________________________
+Double_t
+AliFMDFlowTrueBin::Correction(Double_t& e2, CorType) const 
+{
+  e2 = fResReal.SqVar() / fResReal.N();
+  return fResReal.Average();
+}
+
+//____________________________________________________________________
+void
+AliFMDFlowTrueBin::Print(Option_t*) const 
+{
+  Double_t e2v, e2r;
+  Double_t v   = 100 * Value(e2v, AliFMDFlowBin::none);
+  Double_t r   = 100 * Correction(e2r, AliFMDFlowBin::none);
+
+  std::streamsize         old_prec  = std::cout.precision(3);
+  std::ios_base::fmtflags old_flags = std::cout.setf(std::ios_base::fixed, 
+                                                    std::ios_base::floatfield);
+  std::cout << "  v" << std::setw(1) << fHarmonic.Order() << ":   True: "
+           << std::setw(6) << v << " +/- " 
+           << std::setw(6) << 100*sqrt(e2v) << " ["
+           << std::setw(7) << r << " +/- " 
+           << std::setw(7) << 100*sqrt(e2r) << "]";
+  std::cout << std::endl;
+  std::cout.precision(old_prec);
+  std::cout.setf(old_flags, std::ios_base::floatfield);                               
+}
+
+//====================================================================
+AliFMDFlowTrue1D::AliFMDFlowTrue1D(UShort_t order, const AliFMDFlowAxis& xaxis)
+  : AliFMDFlowBinned1D(order, xaxis)
+{
+  // Delete old flow objects, and make new "true" ones. 
+  for (UInt_t i = 0; i < xaxis.N(); i++) { 
+    delete fBins[i];
+    fBins[i] = new AliFMDFlowTrueBin(order);
+  }
+}
+
+//____________________________________________________________________
+void
+AliFMDFlowTrue1D::SetPsi(Double_t psi) 
+{ 
+  for (UInt_t i = 0; i < fXAxis.N(); i++) 
+    static_cast<AliFMDFlowTrueBin*>(fBins[i])->SetPsi(psi);
+}
+
+//____________________________________________________________________
+void 
+AliFMDFlowTrue1D::Print(Option_t* option) const
+{
+  TString opt(option);
+  opt.ToLower();
+  Bool_t det = opt.Contains("d");
+  Bool_t sum = opt.Contains("s");
+  if (det) AliFMDFlowBinned1D::Print("d");
+  if (sum) { 
+    std::cout << "    x |              Real \n" 
+             << "------+-------------------" << std::endl;
+
+    std::streamsize         old_p = std::cout.precision(2);
+    std::ios_base::fmtflags old_f = std::cout.setf(std::ios_base::fixed, 
+                                                  std::ios_base::floatfield);
+    for (UShort_t i = 0; i < fXAxis.N(); i++) { 
+      Double_t x   = fXAxis.BinCenter(i);
+      Double_t e2v;
+      Double_t v   = fBins[i]->Value(e2v, AliFMDFlowBin::none);
+      std::cout << std::setprecision(2) << std::setw(5) << x << " | " 
+               << std::setprecision(3) 
+               << std::setw(6) << 100 * v << " +/- " 
+               << std::setw(6) << 100 * sqrt(e2v) 
+               << std::endl; 
+    }
+    std::cout.precision(old_p);
+    std::cout.setf(old_f, std::ios_base::floatfield);
+  }
+}
+
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/flow/AliFMDFlowTrue.h b/FMD/flow/AliFMDFlowTrue.h
new file mode 100644 (file)
index 0000000..eef57b7
--- /dev/null
@@ -0,0 +1,82 @@
+// -*- mode: C++ -*-
+/** @file 
+    @brief Declaration and implementation of classes to deal with a
+    well-known event plane @f$ \Psi_R@f$. */
+#ifndef FLOW_TRUE_H
+#define FLOW_TRUE_H
+#include <flow/AliFMDFlowStat.h>
+#include <flow/AliFMDFlowBin.h>
+#include <flow/AliFMDFlowBinned1D.h>
+
+
+/** @defgroup x_true Classes for handling Flow when the event plane is
+    known. */
+/** @class TrueBin flow/True.h <flow/True.h>
+    @brief A specialised Bin of flow in case the event plane is 
+    well-known. 
+    @ingroup x_true */
+struct AliFMDFlowTrueBin : public AliFMDFlowBin
+{
+  /** Constructor */ 
+  AliFMDFlowTrueBin(UShort_t order) : 
+    AliFMDFlowBin(order, 1), 
+    fPsiR(0), 
+    fResReal()
+  {}
+  /** Set the well-known event plane angle @f$ \Psi_R@f$ for this
+      event. 
+      @param psi @f$ \Psi_R@f$ */ 
+  void SetPsi(Double_t psi) { fPsiR = psi; } 
+  /** Should be called at the end of an event */ 
+  virtual void End();
+  /** Add a contribution @f$ \cos(n(\varphi-\Psi_R))@f$ where   
+      @f$ \Psi_R@f$ is the previously set, well-known event plane
+      angle. 
+      @param phi @f$ \varphi@f$ */ 
+  void AddToHarmonic(Double_t phi, Double_t) { fHarmonic.Add(phi, fPsiR); }
+  /** Get the value in this bin 
+      @param t  Which type of correction
+      @return the value of the harmonic */
+  virtual Double_t Value(CorType t=naive) const;
+  /** Get the value in this bin 
+      @param e2 On return, the square error. 
+      @param t  Which type of correction
+      @return the value of the harmonic */
+  Double_t Value(Double_t& e2, CorType t=naive) const;
+  /** Get the value in this bin 
+      @param e2 On return, the square error. 
+      @param t  Which type  of correction
+      @return the value of the harmonic */
+  Double_t Correction(Double_t& e2, CorType t=naive) const;
+  /** Print to standard out. */ 
+  void Print(Option_t* option="") const;
+  /** The well-known event plane */ 
+  Double_t fPsiR; 
+  /** True resolution */ 
+  AliFMDFlowStat fResReal;
+  /** define for ROOT I/O */
+  ClassDef(AliFMDFlowTrueBin,1);
+  
+}; 
+/** @brief A "histogram" of objects of class TrueBin in case the
+    event plane angle @f$ \Psi_R@f$ is well known. 
+    @ingroup x_true */
+struct AliFMDFlowTrue1D : public AliFMDFlowBinned1D
+{
+  /** Constructor */ 
+  AliFMDFlowTrue1D(UShort_t order, const AliFMDFlowAxis& xaxis);
+  /** Set the well-known event plane angle @f$ \Psi_R@f$ for this
+      event. 
+      @param psi @f$ \Psi_R@f$ */ 
+  void SetPsi(Double_t psi);
+  /** Print to standard out */ 
+  virtual void Print(Option_t* option="") const;
+  /** define for ROOT I/O */
+  ClassDef(AliFMDFlowTrue1D,1);
+};
+
+
+#endif
+//
+// EOF
+// 
diff --git a/FMD/flow/AliFMDFlowUtil.h b/FMD/flow/AliFMDFlowUtil.h
new file mode 100644 (file)
index 0000000..81d76f8
--- /dev/null
@@ -0,0 +1,29 @@
+// -*- mode: C++ -*-
+#ifndef FLOW_UTIL_H
+#define FLOW_UTIL_H
+#include <cmath>
+#ifndef M_PI
+# define M_PI 3.14159265358979323846264338327
+#endif
+
+/** @defgroup u_utils Utilities 
+    @brief Group of utility classes and functions */
+//__________________________________________________________________
+/** Normalize the angle @a ang to the interval @f$ [0,2\pi)@f$ 
+    @ingroup u_utils
+    @param ang Angle to normalize 
+    @return the normalised angle */
+inline Double_t 
+NormalizeAngle(Double_t ang) 
+{ 
+  while (ang <  0)      ang += 2 * M_PI;
+  while (ang >= 2*M_PI) ang -= 2 * M_PI;
+  return ang;
+}
+
+#endif
+//____________________________________________________________________
+//
+// EOF
+//
+
diff --git a/FMD/flow/FMDflowLinkDef.h b/FMD/flow/FMDflowLinkDef.h
new file mode 100644 (file)
index 0000000..c6a12d6
--- /dev/null
@@ -0,0 +1,47 @@
+// -*- mode: c++ -*- 
+#ifdef __CINT__
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+/** @file    FMDbaseLinkDef.h
+    @author  Christian Holm Christensen <cholm@nbi.dk>
+    @date    Mon Mar 27 14:18:46 2006
+    @brief   Link specifications for base library 
+*/
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliFMDFlowAxis+;
+#pragma link C++ class AliFMDFlowBin+;
+#pragma link C++ class AliFMDFlowBin;+;
+#pragma link C++ class AliFMDFlowBinned1D+;
+#pragma link C++ class AliFMDFlowBin;+;
+#pragma link C++ class AliFMDFlowBinned2D+;
+#pragma link C++ class AliFMDFlowEventPlane+;
+#pragma link C++ class AliFMDFlowHarmonic+;
+#pragma link C++ class AliFMDFlowResolution+;
+#pragma link C++ class AliFMDFlowResolutionStar+;
+#pragma link C++ class AliFMDFlowResolutionTDR+;
+#pragma link C++ class AliFMDFlowStat+;
+#pragma link C++ class AliFMDFlowTrueBin+;
+#pragma link C++ class AliFMDFlowTrue1D+;
+#pragma link C++ function  NormalizeAngle(Double_t);
+  
+
+#pragma link C++ namespace AliFMDFlowBessel;
+#pragma link C++ function  AliFMDFlowBessel::I01(Double_t,Double_t&,Double_t&,Double_t&,Double_t&)
+#pragma link C++ function  AliFMDFlowBessel::Ihalf(Int_t,Double_t,Double_t*,Double_t*);
+#pragma link C++ function  AliFMDFlowBessel::Iwhole(Int_t,Double_t,Double_t*,Double_t*);
+#pragma link C++ function  AliFMDFlowBessel::Inu(Double_t,Double_t,Double_t,Double_t*,Double_t*);
+#pragma link C++ function  AliFMDFlowBessel::I(Double_t,Double_t);
+#pragma link C++ function  AliFMDFlowBessel::DiffI(Double_t,Double_t);
+  
+
+#else
+# error Not for compilation 
+#endif
+//
+// EOF
+//
diff --git a/FMD/libFMDflow.pkg b/FMD/libFMDflow.pkg
new file mode 100644 (file)
index 0000000..2126513
--- /dev/null
@@ -0,0 +1,26 @@
+#-*- Mode: Makefile -*-
+#
+# $Id$
+
+SRCS           =  flow/AliFMDFlowAxis.cxx              \
+                  flow/AliFMDFlowBessel.cxx            \
+                  flow/AliFMDFlowBin.cxx               \
+                  flow/AliFMDFlowBinned1D.cxx          \
+                  flow/AliFMDFlowBinned2D.cxx          \
+                  flow/AliFMDFlowEfficiency.cxx        \
+                  flow/AliFMDFlowEventPlane.cxx        \
+                  flow/AliFMDFlowHarmonic.cxx          \
+                  flow/AliFMDFlowResolution.cxx        \
+                  flow/AliFMDFlowTrue.cxx              
+
+#                 AliFMDAltroIO.cxx            
+
+HDRS           =  $(SRCS:.cxx=.h)                      \
+                  flow/AliFMDFlowStat.h                \
+                  flow/AliFMDFlowUtil.h                
+
+DHDR           := flow/FMDflowLinkDef.h
+
+#
+# EOF
+#