bug fix
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowBinned1D.cxx
CommitLineData
97e94238 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 */
39eefe19 18/** @file
19 @brief Implementation of a 1-dimensional Flow "histogram" */
97e94238 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
39eefe19 28#include "flow/AliFMDFlowBinned1D.h"
29#include "flow/AliFMDFlowBin.h"
9b98d361 30#include "flow/AliFMDFlowSplitter.h"
31// #include <cmath>
39eefe19 32#include <cstdlib>
33#include <iostream>
34#include <iomanip>
35#include <TBrowser.h>
36#include <TString.h>
824e71c1 37#include <TH1.h>
38#include <TH2.h>
39eefe19 39
40//====================================================================
9b98d361 41AliFMDFlowBinned1D::AliFMDFlowBinned1D()
42 : TNamed("", ""),
43 fXAxis(),
44 fN(0),
45 fBins(0),
46 fSplitter(0)
47{
48 // Default CTOR - do not use
49}
50
51//____________________________________________________________________
52AliFMDFlowBinned1D::AliFMDFlowBinned1D(const char* name,
53 const char* title,
54 UShort_t order,
55 UShort_t k,
39eefe19 56 UShort_t nxbins,
9b98d361 57 Double_t* xbins,
58 AliFMDFlowSplitter* splitter)
59 : TNamed(name,title),
60 fXAxis(nxbins, xbins),
6ce810fc 61 fN(0),
9b98d361 62 fBins(0),
63 fSplitter(splitter)
39eefe19 64{
97e94238 65 // Constructor
66 // Parameters:
67 // Order Order
68 // nxbins Number of bins
69 // xbins Bin borders
39eefe19 70 UShort_t n = fXAxis.N();
9b98d361 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;
39eefe19 75}
76//____________________________________________________________________
9b98d361 77AliFMDFlowBinned1D::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),
6ce810fc 85 fN(0),
9b98d361 86 fBins(0),
87 fSplitter(splitter)
39eefe19 88{
97e94238 89 // Constructor
90 // Parameters:
91 // Order Order
92 // xaxis X axis object
39eefe19 93 UShort_t n = fXAxis.N();
9b98d361 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
39eefe19 99}
100//____________________________________________________________________
101AliFMDFlowBinned1D::AliFMDFlowBinned1D(const AliFMDFlowBinned1D& o)
9b98d361 102 : TNamed(o),
103 TAttLine(o),
104 TAttFill(o),
105 TAttMarker(o),
106 fXAxis(o.fXAxis),
6ce810fc 107 fN(0),
108 fBins(0),
9b98d361 109 fSplitter(0)
39eefe19 110{
97e94238 111 // Copy constructor
112 // Parameters:
113 // o Object to copy from
39eefe19 114 UShort_t n = fXAxis.N();
9b98d361 115 fSplitter = new AliFMDFlowShuffle;
116 fN = n;
117 fBins = new AliFMDFlowBin*[n];
39eefe19 118 for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
119}
120//____________________________________________________________________
121AliFMDFlowBinned1D&
122AliFMDFlowBinned1D::operator=(const AliFMDFlowBinned1D& o)
123{
97e94238 124 // Assignment operator
125 // Parameters:
126 // o Object to assign from
127 //
128 // Returns reference to this object
9b98d361 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());
39eefe19 139 if (fBins) {
140 for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
141 delete [] fBins;
142 }
9b98d361 143 if (fSplitter) delete fSplitter;
39eefe19 144 fXAxis = o.fXAxis;
145 UShort_t n = fXAxis.N();
9b98d361 146 fN = n;
147 fSplitter = new AliFMDFlowShuffle;
148 fBins = new AliFMDFlowBin*[n];
39eefe19 149 for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
150 return *this;
151}
152
153//____________________________________________________________________
154AliFMDFlowBinned1D::~AliFMDFlowBinned1D()
155{
97e94238 156 // Destructor
157 // Parameters:
158 // none
39eefe19 159 if (fBins) {
160 for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
161 delete [] fBins;
162 }
9b98d361 163 if (fSplitter) delete fSplitter;
39eefe19 164}
165
166//____________________________________________________________________
167AliFMDFlowBin*
168AliFMDFlowBinned1D::GetBin(UShort_t i) const
169{
97e94238 170 // Get the ith bin
171 // Parameters:
172 // i Bin number
173 //
174 // Return pointer to bin, or null.
39eefe19 175 if (i >= fXAxis.N()) return 0;
176 return fBins[i];
177}
178//____________________________________________________________________
179AliFMDFlowBin*
180AliFMDFlowBinned1D::GetBin(Double_t x) const
181{
97e94238 182 // Get the bin that contains x
183 // Parameters:
184 // x X axis value
185 //
186 // Return pointer to bin, or null.
39eefe19 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//____________________________________________________________________
194void
195AliFMDFlowBinned1D::Begin()
196{
97e94238 197 // Called at the beginning of an event
198 // Parameters:
199 // none
39eefe19 200 for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Begin();
9b98d361 201 fSplitter->Begin();
39eefe19 202}
203//____________________________________________________________________
204void
205AliFMDFlowBinned1D::End()
206{
97e94238 207 // Called at the end of an event
208 // Parameters:
209 // none
39eefe19 210 for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->End();
9b98d361 211 fSplitter->End();
39eefe19 212}
213//____________________________________________________________________
214void
215AliFMDFlowBinned1D::Finish()
216{
97e94238 217 // Called at the end of an job
218 // Parameters:
219 // none
39eefe19 220 for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Finish();
221}
222//____________________________________________________________________
223Bool_t
97e94238 224AliFMDFlowBinned1D::AddToEventPlane(Double_t x, Double_t phi,
225 Double_t w, Bool_t a)
39eefe19 226{
97e94238 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
39eefe19 235 AliFMDFlowBin* bin = GetBin(x);
236 if (!bin) return kFALSE;
237 bin->AddToEventPlane(phi, w, a);
238 return kTRUE;
239}
240
241//____________________________________________________________________
242Bool_t
9b98d361 243AliFMDFlowBinned1D::AddToHarmonic(Double_t x, Double_t phi,
244 Double_t wp, Double_t wh)
39eefe19 245{
97e94238 246 // Called to add a contribution to the harmonic
247 // Parameters:
9b98d361 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
97e94238 253 // Return false if x falls outside the defined range, true otherwise
39eefe19 254 AliFMDFlowBin* bin = GetBin(x);
255 if (!bin) return kFALSE;
9b98d361 256 bin->AddToHarmonic(phi, wp, wh);
39eefe19 257 return kTRUE;
258}
259
260//____________________________________________________________________
9b98d361 261UShort_t
262AliFMDFlowBinned1D::Order() const
263{
264 return GetBin(UShort_t(0))->Order();
265}
266//____________________________________________________________________
267UShort_t
268AliFMDFlowBinned1D::PsiOrder() const
269{
270 return GetBin(UShort_t(0))->Order();
271}
272
273//____________________________________________________________________
39eefe19 274void
9b98d361 275AliFMDFlowBinned1D::Event(ULong_t n, Double_t* phis, Double_t* xs,
276 Double_t* wp, Double_t* wh)
39eefe19 277{
97e94238 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
39eefe19 284 Begin();
9b98d361 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));
39eefe19 290 End();
291}
292
293//____________________________________________________________________
294void
295AliFMDFlowBinned1D::Browse(TBrowser* b)
296{
97e94238 297 // Browse this object.
39eefe19 298 b->Add(&fXAxis, "xaxis");
9b98d361 299 if (fSplitter) b->Add(fSplitter, "Splitter");
39eefe19 300 for (UInt_t i = 0; i < fXAxis.N(); i++)
301 b->Add(fBins[i], Form("bin_%03d", i));
302}
303
304//____________________________________________________________________
9b98d361 305TH1*
306AliFMDFlowBinned1D::MakeHistogram(UInt_t which, UInt_t what)
824e71c1 307{
9b98d361 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
824e71c1 320 UShort_t nm = 0;
321 Short_t sm = -1;
9b98d361 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
824e71c1 330 TH1* h = 0;
331 if (nm > 1) {
9b98d361 332 // Make 2D histogram
333 h = new TH2D(Form("%s_%s", GetName(), hNames[what]),
334 Form("%s %s", GetTitle(), hTitles[what]),
824e71c1 335 fXAxis.N(), fXAxis.Bins(), nm, 0, nm);
9b98d361 336
337 // Set titles and such
824e71c1 338 h->SetXTitle("x");
339 h->SetYTitle("method");
9b98d361 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 }
824e71c1 348 h->GetYaxis()->SetNdivisions(nm+1, kFALSE);
9b98d361 349
350 // Set Bin labels for the methods
824e71c1 351 UInt_t j = 0;
352 for (UShort_t i = 0; i < 4; i++)
9b98d361 353 if (which & (1 << i)) h->GetYaxis()->SetBinLabel(++j, names[i]);
824e71c1 354 }
355 else {
9b98d361 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());
824e71c1 363 h->SetXTitle("x");
9b98d361 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 }
824e71c1 372 }
373
374 for (UShort_t i = 0; i < fXAxis.N(); i++) {
375 Double_t x = fXAxis.BinCenter(i);
376 AliFMDFlowBin* bin = GetBin(x);
9b98d361 377 Double_t y=0, e2=0, dummy;
378 if (bin->Counts() <= 0) continue;
824e71c1 379 if (nm == 1) {
9b98d361 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);
824e71c1 390 h->SetBinError(i+1, sqrt(e2));
391 continue;
392 }
393 UInt_t j = 0;
394 for (UShort_t k = 0; k < 4; k++) {
9b98d361 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);
824e71c1 406 h->SetBinError(i+1, j+1, sqrt(e2));
407 j++;
408 }
409 }
9b98d361 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//____________________________________________________________________
424void
425AliFMDFlowBinned1D::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());
824e71c1 468}
469
470
471
472//____________________________________________________________________
39eefe19 473void
474AliFMDFlowBinned1D::Print(Option_t* option) const
475{
97e94238 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 //
39eefe19 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);
97e94238 491 std::streamsize oldP = std::cout.precision(3);
492 std::ios_base::fmtflags oldF = std::cout.setf(std::ios_base::fixed,
39eefe19 493 std::ios_base::floatfield);
494 std::cout << "x=" << std::setw(5) << x << std::endl;
495 fBins[i]->Print();
97e94238 496 std::cout.precision(oldP);
497 std::cout.setf(oldF, std::ios_base::floatfield);
39eefe19 498 }
499 }
500
501 if (sum) {
502 UInt_t nType = 4;
503 const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
97e94238 504 AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone,
505 AliFMDFlowBin::kNaive,
506 AliFMDFlowBin::kStar,
507 AliFMDFlowBin::kTdr };
9b98d361 508 std::cout << GetName() << " - " << GetTitle() << "\n"
509 << " x";
39eefe19 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
97e94238 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);
39eefe19 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 }
97e94238 532 std::cout.precision(oldP);
533 std::cout.setf(oldF, std::ios_base::floatfield);
39eefe19 534 }
535}
536
537
538//____________________________________________________________________
539//
540// EOF
541//