FMD/flow/AliFMDFlowEventPlane.cxx

   1 /* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
   2  *
   3  * This library is free software; you can redistribute it and/or
   4  * modify it under the terms of the GNU Lesser General Public License
   5  * as published by the Free Software Foundation; either version 2.1 of
   6  * the License, or (at your option) any later version.
   7  *
   8  * This library is distributed in the hope that it will be useful, but
   9  * WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11  * Lesser General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Lesser General Public
  14  * License along with this library; if not, write to the Free Software
  15  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
  16  * USA
  17  */
  18 /** @file
  19     @brief Implementation of an EventPlane class */
  20 //____________________________________________________________________
  21 //
  22 // Class to determine the event plane
  23 //
  24 // The event plane is calculated as
  25 //
  26 //    Psi_n = 1/n * atan((sum_i(w_i sin(n phi_i)))
  27 //                        sum_i(w_i cos(n phi_i))))
  28 //
  29 // where i runs over all observations of phi in an event, and
  30 // w_i is the weight of the ith observation of phi
  31 #include "flow/AliFMDFlowEventPlane.h"
  32 #include "flow/AliFMDFlowUtil.h"
  33 #include <TMath.h>
  34 #include <TBrowser.h>
  35 // #include <iostream>
  36 // #include <cmath>
  37 #ifndef _GNU_SOURCE
  38 //#warning Using private implementation of sincos
  39 extern "C"
  40 {
  41   /** Function to caculate @f$ \sin(a), \cos(a)@f$ in one go.  Note,
  42       that with GCC, this is a built-in, and we only have this to
  43       serve as a replacement function
  44       @ingroup utils */
  45   inline void sincos(Double_t a, Double_t* sina, Double_t* cosa) {
  46     *sina = sin(a);
  47     *cosa = cos(a);
  48   }
  49 }
  50 #endif
  51
  52 //====================================================================
  53 AliFMDFlowEventPlane::AliFMDFlowEventPlane(UShort_t m)
  54   : fSumSinMPhi(0),
  55     fSumCosMPhi(0),
  56     fOrder(m),
  57     fCache(0),
  58     fSum("sum", Form("#sumw#cos(%d#varphi) vs #sum#sin(%d#varphi)", m, m),
  59          100, -1.1, 1.1, 100, -1.1, 1.1),
  60     fPsi("psi", Form("#Psi_{%d}", m), 80, 0, 2*TMath::Pi()),
  61     fScale(0)
  62 {
  63   Clear();
  64   fSum.SetDirectory(0);
  65   fSum.SetXTitle(Form("#sum_{i}w_{i}cos(%d#varphi_{i})", fOrder));
  66   fSum.SetYTitle(Form("#sum_{i}w_{i}sin(%d#varphi_{i})", fOrder));
  67   fSum.SetMarkerStyle(20);
  68   fPsi.SetDirectory(0);
  69   fPsi.SetXTitle(Form("#Psi_{%d}", fOrder));
  70 }
  71
  72 //____________________________________________________________________
  73 AliFMDFlowEventPlane::AliFMDFlowEventPlane(const AliFMDFlowEventPlane& o)
  74   : TObject(o),
  75     fSumSinMPhi(o.fSumSinMPhi),
  76     fSumCosMPhi(o.fSumCosMPhi),
  77     fOrder(o.fOrder),
  78     fCache(-1),
  79     fSum(o.fSum),
  80     fPsi(o.fPsi),
  81     fScale(o.fScale)
  82 {
  83   // copy cosntructor
  84   // Parameters
  85   //   o  Object to copy from.
  86   fSum.SetDirectory(0);
  87   fSum.SetXTitle(Form("#sum_{i}w_{i}cos(%d#varphi_{i})", fOrder));
  88   fSum.SetYTitle(Form("#sum_{i}w_{i}sin(%d#varphi_{i})", fOrder));
  89   fSum.SetMarkerStyle(20);
  90   fPsi.SetDirectory(0);
  91   fPsi.SetXTitle(Form("#Psi_{%d}", fOrder));
  92 }
  93 //____________________________________________________________________
  94 AliFMDFlowEventPlane&
  95 AliFMDFlowEventPlane::operator=(const AliFMDFlowEventPlane& o)
  96 {
  97   // Assignment operator.
  98   // Parameters:
  99   //  o Object to assign from
 100   fSumSinMPhi = o.fSumSinMPhi;
 101   fSumCosMPhi = o.fSumCosMPhi;
 102   fOrder      = o.fOrder;
 103   fCache      = -1;
 104   fScale      = o.fScale;
 105
 106   fSum.Reset();
 107   fSum.Add(&o.fSum);
 108   fPsi.Reset();
 109   fPsi.Add(&o.fPsi);
 110
 111   return *this;
 112 }
 113
 114 //____________________________________________________________________
 115 void
 116 AliFMDFlowEventPlane::Clear(Option_t*)
 117 {
 118   // clear internal variables.
 119   fSumSinMPhi = 0;
 120   fSumCosMPhi = 0;
 121   fCache      = -1;
 122   fScale      = 0;
 123 }
 124 //____________________________________________________________________
 125 void
 126 AliFMDFlowEventPlane::Add(Double_t phi, Double_t weight)
 127 {
 128   // Add a data point
 129   // Parameters:
 130   //   phi     The angle phi in[0,2pi]
 131   //   weight  The weight
 132   Double_t a = NormalizeAngle(fOrder * phi);
 133   Double_t s, c;
 134   sincos(a, &s, &c);
 135   if (TMath::IsNaN(s) || !TMath::Finite(s) ||
 136       TMath::IsNaN(c) || !TMath::Finite(s)) return;
 137   if (weight == 0) return;
 138   fScale += 1./weight;
 139   fSumSinMPhi += weight * s;
 140   fSumCosMPhi += weight * c;
 141 }
 142
 143 //____________________________________________________________________
 144 void
 145 AliFMDFlowEventPlane::End()
 146 {
 147   Double_t r2 = fSumCosMPhi*fSumCosMPhi + fSumSinMPhi*fSumSinMPhi;
 148   Double_t r  = (r2 < 0 ? 1 : TMath::Sqrt(r2));
 149   fSum.Fill(fSumCosMPhi/(r!=0?r:1),fSumSinMPhi/(r!=0?r:1));
 150   fPsi.Fill(Psi());
 151 }
 152
 153 //____________________________________________________________________
 154 void
 155 AliFMDFlowEventPlane::Browse(TBrowser* b)
 156 {
 157   b->Add(&fSum);
 158   b->Add(&fPsi);
 159 }
 160
 161 //____________________________________________________________________
 162 Double_t
 163 AliFMDFlowEventPlane::Psi() const
 164 {
 165   // Get the event plane
 166   // Parameters:
 167   //   none
 168   if (fCache < 0) {
 169     fCache = DoPsi(fSumSinMPhi, fSumCosMPhi);
 170   }
 171   return fCache;
 172 }
 173 //____________________________________________________________________
 174 Double_t
 175 AliFMDFlowEventPlane::Psi(Double_t phi, Double_t w) const
 176 {
 177   // Get the event plane angle Psi_k disregarding the contribution
 178   // from the observation phi with weight w.  This is to avoid
 179   // auto-correlations
 180   //
 181   // Parameters:
 182   //  phi   The observation  phi
 183   //  w     The weight w of the obervation.
 184   //
 185   // Returns The event plane angle Psi with out the contribution from
 186   // phi_i
 187   Double_t a = NormalizeAngle(fOrder * phi);
 188   Double_t s, c;
 189   sincos(a, &s, &c);
 190   if (TMath::IsNaN(s) || !TMath::Finite(s) ||
 191       TMath::IsNaN(c) || !TMath::Finite(s)) return Psi();
 192   Double_t psi = DoPsi(fSumSinMPhi - w * s, fSumCosMPhi - w * c);
 193   return psi;
 194 }
 195
 196 //____________________________________________________________________
 197 Double_t
 198 AliFMDFlowEventPlane::DoPsi(Double_t sumsin, Double_t sumcos) const
 199 {
 200   // Calculate the event plane
 201   // Parameters:
 202   //   sumsin    Sum of sines
 203   //   sumcos    Sum of cosines
 204   Double_t psi = 0;
 205   // Make sure we get an angle everywhere
 206   if      (sumcos != 0) psi =  atan2(sumsin, sumcos);
 207   else if (sumsin == 0) psi =  0;
 208   else if (sumsin >  0) psi =  M_PI / 2;
 209   else                  psi = -M_PI / 2;
 210   psi =  NormalizeAngle(psi);
 211   psi /= fOrder;
 212   return psi;
 213 }
 214
 215 //____________________________________________________________________
 216 //
 217 // EOF
 218 //
 219