bug fix
[u/mrichter/AliRoot.git] / 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 //