FMD/flow/AliFMDFlowBin.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 a Bin in a Flow histogram */
  20 //____________________________________________________________________
  21 //
  22 // This contains an of class AliFMDFlowHarmonic and an object of
  23 // class AliFMDFlowEventPlane to calculate v_n and \Psi_k.  It contain
  24 // two objects of class AliFMDFlowEventPlane to calculate the
  25 // sub-event event planes Psi_A, \Psi_B.  It also contain 3 objects of
  26 // class AliFMDFlowResolution to calculate the event plane angle
  27 // resolution.
  28 //
  29 #include "flow/AliFMDFlowBin.h"
  30 #include "flow/AliFMDFlowUtil.h"
  31 #include <cmath>
  32 #include <iostream>
  33 #include <iomanip>
  34 #include <TBrowser.h>
  35
  36 //====================================================================
  37 AliFMDFlowBin::AliFMDFlowBin(UShort_t order, UShort_t k)
  38   : fPsi(order / k),
  39     fPsiA(order / k),
  40     fPsiB(order / k),
  41     fRes(order / k),
  42     fResStar(order / k),
  43     fResTdr(order / k),
  44     fHarmonic(order),
  45     fSplit("ab_split", "Relative split in A,B sub-event", 2, 0, 2, 100, 0, 1),
  46     fPhi("phi", "Distribution of #varphi-#Psi", 40, 0, 2*TMath::Pi()),
  47     fNA(0),
  48     fNB(0),
  49     fN(0),
  50     fAB("psiAB", "Distribution of #Psi_{A} vs. #Psi_{B}",
  51         40, 0, 2*TMath::Pi()/fPsiA.Order(), 40, 0, 2*TMath::Pi()/fPsiA.Order())
  52 {
  53   fSplit.SetDirectory(0);
  54   fSplit.GetXaxis()->SetBinLabel(1, "A");
  55   fSplit.GetXaxis()->SetBinLabel(2, "B");
  56   fSplit.SetXTitle("Sub-event");
  57   fSplit.SetYTitle("Fraction");
  58   fPhi.SetDirectory(0);
  59   fPhi.SetXTitle("#varphi");
  60   fAB.SetDirectory(0);
  61   fAB.SetXTitle("#Psi_{A}");
  62   fAB.SetYTitle("#Psi_{B}");
  63 }
  64
  65 //____________________________________________________________________
  66 AliFMDFlowBin::AliFMDFlowBin(const AliFMDFlowBin& o)
  67   : TObject(o),
  68     fPsi(o.fPsi),
  69     fPsiA(o.fPsiA),
  70     fPsiB(o.fPsiB),
  71     fRes(o.fRes),
  72     fResStar(o.fResStar),
  73     fResTdr(o.fResTdr),
  74     fHarmonic(o.fHarmonic),
  75     fSplit(o.fSplit),
  76     fPhi(o.fPhi),
  77     fNA(o.fNA),
  78     fNB(o.fNB),
  79     fN(o.fN),
  80     fAB(o.fAB)
  81 {
  82   // Copy constructor
  83   // Parameters:
  84   //    o Object to copy from
  85   fSplit.SetDirectory(0);
  86   fSplit.GetXaxis()->SetBinLabel(1, "A");
  87   fSplit.GetXaxis()->SetBinLabel(2, "B");
  88   fSplit.SetXTitle("Sub-event");
  89   fSplit.SetYTitle("Fraction");
  90
  91   fPhi.SetDirectory(0);
  92   fPhi.SetXTitle("#varphi");
  93
  94   fAB.SetDirectory(0);
  95   fAB.SetXTitle("#Psi_{A}");
  96   fAB.SetYTitle("#Psi_{B}");
  97 }
  98
  99 //____________________________________________________________________
 100 AliFMDFlowBin&
 101 AliFMDFlowBin::operator=(const AliFMDFlowBin& o)
 102 {
 103   // Assignment operator
 104   // Parameters:
 105   //    o Object to assign from
 106   fPsi      = o.fPsi;
 107   fPsiA     = o.fPsiA;
 108   fPsiB     = o.fPsiB;
 109   fRes      = o.fRes;
 110   fResStar  = o.fResStar;
 111   fResTdr   = o.fResTdr;
 112   fHarmonic = o.fHarmonic;
 113   fNA       = o.fNA;
 114   fNB       = o.fNB;
 115   fN        = o.fN;
 116
 117
 118   fSplit.Reset();
 119   fSplit.Add(&(o.fSplit));
 120
 121   fPhi.Reset();
 122   fPhi.Add(&(o.fPhi));
 123
 124   fAB.Reset();
 125   fAB.Add(&(o.fAB));
 126
 127   return *this;
 128 }
 129
 130 //____________________________________________________________________
 131 void
 132 AliFMDFlowBin::Begin()
 133 {
 134   // Clear event plane calculators
 135   fPsi.Clear();
 136   fPsiA.Clear();
 137   fPsiB.Clear();
 138   fNA = fNB = fN = 0;
 139 }
 140
 141 //____________________________________________________________________
 142 void
 143 AliFMDFlowBin::AddToEventPlane(Double_t phi, Double_t w, Bool_t a)
 144 {
 145   // Called to add a contribution to the event plane
 146   // Parameters
 147   //    phi   The angle phi in [0,2pi]
 148   //    w     Weight
 149   //    a    If true, add to sub-event A, otherwise to sub-event B.
 150   fPsi.Add(phi, w);
 151   if (a) fPsiA.Add(phi, w);
 152   else   fPsiB.Add(phi, w);
 153   if (a) fNA++; else fNB++;
 154   fN++;
 155 }
 156
 157 //____________________________________________________________________
 158 void
 159 AliFMDFlowBin::AddToHarmonic(Double_t phi, Double_t wp, Double_t wh)
 160 {
 161   // Called to add a contribution to the harmonic.
 162   // Parameters:
 163   //   phi   The angle phi in [0,2pi]
 164   //   wp    Weight of phi (only used in the calculation of
 165   //         the event plane).
 166   //   wh    Weight of observation.
 167
 168   // Disregard the obervation of phi from the event plane angle.
 169   Double_t psi   = fPsi.Psi(phi, wp);
 170   fHarmonic.Add(phi, psi, wh);
 171   fPhi.Fill(NormalizeAngle(phi-psi));
 172 }
 173
 174 //____________________________________________________________________
 175 void
 176 AliFMDFlowBin::End()
 177 {
 178   // Should be called at the end of an event
 179   fPsi.End();
 180   fPsiA.End();
 181   fPsiB.End();
 182   Double_t psiA = fPsiA.Psi();
 183   Double_t psiB = fPsiB.Psi();
 184
 185   // Update the resolutions
 186   fRes.Add(psiA, psiB);
 187   fResStar.Add(psiA, psiB);
 188   fResTdr.Add(psiA, psiB);
 189   if (fN != 0) {
 190     fSplit.Fill(.5,  float(fNA)/fN);
 191     fSplit.Fill(1.5, float(fNB)/fN);
 192   }
 193   fAB.Fill(psiA, psiB);
 194 }
 195
 196 //____________________________________________________________________
 197 void
 198 AliFMDFlowBin::Event(UInt_t n, Double_t* phis, Double_t* wp, Double_t* wh)
 199 {
 200   // Analyse events
 201   // Parameters
 202   //  n      Size of phis and possibly ws
 203   //  phis   Array of phi, (phi_1, ..., phi_n)
 204   //  wp     Weights of event plane (optional)
 205   //  wh     Weights of harmonic (optional)
 206   Begin();
 207
 208   // Calculate split.
 209   UInt_t split = n / 2;
 210   // First sub-event.
 211   for (UInt_t i = 0; i < split; i++)
 212     AddToEventPlane(phis[i], (wp ? wp[i] : 1), kTRUE);
 213   // Second sub-event.
 214   for (UInt_t i = split; i < n; i++)
 215     AddToEventPlane(phis[i], (wp ? wp[i] : 1), kFALSE);
 216   // Add contributions to the harmonic.
 217   for (UInt_t i = 0; i < n; i++)
 218     AddToHarmonic(phis[i], (wp ? wp[i] : 1), (wh ? wh[i] : 1));
 219
 220   End();
 221 }
 222
 223 //____________________________________________________________________
 224 Double_t
 225 AliFMDFlowBin::Value(CorType t) const
 226 {
 227   // Get the value in this bin
 228   // Parameters:
 229   //   t  Which type of correction
 230   //
 231   // return the value of the harmonic
 232   Double_t e;
 233   return Value(e, t);
 234 }
 235
 236 //____________________________________________________________________
 237 Double_t
 238 AliFMDFlowBin::EValue(CorType t) const
 239 {
 240   // Get the value in this bin
 241   // Parameters:
 242   //    t    Which type of correction
 243   //
 244   // return the error on the value of the harmonic
 245   Double_t e2;
 246   Value(e2, t);
 247   return sqrt(e2);
 248 }
 249
 250 //____________________________________________________________________
 251 Double_t
 252 AliFMDFlowBin::Value(Double_t& e2, CorType t) const
 253 {
 254   // Get the value in this bin
 255   // Parameters:
 256   //    e2   On return, the square error.
 257   //    t    Which type  of correction
 258   //
 259   // return the value of the harmonic
 260   Double_t r, er2;
 261   r = Correction(er2, t);
 262   return fHarmonic.Value(r, er2, e2);
 263 }
 264
 265 //____________________________________________________________________
 266 Double_t
 267 AliFMDFlowBin::Correction(Double_t& er2, CorType t) const
 268 {
 269   // Get the value in this bin
 270   // Parameters:
 271   //    e2   On return, the square error.
 272   //    t    Which type  of correction
 273   //
 274   // return the value of the Correction
 275   Double_t r = 1;
 276   UShort_t k = fHarmonic.Order()/fRes.Order();
 277   switch (t) {
 278   case kNaive: r = fRes.Correction(k, er2);     break;
 279   case kStar:  r = fResStar.Correction(k, er2); break;
 280   case kTdr:   r = fResTdr.Correction(k, er2);  break;
 281   default:     r = 1; er2 = 0;                  break;
 282   }
 283   return r;
 284 }
 285
 286 //____________________________________________________________________
 287 ULong_t
 288 AliFMDFlowBin::Counts() const
 289 {
 290   // Return the number of counts used in this bin.
 291   // Return:
 292   //  Number of counts that is used in this bin.
 293   return fHarmonic.N();
 294 }
 295
 296 //____________________________________________________________________
 297 void
 298 AliFMDFlowBin::Finish()
 299 {
 300   // Called at the end of the event
 301 }
 302
 303 //____________________________________________________________________
 304 void
 305 AliFMDFlowBin::Browse(TBrowser* b)
 306 {
 307   // Browse this item
 308   b->Add(&fPsi,      "Full event plane");
 309   b->Add(&fPsiA,     "Sub-event A event plane");
 310   b->Add(&fPsiB,     "Sub-event B event plane");
 311   b->Add(&fRes,      "Naive resolution");
 312   b->Add(&fResStar,  "STAR resolution");
 313   b->Add(&fResTdr,   "TDR resolution");
 314   b->Add(&fHarmonic, "Harmonic");
 315   b->Add(&fSplit,    "Split");
 316   b->Add(&fPhi,      "Phi");
 317   b->Add(&fAB,       "AB");
 318 }
 319
 320 //____________________________________________________________________
 321 void
 322 AliFMDFlowBin::Print(Option_t*) const
 323 {
 324   // Print information
 325   Double_t e2v[4], v[4], r[4], e2r[4];
 326   const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
 327   v[0] = 100 * Value(e2v[0], AliFMDFlowBin::kNone);
 328   v[1] = 100 * Value(e2v[1], AliFMDFlowBin::kNaive);
 329   v[2] = 100 * Value(e2v[2], AliFMDFlowBin::kStar);
 330   v[3] = 100 * Value(e2v[3], AliFMDFlowBin::kTdr);
 331   r[0] = 100 * Correction(e2r[0], AliFMDFlowBin::kNone);
 332   r[1] = 100 * Correction(e2r[1], AliFMDFlowBin::kNaive);
 333   r[2] = 100 * Correction(e2r[2], AliFMDFlowBin::kStar);
 334   r[3] = 100 * Correction(e2r[3], AliFMDFlowBin::kTdr);
 335
 336   std::streamsize         oldP = std::cout.precision(3);
 337   std::ios_base::fmtflags oldF = std::cout.setf(std::ios_base::fixed,
 338                                                 std::ios_base::floatfield);
 339   std::cout << "  v" << std::setw(1) << fHarmonic.Order() << ": ";
 340   for (UInt_t i = 0; i < 4; i++)
 341     std::cout << std::setw(6+(i == 0 ? 0 : 6)) << names[i] << ": "
 342               << std::setw(6) << v[i] << " +/- "
 343               << std::setw(6) << 100*sqrt(e2v[i]) << " ["
 344               << std::setw(7) << r[i] << " +/- "
 345               << std::setw(7) << 100*sqrt(e2r[i]) << "]\n";
 346   std::cout << std::flush;
 347   std::cout.precision(oldP);
 348   std::cout.setf(oldF, std::ios_base::floatfield);
 349 }
 350
 351
 352 //____________________________________________________________________
 353 //
 354 // EOF
 355 //