FMD/flow/AliFMDFlowResolution.h

   1 // -*- mode: C++ -*-
   2 /* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Lesser General Public License
   6  * as published by the Free Software Foundation; either version 2.1 of
   7  * the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful, but
  10  * WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Lesser General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Lesser General Public
  15  * License along with this library; if not, write to the Free Software
  16  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
  17  * USA
  18  */
  19 /** @file
  20     @brief Declaration of an Resolution class */
  21 //____________________________________________________________________
  22 //
  23 // Calculate the event plane resolution.
  24 // Input is the observed phis.
  25 // There's a number of implementations of this.
  26 #ifndef ALIFMDFLOWRESOLUTION_H
  27 #define ALIFMDFLOWRESOLUTION_H
  28 #include <flow/AliFMDFlowStat.h>
  29 #include <TH1D.h>
  30 #include <TH2D.h>
  31 class TBrowser;
  32
  33 //______________________________________________________
  34 /** @class AliFMDFlowResolution flow/AliFMDFlowResolution.h <flow/AliFMDFlowResolution.h>
  35     @brief Class to calculate the event plane resolution based on two
  36     sub-events
  37     @ingroup a_basic
  38
  39     This class calculates the event plane resolution based on the
  40     basic formulas given in Phys. Rev. @b C58, 1671.   That is, the
  41     resolution is given by
  42     @f[
  43     R_{k} = \langle\cos(km(\Psi_m-\Psi_R))\rangle
  44     @f]
  45     where @f$ \Psi_R@f$ is the unknown @e true event plane angle, and
  46     @f$ n = km@f$
  47
  48     Using two random sub-events, @f$ A, B@f$ we get that
  49     @f[
  50     \langle\cos(n(\Psi^A_m-\Psi^B_m))\rangle =
  51     \langle\cos(n(\Psi^A_m-\Psi_R))\rangle
  52     \langle\cos(n(\Psi^B_m-\Psi_R))\rangle
  53     @f]
  54     If the sub-events are of equal size, and randomly chosen, then we
  55     get that
  56     @f[
  57     \langle\cos(n(\Psi^A_m-\Psi_R))\rangle =
  58     \sqrt{\langle\cos(n(\Psi^A_m-\Psi^B_m))}
  59     @f]
  60     and it follows that
  61     @f{eqnarray*}
  62     \langle\cos(km(\Psi_m-\Psi_R))\rangle & = &
  63     \sqrt{2}\langle\cos(n(\Psi^A_m-\Psi_R))\rangle\\
  64     \sqrt{2}\sqrt{\langle\cos(n(\Psi^A_m-\Psi^B_m))}
  65     @f}
  66
  67     Hence, the event-plane resolution is simply the square root of the
  68     scaled average distance between the two sub-events, multiplied by
  69     @f$ \sqrt{s}@f$
  70
  71     The error is therefor @f$ \sqrt{s}@f$ times the variance of the
  72     cosine of the distance between the two sub-events.
  73 */
  74 class AliFMDFlowResolution : public AliFMDFlowStat
  75 {
  76 public:
  77   /** Constructor
  78       @param n Harmonic order */
  79   AliFMDFlowResolution(UShort_t n=0);
  80   /** Destructor */
  81   virtual ~AliFMDFlowResolution() {}
  82   /** Copy constructor
  83       @param o Object to copy from */
  84   AliFMDFlowResolution(const AliFMDFlowResolution& o);
  85   /** Assignment operator
  86       @param o Object to copy from */
  87   AliFMDFlowResolution& operator=(const AliFMDFlowResolution& o);
  88
  89   /** @{
  90       @name Processing */
  91   /** add data point
  92       @param psiA A sub-event plane angle @f$ \Psi_A \in[0,2\pi]@f$
  93       @param psiB B sub-event plane angle @f$ \Psi_B \in[0,2\pi]@f$ */
  94   virtual void Add(Double_t psiA, Double_t psiB);
  95   /** @} */
  96
  97   /** @{
  98       @name Information */
  99   /** Get the correction for harmonic strength of order @a k
 100       @param k  The harminic strenght order to get the correction for
 101       @param e2 The square error on the correction
 102       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */
 103   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
 104   /** Get the correction for harmonic strength of order @a k
 105       @param k The harminic strenght order to get the correction for
 106       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */
 107   virtual Double_t Correction(UShort_t k) const;
 108   /** Get the harmnic order */
 109   UShort_t Order() const { return fOrder; }
 110   /** @} */
 111
 112   /** @{
 113       @name Utility */
 114   /** Draw this corrrection function
 115       @param option Options passed to drawing */
 116   virtual void Draw(Option_t* option=""); //*MENU*
 117   /** This is a folder */
 118   Bool_t IsFolder() const { return kTRUE; }
 119   /** Look through this object */
 120   virtual void Browse(TBrowser* b);
 121   /** @} */
 122
 123   /** @{
 124       @name histograms */
 125   /** @return Contrib histogram */
 126   const TH1& ContribHistogram() const { return fContrib; }
 127   /** @} */
 128 protected:
 129   /** Order */
 130   UShort_t fOrder; // Order
 131   /** Contributions */
 132   TH1D fContrib;
 133   /** Define for ROOT I/O */
 134   ClassDef(AliFMDFlowResolution,1);
 135 };
 136
 137 //______________________________________________________
 138 /**
 139     @ingroup a_basic
 140     The event plane angle resolution function is
 141     @f[
 142     R_k(\chi) = \frac{\pi} \chi e^{-\chi^2/4}\left(
 143     I_{\frac{k-1}2}(\chi^2/4) + I_{\frac{k+1}2}(\chi^2/4)\right)
 144     @f]
 145     Where @f$ I_n(x)@f$ is the modified Bessel function of the first
 146     kind.   Identifying
 147     @f[
 148     y = \chi^2/4\quad C=\frac{\sqrt{\pi}e^{-y}}{2\sqrt{2}}
 149     @f]
 150     and the short hands @f$ f2(y) = I_{\frac{k-1}2}@f$ and @f$ f3(y) =
 151     I_{\frac{k+1}2}(y)@f$ we can write this more compact as
 152     @f[
 153     R_k(y) = C y (f2(y) + f3(y))
 154     @f]
 155     The derivative of the resolution function is
 156     @f[
 157     R_k'(y) = \frac{C}{2}\left[4\sqrt{y}\left(f2'(y)-f3'(y)\right)
 158     - (4 y - 2)\left(f2(y) + f3(y)\right)\right]
 159     @f]
 160     Since
 161     @f[
 162     I_\nu'(x) = I_{\nu-1}(x) - \frac{\nu}{x} I_\nu(x)\quad,
 163     @f]
 164     and setting @f$ f1(y) = I_{\frac{k-3}2}(y)@f$, we get
 165     @f[
 166     R_k'(y) = \frac{C}{2}\left[4y f1(y) + (4-2k) f2(y) - (4y+2k)
 167     f3(y)\right]
 168     @f]
 169
 170     In this class, the argument @f$ \chi@f$ is estimated by finding
 171     the minima of @f$ R_k(\chi)@f$ near the average of
 172     @f$\cos(n(\Psi_A-\Psi_B))@f$.  The error @f$ \delta\chi@f$ is
 173     estimated as the largest step size in the minimisation.
 174
 175     The total error on the correction is then
 176     @f[
 177     \delta R_k = R_k'(\chi) \delta\chi
 178     @f]
 179 */
 180
 181 class AliFMDFlowResolutionStar : public AliFMDFlowResolution
 182 {
 183 public:
 184   /** Constructor
 185       @param n Harmonic order */
 186   AliFMDFlowResolutionStar(UShort_t n=0) : AliFMDFlowResolution(n) {}
 187   /** Destructor */
 188   virtual ~AliFMDFlowResolutionStar() {}
 189
 190   /** @{
 191       @name Information */
 192   /** Get the correction for harmonic strength of order @a k
 193       @param k The harminic strenght order to get the correction for
 194       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */
 195   virtual Double_t Correction(UShort_t k) const;
 196   /** Get the correction for harmonic strength of order @a k
 197       @param k  The harminic strenght order to get the correction for
 198       @param e2 The square error on the correction
 199       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */
 200   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
 201   /** Get @f$ \chi@f$
 202       @param res  First shot at the resolution.
 203       @param k    Order
 204       @param delta On return, the last step size in @f$ \chi@f$ -
 205       which is taken to be @f$ \delta\chi@f$
 206       @return @f$\chi@f$ */
 207   virtual Double_t Chi(Double_t res, UShort_t k, Double_t& delta) const;
 208   /** @} */
 209
 210   /** @{
 211       @name Utilities */
 212   /** Draw this corrrection function
 213       @param option Options passed to drawing */
 214   virtual void Draw(Option_t* option=""); //*MENU*
 215   /** @} */
 216 protected:
 217   /** Calculate resolution
 218       @param chi @f$ \chi@f$
 219       @param k   Order factor
 220       @param dr  On return, the derivative of @f$ R(\chi)@f$
 221       @return
 222       @f[
 223       \frac{\sqrt{\pi/2}}{2}\chi e^{-\chi^2/4}
 224       (I_{\frac{(k-1)}{2}}(\chi^2/4)+ I_{\frac{(k+1)}{2}}(\chi^2/4))
 225       @f] */
 226   Double_t Res(Double_t chi, UShort_t k, Double_t& dr) const;
 227   /** Define for ROOT I/O */
 228   ClassDef(AliFMDFlowResolutionStar,1);
 229 };
 230
 231 //______________________________________________________
 232 /**
 233     @ingroup a_basic
 234
 235     For more on the event plane angle resolution function, please
 236     refer to the class description of ResolutionStar.
 237
 238     In this class @f$ \chi@f$ is calculated from the ratio @f$ k/N@f$
 239     of events with @f$ |\Psi_A - \Psi_B| > \pi/2@f$ to the total
 240     number of events.
 241
 242     The pre-print @c nucl-ex/9711003v2 gives the formula
 243     @f[
 244     \frac{k}{N} = \frac{e^{-\chi^2/2}}{2}
 245     @f]
 246     for @f$ \chi@f$.  Note, that this differs from the @f$ \chi@f$
 247     used in ResolutionStar by a factor of @f$ 1/sqrt{2}@f$.
 248
 249     We can isolate @f$ \chi = \mp\sqrt{-2\log(2k/N)}@f$ from the
 250     equation above.
 251
 252     Since @f$ r=k/N@f$ is obviously a efficiency-like ratio, we get
 253     that error @f$ \delta r@f$ is given by Binomial errors
 254     @f[
 255     \delta r =  \sqrt{r\frac{1 - r}{N}}\quad.
 256     @f]
 257     The total error @f$ \delta\chi@f$ then becomes
 258     @f[
 259     \delta^2\chi = \left(\frac{d\chi}{dr}\right)^2 \delta^2r
 260     = \frac{r - 1}{4 k log(2 r)}
 261     @f]
 262
 263     The total error on the correction is, as in
 264     ResolutionStar, then given by
 265     @f[
 266     \delta R_k = R_k'(\chi) \delta\chi
 267     @f]
 268 */
 269 class AliFMDFlowResolutionTDR : public AliFMDFlowResolution
 270 {
 271 public:
 272   /** Constructor
 273       @param n Harmonic order */
 274   AliFMDFlowResolutionTDR(UShort_t n=0);
 275   /** DTOR */
 276   ~AliFMDFlowResolutionTDR() {}
 277   /** Copy constructor
 278       @param o  Object to copy from */
 279   AliFMDFlowResolutionTDR(const AliFMDFlowResolutionTDR& o);
 280   /** Assignment operator
 281       @param o  Object to assign from
 282       @return reference to this */
 283   AliFMDFlowResolutionTDR& operator=(const AliFMDFlowResolutionTDR& o);
 284
 285   /** @{
 286       @name Processing */
 287   virtual void Clear(Option_t* option="");
 288   /** add a data  point */
 289   virtual void Add(Double_t psiA, Double_t psiB);
 290   /** @} */
 291
 292   /** @{
 293       @name Information */
 294   /** Get the correction for harmonic strength of order @a k
 295       @param k The harminic strenght order to get the correction for
 296       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */
 297   virtual Double_t Correction(UShort_t k) const;
 298   /** Get the correction for harmonic strength of order @a k
 299       @param e2 The square error on the correction
 300       @param k  The harminic strenght order to get the correction for
 301       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */
 302   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
 303   /** Get @f$ \chi^2/2@f$
 304       @param e2 The square error on the correction
 305       @return @f$ \chi^2/2@f$ */
 306   virtual Double_t Chi2Over2(Double_t r, Double_t& e2) const;
 307   /** @} */
 308
 309   /** @{
 310       @name Utilities */
 311   /** Draw this corrrection function
 312       @param option Options passed to drawing */
 313   virtual void Draw(Option_t* option=""); //*MENU*
 314   /** Look through this object */
 315   virtual void Browse(TBrowser* b);
 316   /** @} */
 317
 318   /** @{
 319       @name histograms */
 320   /** @return NK histogram */
 321   const TH1& NKHistogram() const { return fNK; }
 322   /** @} */
 323 protected:
 324   /** Calculate resolution
 325       @param k    Order factor
 326       @param y    @f$ \chi^2/2@f$
 327       @param echi @f$\delta\chi@f$
 328       @param dr   On return, the derivative of @f$ R(\chi)@f$
 329       @return
 330       @f[
 331       \frac{\sqrt{\pi/2}}{2}\chi e^{-\chi^2/2}
 332       (I_{\frac{(k-1)}{2}}(\chi^2/2)+ I_{\frac{(k+1)}{2}}(\chi^2/2))
 333       @f] */
 334   Double_t Res(UShort_t k, Double_t y, Double_t echi2, Double_t& e2) const;
 335   /** Number of events with large diviation */
 336   ULong_t fLarge; // Number of events with large angle
 337   /** Histogram of large and all angles */
 338   TH1D fNK;
 339   /** Define for ROOT I/O */
 340   ClassDef(AliFMDFlowResolutionTDR,1);
 341 };
 342
 343
 344 #endif
 345 //
 346 // EOF
 347 //