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 "flow/AliFMDFlowSplitter.h"
  31 // #include <cmath>
  32 #include <cstdlib>
  33 #include <iostream>
  34 #include <iomanip>
  35 #include <TBrowser.h>
  36 #include <TString.h>
  37 #include <TH1.h>
  38 #include <TH2.h>
  39
  40 //====================================================================
  41 AliFMDFlowBinned1D::AliFMDFlowBinned1D()
  42   : TNamed("", ""),
  43     fXAxis(),
  44     fN(0),
  45     fBins(0),
  46     fSplitter(0)
  47 {
  48   // Default CTOR - do not use
  49 }
  50
  51 //____________________________________________________________________
  52 AliFMDFlowBinned1D::AliFMDFlowBinned1D(const char* name,
  53                                        const char* title,
  54                                        UShort_t order,
  55                                        UShort_t k,
  56                                        UShort_t nxbins,
  57                                        Double_t* xbins,
  58                                        AliFMDFlowSplitter* splitter)
  59   : TNamed(name,title),
  60     fXAxis(nxbins, xbins),
  61     fN(0),
  62     fBins(0),
  63     fSplitter(splitter)
  64 {
  65   // Constructor
  66   // Parameters:
  67   //   Order    Order
  68   //   nxbins   Number of bins
  69   //   xbins    Bin borders
  70   UShort_t n = fXAxis.N();
  71   fN         = n;
  72   fBins      = new AliFMDFlowBin*[n];
  73   for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order,k);
  74   if (!fSplitter) fSplitter = new AliFMDFlowShuffle;
  75 }
  76 //____________________________________________________________________
  77 AliFMDFlowBinned1D::AliFMDFlowBinned1D(const char*           name,
  78                                        const char*           title,
  79                                        UShort_t              order,
  80                                        UShort_t              k,
  81                                        const AliFMDFlowAxis& xaxis,
  82                                        AliFMDFlowSplitter*   splitter)
  83   : TNamed(name,title),
  84     fXAxis(xaxis),
  85     fN(0),
  86     fBins(0),
  87     fSplitter(splitter)
  88 {
  89   // Constructor
  90   // Parameters:
  91   //   Order    Order
  92   //   xaxis    X axis object
  93   UShort_t n = fXAxis.N();
  94   fN         = n;
  95   fBins      = new AliFMDFlowBin*[n];
  96   for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order,k);
  97   if (!fSplitter) fSplitter = new AliFMDFlowShuffle;
  98
  99 }
 100 //____________________________________________________________________
 101 AliFMDFlowBinned1D::AliFMDFlowBinned1D(const AliFMDFlowBinned1D& o)
 102   : TNamed(o),
 103     TAttLine(o),
 104     TAttFill(o),
 105     TAttMarker(o),
 106     fXAxis(o.fXAxis),
 107     fN(0),
 108     fBins(0),
 109     fSplitter(0)
 110 {
 111   // Copy constructor
 112   // Parameters:
 113   //   o   Object to copy from
 114   UShort_t n = fXAxis.N();
 115   fSplitter  = new AliFMDFlowShuffle;
 116   fN         = n;
 117   fBins      = new AliFMDFlowBin*[n];
 118   for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
 119 }
 120 //____________________________________________________________________
 121 AliFMDFlowBinned1D&
 122 AliFMDFlowBinned1D::operator=(const AliFMDFlowBinned1D& o)
 123 {
 124   // Assignment operator
 125   // Parameters:
 126   //   o Object to assign from
 127   //
 128   // Returns reference to this object
 129   SetLineColor(o.GetLineColor());
 130   SetLineStyle(o.GetLineStyle());
 131   SetLineWidth(o.GetLineWidth());
 132   SetFillColor(o.GetFillColor());
 133   SetFillStyle(o.GetFillStyle());
 134   SetMarkerColor(o.GetMarkerColor());
 135   SetMarkerStyle(o.GetMarkerStyle());
 136   SetMarkerSize(o.GetMarkerSize());
 137   this->SetName(o.GetName());
 138   this->SetTitle(o.GetTitle());
 139   if (fBins) {
 140     for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
 141     delete [] fBins;
 142   }
 143   if (fSplitter) delete fSplitter;
 144   fXAxis     = o.fXAxis;
 145   UShort_t n = fXAxis.N();
 146   fN         = n;
 147   fSplitter  = new AliFMDFlowShuffle;
 148   fBins      = new AliFMDFlowBin*[n];
 149   for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
 150   return *this;
 151 }
 152
 153 //____________________________________________________________________
 154 AliFMDFlowBinned1D::~AliFMDFlowBinned1D()
 155 {
 156   // Destructor
 157   // Parameters:
 158   //    none
 159   if (fBins) {
 160     for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
 161     delete [] fBins;
 162   }
 163   if (fSplitter) delete fSplitter;
 164 }
 165
 166 //____________________________________________________________________
 167 AliFMDFlowBin*
 168 AliFMDFlowBinned1D::GetBin(UShort_t i) const
 169 {
 170   // Get the ith bin
 171   // Parameters:
 172   //    i  Bin number
 173   //
 174   // Return pointer to bin, or null.
 175   if (i >= fXAxis.N()) return 0;
 176   return fBins[i];
 177 }
 178 //____________________________________________________________________
 179 AliFMDFlowBin*
 180 AliFMDFlowBinned1D::GetBin(Double_t x) const
 181 {
 182   // Get the bin that contains x
 183   // Parameters:
 184   //    x  X axis value
 185   //
 186   // Return pointer to bin, or null.
 187   Int_t i = fXAxis.FindBin(x);
 188   if (i < 0) return 0;
 189   UShort_t j = i;
 190   return GetBin(j);
 191 }
 192
 193 //____________________________________________________________________
 194 void
 195 AliFMDFlowBinned1D::Begin()
 196 {
 197   // Called at the beginning of an event
 198   // Parameters:
 199   //   none
 200   for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Begin();
 201   fSplitter->Begin();
 202 }
 203 //____________________________________________________________________
 204 void
 205 AliFMDFlowBinned1D::End()
 206 {
 207   // Called at the end of an event
 208   // Parameters:
 209   //   none
 210   for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->End();
 211   fSplitter->End();
 212 }
 213 //____________________________________________________________________
 214 void
 215 AliFMDFlowBinned1D::Finish()
 216 {
 217   // Called at the end of an job
 218   // Parameters:
 219   //   none
 220   for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Finish();
 221 }
 222 //____________________________________________________________________
 223 Bool_t
 224 AliFMDFlowBinned1D::AddToEventPlane(Double_t x, Double_t phi,
 225                                     Double_t w, Bool_t a)
 226 {
 227   // Called to add a contribution to the event plane
 228   // Parameters:
 229   //    x   Bin to fill into
 230   //    w   Weight
 231   //    phi The angle phi in radians
 232   //    a   If true, add to sub-event A, otherwise sub-event B
 233   //
 234   // Return false if x falls outside the defined range, true otherwise
 235   AliFMDFlowBin* bin = GetBin(x);
 236   if (!bin) return kFALSE;
 237   bin->AddToEventPlane(phi, w, a);
 238   return kTRUE;
 239 }
 240
 241 //____________________________________________________________________
 242 Bool_t
 243 AliFMDFlowBinned1D::AddToHarmonic(Double_t x,  Double_t phi,
 244                                   Double_t wp, Double_t wh)
 245 {
 246   // Called to add a contribution to the harmonic
 247   // Parameters:
 248   //    x    Bin to fill into
 249   //    phi  The angle phi in radians
 250   //    wp   weight of event plane
 251   //    wh   weight of harmonic
 252
 253   // Return false if x falls outside the defined range, true otherwise
 254   AliFMDFlowBin* bin = GetBin(x);
 255   if (!bin) return kFALSE;
 256   bin->AddToHarmonic(phi, wp, wh);
 257   return kTRUE;
 258 }
 259
 260 //____________________________________________________________________
 261 UShort_t
 262 AliFMDFlowBinned1D::Order() const
 263 {
 264   return GetBin(UShort_t(0))->Order();
 265 }
 266 //____________________________________________________________________
 267 UShort_t
 268 AliFMDFlowBinned1D::PsiOrder() const
 269 {
 270   return GetBin(UShort_t(0))->Order();
 271 }
 272
 273 //____________________________________________________________________
 274 void
 275 AliFMDFlowBinned1D::Event(ULong_t   n,  Double_t* phis, Double_t* xs,
 276                           Double_t* wp, Double_t* wh)
 277 {
 278   // Process a full event.
 279   // Parameters:
 280   //   phis   List of n phi=[0,2pi] angles
 281   //   xs     List of n x values.
 282   //   ws     Weights
 283   //   n      Size of phis and xs
 284   Begin();
 285   fSplitter->Event(phis, xs, n);
 286   for (ULong_t i = 0; i < n; i++)
 287     AddToEventPlane(xs[i], phis[i], (wp ? wp[i] : 1), fSplitter->Select(i));
 288   for (ULong_t i = 0; i < n; i++)
 289     AddToHarmonic(xs[i], phis[i], (wp ? wp[i] : 1), (wh ? wh[i] : 1));
 290   End();
 291 }
 292
 293 //____________________________________________________________________
 294 void
 295 AliFMDFlowBinned1D::Browse(TBrowser* b)
 296 {
 297   // Browse this object.
 298   b->Add(&fXAxis, "xaxis");
 299   if (fSplitter) b->Add(fSplitter, "Splitter");
 300   for (UInt_t i = 0; i < fXAxis.N(); i++)
 301     b->Add(fBins[i], Form("bin_%03d", i));
 302 }
 303
 304 //____________________________________________________________________
 305 TH1*
 306 AliFMDFlowBinned1D::MakeHistogram(UInt_t which, UInt_t what)
 307 {
 308   // Make a histogram of some of the stuff in this object
 309
 310   // Some strings
 311   const char* names[]   = { "Bare",  "Naive", "STAR", "TDR" };
 312   const char* hNames[]  = { "flow",  "resolution", "counts" };
 313   const char* hTitles[] = { "Flow",  "Resolution", "Counts" };
 314   const char* hUnits[]  = { "v_{%d}","<cos(%d(#Psi_{%d}-#Psi_{R}))>", "N" };
 315   AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone,
 316                                      AliFMDFlowBin::kNaive,
 317                                      AliFMDFlowBin::kStar,
 318                                      AliFMDFlowBin::kTdr };
 319   // Figure out how many things to draw
 320   UShort_t nm = 0;
 321   Short_t  sm = -1;
 322   for (UShort_t i = 0; i < 4; i++) {
 323     if ((which & (1 << i)) != 0) {
 324       nm++;
 325       sm = i;
 326     }
 327   }
 328   if (what == AliFMDFlowBin::kCounts) { nm = 1; sm = 0; }
 329
 330   TH1* h = 0;
 331   if (nm > 1) {
 332     // Make 2D histogram
 333     h = new TH2D(Form("%s_%s", GetName(), hNames[what]),
 334                  Form("%s %s", GetTitle(), hTitles[what]),
 335                  fXAxis.N(), fXAxis.Bins(), nm, 0, nm);
 336
 337     // Set titles and such
 338     h->SetXTitle("x");
 339     h->SetYTitle("method");
 340     switch (what) {
 341     case AliFMDFlowBin::kHarmonic:
 342       h->SetZTitle(Form(hUnits[what], Order())); break;
 343     case AliFMDFlowBin::kResolution:
 344       h->SetZTitle(Form(hUnits[what], Order(), PsiOrder())); break;
 345     default:
 346       h->SetZTitle(hUnits[what]); break;
 347     }
 348     h->GetYaxis()->SetNdivisions(nm+1, kFALSE);
 349
 350     // Set Bin labels for the methods
 351     UInt_t j = 0;
 352     for (UShort_t i = 0; i < 4; i++)
 353       if (which & (1 << i)) h->GetYaxis()->SetBinLabel(++j, names[i]);
 354   }
 355   else {
 356     TString name(what == AliFMDFlowBin::kCounts ?
 357                  Form("%s_%s", GetName(), hNames[what]) :
 358                  Form("%s_%s_%s", GetName(), hNames[what], names[sm]));
 359     TString title(what == AliFMDFlowBin::kCounts ?
 360                   Form("%s %s", GetTitle(), hTitles[what]) :
 361                   Form("%s %s %s", GetTitle(), hTitles[what], names[sm]));
 362     h = new TH1D(name.Data(), title.Data(), fXAxis.N(), fXAxis.Bins());
 363     h->SetXTitle("x");
 364     switch (what) {
 365     case AliFMDFlowBin::kHarmonic:
 366       h->SetYTitle(Form(hUnits[what], Order())); break;
 367     case AliFMDFlowBin::kResolution:
 368       h->SetYTitle(Form(hUnits[what], Order(), PsiOrder())); break;
 369     default:
 370       h->SetYTitle(hUnits[what]); break;
 371     }
 372   }
 373
 374   for (UShort_t i = 0; i < fXAxis.N(); i++) {
 375     Double_t       x   = fXAxis.BinCenter(i);
 376     AliFMDFlowBin* bin = GetBin(x);
 377     Double_t       y=0, e2=0, dummy;
 378     if (bin->Counts() <= 0) continue;
 379     if (nm == 1) {
 380       switch (what) {
 381       case AliFMDFlowBin::kHarmonic:
 382         if (bin->Correction(dummy, types[sm]) < .01) continue;
 383         y = bin->Value(e2, types[sm]); break;
 384       case AliFMDFlowBin::kResolution:
 385         y = bin->Correction(e2, types[sm]); break;
 386       case AliFMDFlowBin::kCounts:
 387         y = bin->Counts(); break;
 388       }
 389       h->SetBinContent(i+1, y);
 390       h->SetBinError(i+1, sqrt(e2));
 391       continue;
 392     }
 393     UInt_t j = 0;
 394     for (UShort_t k = 0; k < 4; k++)  {
 395       if (!(which & (1 << k))) continue;
 396       switch (what) {
 397       case AliFMDFlowBin::kHarmonic:
 398         if (bin->Correction(dummy,types[k]) < .01) continue;
 399         y = bin->Value(e2, types[k]); break;
 400       case AliFMDFlowBin::kResolution:
 401         y = bin->Correction(e2, types[k]); break;
 402       case AliFMDFlowBin::kCounts:
 403         y = bin->Counts(); break;
 404       }
 405       h->SetBinContent(i+1, j+1, y);
 406       h->SetBinError(i+1, j+1, sqrt(e2));
 407       j++;
 408     }
 409   }
 410   h->SetLineColor(GetLineColor());
 411   h->SetLineWidth(GetLineWidth());
 412   h->SetLineStyle(GetLineStyle());
 413   h->SetFillColor(GetFillColor());
 414   h->SetFillStyle(GetFillStyle());
 415   h->SetMarkerColor(GetMarkerColor());
 416   h->SetMarkerStyle(GetMarkerStyle());
 417   h->SetMarkerSize(GetMarkerSize());
 418   h->SetDirectory(0);
 419   return h;
 420 }
 421
 422
 423 //____________________________________________________________________
 424 void
 425 AliFMDFlowBinned1D::Draw(Option_t* option)
 426 {
 427   // Draw the distribution of the harmonics or the event plane
 428   // resolution.
 429   // Parameters:
 430   //    option     String of options
 431   //
 432   // Options:
 433   //    b          Draw bare averages of cos(n(phi-Psi))
 434   //    n          Draw harmonic scaled by naive correction
 435   //    s          Draw harmonic scaled by STAR correction
 436   //    t          Draw harmonic scaled by TDR correction (*)
 437   //    r          Draw resolution rather than harmonic
 438   //    c          Draw counts rather than harmonic
 439   //    h          Draw harmonics (*)
 440   //
 441   // (*) This is the default value
 442   //
 443   // If more than one of b, n, s, or t is given, a 2D histogram is
 444   // drawn.
 445   //
 446   // Only one of r, c, h can be specified.  The first specified wins.
 447   //
 448   TString opt(option);
 449   opt.ToLower();
 450   TString dopt = "";
 451   Int_t idx = opt.Index(":");
 452   if (idx != kNPOS) {
 453     dopt = opt(idx+1, opt.Length()-idx-1);
 454     opt.Remove(idx, opt.Length()-idx);
 455   }
 456   UInt_t which = 0;
 457   if (opt.Contains("b")) which |= AliFMDFlowBin::kNone;
 458   if (opt.Contains("n")) which |= AliFMDFlowBin::kNaive;
 459   if (opt.Contains("s")) which |= AliFMDFlowBin::kStar;
 460   if (opt.Contains("t")) which |= AliFMDFlowBin::kTdr;
 461
 462   UInt_t what = AliFMDFlowBin::kHarmonic;
 463   if (opt.Contains("c")) what = AliFMDFlowBin::kCounts;
 464   if (opt.Contains("r")) what = AliFMDFlowBin::kResolution;
 465
 466   TH1* h = MakeHistogram(which, what);
 467   h->Draw(dopt.Data());
 468 }
 469
 470
 471
 472 //____________________________________________________________________
 473 void
 474 AliFMDFlowBinned1D::Print(Option_t* option) const
 475 {
 476   // Print information to standard output.
 477   // Parameters:
 478   //     option    String of options
 479   //
 480   // Options:
 481   //     d         Print details
 482   //     s         Print summary
 483   //
 484   TString opt(option);
 485   opt.ToLower();
 486   Bool_t det = opt.Contains("d");
 487   Bool_t sum = opt.Contains("s");
 488   if (det) {
 489     for (UShort_t i = 0; i < fXAxis.N(); i++) {
 490       Double_t x = fXAxis.BinCenter(i);
 491       std::streamsize         oldP = std::cout.precision(3);
 492       std::ios_base::fmtflags oldF = std::cout.setf(std::ios_base::fixed,
 493                                                      std::ios_base::floatfield);
 494       std::cout << "x=" << std::setw(5) << x << std::endl;
 495       fBins[i]->Print();
 496       std::cout.precision(oldP);
 497       std::cout.setf(oldF, std::ios_base::floatfield);
 498     }
 499   }
 500
 501   if (sum) {
 502     UInt_t       nType = 4;
 503     const char*  names[] = { "Bare",    "Naive",    "STAR",    "TDR" };
 504     AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone,
 505                                        AliFMDFlowBin::kNaive,
 506                                        AliFMDFlowBin::kStar,
 507                                        AliFMDFlowBin::kTdr };
 508     std::cout << GetName() << " - " << GetTitle() << "\n"
 509               << "    x";
 510     for (UInt_t i = 0; i < nType; i++)
 511       std::cout << " | " << std::setw(6+6+5) << names[i];
 512     std::cout << "\n-----" << std::setfill('-');
 513     for (UInt_t i = 0; i < nType; i++)
 514       std::cout << "-+-" <<  std::setw(6+6+5) << "-";
 515     std::cout << std::setfill(' ') << std::endl;
 516
 517     std::streamsize         oldP = std::cout.precision(2);
 518     std::ios_base::fmtflags oldF = std::cout.setf(std::ios_base::fixed,
 519                                                   std::ios_base::floatfield);
 520     for (UShort_t i = 0; i < fXAxis.N(); i++) {
 521       Double_t x = fXAxis.BinCenter(i);
 522       std::cout << std::setprecision(2) << std::setw(5) << x << std::flush;
 523       for (UShort_t j = 0; j < nType; j++) {
 524         Double_t e2v;
 525         Double_t v  = fBins[i]->Value(e2v, types[j]);
 526         std::cout << std::setprecision(3)    << " | "
 527                   << std::setw(6) << 100 * v << " +/- "
 528                   << std::setw(6) << 100 * sqrt(e2v);
 529       }
 530       std::cout << std::endl;
 531     }
 532     std::cout.precision(oldP);
 533     std::cout.setf(oldF, std::ios_base::floatfield);
 534   }
 535 }
 536
 537
 538 //____________________________________________________________________
 539 //
 540 // EOF
 541 //