coverity 15108 fixed
[u/mrichter/AliRoot.git] / 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 //