--- /dev/null
+/** @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
+//
--- /dev/null
+// -*- 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
+//
--- /dev/null
+/** @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
+//
--- /dev/null
+// -*- 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
+//
+
+
--- /dev/null
+/** @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
+//
--- /dev/null
+// -*- 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
+//
--- /dev/null
+/** @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
+//
--- /dev/null
+// -*- 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
+//
--- /dev/null
+/** @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
+//
--- /dev/null
+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
+//
+
--- /dev/null
+#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
+//
+
+
+
+
--- /dev/null
+// -*- 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
+//
+
--- /dev/null
+/** @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
+//
+
--- /dev/null
+// -*- 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
+//
--- /dev/null
+/** @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
+//
--- /dev/null
+// -*- 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
+//
+
--- /dev/null
+/** @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
+//
--- /dev/null
+// -*- 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
+//
--- /dev/null
+// -*- 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
+//
--- /dev/null
+#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
+//
--- /dev/null
+// -*- 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
+//
--- /dev/null
+// -*- 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
+//
+
--- /dev/null
+// -*- 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
+//
--- /dev/null
+#-*- 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
+#
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff