FMD/flow/AliFMDFlowBinned1D.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 1-dimensional Flow "histogram" */
  20 //____________________________________________________________________
  21 //
  22 // A histogram of flow bins.  The axis can by anything
  23 // (pseudo-rapidity, transvers momentum) - there's no assumption on
  24 // what is the basis of the histogram.  The method Event can be used
  25 // to calculate everything in one go.   Alternatively, one can use the
  26 // methods AddToEventPlane and AddToHarmonic.  See also the example
  27 // TestFlow.C
  28 #include "flow/AliFMDFlowBinned1D.h"
  29 #include "flow/AliFMDFlowBin.h"
  30 #include <cmath>
  31 #include <cstdlib>
  32 #include <iostream>
  33 #include <iomanip>
  34 #include <TBrowser.h>
  35 #include <TString.h>
  36 #include <TH1.h>
  37 #include <TH2.h>
  38
  39 //====================================================================
  40 AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order,
  41                                        UShort_t nxbins,
  42                                        Double_t* xbins)
  43   : fXAxis(nxbins, xbins),
  44     fBins(0)
  45 {
  46   // Constructor
  47   // Parameters:
  48   //   Order    Order
  49   //   nxbins   Number of bins
  50   //   xbins    Bin borders
  51   UShort_t n = fXAxis.N();
  52   fBins   = new AliFMDFlowBin*[n];
  53   for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
  54 }
  55 //____________________________________________________________________
  56 AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order,
  57                                        const AliFMDFlowAxis& xaxis)
  58   : fXAxis(xaxis)
  59 {
  60   // Constructor
  61   // Parameters:
  62   //   Order    Order
  63   //   xaxis    X axis object
  64   UShort_t n = fXAxis.N();
  65   fBins   = new AliFMDFlowBin*[n];
  66   for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
  67 }
  68 //____________________________________________________________________
  69 AliFMDFlowBinned1D::AliFMDFlowBinned1D(const AliFMDFlowBinned1D& o)
  70   : TObject(o),
  71     fXAxis(o.fXAxis)
  72 {
  73   // Copy constructor
  74   // Parameters:
  75   //   o   Object to copy from
  76   UShort_t n = fXAxis.N();
  77   fBins   = new AliFMDFlowBin*[n];
  78   for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
  79 }
  80 //____________________________________________________________________
  81 AliFMDFlowBinned1D&
  82 AliFMDFlowBinned1D::operator=(const AliFMDFlowBinned1D& o)
  83 {
  84   // Assignment operator
  85   // Parameters:
  86   //   o Object to assign from
  87   //
  88   // Returns reference to this object
  89   if (fBins) {
  90     for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
  91     delete [] fBins;
  92   }
  93   fXAxis     = o.fXAxis;
  94   UShort_t n = fXAxis.N();
  95   fBins   = new AliFMDFlowBin*[n];
  96   for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
  97   return *this;
  98 }
  99
 100 //____________________________________________________________________
 101 AliFMDFlowBinned1D::~AliFMDFlowBinned1D()
 102 {
 103   // Destructor
 104   // Parameters:
 105   //    none
 106   if (fBins) {
 107     for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
 108     delete [] fBins;
 109   }
 110 }
 111
 112 //____________________________________________________________________
 113 AliFMDFlowBin*
 114 AliFMDFlowBinned1D::GetBin(UShort_t i) const
 115 {
 116   // Get the ith bin
 117   // Parameters:
 118   //    i  Bin number
 119   //
 120   // Return pointer to bin, or null.
 121   if (i >= fXAxis.N()) return 0;
 122   return fBins[i];
 123 }
 124 //____________________________________________________________________
 125 AliFMDFlowBin*
 126 AliFMDFlowBinned1D::GetBin(Double_t x) const
 127 {
 128   // Get the bin that contains x
 129   // Parameters:
 130   //    x  X axis value
 131   //
 132   // Return pointer to bin, or null.
 133   Int_t i = fXAxis.FindBin(x);
 134   if (i < 0) return 0;
 135   UShort_t j = i;
 136   return GetBin(j);
 137 }
 138
 139 //____________________________________________________________________
 140 void
 141 AliFMDFlowBinned1D::Begin()
 142 {
 143   // Called at the beginning of an event
 144   // Parameters:
 145   //   none
 146   for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Begin();
 147 }
 148 //____________________________________________________________________
 149 void
 150 AliFMDFlowBinned1D::End()
 151 {
 152   // Called at the end of an event
 153   // Parameters:
 154   //   none
 155   for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->End();
 156 }
 157 //____________________________________________________________________
 158 void
 159 AliFMDFlowBinned1D::Finish()
 160 {
 161   // Called at the end of an job
 162   // Parameters:
 163   //   none
 164   for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Finish();
 165 }
 166 //____________________________________________________________________
 167 Bool_t
 168 AliFMDFlowBinned1D::AddToEventPlane(Double_t x, Double_t phi,
 169                                     Double_t w, Bool_t a)
 170 {
 171   // Called to add a contribution to the event plane
 172   // Parameters:
 173   //    x   Bin to fill into
 174   //    w   Weight
 175   //    phi The angle phi in radians
 176   //    a   If true, add to sub-event A, otherwise sub-event B
 177   //
 178   // Return false if x falls outside the defined range, true otherwise
 179   AliFMDFlowBin* bin = GetBin(x);
 180   if (!bin) return kFALSE;
 181   bin->AddToEventPlane(phi, w, a);
 182   return kTRUE;
 183 }
 184
 185 //____________________________________________________________________
 186 Bool_t
 187 AliFMDFlowBinned1D::AddToHarmonic(Double_t x, Double_t phi)
 188 {
 189   // Called to add a contribution to the harmonic
 190   // Parameters:
 191   //    x   Bin to fill into
 192   //    phi The angle phi in radians
 193   //
 194   // Return false if x falls outside the defined range, true otherwise
 195   AliFMDFlowBin* bin = GetBin(x);
 196   if (!bin) return kFALSE;
 197   bin->AddToHarmonic(phi);
 198   return kTRUE;
 199 }
 200
 201 //____________________________________________________________________
 202 void
 203 AliFMDFlowBinned1D::Event(Double_t* phis, Double_t* xs, Double_t* ws, ULong_t n)
 204 {
 205   // Process a full event.
 206   // Parameters:
 207   //   phis   List of n phi=[0,2pi] angles
 208   //   xs     List of n x values.
 209   //   ws     Weights
 210   //   n      Size of phis and xs
 211   Begin();
 212   for (UInt_t i = 0; i < n; i++)
 213     AddToEventPlane(xs[i], phis[i], (ws ? ws[i] : 1),
 214                     Float_t(rand()) / RAND_MAX > 0.5);
 215   for (UInt_t i = 0; i < n; i++)
 216     AddToHarmonic(xs[i], phis[i]);
 217   End();
 218 }
 219
 220 //____________________________________________________________________
 221 void
 222 AliFMDFlowBinned1D::Browse(TBrowser* b)
 223 {
 224   // Browse this object.
 225   b->Add(&fXAxis, "xaxis");
 226   for (UInt_t i = 0; i < fXAxis.N(); i++)
 227     b->Add(fBins[i], Form("bin_%03d", i));
 228 }
 229
 230 //____________________________________________________________________
 231 void
 232 AliFMDFlowBinned1D::Draw(Option_t* option)
 233 {
 234   // Draw the distribution of the harmonics or the event plane
 235   // resolution.
 236   // Parameters:
 237   //    option     String of options
 238   //
 239   // Options:
 240   //    b          Draw bare averages of cos(n(phi-Psi))
 241   //    n          Draw harmonic scaled by naive correction
 242   //    s          Draw harmonic scaled by STAR correction
 243   //    t          Draw harmonic scaled by TDR correction
 244   //    r          Draw resolution rather than harmonic
 245   //
 246   // If more than one of b, n, s, or t is given, a 2D histogram is
 247   // drawn.
 248   TString opt(option);
 249   opt.ToLower();
 250   const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
 251   const AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone,
 252                                            AliFMDFlowBin::kNaive,
 253                                            AliFMDFlowBin::kStar,
 254                                            AliFMDFlowBin::kTdr };
 255   Bool_t meths[] = { opt.Contains("b"),
 256                      opt.Contains("n"),
 257                      opt.Contains("s"),
 258                      opt.Contains("t") };
 259   Bool_t res     = opt.Contains("r");
 260   UShort_t nm = 0;
 261   Short_t  sm = -1;
 262   for (UShort_t i = 0; i < 4; i++) { if (meths[i]) { nm++; sm = i; } }
 263   TH1* h = 0;
 264   if (nm > 1) {
 265     h = new TH2D((res ? "res" : "flow"), (res ? "Resolution" : "Flow"),
 266                  fXAxis.N(), fXAxis.Bins(), nm, 0, nm);
 267     h->SetXTitle("x");
 268     h->SetYTitle("method");
 269     h->GetYaxis()->SetNdivisions(nm+1, kFALSE);
 270     h->SetZTitle((res ? "<cos(n(#Psi_{m}-#Psi_{R}))>" : "v"));
 271     UInt_t j = 0;
 272     for (UShort_t i = 0; i < 4; i++)
 273       if (meths[i]) h->GetYaxis()->SetBinLabel(++j, names[i]);
 274   }
 275   else {
 276     h = new TH1D(Form("%s_%s", (res ? "res" : "flow"), names[sm]),
 277                  Form("%s_%s", (res ? "Resolution" : "Flow"), names[sm]),
 278                  fXAxis.N(), fXAxis.Bins());
 279     h->SetXTitle("x");
 280     h->SetYTitle((res ? "<cos(n(#Psi_{m}-#Psi_{R}))>" : "v"));
 281   }
 282
 283   for (UShort_t i = 0; i < fXAxis.N(); i++) {
 284     Double_t       x   = fXAxis.BinCenter(i);
 285     AliFMDFlowBin* bin = GetBin(x);
 286     Double_t       v, e2;
 287     if (nm == 1) {
 288       if (res) v = bin->Correction(e2, types[sm]);
 289       else     v = bin->Value(e2, types[sm]);
 290       h->SetBinContent(i+1, v);
 291       h->SetBinError(i+1, sqrt(e2));
 292       continue;
 293     }
 294     UInt_t j = 0;
 295     for (UShort_t k = 0; k < 4; k++)  {
 296       if (!meths[k]) continue;
 297       if (res) v = bin->Correction(e2, types[k]);
 298       else     v = bin->Value(e2, types[k]);
 299       h->SetBinContent(i+1, j+1, v);
 300       h->SetBinError(i+1, j+1, sqrt(e2));
 301       j++;
 302     }
 303   }
 304   h->Draw(option);
 305 }
 306
 307
 308
 309 //____________________________________________________________________
 310 void
 311 AliFMDFlowBinned1D::Print(Option_t* option) const
 312 {
 313   // Print information to standard output.
 314   // Parameters:
 315   //     option    String of options
 316   //
 317   // Options:
 318   //     d         Print details
 319   //     s         Print summary
 320   //
 321   TString opt(option);
 322   opt.ToLower();
 323   Bool_t det = opt.Contains("d");
 324   Bool_t sum = opt.Contains("s");
 325   if (det) {
 326     for (UShort_t i = 0; i < fXAxis.N(); i++) {
 327       Double_t x = fXAxis.BinCenter(i);
 328       std::streamsize         oldP = std::cout.precision(3);
 329       std::ios_base::fmtflags oldF = std::cout.setf(std::ios_base::fixed,
 330                                                      std::ios_base::floatfield);
 331       std::cout << "x=" << std::setw(5) << x << std::endl;
 332       fBins[i]->Print();
 333       std::cout.precision(oldP);
 334       std::cout.setf(oldF, std::ios_base::floatfield);
 335     }
 336   }
 337
 338   if (sum) {
 339     UInt_t       nType = 4;
 340     const char*  names[] = { "Bare",    "Naive",    "STAR",    "TDR" };
 341     AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone,
 342                                        AliFMDFlowBin::kNaive,
 343                                        AliFMDFlowBin::kStar,
 344                                        AliFMDFlowBin::kTdr };
 345     std::cout << "    x";
 346     for (UInt_t i = 0; i < nType; i++)
 347       std::cout << " | " << std::setw(6+6+5) << names[i];
 348     std::cout << "\n-----" << std::setfill('-');
 349     for (UInt_t i = 0; i < nType; i++)
 350       std::cout << "-+-" <<  std::setw(6+6+5) << "-";
 351     std::cout << std::setfill(' ') << std::endl;
 352
 353     std::streamsize         oldP = std::cout.precision(2);
 354     std::ios_base::fmtflags oldF = std::cout.setf(std::ios_base::fixed,
 355                                                   std::ios_base::floatfield);
 356     for (UShort_t i = 0; i < fXAxis.N(); i++) {
 357       Double_t x = fXAxis.BinCenter(i);
 358       std::cout << std::setprecision(2) << std::setw(5) << x << std::flush;
 359       for (UShort_t j = 0; j < nType; j++) {
 360         Double_t e2v;
 361         Double_t v  = fBins[i]->Value(e2v, types[j]);
 362         std::cout << std::setprecision(3)    << " | "
 363                   << std::setw(6) << 100 * v << " +/- "
 364                   << std::setw(6) << 100 * sqrt(e2v);
 365       }
 366       std::cout << std::endl;
 367     }
 368     std::cout.precision(oldP);
 369     std::cout.setf(oldF, std::ios_base::floatfield);
 370   }
 371 }
 372
 373
 374 //____________________________________________________________________
 375 //
 376 // EOF
 377 //