4a1028f6408c1c8f8fe3784c3f34ec69bc8356d3
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowAnalyser.cxx
1 //////////////////////////////////////////////////////////////////////
2 //
3 // $Id$
4 //
5 // Author: Emanuele Simili
6 //
7 //////////////////////////////////////////////////////////////////////
8 //_____________________________________________________________
9 //
10 // Description: 
11 //         the AliFlowAnalyser provides the methods to perform an Event 
12 // plane flow analysis over AliFlowEvents. 
13 // - The method Init() generates the histograms.
14 // - The method Analyse(AliFlowEvent*) calculates/extracts flow observables,  
15 //  fills some histograms and performs the E.P. analysis.
16 // - The method Resolution() calculates the resolution correction and 
17 //  applies it to observed v_n values.
18 // - The method Finish() saves the histograms and closes the file.
19 // Weights calculation has been moved to the class AliFlowWeighter.
20 //
21 // The most important features of the Flow Analysis are:
22 //  - Reaction Plane calculation: for different harmonics and selections
23 //   with user defined cuts. (The calculation of the R.P. is a method in
24 //   the AliFlowEvent class).
25 //  - Correlation Analysis: particles are related to the R.P. and Vn 
26 //   coefficient are extracted. The AliFlowAnalysisMaker also removes 
27 //   auto-correlations of each particle with respect to the event plane, 
28 //   by subtracting track by track from the Q vector before the correlation 
29 //   study. Observable Flow is plotted.
30 //  - Differential Flow: for each harmonic and selection, flow values 
31 //   (Vn) are plotted versus rapidity/pseudorapidity and pT. 
32 //  - Resulution Estimate: event plane resolution is estimated via the 
33 //   sub-event method (res ~ sqrt(cos(psi1-psi2))), and flow values are 
34 //   then corrected for it. "True" Flow is plotted. 
35 //  - Integrated Flow: for the doubly-integrated (pT & eta) Vn values 
36 //   the error on the event plane resolution is folded in to the error. 
37 //  - Phi and Baiesian weights calculation (from event observables). 
38 //
39 // The AliFlowAnalysisMaker Class is adapted from the original 
40 // StFlowAnalysisMaker, succesfully used to study flow in the 
41 // STAR experiment at RICH.
42 // Original Authors:                 Raimond Snellings & Art Poskanzer
43 //
44 //////////////////////////////////////////////////////////////////////
45
46 #ifndef ALIFLOWANALYSER_CXX
47 #define ALIFLOWANALYSER_CXX
48
49
50 // ROOT things
51 #include <TROOT.h>
52 #include <TMath.h>
53 #include <TFile.h>
54 #include <TTree.h>
55 #include <TChain.h>
56 #include <TString.h>
57 #include <TObject.h>
58 #include <TClonesArray.h>
59 #include <TOrdCollection.h>
60 #include <TParticle.h>
61 #include <TParticlePDG.h>
62 #include <TDatabasePDG.h>
63 #include <TH1.h>
64 #include <TH2.h>
65 #include <TH3.h>
66 #include <TProfile.h>
67 #include <TProfile2D.h>
68 #include <TVector2.h>
69
70 // Flow things
71 #include "AliFlowEvent.h"
72 #include "AliFlowTrack.h"
73 #include "AliFlowV0.h"
74 #include "AliFlowConstants.h"
75 #include "AliFlowSelection.h"
76 #include "AliFlowAnalyser.h"
77
78 // ANSI things
79 #include <math.h>
80 #include <stdio.h>
81 #include <stdlib.h>
82 #include <iostream>
83 using namespace std; //required for resolving the 'cout' symbol
84
85 ClassImp(AliFlowAnalyser) 
86 //-----------------------------------------------------------------------
87 AliFlowAnalyser::AliFlowAnalyser(const AliFlowSelection* flowSelect)
88 {
89  // default constructor (selection given or default selection)
90
91  if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect) ; }
92  else           { fFlowSelect = new AliFlowSelection() ; }
93  
94  // output file (histograms)
95  fHistFileName = "flowAnalysPlot.root" ;
96  fHistFile     = 0 ;
97
98  // for phi weights
99  fPhiWgtFile   = 0 ;
100  fPhiBins = AliFlowConstants::kPhiBins ;     
101  fPhiMin  = 0.;
102  fPhiMax  = 2*TMath::Pi() ; 
103
104  // flags 
105  fTrackLoop     = kTRUE ;  // main loop for charged tracks
106  fV0loop        = kTRUE ;  // correlation analysis is done also for neutral secundary vertex
107  fShuffle       = kFALSE ; // randomly reshuffles tracks
108  fV1Ep1Ep2      = kFALSE ; // disabled by default               
109  fEtaSub        = kFALSE ; // eta subevents
110  fReadPhiWgt    = kFALSE ; // if kTRUE & if flowPhiWgt file is there -> Phi Weights are used                                            
111  fBayWgt        = kFALSE ; // Bayesian Weights for P.Id. used   
112  fRePid         = kFALSE ; // recalculates the P.Id. (becomes kTRUE if new bayesian weights are plugged in)                     
113  fPtWgt         = kFALSE ; // pT as a weight
114  fEtaWgt        = kFALSE ; // eta as a weight
115  fOnePhiWgt     = kTRUE  ; // one/three phi-wgt histogram(s)
116 //  fMaxLabel   = 1000 ;     // for labels histogram
117
118  fPhiWgtHistList = 0 ;
119  fVnResHistList = 0 ;
120 }
121 //-----------------------------------------------------------------------
122 AliFlowAnalyser::~AliFlowAnalyser()
123 {
124  // default distructor (no actions) 
125 }
126 //-----------------------------------------------------------------------
127 Bool_t AliFlowAnalyser::Init() 
128 {
129  // sets some defaults for the analysis
130  
131  cout << "* FlowAnalysis *  -  Init()" << endl ; cout << endl ; 
132  
133  // Open output files (->plots)
134  fHistFile = new TFile(fHistFileName.Data(), "RECREATE");
135  fHistFile->cd() ;
136
137  // counters and pointers to 0
138  fEventNumber = 0 ;
139  fTrackNumber = 0 ;     
140  fV0Number    = 0 ;
141  fPidId       = -1 ;
142  //fNumberOfEvents = 0 ;        
143  fNumberOfTracks = 0 ;  
144  fNumberOfV0s    = 0 ;  
145  fFlowEvent  = 0 ;
146  fFlowTrack  = 0 ;
147  fFlowV0     = 0 ;
148  fFlowTracks = 0 ;
149  fFlowV0s    = 0 ;
150  fSelParts = 0 ;
151  fSelV0s   = 0 ;
152  for(Int_t ii=0;ii<3;ii++) { fVertex[ii] = 0 ; }
153
154  // Check for Reading weights: if weight file is not there, wgt arrays are not used (phiwgt & bayesian)
155  if(!fPhiWgtFile) 
156  { 
157   fReadPhiWgt = kFALSE ; 
158   fBayWgt     = kFALSE ;
159   cout << "  No wgt file. Phi & Bayesian weights will not be used ! " << endl ;
160  }
161  
162  // Histogram settings
163  fLabel = "Pseudorapidity" ;                    if(strlen(fFlowSelect->PidPart()) != 0)            { fLabel = "Rapidity"; }
164  fPtMaxPart = AliFlowConstants::fgPtMaxPart ;   if(fFlowSelect->PtMaxPart())    { fPtMaxPart = fFlowSelect->PtMaxPart() ; }
165  fPtBinsPart = AliFlowConstants::kPtBinsPart ;  if(fFlowSelect->PtBinsPart()) { fPtBinsPart = fFlowSelect->PtBinsPart() ; }
166  fPtMin = AliFlowConstants::fgPtMin ;       
167  fPtMax = AliFlowConstants::fgPtMax ;       
168  fEtaMin = AliFlowConstants::fgEtaMin ;
169  fEtaMax = AliFlowConstants::fgEtaMax ;
170  fPtBins  = AliFlowConstants::kPtBins ;      
171  fEtaBins = AliFlowConstants::kEtaBins ;      
172  // -
173  SetPtRangevEta(1.,0.);   // (AliFlowConstants::fgPtMin , AliFlowConstants::fgPtMax) ;
174  SetEtaRangevPt(1.,0.);   // (AliFlowConstants::fgEtaMin , AliFlowConstants::fgEtaMax) ;
175  // -
176  float triggerMin      =  -1.5;
177  float triggerMax      =  10.5;
178  float chargeMin       =  -2.5;
179  float chargeMax       =   2.5; 
180  float dcaMin          =   0.0;
181  float dcaMax          =  10.0; //  0.5; 
182  float glDcaMax        =  50.0; 
183  float chi2Min         =    0.;
184  float chi2Max         =  500.; 
185  float chi2normMax     =   50.; 
186  float chi2MaxC        =   30.; 
187  float chi2MaxI        =   60.; 
188  float fitPtsMin       =  -0.5;
189  float fitPtsMaxTpc    = 160.5;
190  float fitPtsMaxIts    =   6.5;
191  float fitPtsMaxTrd    = 130.5;
192  float fitPtsMaxTof    =   1.5;
193  float fitOverMaxMin   =    0.;
194  float fitOverMaxMax   =   1.1;
195  float multOverOrigMin =    0.;
196  float multOverOrigMax =    1.;
197  float vertexZMin      =  -10.;
198  float vertexZMax      =   10.; 
199  float vertexXYMin     =  -0.1;
200  float vertexXYMax     =   0.1; 
201  float etaSymMinPart   =   -1.;
202  float etaSymMaxPart   =    1.;
203  float etaSymMin       =   -2.;
204  float etaSymMax       =    2.;
205  float psiMin          =    0.;
206  float psiMax          = 2*TMath::Pi() ; 
207  float multMin         =    0. ;
208  float multMax         =25000. ;
209  float multV0          = 5000. ;
210  float qMin            =    0. ;
211  float qMax            =   20. ;
212  float pidMin          =    0. ;
213  float pidMax          =    1. ;
214  float centMin         =  -0.5 ;
215  float centMax         =   9.5 ;
216  float logpMin         =  -2.5 ;
217  float logpMax         =   2.5 ;
218  float pMin            =    0. ;
219  float pMax            =  100. ;
220  float dEdxMax         = 1000. ;
221  float dEdxMaxTPC      = 1500. ;
222  float tofmin          = 10000. ;
223  float tofmax          = 50000. ;
224  float trdmax          = 2000. ;
225  float lgMin           = -1000. ;
226  float lgMax           = 1000. ;
227  float lgMinV0         =    0. ;
228  float lgMaxV0         =   50. ;
229  float massMin         =  -0.1 ;
230  float massMax         =   2.1 ;
231  float zdcpartMin      =  -0.5 ;
232  float zdcpartMax      = 1000.5 ;
233  float zdceMin         =    0. ;
234  float zdceMax         = 10000. ;
235  // - 
236  enum { kTriggerBins      = 12,
237         kChargeBins       = 5,
238         kDcaBins          = 1000,
239         kChi2Bins         = 150,
240         kFitPtsBinsTpc    = 161,
241         kFitPtsBinsIts    = 7,
242         kFitPtsBinsTrd    = 131,
243         kFitPtsBinsTof    = 2,
244         kFitOverMaxBins   = 55,
245         kMultOverOrigBins =100,
246         kVertexZBins      = 200,
247         kVertexXYBins     = 1000,
248         kEtaSymBins       = 200,
249         kPhi3DBins        = 60,
250         kPsiBins          = 36,
251         kMultBins         = 250,
252         kPidBins          = 100,
253         kCentBins         = 10,
254         kDedxBins         = 1000,
255         kTofBins          = 1000,
256         kTrdBins          = 100,
257         kMomenBins        = 500,
258         kQbins            =  50,
259         kLgBins           = 500,
260         kMassBins         = 2200,
261         kZDCpartBins      = 1001,
262         kZDCeBins         = 200
263       } ;
264
265  // Histograms Booking ...
266  // Trigger
267  fHistTrigger = new TH1F("Flow_Trigger", "Flow_Trigger", kTriggerBins, triggerMin, triggerMax);
268  fHistTrigger->SetXTitle("Trigger: 1 minbias, 2 central, 3 laser, 10 other");
269  fHistTrigger->SetYTitle("Counts");
270
271  // Charge
272  fHistCharge = new TH1F("Flow_Charge", "Flow_Charge", kChargeBins, chargeMin, chargeMax);
273  fHistCharge->SetXTitle("Charge");
274  fHistCharge->SetYTitle("Counts");
275
276  // Reconstructed Momentum (constrained, if so)
277  fHistPtot = new TH1F("Flow_totalP","Flow_totalP", fPtBins, fPtMin, fPtMax);
278  fHistPtot->Sumw2();
279  fHistPtot->SetXTitle("P (GeV/c)");
280  fHistPtot->SetYTitle("Counts");
281
282  // Reconstructed pT (constrained, if so)
283  fHistPt = new TH1F("Flow_Pt","Flow_Pt", kMomenBins, pMin, pMax);
284  fHistPt->Sumw2();
285  fHistPt->SetXTitle("Pt (GeV/c)");
286  fHistPt->SetYTitle("Counts");
287    
288  // Reconstructed P.id vs pT 
289  fHistPidPt = new TH2F("Flow_PidPt","Flow_PidPt", 12, 0., 12., fPtBins, fPtMin, fPtMax);
290  fHistPidPt->Sumw2();
291  fHistPidPt->SetXTitle("e-, e+, mu-, mu+, pi-, pi+, K-, K+, pr-, pr+, d-, d+");
292  fHistPidPt->SetYTitle("Pt (GeV/c)");
293    
294  // Distance of closest approach - Transverse, Unsigned
295  fHistDca = new TH1F("Flow_TransDCA", "Flow_TransverseDCA", kDcaBins, dcaMin, dcaMax);
296  fHistDca->Sumw2();
297  fHistDca->SetXTitle("|Track's Signed dca (cm)|");
298  fHistDca->SetYTitle("Counts");
299
300  // Distance of closest approach - Transverse
301  fHistTransDca = new TH1F("Flow_TransDCA_Signed", "Flow_TransverseDCA_Signed", kDcaBins, -dcaMax, dcaMax);
302  fHistTransDca->Sumw2();
303  fHistTransDca->SetXTitle("Track's Transverse dca (cm)");
304  fHistTransDca->SetYTitle("Counts");
305
306  // Distance of closest approach - 3d
307  fHistDcaGlobal = new TH1F("Flow_3dDcaGlobal", "Flow_3dDcaGlobal", kDcaBins, dcaMin, glDcaMax);
308  fHistDcaGlobal->Sumw2();
309  fHistDcaGlobal->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
310  fHistDcaGlobal->SetYTitle("Counts");
311
312  // Distance of closest approach for particles correlated with the event plane - 3d
313  fHistDcaGlobalPart = new TH1F("Flow_3dDcaGlobalPart", "Flow_3dDcaGlobalPart", kDcaBins, dcaMin, glDcaMax);
314  fHistDcaGlobalPart->Sumw2();
315  fHistDcaGlobalPart->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
316  fHistDcaGlobalPart->SetYTitle("Counts");
317
318  // Distance of closest approach for particles NOT SELECTED for correlation with the event plane - 3d
319  fHistDcaGlobalOut = new TH1F("Flow_3dDcaGlobalOut", "Flow_3dDcaGlobalOut", kDcaBins, dcaMin, glDcaMax);
320  fHistDcaGlobalOut->Sumw2();
321  fHistDcaGlobalOut->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
322  fHistDcaGlobalOut->SetYTitle("Counts");
323
324  // Yield Pt-Phi for constrained tracks
325  fHistPhiPtCon = new TH2D("Flow_PtPhi_Con", "Flow_PtPhi_Con", fPhiBins, fPhiMin, fPhiMax, fPtBins, fPtMin, fPtMax);
326  fHistPhiPtCon->Sumw2();
327  fHistPhiPtCon->SetXTitle("Phi");
328  fHistPhiPtCon->SetYTitle("Pt (GeV/c)");
329  
330  // Yield Pt-Phi for UNconstrained tracks
331  fHistPhiPtUnc = new TH2D("Flow_PtPhi_Unc", "Flow_PtPhi_Unc", fPhiBins, fPhiMin, fPhiMax, fPtBins, fPtMin, fPtMax);
332  fHistPhiPtUnc->Sumw2();
333  fHistPhiPtUnc->SetXTitle("Phi");
334  fHistPhiPtUnc->SetYTitle("Pt (GeV/c)");
335  
336  // Chi2 
337  // at the main vertex
338  fHistChi2 = new TH1F("Flow_Chi2", "Flow_Chi2", kChi2Bins, chi2Min, chi2MaxC);
339  fHistChi2->SetXTitle("Chi square at Main Vertex");
340  fHistChi2->SetYTitle("Counts");
341  // TPC
342  fHistChi2TPC = new TH1F("Flow_Chi2_TPC", "Flow_Chi2_TPC", kChi2Bins, chi2Min, chi2Max);
343  fHistChi2TPC->SetXTitle("Chi square for TPC");
344  fHistChi2TPC->SetYTitle("Counts");
345  // ITS
346  fHistChi2ITS = new TH1F("Flow_Chi2_ITS", "Flow_Chi2_ITS", kChi2Bins, chi2Min, chi2MaxI);
347  fHistChi2ITS->SetXTitle("Chi square for ITS");
348  fHistChi2ITS->SetYTitle("Counts");
349  // TRD
350  fHistChi2TRD = new TH1F("Flow_Chi2_TRD", "Flow_Chi2_TRD", kChi2Bins, chi2Min, chi2Max);
351  fHistChi2TRD->SetXTitle("Chi square for TRD");
352  fHistChi2TRD->SetYTitle("Counts");
353  // TOF
354  fHistChi2TOF = new TH1F("Flow_Chi2_TOF", "Flow_Chi2_TOF", kChi2Bins, chi2Min, chi2Max);
355  fHistChi2TOF->SetXTitle("Chi square for TOF");
356  fHistChi2TOF->SetYTitle("Counts");
357
358  // Chi2 normalized
359  // TPC
360  fHistChi2normTPC = new TH1F("Flow_Chi2norm_TPC", "Flow_Chi2norm_TPC", kChi2Bins, chi2Min, chi2normMax);
361  fHistChi2normTPC->SetXTitle("Normalized #Chi^{2} for TPC");
362  fHistChi2normTPC->SetYTitle("Counts");
363  // ITS
364  fHistChi2normITS = new TH1F("Flow_Chi2norm_ITS", "Flow_Chi2norm_ITS", kChi2Bins, chi2Min, chi2normMax);
365  fHistChi2normITS->SetXTitle("Normalized #Chi^{2} for ITS");
366  fHistChi2normITS->SetYTitle("Counts");
367  // TRD
368  fHistChi2normTRD = new TH1F("Flow_Chi2norm_TRD", "Flow_Chi2norm_TRD", kChi2Bins, chi2Min, chi2normMax);
369  fHistChi2normTRD->SetXTitle("Normalized #Chi^{2} for TRD");
370  fHistChi2normTRD->SetYTitle("Counts");
371  // TOF
372  fHistChi2normTOF = new TH1F("Flow_Chi2norm_TOF", "Flow_Chi2norm_TOF", kChi2Bins, chi2Min, chi2normMax);
373  fHistChi2normTOF->SetXTitle("Normalized #Chi^{2} for TOF");
374  fHistChi2normTOF->SetYTitle("Counts");
375  
376  // FitPts
377  // TPC
378  fHistFitPtsTPC = new TH1F("Flow_FitPts_TPC", "Flow_FitPts_TPC", kFitPtsBinsTpc, fitPtsMin, fitPtsMaxTpc);
379  fHistFitPtsTPC->SetXTitle("Fit Points");
380  fHistFitPtsTPC->SetYTitle("Counts");
381  // ITS
382  fHistFitPtsITS = new TH1F("Flow_HitPts_ITS", "Flow_HitPts_ITS", kFitPtsBinsIts, fitPtsMin, fitPtsMaxIts);
383  fHistFitPtsITS->SetXTitle("Fit Points");
384  fHistFitPtsITS->SetYTitle("Counts");
385  // TRD
386  fHistFitPtsTRD = new TH1F("Flow_FitPts_TRD", "Flow_FitPts_TRD", kFitPtsBinsTrd, fitPtsMin, fitPtsMaxTrd);
387  fHistFitPtsTRD->SetXTitle("Fit Points");
388  fHistFitPtsTRD->SetYTitle("Counts");
389  // TOF
390  fHistFitPtsTOF = new TH1F("Flow_HitPts_TOF", "Flow_HitPts_TOF", kFitPtsBinsTof, fitPtsMin, fitPtsMaxTof);
391  fHistFitPtsTOF->SetXTitle("Fit Points");
392  fHistFitPtsTOF->SetYTitle("Counts");
393
394  // MaxPts 
395  // TPC
396  fHistMaxPtsTPC = new TH1F("Flow_MaxPts_TPC", "Flow_MaxPts_TPC", kFitPtsBinsTpc, fitPtsMin, fitPtsMaxTpc);
397  fHistMaxPtsTPC->SetXTitle("Max Points");
398  fHistMaxPtsTPC->SetYTitle("Counts");
399  // ITS
400  fHistMaxPtsITS = new TH1F("Flow_MaxPts_ITS", "Flow_MaxPts_ITS", kFitPtsBinsIts, fitPtsMin, fitPtsMaxIts);
401  fHistMaxPtsITS->SetXTitle("Max Points");
402  fHistMaxPtsITS->SetYTitle("Counts");
403  // TRD
404  fHistMaxPtsTRD = new TH1F("Flow_MaxPts_TRD", "Flow_MaxPts_TRD", kFitPtsBinsTrd, fitPtsMin, fitPtsMaxTrd);
405  fHistMaxPtsTRD->SetXTitle("Max Points");
406  fHistMaxPtsTRD->SetYTitle("Counts");
407  // TOF
408  fHistMaxPtsTOF = new TH1F("Flow_MaxPts_TOF", "Flow_MaxPts_TOF", kFitPtsBinsTof, fitPtsMin, fitPtsMaxTof);
409  fHistMaxPtsTOF->SetXTitle("Max Points");
410  fHistMaxPtsTOF->SetYTitle("Counts");
411
412  // FitOverMax 
413  // Tpc
414  fHistFitOverMaxTPC = new TH1F("Flow_FitOverMax_TPC", "Flow_FitOverMax_TPC", kFitOverMaxBins, fitOverMaxMin, fitOverMaxMax);
415  fHistFitOverMaxTPC->SetXTitle("(Fit Points - 1) / Max Points");
416  fHistFitOverMaxTPC->SetYTitle("Counts");  
417  // All
418  fHistFitOverMax = new TH1F("Flow_FitOverMax", "Flow_FitOverMax", kFitOverMaxBins, fitOverMaxMin, fitOverMaxMax);
419  fHistFitOverMax->SetXTitle("(Fit Points - 1) / Max Points");
420  fHistFitOverMax->SetYTitle("Counts");  
421
422  // lenght
423  fHistLenght = new TH1F("Flow_TrackLenght", "Flow_TrackLenght", kLgBins, lgMin, lgMax);
424  fHistLenght->SetXTitle("Lenght of the Track (cm)");
425  fHistLenght->SetYTitle("Counts");
426
427   // OrigMult
428  fHistOrigMult = new TH1F("Flow_OrigMult", "Flow_OrigMult", kMultBins, multMin, multMax);
429  fHistOrigMult->SetXTitle("Original Mult");
430  fHistOrigMult->SetYTitle("Counts");
431
432  // MultEta
433  fHistMultEta = new TH1F("Flow_MultEta", "Flow_MultEta", kMultBins, multMin, multMax);
434  fHistMultEta->SetXTitle("Mult for Centrality");
435  fHistMultEta->SetYTitle("Counts");
436    
437  // Mult
438  fHistMult = new TH1F("Flow_Mult", "Flow_Mult", kMultBins, multMin, multMax);
439  fHistMult->SetXTitle("Mult");
440  fHistMult->SetYTitle("Counts");
441
442  // V0s multiplicity
443  fHistV0Mult = new TH1F("FlowV0_Mult","FlowV0_Mult", kMultBins, multMin, multV0);
444  fHistV0Mult->SetXTitle("V0s Multiplicity");
445  fHistV0Mult->SetYTitle("Counts");
446    
447  // MultOverOrig
448  fHistMultOverOrig = new TH1F("Flow_MultOverOrig", "Flow_MultOverOrig", kMultOverOrigBins, multOverOrigMin, multOverOrigMax);
449  fHistMultOverOrig->SetXTitle("Mult / Orig. Mult");
450  fHistMultOverOrig->SetYTitle("Counts");
451    
452  // Mult correlated with the event planes
453  fHistMultPart = new TH1F("Flow_MultPart", "Flow_MultPart", 2*kMultBins, multMin, multMax);
454  fHistMultPart->SetXTitle("Mult of Correlated Particles");
455  fHistMultPart->SetYTitle("Counts");
456    
457  // Mult correlated with the event planes in 1 unit rapidity
458  fHistMultPartUnit = new TH1F("Flow_MultPartUnit", "Flow_MultPartUnit", 2*kMultBins, multMin, multMax/2);
459  fHistMultPartUnit->SetXTitle("Mult of Correlated Particles (-0.5 < eta < 0.5)");
460  fHistMultPartUnit->SetYTitle("Counts");
461    
462  // Mult of V0s correlated with the event planes
463  fHistV0MultPart = new TH1F("FlowV0_MultPart", "FlowV0_MultPart", kMultBins, multMin, multV0);
464  fHistV0MultPart->SetXTitle("Mult of Correlated V0s");
465  fHistV0MultPart->SetYTitle("Counts");
466    
467  // VertexZ
468  fHistVertexZ = new TH1F("Flow_VertexZ", "Flow_VertexZ", kVertexZBins, vertexZMin, vertexZMax);
469  fHistVertexZ->SetXTitle("Vertex Z (cm)");
470  fHistVertexZ->SetYTitle("Counts");
471    
472  // VertexXY
473  fHistVertexXY2D = new TH2F("Flow_VertexXY2D", "Flow_VertexXY2D", kVertexXYBins, vertexXYMin, vertexXYMax, kVertexXYBins, vertexXYMin, vertexXYMax);
474  fHistVertexXY2D->SetXTitle("Vertex X (cm)");
475  fHistVertexXY2D->SetYTitle("Vertex Y (cm)");
476    
477  // EtaSym vs. Vertex Z Tpc
478  fHistEtaSymVerZ2D = new TH2F("Flow_EtaSymVerZ2D", "Flow_EtaSymVerZ2D", kVertexZBins, vertexZMin, vertexZMax, kEtaSymBins, etaSymMin, etaSymMax);
479  fHistEtaSymVerZ2D->SetXTitle("Vertex Z (cm)");
480  fHistEtaSymVerZ2D->SetYTitle("Eta Symmetry TPC");
481
482  // EtaSym Tpc
483  fHistEtaSym = new TH1F("Flow_EtaSym_TPC", "Flow_EtaSym_TPC", kEtaSymBins, etaSymMin, etaSymMax);
484  fHistEtaSym->SetXTitle("Eta Symmetry Ratio TPC");
485  fHistEtaSym->SetYTitle("Counts");
486    
487  // EtaSym vs. Vertex Z Tpc - correlation analysis
488  fHistEtaSymVerZ2DPart = new TH2F("Flow_EtaSymVerZ2D_part", "Flow_EtaSymVerZ2D_part", kVertexZBins, vertexZMin, vertexZMax, kEtaSymBins, etaSymMinPart, etaSymMaxPart);
489  fHistEtaSymVerZ2DPart->SetXTitle("Vertex Z (cm)");
490  fHistEtaSymVerZ2DPart->SetYTitle("Eta Symmetry TPC");
491
492  // EtaSym Tpc - correlation analysis
493  fHistEtaSymPart = new TH1F("Flow_EtaSym_TPC_part", "Flow_EtaSym_TPC_part", kEtaSymBins, etaSymMinPart, etaSymMaxPart);
494  fHistEtaSymPart->SetXTitle("Eta Symmetry Ratio TPC");
495  fHistEtaSymPart->SetYTitle("Counts");
496    
497  // phi , whatever
498  fHistPhi = new TH1F("Flow_xPhi", "Flow_xPhi (all)", fPhiBins, fPhiMin, fPhiMax);
499  fHistPhi->SetXTitle("Phi (rad)");
500  fHistPhi->SetYTitle("Counts");
501
502  // phi constrained
503  fHistPhiCons = new TH1F("Flow_cPhi", "Flow_cPhi", fPhiBins, fPhiMin, fPhiMax);
504  fHistPhiCons->SetXTitle("cPhi (rad)");
505  fHistPhiCons->SetYTitle("Counts");
506    
507  // EtaPtPhi , whatever
508  fHistAllEtaPtPhi3D = new TH3F("Flow_EtaPtPhi3Dall", "Flow_EtaPtPhi3Dall (whatever)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
509  fHistAllEtaPtPhi3D->SetXTitle("Eta");
510  fHistAllEtaPtPhi3D->SetYTitle("Pt (GeV/c)");
511  fHistAllEtaPtPhi3D->SetZTitle("Phi (rad)");
512    
513  // Constrained EtaPtPhi
514  fHistConsEtaPtPhi3D = new TH3F("Flow_consEtaPtPhi3D", "Flow_consEtaPtPhi3D (constrainable)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
515  fHistConsEtaPtPhi3D->SetXTitle("cEta");
516  fHistConsEtaPtPhi3D->SetYTitle("cPt (GeV/c)");
517  fHistConsEtaPtPhi3D->SetZTitle("cPhi (rad)");
518    
519  // Global EtaPtPhi
520  fHistGlobEtaPtPhi3D = new TH3F("Flow_globEtaPtPhi3D", "Flow_globEtaPtPhi3D (constrainable)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
521  fHistGlobEtaPtPhi3D->SetXTitle("gEta");
522  fHistGlobEtaPtPhi3D->SetYTitle("gPt (GeV/c)");
523  fHistGlobEtaPtPhi3D->SetZTitle("gPhi (rad)");
524
525  // UnConstrained EtaPtPhi
526  fHistUncEtaPtPhi3D = new TH3F("Flow_uncEtaPtPhi3D", "Flow_uncEtaPtPhi3D (un-constrainable)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
527  fHistUncEtaPtPhi3D->SetXTitle("gEta");
528  fHistUncEtaPtPhi3D->SetYTitle("gPt (GeV/c)");
529  fHistUncEtaPtPhi3D->SetZTitle("gPhi (rad)");
530    
531  // EtaPtPhi for particles correlated with the event plane
532  fHistEtaPtPhi3DPart = new TH3F("Flow_EtaPtPhi3Dpart", "Flow_EtaPtPhi3Dpart (selected part)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
533  fHistEtaPtPhi3DPart->SetXTitle("Eta");
534  fHistEtaPtPhi3DPart->SetYTitle("Pt (GeV/c)");
535  fHistEtaPtPhi3DPart->SetZTitle("Phi (rad)");
536    
537  // EtaPtPhi for particles NOT SELECTED for correlation with the event plane
538  fHistEtaPtPhi3DOut = new TH3F("Flow_EtaPtPhi3Dout", "Flow_EtaPtPhi3Dout (NOT selected part)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
539  fHistEtaPtPhi3DOut->SetXTitle("Eta");
540  fHistEtaPtPhi3DOut->SetYTitle("Pt (GeV/c)");
541  fHistEtaPtPhi3DOut->SetZTitle("Phi (rad)");
542    
543  // Yield Pt-Phi for all positive
544  fHistPtPhiPos = new TH2D("Flow_PtPhi_Plus", "Flow_PtPhi_Plus", fPhiBins, fPhiMin, fPhiMax, 50, fPtMin, fPtMax);
545  fHistPtPhiPos->Sumw2();
546  fHistPtPhiPos->SetXTitle("Phi");
547  fHistPtPhiPos->SetYTitle("Pt (GeV/c)");
548
549  // Yield Pt-Phi for all negative
550  fHistPtPhiNeg = new TH2D("Flow_PtPhi_Minus", "Flow_PtPhi_Minus", fPhiBins, fPhiMin, fPhiMax, 50, fPtMin, fPtMax);
551  fHistPtPhiNeg->Sumw2();
552  fHistPtPhiNeg->SetXTitle("Phi");
553  fHistPtPhiNeg->SetYTitle("Pt (GeV/c)");
554    
555  // Yield for all particles
556  fHistYieldAll2D = new TH2D("Flow_YieldAll2D", "Flow_YieldAll2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
557  fHistYieldAll2D->Sumw2();
558  fHistYieldAll2D->SetXTitle("Pseudorapidty");
559  fHistYieldAll2D->SetYTitle("Pt (GeV/c)");
560
561  // Yield for constrainable tracks
562  fHistYieldCon2D = new TH2D("Flow_YieldCons2D", "Flow_YieldCons2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
563  fHistYieldCon2D->Sumw2();
564  fHistYieldCon2D->SetXTitle("Pseudorapidty");
565  fHistYieldCon2D->SetYTitle("Pt (GeV/c)");
566
567  // Yield for un-constrainable tracks
568  fHistYieldUnc2D = new TH2D("Flow_YieldUnc2D", "Flow_YieldUnc2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
569  fHistYieldUnc2D->Sumw2();
570  fHistYieldUnc2D->SetXTitle("Pseudorapidty");
571  fHistYieldUnc2D->SetYTitle("Pt (GeV/c)");
572
573  // Yield for particles correlated with the event plane
574  fHistYieldPart2D = new TH2D("Flow_YieldPart2D", "Flow_YieldPart2D (selected part)", fEtaBins, fEtaMin, fEtaMax, fPtBinsPart, fPtMin, fPtMaxPart);
575  fHistYieldPart2D->Sumw2();
576  fHistYieldPart2D->SetXTitle((char*)fLabel.Data());
577  fHistYieldPart2D->SetYTitle("Pt (GeV/c)");
578
579  // Yield for particles NOT SELECTED for correlation with the event plane
580  fHistYieldOut2D = new TH2D("Flow_YieldOut2D", "Flow_YieldOut2D (NOT selected part)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
581  fHistYieldOut2D->Sumw2();
582  fHistYieldOut2D->SetXTitle("Pseudorapidty");
583  fHistYieldOut2D->SetYTitle("Pt (GeV/c)");
584  
585  // invariant Mass for all particles (from TOF)
586  fHistInvMass = new TH1F("Flow_InvMass", "Flow_InvMass (tof)", kMassBins, massMin, massMax);
587  fHistInvMass->SetXTitle("Invariant Mass (GeV)");
588  fHistInvMass->SetYTitle("Counts");
589  
590  // invariant Mass for particles correlated with the event plane (from TOF)
591  fHistInvMassPart = new TH1F("Flow_InvMassPart", "Flow_InvMassPart (tof)", kMassBins, massMin, massMax);
592  fHistInvMassPart->SetXTitle("Invariant Mass (GeV)");
593  fHistInvMassPart->SetYTitle("Counts");
594  
595  // invariant Mass for particles NOT SELECTED for correlation with the event plane (from TOF)
596  fHistInvMassOut = new TH1F("Flow_InvMassOut", "Flow_InvMassOut (tof)", kMassBins, massMin, massMax);
597  fHistInvMassOut->SetXTitle("Invariant Mass (GeV)");
598  fHistInvMassOut->SetYTitle("Counts");
599
600  // Mean Eta in each bin for particles correlated with the event plane
601  fHistBinEta = new TProfile("Flow_Bin_Eta", "Flow_Bin_Eta_part (selected part)", fEtaBins, fEtaMin, fEtaMax, fEtaMin, fEtaMax, "");
602  fHistBinEta->SetXTitle((char*)fLabel.Data());
603  fHistBinEta->SetYTitle((char*)fLabel.Data());
604  
605  // Mean Pt in each bin for particles correlated with the event plane
606  fHistBinPt = new TProfile("Flow_Bin_Pt", "Flow_Bin_Pt_part (selected part)", fPtBinsPart, fPtMin, fPtMaxPart, fPtMin, fPtMaxPart, "");
607  fHistBinPt->SetXTitle("Pt (GeV/c)");
608  fHistBinPt->SetYTitle("<Pt> (GeV/c)");
609  
610  // cos(n*phiLab)
611  fHistCosPhi = new TProfile("Flow_CosPhiLab", "Flow_CosPhiLab", AliFlowConstants::kHars, 0.5, (float)(AliFlowConstants::kHars) + 0.5, -100., 100., "");
612  fHistCosPhi->SetXTitle("Harmonic");
613  fHistCosPhi->SetYTitle("<cos(n*PhiLab)> (%)");
614    
615  // PID pi+
616  fHistPidPiPlus = new TH1F("Flow_PidPiPlus", "Flow_PidPiPlus", kPidBins, pidMin, pidMax);
617  fHistPidPiPlus->SetXTitle("ALICE P.Id.");
618  fHistPidPiPlus->SetYTitle("Counts");
619    
620  // PID pi-
621  fHistPidPiMinus = new TH1F("Flow_PidPiMinus", "Flow_PidPiMinus", kPidBins, pidMin, pidMax);
622  fHistPidPiMinus->SetXTitle("ALICE P.Id.");
623  fHistPidPiMinus->SetYTitle("Counts");
624    
625  // PID proton
626  fHistPidProton = new TH1F("Flow_PidProton", "Flow_PidProton", kPidBins, pidMin, pidMax);
627  fHistPidProton->SetXTitle("ALICE P.Id.");
628  fHistPidProton->SetYTitle("Counts");
629
630  // PID anti proton
631  fHistPidAntiProton = new TH1F("Flow_PidAntiProton", "Flow_PidAntiProton", kPidBins, pidMin, pidMax);
632  fHistPidAntiProton->SetXTitle("ALICE P.Id.");
633  fHistPidAntiProton->SetYTitle("Counts");
634
635  // PID Kplus
636  fHistPidKplus = new TH1F("Flow_PidKplus", "Flow_PidKplus", kPidBins, pidMin, pidMax);
637  fHistPidKplus->SetXTitle("ALICE P.Id.");
638  fHistPidKplus->SetYTitle("Counts");
639
640  // PID Kminus
641  fHistPidKminus = new TH1F("Flow_PidKminus", "Flow_PidKminus", kPidBins, pidMin, pidMax);
642  fHistPidKminus->SetXTitle("ALICE P.Id.");
643  fHistPidKminus->SetYTitle("Counts");
644
645  // PID deuteron
646  fHistPidDeuteron = new TH1F("Flow_PidDeuteron", "Flow_PidDeuteron", kPidBins, pidMin, pidMax);
647  fHistPidDeuteron->SetXTitle("ALICE P.Id.");
648  fHistPidDeuteron->SetYTitle("Counts");
649
650  // PID anti deuteron
651  fHistPidAntiDeuteron = new TH1F("Flow_PidAntiDeuteron", "Flow_PidAntiDeuteron", kPidBins, pidMin, pidMax);
652  fHistPidAntiDeuteron->SetXTitle("ALICE P.Id.");
653  fHistPidAntiDeuteron->SetYTitle("Counts");
654
655  // PID electron
656  fHistPidElectron = new TH1F("Flow_PidElectron", "Flow_PidElectron", kPidBins, pidMin, pidMax);
657  fHistPidElectron->SetXTitle("ALICE P.Id.");
658  fHistPidElectron->SetYTitle("Counts");
659
660  // PID positron
661  fHistPidPositron = new TH1F("Flow_PidPositron", "Flow_PidPositron", kPidBins, pidMin, pidMax);
662  fHistPidPositron->SetXTitle("ALICE P.Id.");
663  fHistPidPositron->SetYTitle("Counts");
664
665  // PID Muon+
666  fHistPidMuonPlus = new TH1F("Flow_PidMuonPlus", "Flow_PidMuonPlus", kPidBins, pidMin, pidMax);
667  fHistPidMuonPlus->SetXTitle("ALICE P.Id.");
668  fHistPidMuonPlus->SetYTitle("Counts");
669
670  // PID Muon-
671  fHistPidMuonMinus = new TH1F("Flow_PidMuonMinus", "Flow_PidMuonMinus", kPidBins, pidMin, pidMax);
672  fHistPidMuonMinus->SetXTitle("ALICE P.Id.");
673  fHistPidMuonMinus->SetYTitle("Counts");
674
675  // PID pi+ selected
676  fHistPidPiPlusPart = new TH1F("Flow_PidPiPlusPart", "Flow_PidPiPlusPart", kPidBins, pidMin, pidMax);
677  fHistPidPiPlusPart->SetXTitle("ALICE P.Id. (pid = pi+)");
678  fHistPidPiPlusPart->SetYTitle("Counts");
679    
680  // PID pi- selected
681  fHistPidPiMinusPart = new TH1F("Flow_PidPiMinusPart", "Flow_PidPiMinusPart", kPidBins, pidMin, pidMax);
682  fHistPidPiMinusPart->SetXTitle("ALICE P.Id. (pid = pi-)");
683  fHistPidPiMinusPart->SetYTitle("Counts");
684    
685  // PID proton selected
686  fHistPidProtonPart = new TH1F("Flow_PidProtonPart", "Flow_PidProtonPart", kPidBins, pidMin, pidMax);
687  fHistPidProtonPart->SetXTitle("ALICE P.Id. (pid = p)");
688  fHistPidProtonPart->SetYTitle("Counts");
689
690  // PID anti proton selected
691  fHistPidAntiProtonPart = new TH1F("Flow_PidAntiProtonPart", "Flow_PidAntiProtonPart", kPidBins, pidMin, pidMax);
692  fHistPidAntiProtonPart->SetXTitle("ALICE P.Id. (pid = p-)");
693  fHistPidAntiProtonPart->SetYTitle("Counts");
694
695  // PID Kplus selected
696  fHistPidKplusPart = new TH1F("Flow_PidKplusPart", "Flow_PidKplusPart", kPidBins, pidMin, pidMax);
697  fHistPidKplusPart->SetXTitle("ALICE P.Id. (pid = K+)");
698  fHistPidKplusPart->SetYTitle("Counts");
699
700  // PID Kminus selected
701  fHistPidKminusPart = new TH1F("Flow_PidKminusPart", "Flow_PidKminusPart", kPidBins, pidMin, pidMax);
702  fHistPidKminusPart->SetXTitle("ALICE P.Id. (pid = K-)");
703  fHistPidKminusPart->SetYTitle("Counts");
704
705  // PID deuteron selected
706  fHistPidDeuteronPart = new TH1F("Flow_PidDeuteronPart", "Flow_PidDeuteronPart", kPidBins, pidMin, pidMax);
707  fHistPidDeuteronPart->SetXTitle("ALICE P.Id. (pid = d)");
708  fHistPidDeuteronPart->SetYTitle("Counts");
709
710  // PID anti deuteron selected
711  fHistPidAntiDeuteronPart = new TH1F("Flow_PidAntiDeuteronPart", "Flow_PidAntiDeuteronPart", kPidBins, pidMin, pidMax);
712  fHistPidAntiDeuteronPart->SetXTitle("ALICE P.Id. (pid = d--)");
713  fHistPidAntiDeuteronPart->SetYTitle("Counts");
714
715  // PID electron selected
716  fHistPidElectronPart = new TH1F("Flow_PidElectronPart", "Flow_PidElectronPart", kPidBins, pidMin, pidMax);
717  fHistPidElectronPart->SetXTitle("ALICE P.Id. (pid = e-)");
718  fHistPidElectronPart->SetYTitle("Counts");
719
720  // PID positron selected
721  fHistPidPositronPart = new TH1F("Flow_PidPositronPart", "Flow_PidPositronPart", kPidBins, pidMin, pidMax);
722  fHistPidPositronPart->SetXTitle("ALICE P.Id. (pid = e+)");
723  fHistPidPositronPart->SetYTitle("Counts");
724
725  // PID Muon+ selected
726  fHistPidMuonPlusPart = new TH1F("Flow_PidMuonPlusPart", "Flow_PidMuonPlusPart", kPidBins, pidMin, pidMax);
727  fHistPidMuonPlusPart->SetXTitle("ALICE P.Id. (pid = mu+)");
728  fHistPidMuonPlusPart->SetYTitle("Counts");
729
730  // PID Muon- selected
731  fHistPidMuonMinusPart = new TH1F("Flow_PidMuonMinusPart", "Flow_PidMuonMinusPart", kPidBins, pidMin, pidMax);
732  fHistPidMuonMinusPart->SetXTitle("ALICE P.Id. (pid = mu-)");
733  fHistPidMuonMinusPart->SetYTitle("Counts");
734
735  // PID multiplicities (all)
736  fHistPidMult = new TProfile("Flow_PidMult", "Flow_PidMult", 15, 0.5, 15.5, "");
737  fHistPidMult->SetXTitle("all, h+, h-, pi+, pi-, pr+, pr-, K+, K-, d+, d-, e-, e+, mu-, mu+");
738  fHistPidMult->SetYTitle("Multiplicity");
739
740  // PID for all tracks
741  fHistBayPidMult = new TH1F("Flow_BayPidMult","Flow_BayPidMult",AliFlowConstants::kPid,-0.5,((float)AliFlowConstants::kPid-0.5));
742  fHistBayPidMult->Sumw2() ;
743  fHistBayPidMult->SetXTitle("e+/- , mu+/- , pi+/- , K+/- , p+/- , d+/- ");
744  fHistBayPidMult->SetYTitle("Counts");
745    
746  // PID for particles correlated with the event plane
747  fHistBayPidMultPart = new TH1F("Flow_BayPidMult_Part","Flow_BayPidMult_Part",AliFlowConstants::kPid,-0.5,((float)AliFlowConstants::kPid-0.5));
748  fHistBayPidMultPart->Sumw2() ;
749  fHistBayPidMultPart->SetXTitle("e+/- , mu+/- , pi+/- , K+/- , p+/- , d+/- ");
750  fHistBayPidMultPart->SetYTitle("Counts");
751
752  // Centrality
753  fHistCent = new TH1F("Flow_Cent", "Flow_Cent", kCentBins, centMin, centMax);
754  fHistCent->SetXTitle("Centrality Bin");
755  fHistCent->SetYTitle("Counts");
756    
757  // Deposited Energy in ZDC
758  fHistEnergyZDC = new TH2F("Flow_ZDC_E", "Flow_ZDC_E", 3, 0., 3., kZDCeBins, zdceMin, zdceMax);
759  fHistEnergyZDC->SetXTitle("neutron , proton , e.m.");
760  fHistEnergyZDC->SetYTitle("ZDC energy");
761
762  // Part. versus Energy in ZDC
763  fHistPartZDC = new TH1F("Flow_ZDC_Participants", "Flow_ZDC_Participants", kZDCpartBins, zdcpartMin, zdcpartMax);
764  fHistPartZDC->SetXTitle("ZDC part");
765  fHistPartZDC->SetYTitle("Counts");
766
767  // MeanDedxPos TPC
768  fHistMeanDedxPos2D = new TH2F("Flow_MeanDedxPos2D_TPC","Flow_MeanDedxPos2D_TPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
769  fHistMeanDedxPos2D->SetXTitle("log(momentum) (GeV/c)");
770  fHistMeanDedxPos2D->SetYTitle("mean dEdx (+)");
771  
772  // MeanDedxNeg TPC
773  fHistMeanDedxNeg2D = new TH2F("Flow_MeanDedxNeg2D_TPC","Flow_MeanDedxNeg2D_TPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
774  fHistMeanDedxNeg2D->SetXTitle("log(momentum) (GeV/c)");
775  fHistMeanDedxNeg2D->SetYTitle("mean dEdx (-)");
776
777  // MeanDedxPos ITS
778  fHistMeanDedxPos2DITS = new TH2F("Flow_MeanDedxPos2D_ITS","Flow_MeanDedxPos2D_ITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
779  fHistMeanDedxPos2DITS->SetXTitle("log(momentum) (GeV/c)");
780  fHistMeanDedxPos2DITS->SetYTitle("mean dEdx (+)");
781  
782  // MeanDedxNeg ITS
783  fHistMeanDedxNeg2DITS = new TH2F("Flow_MeanDedxNeg2D_ITS","Flow_MeanDedxNeg2D_ITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
784  fHistMeanDedxNeg2DITS->SetXTitle("log(momentum) (GeV/c)");
785  fHistMeanDedxNeg2DITS->SetYTitle("mean dEdx (-)");
786
787  // MeanDedxPos TPC 3D selected Part
788  fHistMeanDedxPos3DPart = new TH3F("Flow_MeanDedxPos3DPart_TPC","Flow_MeanDedxPos3DPart_TPC", kMomenBins, logpMin, logpMax, kDedxBins, 0., dEdxMaxTPC, AliFlowConstants::kPid, 0., AliFlowConstants::kPid);
789  fHistMeanDedxPos3DPart->SetXTitle("log(momentum) (GeV/c)");
790  fHistMeanDedxPos3DPart->SetYTitle("mean dEdx (+)");
791  fHistMeanDedxPos3DPart->SetZTitle("e , mu , pi , k , p , d");
792  
793  // MeanDedxNeg TPC 3D selected Part
794  fHistMeanDedxNeg3DPart = new TH3F("Flow_MeanDedxNeg3DPart_TPC","Flow_MeanDedxNeg3DPart_TPC", kMomenBins, logpMin, logpMax, kDedxBins, 0., dEdxMaxTPC, AliFlowConstants::kPid, 0., AliFlowConstants::kPid);
795  fHistMeanDedxNeg3DPart->SetXTitle("log(momentum) (GeV/c)");
796  fHistMeanDedxNeg3DPart->SetYTitle("mean dEdx (-)");
797  fHistMeanDedxNeg3DPart->SetZTitle("e , mu , pi , k , p , d");
798
799  // MeanDedxPos ITS 3D selected Part
800  fHistMeanDedxPos3DPartITS = new TH3F("Flow_MeanDedxPos3DPart_ITS","Flow_MeanDedxPos3DPart_ITS", kMomenBins, logpMin, logpMax, kDedxBins, 0., dEdxMax, AliFlowConstants::kPid, 0., AliFlowConstants::kPid);
801  fHistMeanDedxPos3DPartITS->SetXTitle("log(momentum) (GeV/c)");
802  fHistMeanDedxPos3DPartITS->SetYTitle("mean dEdx (+)");
803  fHistMeanDedxPos3DPartITS->SetZTitle("e , mu , pi , k , p , d");
804  
805  // MeanDedxNeg ITS 3D selected Part
806  fHistMeanDedxNeg3DPartITS = new TH3F("Flow_MeanDedxNeg3DPart_ITS","Flow_MeanDedxNeg3DPart_ITS", kMomenBins, logpMin, logpMax, kDedxBins, 0., dEdxMax, AliFlowConstants::kPid, 0., AliFlowConstants::kPid);
807  fHistMeanDedxNeg3DPartITS->SetXTitle("log(momentum) (GeV/c)");
808  fHistMeanDedxNeg3DPartITS->SetYTitle("mean dEdx (-)");
809  fHistMeanDedxNeg3DPartITS->SetZTitle("e , mu , pi , k , p , d");
810
811  // MeanSignalPos TRD
812  fHistMeanDedxPos2DTRD = new TH2F("Flow_MeanSignalPos2D_TRD","Flow_MeanSignalPos2D_TRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
813  fHistMeanDedxPos2DTRD->SetXTitle("log(momentum) (GeV/c)");
814  fHistMeanDedxPos2DTRD->SetYTitle("mean TRD (+)");
815  
816  // MeanSignalNeg TRD
817  fHistMeanDedxNeg2DTRD = new TH2F("Flow_MeanSignalNeg2D_TRD","Flow_MeanSignalNeg2D_TRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
818  fHistMeanDedxNeg2DTRD->SetXTitle("log(momentum) (GeV/c)");
819  fHistMeanDedxNeg2DTRD->SetYTitle("mean TRD (-)");
820
821  // MeanTimePos TOF
822  fHistMeanDedxPos2DTOF = new TH2F("Flow_MeanTimePos2D_TOF","Flow_MeanTimePos2D_TOF", kMomenBins, logpMin, logpMax, kTofBins, tofmin, tofmax);
823  fHistMeanDedxPos2DTOF->SetXTitle("log(momentum) (GeV/c)");
824  fHistMeanDedxPos2DTOF->SetYTitle("mean Time of Flight (+)");
825  
826  // MeanTimeNeg TOF
827  fHistMeanDedxNeg2DTOF = new TH2F("Flow_MeanTimeNeg2D_TOF","Flow_MeanTimeNeg2D_TOF", kMomenBins, logpMin, logpMax, kTofBins, tofmin, tofmax);
828  fHistMeanDedxNeg2DTOF->SetXTitle("log(momentum) (GeV/c)");
829  fHistMeanDedxNeg2DTOF->SetYTitle("mean Time of Flight (-)");
830  
831  // Detector response for each particle { ... 
832  // MeanDedx PiPlus - TPC
833  fHistMeanTPCPiPlus = new TH2F("Flow_MeanDedxPiPlusTPC", "Flow_MeanDedxPiPlusTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
834  fHistMeanTPCPiPlus->SetXTitle("Log10(p) (GeV/c)");
835  fHistMeanTPCPiPlus->SetYTitle("mean dEdx (pi+)");
836  // MeanDedx PiMinus
837  fHistMeanTPCPiMinus = new TH2F("Flow_MeanDedxPiMinusTPC", "Flow_MeanDedxPiMinusTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
838  fHistMeanTPCPiMinus->SetXTitle("Log10(p) (GeV/c)");
839  fHistMeanTPCPiMinus->SetYTitle("mean dEdx (pi-)");
840  // MeanDedx Proton
841  fHistMeanTPCProton = new TH2F("Flow_MeanDedxProtonTPC", "Flow_MeanDedxProtonTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
842  fHistMeanTPCProton->SetXTitle("Log10(p) (GeV/c)");
843  fHistMeanTPCProton->SetYTitle("mean dEdx (pr+)");
844  // MeanDedx Pbar
845  fHistMeanTPCPbar = new TH2F("Flow_MeanDedxPbarTPC", "Flow_MeanDedxPbarTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
846  fHistMeanTPCPbar->SetXTitle("Log10(p) (GeV/c)");
847  fHistMeanTPCPbar->SetYTitle("mean dEdx (pr-)");
848  // MeanDedx Kplus
849  fHistMeanTPCKplus = new TH2F("Flow_MeanDedxKplusTPC", "Flow_MeanDedxKplusTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
850  fHistMeanTPCKplus->SetXTitle("Log10(p) (GeV/c)");
851  fHistMeanTPCKplus->SetYTitle("mean dEdx (K+)");
852  // MeanDedx Kminus
853  fHistMeanTPCKminus = new TH2F("Flow_MeanDedxKminusTPC", "Flow_MeanDedxKminusTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
854  fHistMeanTPCKminus->SetXTitle("Log10(p) (GeV/c)");
855  fHistMeanTPCKminus->SetYTitle("mean dEdx (K-)");
856  // MeanDedx Deuteron
857  fHistMeanTPCDeuteron = new TH2F("Flow_MeanDedxDeuteronTPC", "Flow_MeanDedxDeuteronTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
858  fHistMeanTPCDeuteron->SetXTitle("Log10(p) (GeV/c)");
859  fHistMeanTPCDeuteron->SetYTitle("mean dEdx (d+)");
860  // MeanDedx AntiDeuteron
861  fHistMeanTPCAntiDeuteron = new TH2F("Flow_MeanDedxAntiDeuteronTPC", "Flow_MeanDedxAntiDeuteronTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
862  fHistMeanTPCAntiDeuteron->SetXTitle("Log10(p) (GeV/c)");
863  fHistMeanTPCAntiDeuteron->SetYTitle("mean dEdx (d-)");
864  // MeanDedxElectron
865  fHistMeanTPCElectron = new TH2F("Flow_MeanDedxElectronTPC", "Flow_MeanDedxElectronTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
866  fHistMeanTPCElectron->SetXTitle("Log10(p) (GeV/c)");
867  fHistMeanTPCElectron->SetYTitle("mean dEdx (e-)");
868  // MeanDedx Positron
869  fHistMeanTPCPositron = new TH2F("Flow_MeanDedxPositronTPC", "Flow_MeanDedxPositronTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
870  fHistMeanTPCPositron->SetXTitle("Log10(p) (GeV/c)");
871  fHistMeanTPCPositron->SetYTitle("mean dEdx (e+)");
872  // MeanDedx MuonPlus
873  fHistMeanTPCMuonPlus = new TH2F("Flow_MeanDedxMuonPlusTPC", "Flow_MeanDedxMuonPlusTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
874  fHistMeanTPCMuonPlus->SetXTitle("Log10(p) (GeV/c)");
875  fHistMeanTPCMuonPlus->SetYTitle("mean dEdx (mu+)");
876  // MeanDedx MuonMinus
877  fHistMeanTPCMuonMinus = new TH2F("Flow_MeanDedxMuonMinusTPC", "Flow_MeanDedxMuonMinusTPC", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMaxTPC);
878  fHistMeanTPCMuonMinus->SetXTitle("Log10(p) (GeV/c)");
879  fHistMeanTPCMuonMinus->SetYTitle("mean dEdx (mu-)");
880  
881  // MeanDedx PiPlus - ITS
882  fHistMeanITSPiPlus = new TH2F("Flow_MeanDedxPiPlusITS", "Flow_MeanDedxPiPlusITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
883  fHistMeanITSPiPlus->SetXTitle("Log10(p) (GeV/c)");
884  fHistMeanITSPiPlus->SetYTitle("mean dEdx (pi+)");
885  // MeanDedx PiMinus
886  fHistMeanITSPiMinus = new TH2F("Flow_MeanDedxPiMinusITS", "Flow_MeanDedxPiMinusITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
887  fHistMeanITSPiMinus->SetXTitle("Log10(p) (GeV/c)");
888  fHistMeanITSPiMinus->SetYTitle("mean dEdx (pi-)");
889  // MeanDedx Proton
890  fHistMeanITSProton = new TH2F("Flow_MeanDedxProtonITS", "Flow_MeanDedxProtonITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
891  fHistMeanITSProton->SetXTitle("Log10(p) (GeV/c)");
892  fHistMeanITSProton->SetYTitle("mean dEdx (pr+)");
893  // MeanDedx Pbar
894  fHistMeanITSPbar = new TH2F("Flow_MeanDedxPbarITS", "Flow_MeanDedxPbarITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
895  fHistMeanITSPbar->SetXTitle("Log10(p) (GeV/c)");
896  fHistMeanITSPbar->SetYTitle("mean dEdx (pr-)");
897  // MeanDedx Kplus
898  fHistMeanITSKplus = new TH2F("Flow_MeanDedxKplusITS", "Flow_MeanDedxKplusITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
899  fHistMeanITSKplus->SetXTitle("Log10(p) (GeV/c)");
900  fHistMeanITSKplus->SetYTitle("mean dEdx (K+)");
901  // MeanDedx Kminus
902  fHistMeanITSKminus = new TH2F("Flow_MeanDedxKminusITS", "Flow_MeanDedxKminusITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
903  fHistMeanITSKminus->SetXTitle("Log10(p) (GeV/c)");
904  fHistMeanITSKminus->SetYTitle("mean dEdx (K-)");
905  // MeanDedx Deuteron
906  fHistMeanITSDeuteron = new TH2F("Flow_MeanDedxDeuteronITS", "Flow_MeanDedxDeuteronITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
907  fHistMeanITSDeuteron->SetXTitle("Log10(p) (GeV/c)");
908  fHistMeanITSDeuteron->SetYTitle("mean dEdx (d+)");
909  // MeanDedx AntiDeuteron
910  fHistMeanITSAntiDeuteron = new TH2F("Flow_MeanDedxAntiDeuteronITS", "Flow_MeanDedxAntiDeuteronITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
911  fHistMeanITSAntiDeuteron->SetXTitle("Log10(p) (GeV/c)");
912  fHistMeanITSAntiDeuteron->SetYTitle("mean dEdx (d-)");
913  // MeanDedx Electron
914  fHistMeanITSElectron = new TH2F("Flow_MeanDedxElectronITS", "Flow_MeanDedxElectronITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
915  fHistMeanITSElectron->SetXTitle("Log10(p) (GeV/c)");
916  fHistMeanITSElectron->SetYTitle("mean dEdx (e-)");
917  // MeanDedx Positron
918  fHistMeanITSPositron = new TH2F("Flow_MeanDedxPositronITS", "Flow_MeanDedxPositronITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
919  fHistMeanITSPositron->SetXTitle("Log10(p) (GeV/c)");
920  fHistMeanITSPositron->SetYTitle("mean dEdx (e+)");
921  // MeanDedx MuonPlus
922  fHistMeanITSMuonPlus = new TH2F("Flow_MeanDedxMuonPlusITS", "Flow_MeanDedxMuonPlusITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
923  fHistMeanITSMuonPlus->SetXTitle("Log10(p) (GeV/c)");
924  fHistMeanITSMuonPlus->SetYTitle("mean dEdx (mu+)");
925  // MeanDedx MuonMinus
926  fHistMeanITSMuonMinus = new TH2F("Flow_MeanDedxMuonMinusITS", "Flow_MeanDedxMuonMinusITS", kMomenBins, logpMin, logpMax, kDedxBins, 0, dEdxMax);
927  fHistMeanITSMuonMinus->SetXTitle("Log10(p) (GeV/c)");
928  fHistMeanITSMuonMinus->SetYTitle("mean dEdx (mu-)");
929  
930  // MeanTrd PiPlus - TRD
931  fHistMeanTRDPiPlus = new TH2F("Flow_MeanTrdPiPlusTRD", "Flow_MeanTrdPiPlusTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
932  fHistMeanTRDPiPlus->SetXTitle("Log10(p) (GeV/c)");
933  fHistMeanTRDPiPlus->SetYTitle("mean t.r.[] (pi+)");
934  // MeanTrd PiMinus
935  fHistMeanTRDPiMinus = new TH2F("Flow_MeanTrdPiMinusTRD", "Flow_MeanTrdPiMinusTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
936  fHistMeanTRDPiMinus->SetXTitle("Log10(p) (GeV/c)");
937  fHistMeanTRDPiMinus->SetYTitle("mean t.r.[] (pi-)");
938  // MeanTrd Proton
939  fHistMeanTRDProton = new TH2F("Flow_MeanTrdProtonTRD", "Flow_MeanTrdProtonTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
940  fHistMeanTRDProton->SetXTitle("Log10(p) (GeV/c)");
941  fHistMeanTRDProton->SetYTitle("mean t.r.[] (pr+)");
942  // MeanTrd Pbar
943  fHistMeanTRDPbar = new TH2F("Flow_MeanTrdPbarTRD", "Flow_MeanTrdPbarTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
944  fHistMeanTRDPbar->SetXTitle("Log10(p) (GeV/c)");
945  fHistMeanTRDPbar->SetYTitle("mean t.r.[] (pr-)");
946  // MeanTrd Kplus
947  fHistMeanTRDKplus = new TH2F("Flow_MeanTrdKplusTRD", "Flow_MeanTrdKplusTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
948  fHistMeanTRDKplus->SetXTitle("Log10(p) (GeV/c)");
949  fHistMeanTRDKplus->SetYTitle("mean t.r.[] (K+)");
950  // MeanTrd Kminus
951  fHistMeanTRDKminus = new TH2F("Flow_MeanTrdKminusTRD", "Flow_MeanTrdKminusTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
952  fHistMeanTRDKminus->SetXTitle("Log10(p) (GeV/c)");
953  fHistMeanTRDKminus->SetYTitle("mean t.r.[] (K-)");
954  // MeanTrd Deuteron
955  fHistMeanTRDDeuteron = new TH2F("Flow_MeanTrdDeuteronTRD", "Flow_MeanTrdDeuteronTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
956  fHistMeanTRDDeuteron->SetXTitle("Log10(p) (GeV/c)");
957  fHistMeanTRDDeuteron->SetYTitle("mean t.r.[] (d+)");
958  // MeanTrd AntiDeuteron
959  fHistMeanTRDAntiDeuteron = new TH2F("Flow_MeanTrdAntiDeuteronTRD", "Flow_MeanTrdAntiDeuteronTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
960  fHistMeanTRDAntiDeuteron->SetXTitle("Log10(p) (GeV/c)");
961  fHistMeanTRDAntiDeuteron->SetYTitle("mean t.r.[] (d-)");
962  // MeanTrd Electron
963  fHistMeanTRDElectron = new TH2F("Flow_MeanTrdElectronTRD", "Flow_MeanTrdElectronTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
964  fHistMeanTRDElectron->SetXTitle("Log10(p) (GeV/c)");
965  fHistMeanTRDElectron->SetYTitle("mean t.r.[] (e-)");
966  // MeanTrd Positron
967  fHistMeanTRDPositron = new TH2F("Flow_MeanTrdPositronTRD", "Flow_MeanTrdPositronTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
968  fHistMeanTRDPositron->SetXTitle("Log10(p) (GeV/c)");
969  fHistMeanTRDPositron->SetYTitle("mean t.r.[] (e+)");
970  // MeanTrd MuonPlus
971  fHistMeanTRDMuonPlus = new TH2F("Flow_MeanTrdMuonPlusTRD", "Flow_MeanTrdMuonPlusTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
972  fHistMeanTRDMuonPlus->SetXTitle("Log10(p) (GeV/c)");
973  fHistMeanTRDMuonPlus->SetYTitle("mean t.r.[] (mu+)");
974  // MeanTrd MuonMinus
975  fHistMeanTRDMuonMinus = new TH2F("Flow_MeanTrdMuonMinusTRD", "Flow_MeanTrdMuonMinusTRD", kMomenBins, logpMin, logpMax, kTrdBins, 0, trdmax);
976  fHistMeanTRDMuonMinus->SetXTitle("Log10(p) (GeV/c)");
977  fHistMeanTRDMuonMinus->SetYTitle("mean t.r.[] (mu-)");
978  
979  // T.O.F. PiPlus - TOF
980  fHistMeanTOFPiPlus = new TH2F("Flow_MeanTofPiPlusTOF", "Flow_MeanTofPiPlusTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
981  fHistMeanTOFPiPlus->SetXTitle("invariant mass (GeV)");
982  fHistMeanTOFPiPlus->SetYTitle("mean t.o.f.[psec] (pi+)");
983  // MeanTof PiMinus
984  fHistMeanTOFPiMinus = new TH2F("Flow_MeanTofPiMinusTOF", "Flow_MeanTofPiMinusTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
985  fHistMeanTOFPiMinus->SetXTitle("invariant mass (GeV)");
986  fHistMeanTOFPiMinus->SetYTitle("mean t.o.f.[psec] (pi-)");
987  // MeanTof Proton
988  fHistMeanTOFProton = new TH2F("Flow_MeanTofProtonTOF", "Flow_MeanTofProtonTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
989  fHistMeanTOFProton->SetXTitle("invariant mass (GeV)");
990  fHistMeanTOFProton->SetYTitle("mean t.o.f.[psec] (pr+)");
991  // Mean TofPbar
992  fHistMeanTOFPbar = new TH2F("Flow_MeanTofPbarTOF", "Flow_MeanTofPbarTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
993  fHistMeanTOFPbar->SetXTitle("invariant mass (GeV)");
994  fHistMeanTOFPbar->SetYTitle("mean t.o.f.[psec] (pr-)");
995  // mean t.o.f.[psec]Kplus
996  fHistMeanTOFKplus = new TH2F("Flow_MeanTofKplusTOF", "Flow_MeanTofKplusTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
997  fHistMeanTOFKplus->SetXTitle("invariant mass (GeV)");
998  fHistMeanTOFKplus->SetYTitle("mean t.o.f.[psec] (K+)");
999  // mean t.o.f.[psec]Kminus
1000  fHistMeanTOFKminus = new TH2F("Flow_MeanTofKminusTOF", "Flow_MeanTofKminusTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
1001  fHistMeanTOFKminus->SetXTitle("invariant mass (GeV)");
1002  fHistMeanTOFKminus->SetYTitle("mean t.o.f.[psec] (K-)");
1003  // MeanTof Deuteron
1004  fHistMeanTOFDeuteron = new TH2F("Flow_MeanTofDeuteronTOF", "Flow_MeanTofDeuteronTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
1005  fHistMeanTOFDeuteron->SetXTitle("invariant mass (GeV)");
1006  fHistMeanTOFDeuteron->SetYTitle("mean t.o.f.[psec] (d+)");
1007  // MeanTof AntiDeuteron
1008  fHistMeanTOFAntiDeuteron = new TH2F("Flow_MeanTofAntiDeuteronTOF", "Flow_MeanTofAntiDeuteronTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
1009  fHistMeanTOFAntiDeuteron->SetXTitle("invariant mass (GeV)");
1010  fHistMeanTOFAntiDeuteron->SetYTitle("mean t.o.f.[psec] (d-)");
1011  // MeanTof Electron
1012  fHistMeanTOFElectron = new TH2F("Flow_MeanTofElectronTOF", "Flow_MeanTofElectronTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
1013  fHistMeanTOFElectron->SetXTitle("invariant mass (GeV)");
1014  fHistMeanTOFElectron->SetYTitle("mean t.o.f.[psec] (e-)");
1015  // MeanTof Positron
1016  fHistMeanTOFPositron = new TH2F("Flow_MeanTofPositronTOF", "Flow_MeanTofPositronTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
1017  fHistMeanTOFPositron->SetXTitle("invariant mass (GeV)");
1018  fHistMeanTOFPositron->SetYTitle("mean t.o.f.[psec] (e+)");
1019  // MeanTof MuonPlus
1020  fHistMeanTOFMuonPlus = new TH2F("Flow_MeanTofMuonPlusTOF", "Flow_MeanTofMuonPlusTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
1021  fHistMeanTOFMuonPlus->SetXTitle("invariant mass (GeV)");
1022  fHistMeanTOFMuonPlus->SetYTitle("mean t.o.f.[psec] (mu+)");
1023  // MeanTof MuonMinus
1024  fHistMeanTOFMuonMinus = new TH2F("Flow_MeanTofMuonMinusTOF", "Flow_MeanTofMuonMinusTOF", kMassBins, massMin, massMax, kTofBins, tofmin, tofmax);
1025  fHistMeanTOFMuonMinus->SetXTitle("invariant mass (GeV)");
1026  fHistMeanTOFMuonMinus->SetYTitle("mean t.o.f.[psec] (mu-)");
1027  // ... }
1028
1029  TString* histTitle;
1030  for (int n = 0; n < AliFlowConstants::kSubs; n++) // for sub-events
1031  {
1032   for (int k = 0; k < AliFlowConstants::kSels; k++) 
1033   {
1034    for (int j = 0; j < AliFlowConstants::kHars; j++) 
1035    {
1036     float order = (float)(j + 1);
1037     int i = AliFlowConstants::kSubs * k + n ;
1038
1039     // event planes
1040     histTitle = new TString("Flow_Psi_Sub");
1041     *histTitle += n+1;
1042     histTitle->Append("_Sel");
1043     *histTitle += k+1;
1044     histTitle->Append("_Har");
1045     *histTitle += j+1;
1046     fHistSub[i].fHistSubHar[j].fHistPsiSubs = new TH1F(histTitle->Data(),histTitle->Data(), kPsiBins, psiMin, (psiMax / order));
1047     fHistSub[i].fHistSubHar[j].fHistPsiSubs->SetXTitle("Event Plane Angle (rad)");
1048     fHistSub[i].fHistSubHar[j].fHistPsiSubs->SetYTitle("Counts");
1049     delete histTitle;
1050    }
1051   }
1052  }
1053  
1054  if(fV0loop)            // All V0s (if there, if flag on)
1055  {
1056  // Mass
1057   fHistV0Mass = new TH1F("FlowV0_InvMass", "FlowV0_InvMass", kMassBins, massMin, massMax);
1058   fHistV0Mass->SetXTitle("Invariant Mass (GeV)");
1059   fHistV0Mass->SetYTitle("Counts");
1060  // Distance of closest approach
1061   fHistV0Dca = new TH1F("FlowV0_Dca", "FlowV0_Dca", kDcaBins, dcaMin, glDcaMax);
1062   fHistV0Dca->SetXTitle("dca between tracks (cm)");
1063   fHistV0Dca->SetYTitle("Counts");
1064  // lenght
1065   fHistV0Lenght = new TH1F("FlowV0_Lenght", "FlowV0_Lenght", kLgBins, lgMinV0, lgMaxV0);
1066   fHistV0Lenght->SetXTitle("Distance of V0s (cm)");
1067   fHistV0Lenght->SetYTitle("Counts");
1068  // Sigma for all particles
1069   fHistV0Sigma = new TH1F("FlowV0_Sigma", "FlowV0_Sigma", kLgBins, lgMinV0, lgMaxV0 );
1070   fHistV0Sigma->SetXTitle("Sigma");
1071   fHistV0Sigma->SetYTitle("Counts");
1072  // Chi2 
1073   fHistV0Chi2 = new TH1F("FlowV0_Chi2", "FlowV0_Chi2", kChi2Bins, chi2Min, chi2MaxC);
1074   fHistV0Chi2->SetXTitle("Chi square at Main Vertex");
1075   fHistV0Chi2->SetYTitle("Counts");
1076  // EtaPtPhi
1077   fHistV0EtaPtPhi3D = new TH3F("FlowV0_EtaPtPhi3D", "FlowV0_EtaPtPhi3D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
1078   fHistV0EtaPtPhi3D->SetXTitle("Eta");
1079   fHistV0EtaPtPhi3D->SetYTitle("Pt (GeV/c)");
1080   fHistV0EtaPtPhi3D->SetZTitle("Phi (rad)");
1081  // Yield for all v0s
1082   fHistV0YieldAll2D = new TH2D("FlowV0_YieldAll2D", "FlowV0_YieldAll2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1083   fHistV0YieldAll2D->Sumw2();
1084   fHistV0YieldAll2D->SetXTitle("Pseudorapidty");
1085   fHistV0YieldAll2D->SetYTitle("Pt (GeV/c)");
1086  // Mass slices on pT
1087   fHistV0MassPtSlices = new TH2D("FlowV0_MassPtSlices", "FlowV0_MassPtSlices", kMassBins, massMin, massMax, fPtBins, fPtMin, fPtMax);
1088   fHistV0MassPtSlices->Sumw2();
1089   fHistV0MassPtSlices->SetXTitle("Invariant Mass (GeV)");
1090   fHistV0MassPtSlices->SetYTitle("Pt (GeV/c)");
1091  // Yield v0s (total P and Rapidity)
1092   fHistV0PYall2D = new TH2D("FlowV0_PYall2D", "FlowV0_PYall2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1093   fHistV0PYall2D->Sumw2();
1094   fHistV0PYall2D->SetXTitle("Rapidty");
1095   fHistV0PYall2D->SetYTitle("P (GeV/c)");
1096
1097  // Selected V0s ...  
1098  // Yield 
1099   fHistV0YieldPart2D = new TH2D("FlowV0_Yield2Dsel", "FlowV0_Yield2Dsel", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1100   fHistV0YieldPart2D->Sumw2();
1101   fHistV0YieldPart2D->SetXTitle("Pseudorapidty");
1102   fHistV0YieldPart2D->SetYTitle("Pt (GeV/c)");
1103  // Mass Window
1104   fHistV0MassWin = new TH1F("FlowV0_MassWinPart", "FlowV0_MassWinPart", kMassBins, massMin, massMax);
1105   fHistV0MassWin->SetXTitle("Invariant Mass (GeV)");
1106   fHistV0MassWin->SetYTitle("Counts");
1107  // EtaPtPhi
1108   fHistV0EtaPtPhi3DPart = new TH3F("FlowV0_EtaPtPhi3Dpart", "FlowV0_EtaPtPhi3Dpart", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
1109   fHistV0EtaPtPhi3DPart->SetXTitle("Eta");
1110   fHistV0EtaPtPhi3DPart->SetYTitle("Pt (GeV/c)");
1111   fHistV0EtaPtPhi3DPart->SetZTitle("Phi (rad)");
1112  // Distance of closest approach
1113   fHistV0DcaPart = new TH1F("FlowV0_DcaPart", "FlowV0_DcaPart", kDcaBins, dcaMin, dcaMax);
1114   fHistV0DcaPart->Sumw2();
1115   fHistV0DcaPart->SetXTitle("dca between tracks (cm)");
1116   fHistV0DcaPart->SetYTitle("Counts");
1117  // lenght
1118   fHistV0LenghtPart = new TH1F("FlowV0_LenghtPart", "FlowV0_LenghtPart", kLgBins, lgMinV0, lgMaxV0);
1119   fHistV0LenghtPart->SetXTitle("Distance of V0s (cm)");
1120   fHistV0LenghtPart->SetYTitle("Counts");
1121  // SideBand Mass (sidebands)
1122   fHistV0sbMassSide = new TH1F("FlowV0sb_MassWinSideBands", "FlowV0sb_MassWinSideBands", kMassBins, massMin, massMax);
1123   fHistV0sbMassSide->SetXTitle("Invariant Mass (GeV)");
1124   fHistV0sbMassSide->SetYTitle("Counts");
1125  // EtaPtPhi (sidebands)
1126   fHistV0sbEtaPtPhi3DPart = new TH3F("FlowV0sb_EtaPtPhi3D", "FlowV0sb_EtaPtPhi3D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
1127   fHistV0sbEtaPtPhi3DPart->SetXTitle("Eta");
1128   fHistV0sbEtaPtPhi3DPart->SetYTitle("Pt (GeV/c)");
1129   fHistV0sbEtaPtPhi3DPart->SetZTitle("Phi (rad)");
1130  // Distance of closest approach (sidebands)
1131   fHistV0sbDcaPart = new TH1F("FlowV0sb_Dca", "FlowV0sb_Dca", kDcaBins, dcaMin, dcaMax);
1132   fHistV0sbDcaPart->Sumw2();
1133   fHistV0sbDcaPart->SetXTitle("dca between tracks (cm)");
1134   fHistV0sbDcaPart->SetYTitle("Counts");
1135  // lenght (sidebands)
1136   fHistV0sbLenghtPart = new TH1F("FlowV0sb_Lenght", "FlowV0sb_Lenght", kLgBins, lgMinV0, lgMaxV0);
1137   fHistV0sbLenghtPart->SetXTitle("Distance of V0s (cm)");
1138   fHistV0sbLenghtPart->SetYTitle("Counts");
1139
1140  // Mean Eta in each bin
1141   fHistV0BinEta = new TProfile("FlowV0_Bin_Eta", "FlowV0_Bin_Eta", fEtaBins, fEtaMin, fEtaMax, fEtaMin, fEtaMax, "");
1142   fHistV0BinEta->SetXTitle((char*)fLabel.Data());
1143   fHistV0BinEta->SetYTitle("<Eta>");
1144  // Mean Pt in each bin
1145   fHistV0BinPt = new TProfile("FlowV0_Bin_Pt", "FlowV0_Bin_Pt", fPtBinsPart, fPtMin, fPtMaxPart, fPtMin, fPtMaxPart, "");
1146   fHistV0BinPt->SetXTitle("Pt (GeV/c)");
1147   fHistV0BinPt->SetYTitle("<Pt> (GeV/c)");
1148  // Mean Eta in each bin (sidebands)
1149   fHistV0sbBinEta = new TProfile("FlowV0sb_Bin_Eta", "FlowV0sb_Bin_Eta", fEtaBins, fEtaMin, fEtaMax, fEtaMin, fEtaMax, "");
1150   fHistV0sbBinEta->SetXTitle((char*)fLabel.Data());
1151   fHistV0sbBinEta->SetYTitle("<Eta>");
1152  // Mean Pt in each bin (sidebands)
1153   fHistV0sbBinPt = new TProfile("FlowV0sb_Bin_Pt", "FlowV0sb_Bin_Pt", fPtBinsPart, fPtMin, fPtMaxPart, fPtMin, fPtMaxPart, "");
1154   fHistV0sbBinPt->SetXTitle("Pt (GeV/c)");
1155   fHistV0sbBinPt->SetYTitle("<Pt> (GeV/c)");
1156  }
1157
1158  for (int k = 0; k < AliFlowConstants::kSels; k++) // for each selection
1159  {
1160   // cos(n*delta_Psi)
1161   histTitle = new TString("Flow_Cos_Sel");
1162   *histTitle += k+1;
1163   fHistFull[k].fHistCos = new TProfile(histTitle->Data(), histTitle->Data(), AliFlowConstants::kHars, 0.5, (float)(AliFlowConstants::kHars) + 0.5, -1., 1., "");
1164   fHistFull[k].fHistCos->SetXTitle("Harmonic");
1165   fHistFull[k].fHistCos->SetYTitle("<cos(n*delta_Psi)>");
1166   delete histTitle;
1167    
1168   // resolution
1169   histTitle = new TString("Flow_Res_Sel");
1170   *histTitle += k+1;
1171   fHistFull[k].fHistRes = new TH1F(histTitle->Data(), histTitle->Data(), AliFlowConstants::kHars, 0.5, (float)(AliFlowConstants::kHars) + 0.5);
1172   fHistFull[k].fHistRes->SetXTitle("Harmonic");
1173   fHistFull[k].fHistRes->SetYTitle("Resolution");
1174   delete histTitle;
1175
1176   // vObs
1177   histTitle = new TString("Flow_vObs_Sel");
1178   *histTitle += k+1;
1179   fHistFull[k].fHistvObs = new TProfile(histTitle->Data(), histTitle->Data(), AliFlowConstants::kHars, 0.5, (float)(AliFlowConstants::kHars) + 0.5, -100., 100., "");
1180   fHistFull[k].fHistvObs->SetXTitle("Harmonic");
1181   fHistFull[k].fHistvObs->SetYTitle("vObs (%)");
1182   delete histTitle;
1183
1184   // vObs V0
1185   histTitle = new TString("FlowV0_vObs_Sel");
1186   *histTitle += k+1;
1187   fHistFull[k].fHistV0vObs = new TProfile(histTitle->Data(), histTitle->Data(), AliFlowConstants::kHars, 0.5, (float)(AliFlowConstants::kHars) + 0.5, -100., 100., "");
1188   fHistFull[k].fHistV0vObs->SetXTitle("Harmonic");
1189   fHistFull[k].fHistV0vObs->SetYTitle("vObs (%)");
1190   delete histTitle;
1191    
1192   // vObs V0 sideband SX
1193   histTitle = new TString("FlowV0sb_vObs_sx_Sel");
1194   *histTitle += k+1;
1195   fHistFull[k].fHistV0sbvObsSx = new TProfile(histTitle->Data(), histTitle->Data(), AliFlowConstants::kHars, 0.5, (float)(AliFlowConstants::kHars) + 0.5, -100., 100., "");
1196   fHistFull[k].fHistV0sbvObsSx->SetXTitle("Harmonic");
1197   fHistFull[k].fHistV0sbvObsSx->SetYTitle("vObs (%)");
1198   delete histTitle;
1199    
1200   // vObs V0 sideband DX
1201   histTitle = new TString("FlowV0sb_vObs_dx_Sel");
1202   *histTitle += k+1;
1203   fHistFull[k].fHistV0sbvObsDx = new TProfile(histTitle->Data(), histTitle->Data(), AliFlowConstants::kHars, 0.5, (float)(AliFlowConstants::kHars) + 0.5, -100., 100., "");
1204   fHistFull[k].fHistV0sbvObsDx->SetXTitle("Harmonic");
1205   fHistFull[k].fHistV0sbvObsDx->SetYTitle("vObs (%)");
1206   delete histTitle;
1207    
1208   // PID for tracks used in R.P.
1209   histTitle = new TString("Flow_BayPidMult_Sel");
1210   *histTitle += k+1;
1211   fHistFull[k].fHistBayPidMult = new TH1F(histTitle->Data(), histTitle->Data(),AliFlowConstants::kPid,-0.5,((float)AliFlowConstants::kPid-0.5));
1212   fHistFull[k].fHistBayPidMult->Sumw2() ;
1213   fHistFull[k].fHistBayPidMult->SetXTitle("e+/-  ,  mu+/-  ,  pi+/-  ,  K+/-  ,  p+/-  ,  d+/- ");
1214   fHistFull[k].fHistBayPidMult->SetYTitle("Counts");
1215   delete histTitle;
1216
1217   for (int j = 0; j < AliFlowConstants::kHars; j++)   // for each harmonic
1218   {
1219    float order  = (float)(j+1);
1220
1221    // multiplicity
1222    histTitle = new TString("Flow_Mul_Sel");
1223    *histTitle += k+1;
1224    histTitle->Append("_Har");
1225    *histTitle += j+1;
1226    fHistFull[k].fHistFullHar[j].fHistMult = new TH1F(histTitle->Data(),histTitle->Data(), kMultBins, multMin, multMax);
1227    fHistFull[k].fHistFullHar[j].fHistMult->SetXTitle("Multiplicity");
1228    fHistFull[k].fHistFullHar[j].fHistMult->SetYTitle("Counts");
1229    delete histTitle;
1230
1231    // event plane
1232    histTitle = new TString("Flow_Psi_Sel");
1233    *histTitle += k+1;
1234    histTitle->Append("_Har");
1235    *histTitle += j+1;
1236    fHistFull[k].fHistFullHar[j].fHistPsi = new TH1F(histTitle->Data(), histTitle->Data(), kPsiBins, psiMin, psiMax / order);
1237    fHistFull[k].fHistFullHar[j].fHistPsi->SetXTitle("Event Plane Angle (rad)");
1238    fHistFull[k].fHistFullHar[j].fHistPsi->SetYTitle("Counts");
1239    delete histTitle;
1240      
1241    // event plane difference of two selections
1242    histTitle = new TString("Flow_Psi_Diff_Sel");
1243    *histTitle += k+1;
1244    histTitle->Append("_Har");
1245    *histTitle += j+1;
1246    if (k == 0 ) 
1247    {
1248     Int_t myOrder = 1;
1249     if (j == 1) { myOrder = 2 ; }
1250     fHistFull[k].fHistFullHar[j].fHistPsiDiff = new TH1F(histTitle->Data(), histTitle->Data(), kPsiBins, -psiMax/myOrder/2., psiMax/myOrder/2.);
1251    } 
1252    else 
1253    {
1254     fHistFull[k].fHistFullHar[j].fHistPsiDiff = new TH1F(histTitle->Data(), histTitle->Data(), kPsiBins, -psiMax/2., psiMax/2.);
1255    }
1256    if (k == 0) 
1257    {
1258     if (j == 0) 
1259     {
1260      fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{1,Sel2}(rad)");
1261     } 
1262     else if (j == 1) 
1263     {
1264      fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{2,Sel1} - #Psi_{2,Sel2}(rad)");
1265     }
1266    } 
1267    else if (k == 1) 
1268    {
1269     if (j == 0) 
1270     {  
1271      fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{2,Sel2}(rad)");
1272     } 
1273     else if (j == 1) 
1274     {
1275      fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{2,Sel1}(rad)");
1276     }
1277    }
1278    fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetYTitle("Counts");
1279    delete histTitle;
1280
1281    // correlation of sub-event planes
1282    histTitle = new TString("Flow_Psi_Sub_Corr_Sel");
1283    *histTitle += k+1;
1284    histTitle->Append("_Har");
1285    *histTitle += j+1;
1286    fHistFull[k].fHistFullHar[j].fHistPsiSubCorr = new TH1F(histTitle->Data(), histTitle->Data(), kPsiBins, psiMin, psiMax / order);
1287    fHistFull[k].fHistFullHar[j].fHistPsiSubCorr->Sumw2();
1288    fHistFull[k].fHistFullHar[j].fHistPsiSubCorr->SetXTitle("Sub-Event Correlation (rad)");
1289    fHistFull[k].fHistFullHar[j].fHistPsiSubCorr->SetYTitle("Counts");
1290    delete histTitle;
1291      
1292    // correlation of sub-event planes of different order
1293    histTitle = new TString("Flow_Psi_Sub_Corr_Diff_Sel");
1294    *histTitle += k+1;
1295    histTitle->Append("_Har");
1296    *histTitle += j+1;
1297    fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff = new TH1F(histTitle->Data(), histTitle->Data(), kPsiBins, psiMin, psiMax / (order+1.));
1298    fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->Sumw2();
1299    fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->SetXTitle("Sub-Event Correlation (rad)");
1300    fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->SetYTitle("Counts");
1301    delete histTitle;
1302      
1303    // q
1304    histTitle = new TString("Flow_NormQ_Sel");
1305    *histTitle += k+1;
1306    histTitle->Append("_Har");
1307    *histTitle += j+1;
1308    fHistFull[k].fHistFullHar[j].fHistQnorm = new TH1F(histTitle->Data(), histTitle->Data(), kQbins, qMin, qMax);
1309    fHistFull[k].fHistFullHar[j].fHistQnorm->Sumw2();
1310    fHistFull[k].fHistFullHar[j].fHistQnorm->SetXTitle("q = |Q|/sqrt(Mult)");
1311    fHistFull[k].fHistFullHar[j].fHistQnorm->SetYTitle("Counts");
1312    delete histTitle;
1313
1314     // particle-plane azimuthal correlation
1315    histTitle = new TString("Flow_Phi_Corr_Sel");
1316    *histTitle += k+1;
1317    histTitle->Append("_Har");
1318    *histTitle += j+1;
1319    fHistFull[k].fHistFullHar[j].fHistPhiCorr = new TH1F(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax / order);
1320    fHistFull[k].fHistFullHar[j].fHistPhiCorr->Sumw2();
1321    fHistFull[k].fHistFullHar[j].fHistPhiCorr->SetXTitle("Particle-Plane Correlation (rad)");
1322    fHistFull[k].fHistFullHar[j].fHistPhiCorr->SetYTitle("Counts");
1323    delete histTitle;
1324      
1325    // neutral particle-plane azimuthal correlation
1326    histTitle = new TString("FlowV0_Phi_Corr_Sel");
1327    *histTitle += k+1;
1328    histTitle->Append("_Har");
1329    *histTitle += j+1;
1330    fHistFull[k].fHistFullHar[j].fHistV0PhiCorr = new TH1F(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax / order);
1331    fHistFull[k].fHistFullHar[j].fHistV0PhiCorr->Sumw2();
1332    fHistFull[k].fHistFullHar[j].fHistV0PhiCorr->SetXTitle("V0-Plane Correlation (rad)");
1333    fHistFull[k].fHistFullHar[j].fHistV0PhiCorr->SetYTitle("Counts");
1334    delete histTitle;
1335
1336    // neutral sidebands-plane azimuthal correlation
1337    histTitle = new TString("FlowV0sb_Phi_Corr_Sel");
1338    *histTitle += k+1;
1339    histTitle->Append("_Har");
1340    *histTitle += j+1;
1341    fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr = new TH1F(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax / order);
1342    fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr->Sumw2();
1343    fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr->SetXTitle("V0sideBands-Plane Correlation (rad)");
1344    fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr->SetYTitle("Counts");
1345    delete histTitle;
1346
1347    // Yield(pt)
1348    histTitle = new TString("Flow_YieldPt_Sel");
1349    *histTitle += k+1;
1350    histTitle->Append("_Har");
1351    *histTitle += j+1;
1352    fHistFull[k].fHistFullHar[j].fHistYieldPt = new TH1F(histTitle->Data(), histTitle->Data(), fPtBins, fPtMin, fPtMax);
1353    fHistFull[k].fHistFullHar[j].fHistYieldPt->Sumw2();
1354    fHistFull[k].fHistFullHar[j].fHistYieldPt->SetXTitle("Pt (GeV/c)");
1355    fHistFull[k].fHistFullHar[j].fHistYieldPt->SetYTitle("Yield");
1356    delete histTitle;
1357      
1358    // EtaPtPhi 
1359    histTitle = new TString("Flow_EtaPtPhi3D_Sel");
1360    *histTitle += k+1;
1361    histTitle->Append("_Har");
1362    *histTitle += j+1;
1363    fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D = new TH3F(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
1364    fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D->SetXTitle("Eta");
1365    fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D->SetYTitle("Pt (GeV/c)");
1366    fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D->SetZTitle("Phi (rad)");
1367    delete histTitle;
1368
1369    // Yield(eta,pt)
1370    histTitle = new TString("Flow_Yield2D_Sel");
1371    *histTitle += k+1;
1372    histTitle->Append("_Har");
1373    *histTitle += j+1;
1374    fHistFull[k].fHistFullHar[j].fHistYield2D = new TH2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1375    fHistFull[k].fHistFullHar[j].fHistYield2D->Sumw2();
1376    fHistFull[k].fHistFullHar[j].fHistYield2D->SetXTitle("Pseudorapidty");
1377    fHistFull[k].fHistFullHar[j].fHistYield2D->SetYTitle("Pt (GeV/c)");
1378    delete histTitle;
1379
1380    // Dca - 3D
1381    histTitle = new TString("Flow_3dDca_Sel") ;
1382    *histTitle += k+1;
1383    histTitle->Append("_Har");
1384    *histTitle += j+1;
1385    fHistFull[k].fHistFullHar[j].fHistDcaGlob = new TH1F(histTitle->Data(), histTitle->Data(), kDcaBins, dcaMin, glDcaMax);
1386    fHistFull[k].fHistFullHar[j].fHistDcaGlob->Sumw2();
1387    fHistFull[k].fHistFullHar[j].fHistDcaGlob->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
1388    delete histTitle;
1389
1390    // Yield(pt) - excluded from R.P.
1391    histTitle = new TString("Flow_YieldPt_outSel");
1392    *histTitle += k+1;
1393    histTitle->Append("_Har");
1394    *histTitle += j+1;
1395    fHistFull[k].fHistFullHar[j].fHistYieldPtout = new TH1F(histTitle->Data(), histTitle->Data(), fPtBins, fPtMin, fPtMax);
1396    fHistFull[k].fHistFullHar[j].fHistYieldPtout->Sumw2();
1397    fHistFull[k].fHistFullHar[j].fHistYieldPtout->SetXTitle("Pt (GeV/c)");
1398    fHistFull[k].fHistFullHar[j].fHistYieldPtout->SetYTitle("Yield");
1399    delete histTitle;
1400      
1401    // EtaPtPhi - excluded from R.P. 
1402    histTitle = new TString("Flow_EtaPtPhi3D_outSel");
1403    *histTitle += k+1;
1404    histTitle->Append("_Har");
1405    *histTitle += j+1;
1406    fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout = new TH3F(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, kPhi3DBins, fPhiMin, fPhiMax);
1407    fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout->SetXTitle("Eta");
1408    fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout->SetYTitle("Pt (GeV/c)");
1409    fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout->SetZTitle("Phi (rad)");
1410    delete histTitle;
1411
1412    // Yield(eta,pt) - excluded from R.P.
1413    histTitle = new TString("Flow_Yield2D_outSel");
1414    *histTitle += k+1;
1415    histTitle->Append("_Har");
1416    *histTitle += j+1;
1417    fHistFull[k].fHistFullHar[j].fHistYield2Dout = new TH2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1418    fHistFull[k].fHistFullHar[j].fHistYield2Dout->Sumw2();
1419    fHistFull[k].fHistFullHar[j].fHistYield2Dout->SetXTitle("Pseudorapidty");
1420    fHistFull[k].fHistFullHar[j].fHistYield2Dout->SetYTitle("Pt (GeV/c)");
1421    delete histTitle;
1422
1423    // Dca - 3D - excluded from R.P.
1424    histTitle = new TString("Flow_3dDca_outSel") ;
1425    *histTitle += k+1;
1426    histTitle->Append("_Har");
1427    *histTitle += j+1;
1428    fHistFull[k].fHistFullHar[j].fHistDcaGlobout = new TH1F(histTitle->Data(), histTitle->Data(), kDcaBins, dcaMin, glDcaMax);
1429    fHistFull[k].fHistFullHar[j].fHistDcaGlobout->Sumw2();
1430    fHistFull[k].fHistFullHar[j].fHistDcaGlobout->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
1431    delete histTitle;
1432
1433    // Flow observed - v_obs,Pt,Eta
1434    histTitle = new TString("Flow_vObs2D_Sel");
1435    *histTitle += k+1;
1436    histTitle->Append("_Har");
1437    *histTitle += j+1;
1438    fHistFull[k].fHistFullHar[j].fHistvObs2D = new TProfile2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1439    fHistFull[k].fHistFullHar[j].fHistvObs2D->SetXTitle((char*)fLabel.Data());
1440    fHistFull[k].fHistFullHar[j].fHistvObs2D->SetYTitle("Pt (GeV/c)");
1441    delete histTitle;
1442
1443    // v_obs,Eta
1444    histTitle = new TString("Flow_vObsEta_Sel");
1445    *histTitle += k+1;
1446    histTitle->Append("_Har");
1447    *histTitle += j+1;
1448    fHistFull[k].fHistFullHar[j].fHistvObsEta = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1449    fHistFull[k].fHistFullHar[j].fHistvObsEta->SetXTitle((char*)fLabel.Data());
1450    fHistFull[k].fHistFullHar[j].fHistvObsEta->SetYTitle("v (%)");
1451    delete histTitle;
1452
1453    // v_obs,Pt
1454    histTitle = new TString("Flow_vObsPt_Sel");
1455    *histTitle += k+1;
1456    histTitle->Append("_Har");
1457    *histTitle += j+1;
1458    fHistFull[k].fHistFullHar[j].fHistvObsPt = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1459    fHistFull[k].fHistFullHar[j].fHistvObsPt->SetXTitle("Pt (GeV/c)");
1460    fHistFull[k].fHistFullHar[j].fHistvObsPt->SetYTitle("v (%)");
1461    delete histTitle;
1462
1463    // neutral Flow observed - Pt,Eta
1464    histTitle = new TString("FlowV0_vObs2D_Sel");
1465    *histTitle += k+1;
1466    histTitle->Append("_Har");
1467    *histTitle += j+1;
1468    fHistFull[k].fHistFullHar[j].fHistV0vObs2D = new TProfile2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1469    fHistFull[k].fHistFullHar[j].fHistV0vObs2D->SetXTitle((char*)fLabel.Data());
1470    fHistFull[k].fHistFullHar[j].fHistV0vObs2D->SetYTitle("Pt (GeV/c)");
1471    delete histTitle;
1472
1473    // neutral Flow observed - Eta
1474    histTitle = new TString("FlowV0_vObsEta_Sel");
1475    *histTitle += k+1;
1476    histTitle->Append("_Har");
1477    *histTitle += j+1;
1478    fHistFull[k].fHistFullHar[j].fHistV0vObsEta = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1479    fHistFull[k].fHistFullHar[j].fHistV0vObsEta->SetXTitle((char*)fLabel.Data());
1480    fHistFull[k].fHistFullHar[j].fHistV0vObsEta->SetYTitle("v (%)");
1481    delete histTitle;
1482
1483    // neutral Flow observed - Pt
1484    histTitle = new TString("FlowV0_vObsPt_Sel");
1485    *histTitle += k+1;
1486    histTitle->Append("_Har");
1487    *histTitle += j+1;
1488    fHistFull[k].fHistFullHar[j].fHistV0vObsPt = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1489    fHistFull[k].fHistFullHar[j].fHistV0vObsPt->SetXTitle("Pt (GeV/c)");
1490    fHistFull[k].fHistFullHar[j].fHistV0vObsPt->SetYTitle("v (%)");
1491    delete histTitle;
1492
1493    // neutral sidebands Flow observed - Pt,Eta
1494    histTitle = new TString("FlowV0sb_vObs2D_Sel");
1495    *histTitle += k+1;
1496    histTitle->Append("_Har");
1497    *histTitle += j+1;
1498    fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D = new TProfile2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1499    fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D->SetXTitle((char*)fLabel.Data());
1500    fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D->SetYTitle("Pt (GeV/c)");
1501    delete histTitle;
1502
1503    // neutral sidebands Flow observed - Eta
1504    histTitle = new TString("FlowV0sb_vObsEta_Sel");
1505    *histTitle += k+1;
1506    histTitle->Append("_Har");
1507    *histTitle += j+1;
1508    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1509    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->SetXTitle((char*)fLabel.Data());
1510    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->SetYTitle("v (%)");
1511    delete histTitle;
1512
1513    // neutral sidebands Flow observed - Pt
1514    histTitle = new TString("FlowV0sb_vObsPt_Sel");
1515    *histTitle += k+1;
1516    histTitle->Append("_Har");
1517    *histTitle += j+1;
1518    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1519    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->SetXTitle("Pt (GeV/c)");
1520    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->SetYTitle("v (%)");
1521    delete histTitle;
1522
1523    // SX neutral sidebands Flow observed - Eta
1524    histTitle = new TString("FlowV0sb_vObsEta_sx_Sel");
1525    *histTitle += k+1;
1526    histTitle->Append("_Har");
1527    *histTitle += j+1;
1528    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1529    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx->SetXTitle((char*)fLabel.Data());
1530    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx->SetYTitle("v (%)");
1531    delete histTitle;
1532
1533    // SX neutral sidebands Flow observed - Pt
1534    histTitle = new TString("FlowV0sb_vObsPt_sx_Sel");
1535    *histTitle += k+1;
1536    histTitle->Append("_Har");
1537    *histTitle += j+1;
1538    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1539    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx->SetXTitle("Pt (GeV/c)");
1540    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx->SetYTitle("v (%)");
1541    delete histTitle;
1542
1543    // DX neutral sidebands Flow observed - Eta
1544    histTitle = new TString("FlowV0sb_vObsEta_dx_Sel");
1545    *histTitle += k+1;
1546    histTitle->Append("_Har");
1547    *histTitle += j+1;
1548    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1549    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx->SetXTitle((char*)fLabel.Data());
1550    fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx->SetYTitle("v (%)");
1551    delete histTitle;
1552
1553    // DX neutral sidebands Flow observed - Pt
1554    histTitle = new TString("FlowV0sb_vObsPt_dx_Sel");
1555    *histTitle += k+1;
1556    histTitle->Append("_Har");
1557    *histTitle += j+1;
1558    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1559    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx->SetXTitle("Pt (GeV/c)");
1560    fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx->SetYTitle("v (%)");
1561    delete histTitle;
1562
1563    // Phi lab
1564    // Tpc (plus)
1565    histTitle = new TString("Flow_Phi_TPCplus_Sel");
1566    *histTitle += k+1;
1567    histTitle->Append("_Har");
1568    *histTitle += j+1;
1569    fHistFull[k].fHistFullHar[j].fHistPhiPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1570    fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetXTitle("Azimuthal Angles (rad)");
1571    fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetYTitle("Counts");
1572    delete histTitle;
1573
1574    // Tpc (minus)
1575    histTitle = new TString("Flow_Phi_TPCminus_Sel");
1576    *histTitle += k+1;
1577    histTitle->Append("_Har");
1578    *histTitle += j+1;
1579    fHistFull[k].fHistFullHar[j].fHistPhiMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1580    fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetXTitle("Azimuthal Angles (rad)");
1581    fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetYTitle("Counts");
1582    delete histTitle;
1583   
1584    // Tpc (cross)
1585    histTitle = new TString("Flow_Phi_TPCcross_Sel");
1586    *histTitle += k+1;
1587    histTitle->Append("_Har");
1588    *histTitle += j+1;
1589    fHistFull[k].fHistFullHar[j].fHistPhiAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1590    fHistFull[k].fHistFullHar[j].fHistPhiAll->SetXTitle("Azimuthal Angles (rad)");
1591    fHistFull[k].fHistFullHar[j].fHistPhiAll->SetYTitle("Counts");
1592    delete histTitle;
1593   
1594    // Tpc
1595    histTitle = new TString("Flow_Phi_TPC_Sel");
1596    *histTitle += k+1;
1597    histTitle->Append("_Har");
1598    *histTitle += j+1;
1599    fHistFull[k].fHistFullHar[j].fHistPhi = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1600    fHistFull[k].fHistFullHar[j].fHistPhi->SetXTitle("Azimuthal Angles (rad)");
1601    fHistFull[k].fHistFullHar[j].fHistPhi->SetYTitle("Counts");
1602    delete histTitle;
1603   
1604    // Phi lab flattened
1605    // Tpc (Plus)
1606    histTitle = new TString("Flow_Phi_Flat_TPCplus_Sel");
1607    *histTitle += k+1;
1608    histTitle->Append("_Har");
1609    *histTitle += j+1;
1610    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1611    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetXTitle("Azimuthal Angles (rad)");
1612    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetYTitle("Counts");
1613    delete histTitle;
1614
1615    // Tpc (Minus)
1616    histTitle = new TString("Flow_Phi_Flat_TPCminus_Sel");
1617    *histTitle += k+1;
1618    histTitle->Append("_Har");
1619    *histTitle += j+1;
1620    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1621    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetXTitle("Azimuthal Angles (rad)");
1622    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetYTitle("Counts");
1623    delete histTitle;
1624
1625    // Tpc (cross)
1626    histTitle = new TString("Flow_Phi_Flat_TPCcross_Sel");
1627    *histTitle += k+1;
1628    histTitle->Append("_Har");
1629    *histTitle += j+1;
1630    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1631    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetXTitle("Azimuthal Angles (rad)");
1632    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetYTitle("Counts");
1633    delete histTitle;
1634
1635    // Tpc
1636    histTitle = new TString("Flow_Phi_Flat_TPC_Sel");
1637    *histTitle += k+1;
1638    histTitle->Append("_Har");
1639    *histTitle += j+1;
1640    fHistFull[k].fHistFullHar[j].fHistPhiFlat = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1641    fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetXTitle("Azimuthal Angles (rad)");
1642    fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetYTitle("Counts");
1643    delete histTitle;
1644   }
1645  }
1646  cout << "    Init()  -  Histograms booked" << endl ; 
1647
1648  return kTRUE ;
1649 }
1650 //-----------------------------------------------------------------------
1651
1652 //-----------------------------------------------------------------------
1653 Bool_t AliFlowAnalyser::Finish() 
1654 {
1655  // Close the analysis and saves the histograms on the histFile .
1656  
1657  cout << "* FlowAnalysis *  -  Finish()" << endl ; cout << endl ;
1658
1659  // Write all histograms
1660  fHistFile->cd() ; 
1661  fHistFile->Write() ; 
1662  
1663  // Write Resolution corrected histograms
1664  if(fVnResHistList)
1665  { 
1666   fVnResHistList->Write(); 
1667   delete fVnResHistList ;
1668  }
1669  else { cout << "  E.P. resolution has not been calculated. No v_n histograms!" << endl ; } 
1670
1671  // Write PhiWgt histograms
1672  if(fPhiWgtHistList)
1673  { 
1674   fPhiWgtHistList->Write();
1675   delete fPhiWgtHistList ; 
1676  }
1677  
1678  fFlowSelect->Write(); 
1679  // delete fFlowSelect ;
1680
1681  fHistFile->Close() ;
1682
1683  cout << "    Finish()  -  Histograms saved : " << fHistFileName.Data() << endl ; cout << endl ; 
1684
1685  return kTRUE ;
1686 }
1687 //-----------------------------------------------------------------------
1688 // ###
1689 //----------------------------------------------------------------------
1690 Float_t AliFlowAnalyser::GetRunBayesian(Int_t nPid, Int_t selN) const 
1691 {
1692  // Returns the normalized particle abundance of "e","mu","pi","k","p","d"
1693  // in all the analysed events (in selection selN).
1694  // Call at the end of the analysis.
1695
1696  if(selN>AliFlowConstants::kSels) { selN = 0 ; }
1697  Double_t totCount = (fHistFull[selN].fHistBayPidMult)->GetSumOfWeights() ;
1698  if(totCount) { return (fHistFull[selN].fHistBayPidMult->GetBinContent(nPid+1) / totCount) ; }
1699  else         { return 1. ; }
1700 }
1701 //-----------------------------------------------------------------------
1702 void AliFlowAnalyser::PrintRunBayesian(Int_t selN) const
1703 {
1704  // Prints the normalized particle abundance of all the analysed events 
1705  // (in selection selN).
1706
1707  if(selN>AliFlowConstants::kSels) { selN = 0 ; } 
1708  Char_t* names[AliFlowConstants::kPid] = {"e","mu","pi","k","p","d"} ;
1709  Double_t bayes = 0. ;
1710  cout << " selN = " << selN << " particles normalized abundance : " ;
1711  for(int i=0;i<AliFlowConstants::kPid;i++)
1712  {
1713   bayes = GetRunBayesian(i, selN) ;
1714   cout << bayes << "_" << names[i] << " ; " ;
1715  }
1716  cout << endl ;
1717  
1718  return ;
1719 }
1720 //-----------------------------------------------------------------------
1721 void AliFlowAnalyser::FillWgtArrays(TFile* wgtFile)
1722 {
1723  // Loads PhiWeights & Bayesian particles' abundance from file (default: flowPhiWgt.hist.root). 
1724  // Weights are stored in a static TH1D* data member, ready to be plugged into the AliFlowEvent.
1725  // The plugging is done by the method ::FillEvtPhiWgt() (if wgt file is there).
1726
1727  fPhiWgtFile = wgtFile ;
1728
1729  TString* histTitle ;
1730  TH1D* hTPCall ; TH1D* hTPCplus ; TH1D* hTPCminus ; TH1D* hTPCcross ;
1731  TH1D* hPIDbay ;
1732  for(int k=0;k<AliFlowConstants::kSels;k++)
1733  {
1734   for(int j=0;j<AliFlowConstants::kHars;j++) 
1735   {
1736   // Tpc (plus)
1737    histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
1738    *histTitle += k+1;
1739    histTitle->Append("_Har");
1740    *histTitle += j+1;
1741    hTPCplus = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1742    delete histTitle;
1743    // Tpc (minus)
1744    histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
1745    *histTitle += k+1;
1746    histTitle->Append("_Har");
1747    *histTitle += j+1;
1748    hTPCminus = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1749    delete histTitle;
1750    // Tpc (cross)
1751    histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
1752    *histTitle += k+1;
1753    histTitle->Append("_Har");
1754    *histTitle += j+1;
1755    hTPCcross = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1756    delete histTitle;
1757
1758    // Tpc
1759    histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
1760    *histTitle += k+1;
1761    histTitle->Append("_Har");
1762    *histTitle += j+1;
1763    hTPCall = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1764    delete histTitle;
1765
1766    for(int n=0;n<fPhiBins;n++) 
1767    {
1768     fPhiWgtPlus[k][j][n]  = hTPCplus->GetBinContent(n+1) ;
1769     fPhiWgtMinus[k][j][n] = hTPCminus->GetBinContent(n+1) ;
1770     fPhiWgtCross[k][j][n] = hTPCcross->GetBinContent(n+1) ;
1771     fPhiWgt[k][j][n]      = hTPCall->GetBinContent(n+1) ;
1772     // cout << "  Weights: " << fPhiWgt[k][j][n] << " ; " << fPhiWgtPlus[k][j][n] << " | " << fPhiWgtMinus[k][j][n] << " | " << fPhiWgtCross[k][j][n] << endl ; 
1773    }
1774   }
1775   
1776   // Bayesian weights
1777   histTitle = new TString("Flow_BayPidMult_Sel");
1778   *histTitle += k+1;
1779   hPIDbay = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1780   delete histTitle;
1781   Double_t totCount = hPIDbay->GetSumOfWeights() ;
1782   for (int n=0;n<AliFlowConstants::kPid;n++) 
1783   { 
1784    if(totCount) { fBayesianWgt[k][n] = hPIDbay->GetBinContent(n+1) / totCount ; }
1785    else         { fBayesianWgt[k][n] = 1. ; }
1786    // cout << "  Bayesian Weights (" << n << ") : " << fBayesianWgt[k][n] << endl ; 
1787   }
1788  }
1789
1790  delete hTPCall ; delete hTPCplus ; delete hTPCminus ; delete hTPCcross ;
1791  delete hPIDbay ;
1792
1793  return ;
1794 }
1795 //-----------------------------------------------------------------------
1796 void AliFlowAnalyser::FillEvtPhiWgt(AliFlowEvent* fFlowEvent)
1797 {
1798  // Plugs phi weights into the static dwgt data member of the AliFlowEvent class.
1799  // Weights are given in special AliFlowConstants::PhiWgt_t arrays (see AliFlowConstants), 
1800  // which are read from the wgt histograms by the method FillWgtArrays(...).
1801  
1802  fFlowEvent->SetPhiWeight(fPhiWgt);
1803  fFlowEvent->SetPhiWeightPlus(fPhiWgtPlus);
1804  fFlowEvent->SetPhiWeightMinus(fPhiWgtMinus);
1805  fFlowEvent->SetPhiWeightCross(fPhiWgtCross); 
1806  
1807  return ;
1808 }
1809 //-----------------------------------------------------------------------
1810 void AliFlowAnalyser::FillBayesianWgt(AliFlowEvent* fFlowEvent)
1811 {
1812  // Plugs Bayesian particle abundance into the current AliFlowEvent.
1813  // A bayesian vector should be used for the PId of any different selection 
1814  // (different sets of cuts -> different particle abundance), but for now  
1815  // just the Selection n.0 (with no cuts) is used .
1816  // (AliFlowEvent::mBayesianCs[6] is a 1-dimensional array, change that first!).
1817
1818  Double_t bayes[AliFlowConstants::kPid] ; 
1819  Double_t bayCheck = 0. ;
1820  for (int n=0;n<AliFlowConstants::kPid;n++) 
1821  {
1822   bayes[n] = fBayesianWgt[0][n] ;
1823   bayCheck += bayes[n] ;   
1824   // cout << "Bayesian V[" << n << "]  =  " << fBayesianWgt[0][n] << endl ; 
1825  }
1826  if(bayCheck)   { fFlowEvent->SetBayesian(bayes) ; fRePid = kTRUE ; }
1827  else           { cout << "An empty bayesian vector is stored !!! - Bayesian weights = {1,1,1,1,1,1} " << endl ; }
1828
1829  return ;
1830 }
1831 //-----------------------------------------------------------------------
1832 void AliFlowAnalyser::Weightening()
1833 {
1834  // Calculates weights, and fills PhiWgt histograms .
1835  // This is called at the end of the event analysis.
1836  
1837  cout << " AliFlowAnalyser::Weightening() " << endl ; cout << endl ;
1838  
1839  // PhiWgt histogram collection
1840  fPhiWgtHistList = new TOrdCollection(4*AliFlowConstants::kSels*AliFlowConstants::kHars) ;
1841  
1842  // Creates PhiWgt Histograms
1843  TString* histTitle ;
1844  for(int k = 0; k < AliFlowConstants::kSels; k++)
1845  {
1846   for(int j = 0; j < AliFlowConstants::kHars; j++) 
1847   {
1848    // Tpc (plus)
1849    histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
1850    *histTitle += k+1;
1851    histTitle->Append("_Har");
1852    *histTitle += j+1;
1853    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1854    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Sumw2();
1855    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetXTitle("Azimuthal Angles (rad)");
1856    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetYTitle("PhiWgt");
1857    delete histTitle;
1858    // Tpc (minus)
1859    histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
1860    *histTitle += k+1;
1861    histTitle->Append("_Har");
1862    *histTitle += j+1;
1863    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1864    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Sumw2();
1865    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetXTitle("Azimuthal Angles (rad)");
1866    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetYTitle("PhiWgt");
1867    delete histTitle;
1868    // Tpc (cross)
1869    histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
1870    *histTitle += k+1;
1871    histTitle->Append("_Har");
1872    *histTitle += j+1;
1873    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1874    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Sumw2();
1875    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetXTitle("Azimuthal Angles (rad)");
1876    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetYTitle("PhiWgt");
1877    delete histTitle;
1878    // Tpc
1879    histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
1880    *histTitle += k+1;
1881    histTitle->Append("_Har");
1882    *histTitle += j+1;
1883    fHistFull[k].fHistFullHar[j].fHistPhiWgt = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1884    fHistFull[k].fHistFullHar[j].fHistPhiWgt->Sumw2();
1885    fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetXTitle("Azimuthal Angles (rad)");
1886    fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetYTitle("PhiWgt");
1887    delete histTitle;
1888
1889    // Calculate PhiWgt
1890    double meanPlus  = fHistFull[k].fHistFullHar[j].fHistPhiPlus->Integral() / (double)fPhiBins ;
1891    double meanMinus = fHistFull[k].fHistFullHar[j].fHistPhiMinus->Integral() / (double)fPhiBins ;
1892    double meanCross = fHistFull[k].fHistFullHar[j].fHistPhiAll->Integral() / (double)fPhiBins ;
1893    double meanTPC = fHistFull[k].fHistFullHar[j].fHistPhi->Integral() / (double)fPhiBins ;
1894
1895    // Tpc
1896    for (int i=0;i<fPhiBins;i++) 
1897    {
1898     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinContent(i+1,meanPlus);
1899     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinError(i+1, 0.);
1900     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinContent(i+1,meanMinus);
1901     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinError(i+1, 0.);
1902     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinContent(i+1,meanCross);
1903     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinError(i+1, 0.);
1904     fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinContent(i+1,meanTPC);
1905     fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinError(i+1, 0.);
1906    }
1907    
1908    if(meanTPC==0) { cout << " Sel." << k << " , Har." << j << " :  empty phi histogram ! " << endl ; }
1909    else 
1910    {
1911     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiPlus);
1912     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiMinus);
1913     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Divide(fHistFull[k].fHistFullHar[j].fHistPhiAll);
1914     fHistFull[k].fHistFullHar[j].fHistPhiWgt->Divide(fHistFull[k].fHistFullHar[j].fHistPhi);
1915    }
1916    
1917    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus);
1918    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus);
1919    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtAll);
1920    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgt);
1921   }
1922  }
1923
1924  return ;
1925 }
1926 //-----------------------------------------------------------------------
1927 // ###
1928 //-----------------------------------------------------------------------
1929 Bool_t AliFlowAnalyser::Analyse(AliFlowEvent* flowEvent)         
1930 {
1931  // Runs the analysis on the AliFlowEvent (* fFlowEvent). 
1932  // This method can be inserted in a loop over a collection of 
1933  // AliFlowEvents or for on-fly analysis on AliESDs if used toghether 
1934  // with the AliFlowMaker.
1935  
1936  cout << " AliFlowAnalyser::Analyze(" << fFlowEvent << " )   -   " << fEventNumber << endl ;
1937  if(!flowEvent) { return kFALSE ; }
1938  else  { fFlowEvent = flowEvent ; }
1939
1940  if(fFlowSelect->Select(fFlowEvent))     // event selected - here below the ANALYSIS FLAGS are setted -
1941  {
1942   cout << " * 1 . Load event (track & v0s) and set flags . " << endl ;
1943   fFlowTracks = fFlowEvent->TrackCollection() ; 
1944   fNumberOfTracks = fFlowTracks->GetEntries() ;
1945   fFlowV0s = fFlowEvent->V0Collection() ; 
1946   fNumberOfV0s = fFlowV0s->GetEntries() ;
1947   cout << "       event ID = " << fFlowEvent->EventID() << " :  found " << fNumberOfTracks << " AliFlowTracks, and " << fNumberOfV0s << " AliFlowV0s . " << endl ;  
1948                                                             
1949   if(fReadPhiWgt) { FillEvtPhiWgt(fFlowEvent) ; }           // phi and bayesian weights are filled previous to the loop (FillWgtArrays(TFile*))
1950   else            { fFlowEvent->SetNoWgt() ; }              // phi weights can be used or not , this plugs them into the event
1951
1952   if(fBayWgt)     { FillBayesianWgt(fFlowEvent) ; }         // bayesian weights can be used or not , this plugs them into the event
1953   if(fRePid)      { fFlowEvent->SetPids() ; }               // re-calculate all p.id. hypotesis with the (new) bayesian array
1954
1955   if(fOnePhiWgt)  { fFlowEvent->SetOnePhiWgt() ; }          // one phi-wgt histogram
1956   else            { fFlowEvent->SetFirstLastPhiWgt() ; }    // three phi-wgt histogram
1957   if(fPtWgt)      { fFlowEvent->SetPtWgt(); ; }             // pT as a weight
1958   if(fEtaWgt)     { fFlowEvent->SetEtaWgt() ; }             // eta as a weight
1959
1960   if(fShuffle)    { fFlowEvent->RandomShuffle() ; }         // tracks re-shuffling
1961
1962   fFlowEvent->SetSelections(fFlowSelect) ;                  // does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
1963   fFlowEvent->SetEtaSubs(fEtaSub) ;                         // setting for the subevents (eta or random)
1964   fFlowEvent->MakeSubEvents() ;                             // makes the subevent, eta or random basing on the previous flag
1965
1966   cout << " * 2 . Calculating event quantities all in one shoot . " << endl ; 
1967   fFlowEvent->MakeAll() ; 
1968   
1969   if(FillFromFlowEvent(fFlowEvent))                         // calculates event quantities
1970   {
1971    cout << " * 3 . Event Histograms and Particles loop . " << endl ;
1972    FillEventHistograms(fFlowEvent);                         // fill histograms from AliFlowEvents
1973    if(fTrackLoop) { FillParticleHistograms(fFlowTracks) ; } // fill histograms from AliFlowTracks
1974    if(fV0loop)    { FillV0Histograms(fFlowV0s) ; }          // fill histograms from AliFlowV0s 
1975    //FillLabels() ;                                         // fill the histogram of MC labels (from the simulation)
1976   }
1977   else 
1978   {
1979    cout << " * 3 . Event psi = 0   -   Skipping! " << endl ; 
1980    return kFALSE ;
1981   }
1982  }
1983  else 
1984  {
1985   cout << " * 0 . Event " << fEventNumber << " (event ID = " << fFlowEvent->EventID() << ") discarded . " << endl ; 
1986   // delete fFlowEvent  ; fFlowEvent = 0 ; 
1987   return kFALSE ;
1988  }
1989  fEventNumber++ ;
1990  
1991  return kTRUE ;
1992 }
1993 //-----------------------------------------------------------------------
1994 Bool_t AliFlowAnalyser::FillFromFlowEvent(AliFlowEvent* fFlowEvent)
1995 {
1996  // gets event quantities, returns kFALSE if Q vector is always 0 
1997
1998  Int_t selCheck = 0 ;
1999  for(Int_t k = 0; k < AliFlowConstants::kSels; k++) 
2000  {
2001   fFlowSelect->SetSelection(k) ;
2002   for(Int_t j = 0; j < AliFlowConstants::kHars; j++) 
2003   {
2004    fFlowSelect->SetHarmonic(j) ;
2005    for(Int_t n = 0; n < AliFlowConstants::kSubs; n++) 
2006    {
2007     fFlowSelect->SetSubevent(n) ;
2008     fPsiSub[n][k][j] = fFlowEvent->Psi(fFlowSelect) ;   // sub-event quantities
2009     fMultSub[n][k][j] = fFlowEvent->Mult(fFlowSelect) ;
2010    }
2011    fFlowSelect->SetSubevent(-1);
2012    fQ[k][j]    = fFlowEvent->Q(fFlowSelect) ;           // full event quantities
2013    fPsi[k][j]  = fFlowEvent->Psi(fFlowSelect) ;
2014    fQnorm[k][j]   = fFlowEvent->NormQ(fFlowSelect).Mod() ; // was: fFlowEvent->q(fFlowSelect) ; // but the normalization was bad (no pT,eta weight)
2015    fMult[k][j] = fFlowEvent->Mult(fFlowSelect) ;
2016    selCheck += fMult[k][j] ; 
2017   }
2018  }
2019  if(!selCheck) { return kFALSE ; }                      // if there are no particles in the selection -> skip the event
2020
2021  return kTRUE ;
2022 }
2023 //-----------------------------------------------------------------------
2024 void AliFlowAnalyser::FillEventHistograms(AliFlowEvent* fFlowEvent)    
2025 {
2026  // fill event histograms
2027
2028  cout << " Fill Event Histograms ... " << endl ; 
2029
2030  float trigger = (float)fFlowEvent->L0TriggerWord() ;
2031  fHistTrigger->Fill(trigger);
2032
2033  // no selections: OrigMult, Centrality, Mult, MultOverOrig, VertexZ, VertexXY
2034  int origMult = fFlowEvent->OrigMult();
2035  fHistOrigMult->Fill((float)origMult);
2036  fHistMultEta->Fill((float)fFlowEvent->MultEta());
2037
2038  int cent = fFlowEvent->Centrality();
2039  fHistCent->Fill((float)cent);
2040
2041  fHistMult->Fill((float)fNumberOfTracks) ;
2042  fHistV0Mult->Fill((float)fNumberOfV0s) ;
2043  if(origMult) { fHistMultOverOrig->Fill((float)fNumberOfTracks/(float)origMult) ; }
2044
2045  fFlowEvent->VertexPos(fVertex) ;
2046  fHistVertexZ->Fill(fVertex[2]) ;
2047  fHistVertexXY2D->Fill(fVertex[0],fVertex[1]) ;
2048
2049  // ZDC info
2050  fHistPartZDC->Fill(fFlowEvent->ZDCpart()) ;
2051  for(int ii=0;ii<3;ii++) { fHistEnergyZDC->Fill(ii,fFlowEvent->ZDCenergy(ii)) ; }
2052
2053  // sub-event Psi_Subs
2054  for(int k = 0; k < AliFlowConstants::kSels; k++) 
2055  {
2056   for(int j = 0; j < AliFlowConstants::kHars; j++) 
2057   {
2058    for(int n = 0; n < AliFlowConstants::kSubs; n++) 
2059    {
2060     int iii = AliFlowConstants::kSubs * k + n ;    //cout << "  " << k << j << n << " , " << iii << endl ;
2061     fHistSub[iii].fHistSubHar[j].fHistPsiSubs->Fill(fPsiSub[n][k][j]) ;
2062    }
2063   }
2064  }
2065
2066  // full event Psi, PsiSubCorr, PsiSubCorrDiff, cos, mult, q
2067  for(int k = 0; k < AliFlowConstants::kSels; k++) 
2068  {
2069   for(int j = 0; j < AliFlowConstants::kHars; j++) 
2070   {
2071    float order = (float)(j+1);
2072    fHistFull[k].fHistFullHar[j].fHistPsi->Fill(fPsi[k][j]);
2073    if(k<2 && j<2)
2074    {
2075     if(k==0) 
2076     { 
2077      float psi1 = fPsi[0][j] ; 
2078      float psi2 = fPsi[1][j] ;
2079      float diff = psi1 - psi2 ;
2080      if(diff < -TMath::Pi()/(j+1))      { diff += 2*TMath::Pi()/(j+1) ; } 
2081      else if(diff > +TMath::Pi()/(j+1)) { diff -= 2*TMath::Pi()/(j+1) ; }
2082      fHistFull[k].fHistFullHar[j].fHistPsiDiff->Fill(diff) ; // k=0
2083     } 
2084     else if(k==1) 
2085     {
2086      float psi1 = 0. ; float psi2 = 0. ;
2087      if (j==0)     { psi1 = fPsi[0][0] ; psi2 = fPsi[1][1] ; }  
2088      else if(j==1) { psi1 = fPsi[0][0] ; psi2 = fPsi[0][1] ; }
2089      float diff = psi1 - psi2 ;
2090      diff = (TMath::Abs(diff) > TMath::Pi()) ? ((diff > 0.) ? -(2*TMath::Pi()-diff) : -(diff+2*TMath::Pi())) : diff ;
2091      fHistFull[k].fHistFullHar[j].fHistPsiDiff->Fill(diff) ; // k=1
2092     }      
2093    }
2094
2095    if(fPsiSub[0][k][j] != 0. && fPsiSub[1][k][j] != 0.)
2096    {
2097     float psiSubCorr;                           // this is:  delta_Psi
2098     if(fV1Ep1Ep2 == kFALSE || order != 1) 
2099     {
2100      psiSubCorr = fPsiSub[0][k][j] - fPsiSub[1][k][j];
2101     }
2102     else // i.e. (fV1Ep1Ep2 == kTRUE && order == 1)
2103     { 
2104      psiSubCorr = fPsiSub[0][k][0] + fPsiSub[1][k][0] - 2*fPsi[k][1];
2105     }
2106     fHistFull[k].fHistCos->Fill(order, (float)cos(order * psiSubCorr)) ;
2107     if(psiSubCorr < 0.)                  { psiSubCorr += 2*TMath::Pi()/order ; }
2108     if(psiSubCorr > 2*TMath::Pi()/order) { psiSubCorr -= 2*TMath::Pi()/order ; } // for v1Ep1Ep2 which gives -2*TMath::Pi() < psiSubCorr < 2*2*TMath::Pi()
2109     fHistFull[k].fHistFullHar[j].fHistPsiSubCorr->Fill(psiSubCorr);
2110    }
2111
2112    if(j < AliFlowConstants::kHars - 1) // subevents of different harmonics
2113    {
2114     int j1 = 0 ; int j2 = 0 ;
2115     float psiSubCorrDiff;
2116     if(j==0)      { j1 = 1, j2 = 2 ; } 
2117     else if(j==1) { j1 = 1, j2 = 3 ; } 
2118     else if(j==2) { j1 = 2, j2 = 4 ; }
2119     psiSubCorrDiff = fmod((double)fPsiSub[0][k][j1-1],2*TMath::Pi()/(double)j2)-fmod((double)fPsiSub[1][k][j2-1],2*TMath::Pi()/(double)j2) ;
2120     if(psiSubCorrDiff < 0.) { psiSubCorrDiff += 2*TMath::Pi()/(float)j2 ; }
2121     fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->Fill(psiSubCorrDiff) ;
2122     psiSubCorrDiff = fmod((double)fPsiSub[0][k][j2-1],2*TMath::Pi()/(double)j2)-fmod((double)fPsiSub[1][k][j1-1],2*TMath::Pi()/(double)j2) ;
2123     if(psiSubCorrDiff < 0.) { psiSubCorrDiff += 2*TMath::Pi()/(float)j2 ; }
2124     fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->Fill(psiSubCorrDiff) ;
2125    }
2126    
2127    fHistFull[k].fHistFullHar[j].fHistMult->Fill((float)fMult[k][j]) ;
2128    fHistFull[k].fHistFullHar[j].fHistQnorm->Fill(fQnorm[k][j]) ;
2129   }
2130  }
2131  
2132  return ;
2133 }
2134 //-----------------------------------------------------------------------
2135 void AliFlowAnalyser::FillParticleHistograms(TClonesArray* fFlowTracks) 
2136 {
2137  // fills tracks histograms
2138
2139  cout << " Tracks Loop . " << endl ; 
2140
2141  float corrMultUnit   = 0. ;
2142  float corrMultN      = 0. ;
2143  float etaSymPosTpcN  = 0. ;
2144  float etaSymNegTpcN  = 0. ;
2145  float etaSymPosTpcNpart = 0. ;
2146  float etaSymNegTpcNpart = 0. ;
2147  float hPlusN         = 0. ;
2148  float hMinusN        = 0. ;
2149  float piPlusN        = 0. ;
2150  float piMinusN       = 0. ;
2151  float protonN        = 0. ;
2152  float pbarN          = 0. ;
2153  float kMinusN        = 0. ;
2154  float kPlusN         = 0. ;
2155  float deuteronN      = 0. ;
2156  float dbarN          = 0. ;
2157  float electronN      = 0. ;
2158  float positronN      = 0. ;
2159  float muonMinusN     = 0. ;
2160  float muonPlusN      = 0. ;
2161
2162  fSelParts = 0 ;
2163  fSelV0s   = 0 ;
2164  for(fTrackNumber=0;fTrackNumber<fNumberOfTracks;fTrackNumber++) 
2165  {
2166   //fFlowTrack = (AliFlowTrack*)fFlowTracks->At(fTrackNumber) ;
2167   
2168   fFlowTrack = (AliFlowTrack*)(fFlowTracks->AddrAt(fTrackNumber)) ;
2169   //cout << "Track n. " << fTrackNumber << endl ; // fFlowTrack->Dump() ; 
2170   
2171   bool constrainable = fFlowTrack->IsConstrainable() ;
2172   // int label = fFlowTrack->Label() ;                  
2173   float dcaGlobal = TMath::Abs(fFlowTrack->Dca()) ;     
2174   float dcaSigned = fFlowTrack->TransDcaSigned() ;      
2175   float dcaTrans = TMath::Abs(dcaSigned) ;              
2176   float eta = fFlowTrack->Eta() ;                       
2177   float phi = fFlowTrack->Phi() ;                       
2178   float pt = fFlowTrack->Pt() ;                         
2179   float etaGlob = fFlowTrack->EtaGlobal() ;                                                                     
2180   float phiGlob = fFlowTrack->PhiGlobal() ;             
2181   float ptGlob = fFlowTrack->PtGlobal() ;               
2182   float totalp = fFlowTrack->P() ;                      
2183   // float logp = TMath::Log10(totalp) ;
2184   // float zFirstPoint = fFlowTrack->ZFirstPoint() ;    
2185   // float zLastPoint = fFlowTrack->ZLastPoint() ;      
2186   float lenght = fFlowTrack->TrackLength() ;            
2187   int   charge = fFlowTrack->Charge() ;                 
2188   float chi2 = fFlowTrack->Chi2() ;                     
2189   int   fitPtsTPC = (int)((float)fFlowTrack->FitPtsTPC()) ;
2190   int   maxPtsTPC = fFlowTrack->MaxPtsTPC() ;           
2191   float chi2TPC = fFlowTrack->Chi2TPC() ;               
2192   int   fitPtsITS = fFlowTrack->FitPtsITS() ;           
2193   int   maxPtsITS = fFlowTrack->MaxPtsITS() ;           
2194   float chi2ITS = fFlowTrack->Chi2ITS() ;               
2195   int   fitPtsTRD = fFlowTrack->NhitsTRD() ;            
2196   int   maxPtsTRD = fFlowTrack->MaxPtsTRD() ;           
2197   float chi2TRD = fFlowTrack->Chi2TRD() ;               
2198   int   fitPtsTOF = fFlowTrack->NhitsTOF() ;            
2199   int   maxPtsTOF = fFlowTrack->MaxPtsTOF() ;           
2200   float chi2TOF = fFlowTrack->Chi2TOF() ;               
2201   float dEdx = fFlowTrack->DedxTPC() ;                  
2202   float its = fFlowTrack->DedxITS() ;                   
2203   float trd = fFlowTrack->SigTRD() ;                    
2204   float tof = fFlowTrack->TofTOF() ;
2205   float lpTPC = 0 ; if(fFlowTrack->PatTPC()>0) { lpTPC = TMath::Log10(fFlowTrack->PatTPC()) ; }
2206   float lpITS = 0 ; if(fFlowTrack->PatITS()>0) { lpITS = TMath::Log10(fFlowTrack->PatITS()) ; }
2207   float lpTRD = 0 ; if(fFlowTrack->PatTRD()>0) { lpTRD = TMath::Log10(fFlowTrack->PatTRD()) ; }
2208   float lpTOF = 0 ; if(fFlowTrack->PatTOF()>0) { lpTOF = TMath::Log10(fFlowTrack->PatTOF()) ; }  
2209   float invMass = fFlowTrack->InvMass() ;              
2210   Char_t pid[10]="0" ; strcpy(pid,fFlowTrack->Pid()) ; 
2211   fPidId        = -1 ; // assigned later
2212
2213   // no selections: Charge, Dca, Chi2, FitPts, MaxPts, FitOverMax, PID
2214   fHistCharge->Fill((float)charge);
2215   fHistDcaGlobal->Fill(dcaGlobal);
2216   fHistDca->Fill(dcaTrans) ;
2217   fHistTransDca->Fill(dcaSigned);
2218   fHistChi2->Fill(chi2);
2219   
2220   // - here TPC   (chi2 & nHits)
2221   fHistChi2TPC->Fill(chi2TPC);
2222   if(fitPtsTPC>0) { fHistChi2normTPC->Fill(chi2TPC/((float)fitPtsTPC)) ; }
2223   fHistFitPtsTPC->Fill((float)fitPtsTPC);
2224   fHistMaxPtsTPC->Fill((float)maxPtsTPC);
2225   if(maxPtsTPC>0) { fHistFitOverMaxTPC->Fill((float)(fitPtsTPC)/(float)maxPtsTPC) ; }
2226
2227   // - here ITS   (chi2 & nHits)
2228   fHistChi2ITS->Fill(chi2ITS);
2229   if(fitPtsITS>0) { fHistChi2normITS->Fill(chi2ITS/((float)fitPtsITS)) ; }
2230   fHistFitPtsITS->Fill((float)fitPtsITS);
2231   fHistMaxPtsITS->Fill((float)maxPtsITS);
2232
2233   // - here TRD   (chi2 & nHits)
2234   fHistChi2TRD->Fill(chi2TRD);
2235   if(fitPtsTRD>0) { fHistChi2normTRD->Fill(chi2TRD/((float)fitPtsTRD)) ; }
2236   fHistFitPtsTRD->Fill((float)fitPtsTRD);
2237   fHistMaxPtsTRD->Fill((float)maxPtsTRD);
2238
2239   // - here TOF   (chi2 & nHits)
2240   fHistChi2TOF->Fill(chi2TOF);
2241   if(fitPtsTOF>0) { fHistChi2normTOF->Fill(chi2TOF/((float)fitPtsTOF)) ; }
2242   fHistFitPtsTOF->Fill((float)fitPtsTOF);
2243   fHistMaxPtsTOF->Fill((float)maxPtsTOF);
2244   
2245   // fit over max (all)
2246   int maxPts = maxPtsITS + maxPtsTPC + maxPtsTRD + maxPtsTOF ;
2247   int fitPts = fitPtsITS + fitPtsTPC + fitPtsTRD + fitPtsTOF ;
2248   if(maxPts>0)    { fHistFitOverMax->Fill((float)(fitPts)/(float)maxPts) ; }
2249   
2250   // lenght
2251   fHistLenght->Fill(lenght) ;
2252   fHistInvMass->Fill(invMass) ;
2253
2254   // PID histograms & multiplicity count (for bayesian histogram)
2255   if(charge == 1) 
2256   {
2257    hPlusN++ ;
2258    fHistMeanDedxPos2D->Fill(lpTPC, dEdx) ;
2259    fHistMeanDedxPos2DITS->Fill(lpITS, its) ;
2260    fHistMeanDedxPos2DTRD->Fill(lpTRD, trd) ;
2261    fHistMeanDedxPos2DTOF->Fill(lpTOF, tof) ;
2262    //
2263    float positron  = fFlowTrack->ElectronPositronProb() ;
2264    fHistPidPositron->Fill(positron) ;
2265    if(strcmp(pid, "e+") == 0) 
2266    {
2267     fPidId = 0 ; positronN++ ; fHistPidPt->Fill(1.5,pt) ;
2268     fHistMeanTPCPositron->Fill(lpTPC, dEdx) ;
2269     fHistMeanITSPositron->Fill(lpITS, its);
2270     fHistMeanTRDPositron->Fill(lpTRD, trd);
2271     fHistMeanTOFPositron->Fill(invMass, tof);
2272     fHistPidPositronPart->Fill(positron) ;
2273    }
2274    float muonPlus  = fFlowTrack->MuonPlusMinusProb() ;
2275    fHistPidMuonPlus->Fill(muonPlus) ;
2276    if(strcmp(pid, "mu+") == 0) 
2277    {
2278     fPidId = 1 ; muonPlusN++ ; fHistPidPt->Fill(3.5,pt) ; 
2279     fHistMeanTPCMuonPlus->Fill(lpTPC, dEdx) ;
2280     fHistMeanITSMuonPlus->Fill(lpITS, its);
2281     fHistMeanTRDMuonPlus->Fill(lpTRD, trd);
2282     fHistMeanTOFMuonPlus->Fill(invMass, tof);
2283     fHistPidMuonPlusPart->Fill(muonPlus) ;
2284    }
2285    float piPlus = fFlowTrack->PionPlusMinusProb() ;
2286    fHistPidPiPlus->Fill(piPlus) ;
2287    if(strcmp(pid, "pi+") == 0) 
2288    { 
2289     fPidId = 2 ; piPlusN++ ; fHistPidPt->Fill(5.5,pt) ;
2290     fHistMeanTPCPiPlus->Fill(lpTPC, dEdx) ;
2291     fHistMeanITSPiPlus->Fill(lpITS, its);
2292     fHistMeanTRDPiPlus->Fill(lpTRD, trd);
2293     fHistMeanTOFPiPlus->Fill(invMass, tof);
2294     fHistPidPiPlusPart->Fill(piPlus) ;
2295    }
2296    float kplus  = fFlowTrack->KaonPlusMinusProb() ;
2297    fHistPidKplus->Fill(kplus) ;
2298    if(strcmp(pid, "k+") == 0) 
2299    {
2300     fPidId = 3 ; kPlusN++ ; fHistPidPt->Fill(7.5,pt) ; 
2301     fHistMeanTPCKplus->Fill(lpTPC, dEdx) ;
2302     fHistMeanITSKplus->Fill(lpITS, its);
2303     fHistMeanTRDKplus->Fill(lpTRD, trd);
2304     fHistMeanTOFKplus->Fill(invMass, tof);
2305     fHistPidKplusPart->Fill(kplus) ;
2306    }
2307    float proton  = fFlowTrack->ProtonPbarProb() ;
2308    fHistPidProton->Fill(proton) ;
2309    if(strcmp(pid, "pr+") == 0) 
2310    {
2311     fPidId = 4 ; protonN++ ; fHistPidPt->Fill(9.5,pt) ;  
2312     fHistMeanTPCProton->Fill(lpTPC, dEdx) ;
2313     fHistMeanITSProton->Fill(lpITS, its);
2314     fHistMeanTRDProton->Fill(lpTRD, trd);
2315     fHistMeanTOFProton->Fill(invMass, tof);
2316     fHistPidProtonPart->Fill(proton) ;
2317    }
2318    float deuteron  = fFlowTrack->DeuteriumAntiDeuteriumProb() ;
2319    fHistPidDeuteron->Fill(deuteron) ;
2320    if(strcmp(pid, "d+") == 0) 
2321    {
2322     fPidId = 6 ; deuteronN++ ; fHistPidPt->Fill(11.5,pt) ; 
2323     fHistMeanTPCDeuteron->Fill(lpTPC, dEdx) ;
2324     fHistMeanITSDeuteron->Fill(lpITS, its);
2325     fHistMeanTRDDeuteron->Fill(lpTRD, trd);
2326     fHistMeanTOFDeuteron->Fill(invMass, tof);
2327     fHistPidDeuteronPart->Fill(deuteron) ;
2328    }
2329   }
2330   else if(charge == -1) 
2331   {
2332    hMinusN++ ;
2333    fHistMeanDedxNeg2D->Fill(lpTPC, dEdx) ;
2334    fHistMeanDedxNeg2DITS->Fill(lpITS, its) ;
2335    fHistMeanDedxNeg2DTRD->Fill(lpTRD, trd) ;
2336    fHistMeanDedxNeg2DTOF->Fill(lpTOF, tof) ;
2337    //
2338    float electron  = fFlowTrack->ElectronPositronProb() ;
2339    fHistPidElectron->Fill(electron);
2340    if(strcmp(pid, "e-") == 0) 
2341    {
2342     fPidId = 0 ; electronN++ ; fHistPidPt->Fill(0.5,pt) ; 
2343     fHistMeanTPCElectron->Fill(lpTPC, dEdx);
2344     fHistMeanITSElectron->Fill(lpITS, its);
2345     fHistMeanTRDElectron->Fill(lpTRD, trd);
2346     fHistMeanTOFElectron->Fill(invMass, tof);
2347     fHistPidElectronPart->Fill(electron);
2348    }
2349    float muonMinus  = fFlowTrack->MuonPlusMinusProb() ;
2350    fHistPidMuonMinus->Fill(muonMinus) ;
2351    if(strcmp(pid, "mu-") == 0) 
2352    {
2353     fPidId = 1 ; muonMinusN++ ; fHistPidPt->Fill(2.5,pt) ;
2354     fHistMeanTPCMuonMinus->Fill(lpTPC, dEdx) ;
2355     fHistMeanITSMuonMinus->Fill(lpITS, its);
2356     fHistMeanTRDMuonMinus->Fill(lpTRD, trd);
2357     fHistMeanTOFMuonMinus->Fill(invMass, tof);
2358     fHistPidMuonMinusPart->Fill(muonMinus) ;
2359    }
2360    float piMinus = fFlowTrack->PionPlusMinusProb() ;
2361    fHistPidPiMinus->Fill(piMinus) ;
2362    if(strcmp(pid, "pi-") == 0) 
2363    {
2364     fPidId = 2 ; piMinusN++ ; fHistPidPt->Fill(4.5,pt) ;
2365     fHistMeanTPCPiMinus->Fill(lpTPC, dEdx);
2366     fHistMeanITSPiMinus->Fill(lpITS, its);
2367     fHistMeanTRDPiMinus->Fill(lpTRD, trd);
2368     fHistMeanTOFPiMinus->Fill(invMass, tof);
2369     fHistPidPiMinusPart->Fill(piMinus);
2370    }
2371    float kminus  = fFlowTrack->KaonPlusMinusProb() ;
2372    fHistPidKminus->Fill(kminus);
2373    if(strcmp(pid, "k-") == 0) 
2374    {
2375     fPidId = 3 ; kMinusN++ ; fHistPidPt->Fill(6.5,pt) ;
2376     fHistMeanTPCKminus->Fill(lpTPC, dEdx);
2377     fHistMeanITSKminus->Fill(lpITS, its);
2378     fHistMeanTRDKminus->Fill(lpTRD, trd);
2379     fHistMeanTOFKminus->Fill(invMass, tof);
2380     fHistPidKminusPart->Fill(kminus);
2381    }
2382    float antiproton  = fFlowTrack->ProtonPbarProb() ;
2383    fHistPidAntiProton->Fill(antiproton);
2384    if(strcmp(pid, "pr-") == 0) 
2385    {
2386     fPidId = 4 ; pbarN++ ; fHistPidPt->Fill(8.5,pt) ;
2387     fHistMeanTPCPbar->Fill(lpTPC, dEdx);
2388     fHistMeanITSPbar->Fill(lpITS, its);
2389     fHistMeanTRDPbar->Fill(lpTRD, trd);
2390     fHistMeanTOFPbar->Fill(invMass, tof);
2391     fHistPidAntiProtonPart->Fill(antiproton);
2392    }
2393    float antideuteron  = fFlowTrack->DeuteriumAntiDeuteriumProb() ;
2394    fHistPidAntiDeuteron->Fill(antideuteron);
2395    if(strcmp(pid, "d-") == 0) 
2396    {
2397     fPidId = 6 ; dbarN++ ; fHistPidPt->Fill(10.5,pt) ;
2398     fHistMeanTPCAntiDeuteron->Fill(lpTPC, dEdx);
2399     fHistMeanITSAntiDeuteron->Fill(lpITS, its);
2400     fHistMeanTRDAntiDeuteron->Fill(lpTRD, trd);
2401     fHistMeanTOFAntiDeuteron->Fill(invMass, tof);
2402     fHistPidAntiDeuteronPart->Fill(antideuteron);
2403    }
2404   }
2405   
2406   // Yield3D, Yield2D, Eta, Pt, Phi, bayP.Id.
2407   fHistPtot->Fill(totalp) ;
2408   fHistPt->Fill(pt) ;
2409   fHistPhi->Fill(phi);  
2410   fHistAllEtaPtPhi3D->Fill(eta, pt, phi) ;
2411   fHistYieldAll2D->Fill(eta, pt) ;
2412   fHistBayPidMult->Fill(fPidId) ;
2413   if(constrainable) 
2414   { 
2415    fHistPhiCons->Fill(phi);  
2416    fHistPhiPtCon->Fill(phi, pt);  
2417    fHistYieldCon2D->Fill(eta, pt) ;
2418    fHistConsEtaPtPhi3D->Fill(eta, pt, phi) ;
2419    fHistGlobEtaPtPhi3D->Fill(etaGlob, ptGlob, phiGlob) ;
2420   }
2421   else 
2422   { 
2423    fHistYieldUnc2D->Fill(etaGlob, ptGlob) ;
2424    fHistUncEtaPtPhi3D->Fill(etaGlob, ptGlob, phiGlob) ;
2425    fHistPhiPtUnc->Fill(phiGlob, ptGlob) ; 
2426   }
2427   if(fFlowTrack->Charge()>0)       { fHistPtPhiPos->Fill(phi, pt); }
2428   else if(fFlowTrack->Charge()<0)  { fHistPtPhiNeg->Fill(phi, pt); }
2429
2430   // fills selected part histograms
2431   if(fFlowSelect->SelectPart(fFlowTrack)) 
2432   {
2433    fSelParts++ ;
2434   
2435    if(strlen(fFlowSelect->PidPart()) != 0) 
2436    {
2437     float rapidity = fFlowTrack->Y();
2438     fHistBinEta->Fill(rapidity, rapidity);
2439     fHistYieldPart2D->Fill(rapidity, pt);
2440    }
2441    else 
2442    {
2443     fHistBinEta->Fill(eta, eta) ;
2444     fHistYieldPart2D->Fill(eta, pt) ;
2445    }
2446    fHistBayPidMultPart->Fill(fPidId) ;
2447    fHistBinPt->Fill(pt, pt) ;
2448    fHistEtaPtPhi3DPart->Fill(eta,pt,phi) ;
2449    fHistDcaGlobalPart->Fill(dcaGlobal) ;
2450    fHistInvMassPart->Fill(invMass) ;
2451    if(charge == 1) 
2452    {
2453     fHistMeanDedxPos3DPart->Fill(lpITS, dEdx, fPidId) ;
2454     fHistMeanDedxPos3DPartITS->Fill(lpITS, its, fPidId) ;
2455    }
2456    else if(charge == -1) 
2457    {
2458     fHistMeanDedxNeg3DPart->Fill(lpITS, dEdx, fPidId) ;
2459     fHistMeanDedxNeg3DPartITS->Fill(lpITS, its, fPidId) ;
2460    }
2461    if(eta > 0.) { etaSymPosTpcNpart++ ; }  // for fHistEtaSymPart 
2462    else         { etaSymNegTpcNpart++ ; }
2463   }
2464   else         
2465   {
2466    fHistEtaPtPhi3DOut->Fill(eta,pt,phi) ;
2467    fHistYieldOut2D->Fill(eta, pt) ;
2468    fHistDcaGlobalOut->Fill(dcaGlobal) ;
2469    fHistInvMassOut->Fill(invMass) ;
2470   }
2471
2472   // cos(n*phiLab)
2473   for(int j = 0; j < AliFlowConstants::kHars; j++) 
2474   {
2475    bool oddHar = (j+1) % 2 ;
2476    float order = (float)(j+1) ;
2477    float vIn   = 100 * cos((double)order * phi) ;
2478    if(eta < 0 && oddHar) { vIn *= -1 ; }
2479    fHistCosPhi->Fill(order, vIn);
2480   }
2481
2482   //For Eta symmetry TPC
2483   if(fFlowTrack->FitPtsTPC())
2484   {
2485    if(eta > 0.) { etaSymPosTpcN++ ; } // for fHistEtaSym 
2486    else         { etaSymNegTpcN++ ; }
2487   }
2488   
2489   // not to call it twice ...
2490   int nEtaS = HarmonicsLoop(fFlowTrack) ;
2491
2492   //For Correlated Multiplicity 
2493   corrMultN += (float)nEtaS ;
2494
2495   //For Correlated Multiplicity in 1 unit rapidity
2496   if(TMath::Abs(eta) <= 0.5) { corrMultUnit += (float)nEtaS ; }
2497
2498   //delete pointer
2499   fFlowTrack = 0 ;
2500  }                                                                                               // end of tracks loop
2501                                                      
2502  // EtaSym
2503  float etaSymTpc = 0 ; 
2504  if(etaSymPosTpcN || etaSymNegTpcN)         { etaSymTpc  = (etaSymPosTpcN  - etaSymNegTpcN)  / (etaSymPosTpcN  + etaSymNegTpcN); }
2505  float etaSymTpcPart = 0 ; 
2506  if(etaSymPosTpcNpart || etaSymNegTpcNpart) { etaSymTpcPart  = (etaSymPosTpcNpart  - etaSymNegTpcNpart)  / (etaSymPosTpcNpart  + etaSymNegTpcNpart) ; }
2507  Float_t vertexZ = fVertex[2] ;
2508  // all
2509  fHistEtaSym->Fill(etaSymTpc);
2510  fHistEtaSymVerZ2D->Fill(vertexZ , etaSymTpc);
2511  // selected
2512  fHistEtaSymPart->Fill(etaSymTpc);
2513  fHistEtaSymVerZ2DPart->Fill(vertexZ , etaSymTpcPart);
2514
2515  // PID multiplicities
2516  float totalMult = (float)fFlowTracks->GetEntries() ;
2517  fHistPidMult->Fill(1., totalMult);
2518  fHistPidMult->Fill(2., hPlusN);
2519  fHistPidMult->Fill(3., hMinusN);
2520  fHistPidMult->Fill(4., piPlusN);
2521  fHistPidMult->Fill(5., piMinusN);
2522  fHistPidMult->Fill(6., protonN);
2523  fHistPidMult->Fill(7., pbarN);
2524  fHistPidMult->Fill(8., kPlusN);
2525  fHistPidMult->Fill(9., kMinusN);
2526  fHistPidMult->Fill(10., deuteronN);
2527  fHistPidMult->Fill(11., dbarN);
2528  fHistPidMult->Fill(12., electronN);
2529  fHistPidMult->Fill(13., positronN);
2530  fHistPidMult->Fill(14., muonMinusN);
2531  fHistPidMult->Fill(15., muonPlusN);
2532
2533  // Multiplicity of particles correlated with the event planes
2534  corrMultN /= (float)(AliFlowConstants::kHars * AliFlowConstants::kSels) ; 
2535  if(fSelParts == corrMultN) { fHistMultPart->Fill(corrMultN) ; } // crosscheck
2536  else                       { fHistMultPart->Fill(-1) ; }
2537  // ...in one unit rapidity
2538  corrMultUnit /= (float)(AliFlowConstants::kHars * AliFlowConstants::kSels) ; 
2539  fHistMultPartUnit->Fill(corrMultUnit) ;
2540
2541  return ;
2542 }
2543 //-----------------------------------------------------------------------
2544 Int_t AliFlowAnalyser::HarmonicsLoop(AliFlowTrack* fFlowTrack)
2545 {
2546  // HarmonicsLoop, does the correlation analysis, fills histograms
2547  // (internally called by ::FillParticleHistograms(..) loop) .
2548  
2549  float eta = fFlowTrack->Eta() ;                        
2550  float phi = fFlowTrack->Phi() ;                        
2551  float pt = fFlowTrack->Pt() ;                  
2552  float zFirstPoint = fFlowTrack->ZFirstPoint() ; 
2553  // float zLastPoint = fFlowTrack->ZLastPoint() ;
2554
2555  // Looping over Selections and Harmonics
2556  int corrMultN = 0 ;
2557  for (int k = 0; k < AliFlowConstants::kSels; k++) 
2558  {
2559   fFlowSelect->SetSelection(k) ;
2560   for (int j = 0; j < AliFlowConstants::kHars; j++) 
2561   {
2562    bool oddHar = (j+1) % 2;
2563    fFlowSelect->SetHarmonic(j);
2564    double order  = (double)(j+1);
2565    float psii = 0. ; float psi2 = 0. ;
2566    if(fFlowEvent->EtaSubs())           // particles with the opposite subevent
2567    {
2568     if(eta > 0) { psii = fPsiSub[1][k][j] ; }     //check
2569     else        { psii = fPsiSub[0][k][j] ; }
2570    } 
2571    else if(order > 3. && !oddHar) 
2572    {
2573     psii = fPsi[k][1];  // 2nd harmomic event plane
2574     if(psii > 2*TMath::Pi()/order) { psii -= 2*TMath::Pi()/order ; } 
2575     if(psii > 2*TMath::Pi()/order) { psii -= 2*TMath::Pi()/order ; }
2576    } 
2577    else  // random subevents
2578    {
2579     psii = fPsi[k][j] ;
2580    }
2581
2582    if(fFlowSelect->Select(fFlowTrack)) // Get detID
2583    {
2584     Bool_t kTpcPlus  = kFALSE ;
2585     Bool_t kTpcMinus = kFALSE ;
2586     Bool_t kTpcAll   = kFALSE ;
2587
2588     fHistFull[k].fHistFullHar[j].fHistYieldPt->Fill(pt);
2589     fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D->Fill(eta, pt, phi);
2590     fHistFull[k].fHistFullHar[j].fHistYield2D->Fill(eta, pt);
2591     fHistFull[k].fHistFullHar[j].fHistDcaGlob->Fill(TMath::Abs(fFlowTrack->Dca()));
2592
2593     // Set Tpc (+ and -)
2594     if(fFlowTrack->FitPtsTPC())    //OR*: AliESDtrack:: "TBits& GetTPCClusterMap()" or "Int_t GetTPCclusters(Int_t* idx)" ...
2595     {
2596      if(zFirstPoint >= 0. && eta > 0.)      { kTpcPlus  = kTRUE ; } 
2597      else if(zFirstPoint <= 0. && eta < 0.) { kTpcMinus = kTRUE ; }
2598      else                                   { kTpcAll   = kTRUE ; }       
2599     }
2600
2601     // PID Multiplicities (particle for R.P.) - done just one time for each selection
2602     if(j==0) { fHistFull[k].fHistBayPidMult->Fill(fPidId) ; }
2603
2604     // Calculate weights for filling histograms   
2605     float wt = 1. ; // TMath::Abs(fFlowEvent->Weight(k, j, fFlowTrack)) ; 
2606
2607     // Fill histograms with selections
2608     if(kTpcPlus)       { fHistFull[k].fHistFullHar[j].fHistPhiPlus->Fill(phi,wt) ;  } 
2609     else if(kTpcMinus) { fHistFull[k].fHistFullHar[j].fHistPhiMinus->Fill(phi,wt) ; } 
2610     else if(kTpcAll)   { fHistFull[k].fHistFullHar[j].fHistPhiAll->Fill(phi,wt) ;   } 
2611     fHistFull[k].fHistFullHar[j].fHistPhi->Fill(phi,wt) ;
2612
2613     // Get phiWgt from file
2614     double phiWgt;
2615     if(order > 3. && !oddHar) { phiWgt = fFlowEvent->PhiWeight(k, 1, fFlowTrack) ; } 
2616     else                      { phiWgt = fFlowEvent->PhiWeight(k, j, fFlowTrack) ; }
2617     if(oddHar && eta<0.)      { phiWgt /= -1. ; } // only for flat hists
2618     // cout << " Weight [" << k << "][" << j << "] (" << fFlowTrack->GetName() << ") = " << phiWgt << " or " << wt << " for phi = " << phi << endl ; // chk
2619
2620     // Fill Flat histograms
2621     if(kTpcPlus)       { fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->Fill(phi, phiWgt) ; } 
2622     else if(kTpcMinus) { fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->Fill(phi, phiWgt) ; } 
2623     else if(kTpcAll)   { fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->Fill(phi, phiWgt) ; } 
2624     fHistFull[k].fHistFullHar[j].fHistPhiFlat->Fill(phi,phiWgt) ;
2625     
2626     if(oddHar && eta<0.)  { phiWgt *= -1. ; } // restore value
2627
2628     // Remove autocorrelations
2629     TVector2 qi;
2630     if(!fFlowEvent->EtaSubs())   // random subevents
2631     {
2632      if(order > 3. && !oddHar) // 2nd harmonic event plane
2633      { 
2634       qi.Set(phiWgt * cos(phi * 2), phiWgt * sin(phi * 2));
2635       TVector2 mQi = fQ[k][1] - qi;
2636       psii = mQi.Phi() / 2;
2637       if(psii < 0.) { psii += TMath::Pi() ; } 
2638      } 
2639      else 
2640      {
2641       qi.Set(phiWgt * cos(phi * order), phiWgt * sin(phi * order));
2642       TVector2 mQi = fQ[k][j] - qi;
2643       psii = mQi.Phi() / order;
2644       if(psii < 0.) { psii += 2*TMath::Pi()/order ; }
2645      }
2646     }
2647           
2648     // Remove autocorrelations of the second order 'particles' which are used for v1{EP1,EP2}.
2649     if (fV1Ep1Ep2 == kTRUE && order == 1) 
2650     {
2651      AliFlowSelection usedForPsi2 = *fFlowSelect ;
2652      usedForPsi2.SetHarmonic(1);
2653      if(usedForPsi2.Select(fFlowTrack))  // particle was used for Psi2
2654      {
2655       qi.Set(phiWgt * cos(phi * 2), phiWgt * sin(phi * 2));
2656       TVector2 mQi = fQ[k][1] - qi;
2657       psi2 = mQi.Phi() / 2;
2658       if(psi2 < 0.) { psi2 += TMath::Pi() ; }
2659      }
2660      else                                // particle was not used for Psi2
2661      { 
2662       psi2 = fPsi[k][1];
2663      }
2664     }
2665    }
2666    else
2667    {
2668     fHistFull[k].fHistFullHar[j].fHistYieldPtout->Fill(pt);
2669     fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout->Fill(eta, pt, phi);
2670     fHistFull[k].fHistFullHar[j].fHistYield2Dout->Fill(eta, pt);
2671     fHistFull[k].fHistFullHar[j].fHistDcaGlobout->Fill(TMath::Abs(fFlowTrack->Dca()));
2672    }
2673
2674    // Caculate v for all particles selected for correlation analysis
2675    if(fFlowSelect->SelectPart(fFlowTrack)) 
2676    {
2677     corrMultN++;
2678     
2679     float v;
2680     if (fV1Ep1Ep2 == kFALSE || order != 1) 
2681     {
2682      v = 100 * cos(order * (phi - psii)) ;
2683     }
2684     else // i.e. (fV1Ep1Ep2 == kTRUE && order == 1)
2685     {
2686      v = 100 * cos(phi + psii - 2*psi2) ;
2687     }
2688     
2689     float vFlip = v;
2690     if(eta < 0 && oddHar) { vFlip *= -1 ; }
2691     if(strlen(fFlowSelect->PidPart()) != 0) // pid, fill rapidity 
2692     {
2693      float rapidity = fFlowTrack->Y();
2694      fHistFull[k].fHistFullHar[j].fHistvObs2D->Fill(rapidity, pt, v);
2695    
2696      if(fPtRangevEta[1] > fPtRangevEta[0])  // cut is used
2697      {
2698       if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2699       {
2700        fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(rapidity, v);
2701       }
2702      }
2703      else // cut is not used, fill in any case
2704      { 
2705       fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(rapidity, v);
2706      }
2707     } 
2708     else  // no pid, fill eta
2709     {
2710      fHistFull[k].fHistFullHar[j].fHistvObs2D->Fill(eta, pt, v);
2711
2712      if(fPtRangevEta[1] > fPtRangevEta[0]) // cut is used
2713      {
2714       if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2715       {
2716        fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(eta, v);
2717       }
2718      }
2719      else // cut is not used, fill in any case
2720      { 
2721       fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(eta, v);
2722      }
2723     }
2724
2725     if(fEtaRangevPt[1] > fEtaRangevPt[0]) // cut is used
2726     { 
2727      if(TMath::Abs(eta) < fEtaRangevPt[1] && TMath::Abs(eta) >= fEtaRangevPt[0]) // check cut range, fill if in range
2728      {
2729       fHistFull[k].fHistFullHar[j].fHistvObsPt->Fill(pt, vFlip);  // for odd harmonis /-/
2730      }
2731     }
2732     else  // cut is not used, fill in any case
2733     {
2734      fHistFull[k].fHistFullHar[j].fHistvObsPt->Fill(pt, vFlip);
2735     }
2736
2737     // v_
2738     Bool_t etaPtNoCut = kTRUE;
2739     if(fPtRangevEta[1] > fPtRangevEta[0] && (pt < fPtRangevEta[0] || pt >= fPtRangevEta[1])) 
2740     {
2741      etaPtNoCut = kFALSE;
2742     }
2743     if(fEtaRangevPt[1] > fEtaRangevPt[0] && (TMath::Abs(eta) < fEtaRangevPt[0] || TMath::Abs(eta) >= fEtaRangevPt[1]))
2744     {
2745      etaPtNoCut = kFALSE;
2746     }
2747     if(etaPtNoCut) { fHistFull[k].fHistvObs->Fill(order, vFlip) ; }
2748  
2749     // Correlation of Phi of selected particles with Psi
2750     float phii = phi;
2751     if(eta < 0 && oddHar) 
2752     {
2753      phii += TMath::Pi() ; // backward particle and odd harmonic
2754      if(phii > 2*TMath::Pi()) { phii -= 2*TMath::Pi() ; }
2755     }
2756     float dPhi = phii - psii;
2757     if(dPhi < 0.)              { dPhi += 2*TMath::Pi() ; }
2758     fHistFull[k].fHistFullHar[j].fHistPhiCorr->Fill(fmod((double)dPhi, 2*TMath::Pi() / order));
2759    }
2760   }
2761  }
2762  
2763  return corrMultN ;
2764 }
2765 //-----------------------------------------------------------------------
2766 void AliFlowAnalyser::FillV0Histograms(TClonesArray* fFlowV0s)
2767 {
2768  // v0s histograms
2769
2770  cout << " V0s Loop . " << endl ;
2771
2772  fSelV0s = 0 ;
2773  for(fV0Number=0;fV0Number<fNumberOfV0s;fV0Number++) 
2774  {
2775   //fFlowV0 = (AliFlowV0*)fFlowV0s->At(fV0Number) ;       
2776   
2777   fFlowV0 = (AliFlowV0*)(fFlowV0s->AddrAt(fV0Number)) ;
2778   // cout << " looping v0 n. " << fV0Number << " * " << fFlowV0 << endl ;  
2779
2780   float mass = fFlowV0->Mass() ;                    
2781   float eta = fFlowV0->Eta() ;                      
2782   float rapidity = fFlowV0->Y() ;                   
2783   float phi = fFlowV0->Phi() ;                      
2784   float pt = fFlowV0->Pt() ;                        
2785   float totalp = fFlowV0->P() ;                     
2786   //int charge = fFlowV0->Charge() ;                
2787   float dca = fFlowV0->Dca() ;                      
2788   float lenght = fFlowV0->V0Lenght() ;              
2789   float sigma = fFlowV0->Sigma() ;                  
2790   float chi2 = fFlowV0->Chi2() ;        
2791   Char_t pid[10] ; strcpy(pid, fFlowV0->Pid()) ;    
2792   int labelPlus =  -1 ; 
2793   int labelMinus = -1 ;             
2794   AliFlowTrack* daughterPlus = fFlowV0->DaughterP() ;  
2795   AliFlowTrack* daughterMinus = fFlowV0->DaughterN() ; 
2796   if(daughterPlus)  { labelPlus = daughterPlus->Label() ; }                 
2797   if(daughterMinus) { labelMinus = daughterMinus->Label() ; }                       
2798
2799   fHistV0Mass->Fill(mass) ;
2800   fHistV0EtaPtPhi3D->Fill(eta, pt, phi) ;
2801   fHistV0YieldAll2D->Fill(eta, pt) ;
2802   fHistV0Dca->Fill(dca);
2803   fHistV0Chi2->Fill(chi2);
2804   fHistV0Lenght->Fill(lenght);
2805   fHistV0Sigma->Fill(sigma);
2806   fHistV0MassPtSlices->Fill(mass,pt);
2807   fHistV0PYall2D->Fill(rapidity, totalp) ;
2808
2809   if(fFlowSelect->SelectPart(fFlowV0))
2810   {
2811    bool inWin = fFlowSelect->SelectV0Part(fFlowV0) ;
2812    bool sx = fFlowSelect->SelectV0sxSide(fFlowV0) ;
2813    bool dx = fFlowSelect->SelectV0dxSide(fFlowV0) ;
2814    fSelV0s++ ;
2815
2816    fHistV0YieldPart2D->Fill(eta, pt) ;
2817    if(inWin) 
2818    { 
2819     fHistV0EtaPtPhi3DPart->Fill(eta, pt, phi) ;
2820     fHistV0LenghtPart->Fill(lenght);
2821     fHistV0DcaPart->Fill(dca);
2822     fHistV0MassWin->Fill(mass) ; 
2823     fHistV0BinEta->Fill(eta, eta);
2824     fHistV0BinPt->Fill(pt, pt);   
2825    }
2826    else      
2827    { 
2828     fHistV0sbEtaPtPhi3DPart->Fill(eta, pt, phi) ;
2829     fHistV0sbLenghtPart->Fill(lenght);
2830     fHistV0sbDcaPart->Fill(dca);
2831     fHistV0sbMassSide->Fill(mass) ; 
2832     fHistV0sbBinEta->Fill(eta, eta);
2833     fHistV0sbBinPt->Fill(pt, pt);
2834    }
2835
2836    for(int k = 0; k < AliFlowConstants::kSels; k++)  // sort of HarmonicsLoop - selection number used
2837    {
2838     fFlowSelect->SetSelection(k) ;
2839     for(Int_t j=0;j<AliFlowConstants::kHars;j++) 
2840     {
2841      Bool_t oddHar = (j+1) % 2 ;
2842      Float_t order = (Float_t)(j+1) ;
2843      fFlowSelect->SetHarmonic(j);       
2844
2845      // Remove autocorrelations
2846      Float_t psii, psi2 ;
2847      TVector2 q1, q2 ;
2848      Float_t phiDaughter1 = 0. ; Float_t phiDaughter2 = 0. ;
2849      Double_t phiWgt1 = 0. ; Double_t phiWgt2 = 0. ;
2850      // -
2851      if(daughterPlus)
2852      { // 1
2853       fFlowTrack = daughterPlus ;
2854       phiDaughter1 = fFlowTrack->Phi() ;                  
2855       if(fFlowSelect->Select(fFlowTrack)) // Get phiWgt from file
2856       {
2857        if(order > 3. && !oddHar) { phiWgt1 = fFlowEvent->PhiWeight(k, 1, fFlowTrack) ; } 
2858        else                      { phiWgt1 = fFlowEvent->PhiWeight(k, j, fFlowTrack) ; }
2859       }
2860      }
2861      if(daughterMinus)
2862      { // 2
2863       fFlowTrack = daughterMinus ;
2864       phiDaughter2 = fFlowTrack->Phi() ;
2865       if(fFlowSelect->Select(fFlowTrack)) // Get phiWgt from file
2866       {
2867        if(order > 3. && !oddHar) { phiWgt2 = fFlowEvent->PhiWeight(k, 1, fFlowTrack) ; } 
2868        else                      { phiWgt2 = fFlowEvent->PhiWeight(k, j, fFlowTrack) ; }
2869       } 
2870      }
2871
2872      // psi2
2873      q1.Set(phiWgt1 * cos(phiDaughter1 * 2), phiWgt1 * sin(phiDaughter1 * 2));
2874      q2.Set(phiWgt2 * cos(phiDaughter2 * 2), phiWgt2 * sin(phiDaughter2 * 2));
2875      TVector2 mQi = fQ[k][1] ; mQi -= q1 ; mQi -= q2 ; 
2876      psi2 = mQi.Phi() / 2 ;       
2877      if(psi2 < 0.)           { psi2 += TMath::Pi() ; }
2878      if(psi2 > TMath::Pi())  { psi2 -= TMath::Pi() ; } 
2879
2880      // psii
2881      if(order > 3. && !oddHar) { psii = psi2 ; } 
2882      else 
2883      {
2884       q1.Set(phiWgt1 * cos(phiDaughter1 * order), phiWgt1 * sin(phiDaughter1 * order));
2885       q2.Set(phiWgt2 * cos(phiDaughter2 * order), phiWgt2 * sin(phiDaughter2 * order));
2886       TVector2 mQi = fQ[k][j] ; mQi -= q1 ; mQi -= q2 ;
2887       psii = mQi.Phi()/order ; 
2888       if(psii < 0.)                  { psii += 2*TMath::Pi()/order ; }  
2889       if(psii > 2*TMath::Pi()/order) { psii -= 2*TMath::Pi()/order ; } 
2890      }
2891
2892     // Caculate v for all V0s selected for correlation analysis
2893      float v ;
2894      if(fV1Ep1Ep2 == kFALSE || order != 1) { v = 100 * cos(order * (phi - psii)) ; }
2895      else { v = 100 * cos(phi + psii - 2*psi2) ; }
2896      float vFlip = v ; if(eta < 0 && oddHar) { vFlip *= -1 ; }
2897
2898      // invariant mass windows & sidebands
2899      if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObs2D->Fill(eta, pt, v) ; }
2900      else      { fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D->Fill(eta, pt, v) ; }
2901
2902      if(fPtRangevEta[1] > fPtRangevEta[0]) // cut is used
2903      {
2904       if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2905       {
2906        if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsEta->Fill(eta, v) ; }
2907        else      
2908        { 
2909         fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->Fill(eta, v) ; 
2910         if(sx)      { fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx->Fill(eta, v) ; } 
2911         else if(dx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx->Fill(eta, v) ; } 
2912        }
2913       }
2914      }
2915      else // cut is not used, fill in any case
2916      { 
2917       if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsEta->Fill(eta, v) ; }
2918       else      
2919       { 
2920        fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->Fill(eta, v) ; 
2921        if(sx)      { fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx->Fill(eta, v) ; } 
2922        else if(dx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx->Fill(eta, v) ; } 
2923       }
2924      } 
2925      if(fEtaRangevPt[1] > fEtaRangevPt[0]) // cut is used
2926      {
2927       if(TMath::Abs(eta) < fEtaRangevPt[1] && TMath::Abs(eta) >= fEtaRangevPt[0]) // check cut range, fill if in range
2928       {          
2929        if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsPt->Fill(pt, vFlip) ; }      // for odd harmonis /-/
2930        else      
2931        { 
2932         fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->Fill(pt, vFlip) ; 
2933         if(sx)      { fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx->Fill(eta, v) ; } 
2934         else if(dx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx->Fill(eta, v) ; }  
2935        }
2936       }
2937      }
2938      else  // cut is not used, fill in any case
2939      {
2940       if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsPt->Fill(pt, vFlip) ; } 
2941       else      
2942       { 
2943        fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->Fill(pt, vFlip) ; 
2944        if(sx)      { fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx->Fill(eta, v) ; } 
2945        else if(dx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx->Fill(eta, v) ; }  
2946       }
2947      }
2948      // v_
2949      Bool_t etaPtNoCut = kTRUE;
2950      if(fPtRangevEta[1] > fPtRangevEta[0] && (pt < fPtRangevEta[0] || pt >= fPtRangevEta[1])) 
2951      {
2952       etaPtNoCut = kFALSE;
2953      }
2954      if(fEtaRangevPt[1] > fEtaRangevPt[0] && (TMath::Abs(eta) < fEtaRangevPt[0] || TMath::Abs(eta) >= fEtaRangevPt[1]))
2955      {
2956       etaPtNoCut = kFALSE;
2957      }
2958      if(etaPtNoCut) 
2959      { 
2960       if(inWin) { fHistFull[k].fHistV0vObs->Fill(order, vFlip) ; }
2961       else      
2962       { 
2963        if(sx)      { fHistFull[k].fHistV0sbvObsSx->Fill(order, vFlip) ; } 
2964        else if(dx) { fHistFull[k].fHistV0sbvObsDx->Fill(order, vFlip) ; }  
2965       }
2966      }
2967    
2968      // Correlation of Phi of selected v0s with Psi
2969      float phii = phi;
2970      if(eta < 0 && oddHar)
2971      {
2972       phii += TMath::Pi() ; // backward particle and odd harmonic
2973       if(phii > 2*TMath::Pi()) { phii -= 2*TMath::Pi() ; }
2974      }
2975      float dPhi = phii - psii;
2976      if(dPhi < 0.)  { dPhi += 2*TMath::Pi() ; }
2977
2978      if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0PhiCorr->Fill(fmod((double)dPhi, 2*TMath::Pi() / order)) ; } 
2979      else      { fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr->Fill(fmod((double)dPhi, 2*TMath::Pi() / order)) ; }
2980     }
2981    }
2982   }
2983   //delete fFlowV0 ; fFlowV0 = 0 ;
2984   //delete daughterPlus ; daughterPlus = 0 ;  
2985   //delete daughterMinus ; daughterMinus = 0 ;  
2986  }
2987  fHistV0MultPart->Fill(fSelV0s) ;
2988
2989  return ;
2990 }
2991 //-----------------------------------------------------------------------
2992 // void AliFlowAnalyser::FillLabels()  
2993 // {
2994 //  // Reads the tracks label (i.e. the link from ESD tracking to KineTree)
2995 //  // and ...
2996 //  
2997 //  cout << " Fill Labels . " << endl ;  
2998 // 
2999 //  fLabHist ;
3000 //  return ;
3001 // }
3002 //-----------------------------------------------------------------------
3003 void AliFlowAnalyser::PrintEventQuantities() const
3004 {
3005  // prints event by event calculated quantities
3006
3007  cout << endl ; cout << " # " << fEventNumber << " - Event quantities : " << endl ; cout << endl ; 
3008  
3009  cout << "fQ[k][j] = " ; 
3010  for(int k=0;k<AliFlowConstants::kSels;k++) 
3011  {
3012   for(int j=0;j<AliFlowConstants::kHars;j++) 
3013   { 
3014    cout << (Float_t)fQ[k][j].X() << "," << (Float_t)fQ[k][j].Y() << " ; " ; 
3015   }
3016   cout << endl ;
3017   cout << "            " ; 
3018  }
3019  cout << endl ; cout << "fPsi[k][j] = " ; 
3020  for(int k=0;k<AliFlowConstants::kSels;k++) 
3021  {
3022   for(int j=0;j<AliFlowConstants::kHars;j++) { Float_t aaa = (Float_t)fPsi[k][j] ; cout << aaa << " , " ; }
3023   cout << endl ;
3024   cout << "              " ; 
3025  }
3026  cout << endl ; cout << "fQnorm[k][j] = " ; 
3027  for(int k=0;k<AliFlowConstants::kSels;k++) 
3028  {
3029   for(int j=0;j<AliFlowConstants::kHars;j++) { Float_t aaa = (Float_t)fQnorm[k][j] ; cout << aaa << " , " ; }
3030   cout << endl ;
3031   cout << "             " ; 
3032  }
3033  cout << endl ; cout << "fMult[k][j] = " ; 
3034  for(int k=0;k<AliFlowConstants::kSels;k++) 
3035  {
3036   for(int j=0;j<AliFlowConstants::kHars;j++) { Float_t aaa = (Float_t)fMult[k][j] ; cout << aaa << " , " ; } 
3037   cout << endl ;
3038   cout << "               " ; 
3039  }
3040  cout << endl ; cout << "fMultSub[n][k][j] = " ; 
3041  for(int n=0;n<AliFlowConstants::kSubs;n++) 
3042  {
3043   for(int k=0;k<AliFlowConstants::kSels;k++) 
3044   {
3045    for(int j=0;j<AliFlowConstants::kHars;j++) { Float_t aaa = fMultSub[n][k][j] ; cout << aaa << " , " ; } 
3046    cout << endl ;
3047    cout << "                 " ; 
3048   }
3049  }
3050  cout << endl ; cout << "fPsiSub[n][k][j] = " ; 
3051  for(int n=0;n<AliFlowConstants::kSubs;n++) 
3052  {
3053   for(int k=0;k<AliFlowConstants::kSels;k++) 
3054   {
3055    for(int j=0;j<AliFlowConstants::kHars;j++) { Float_t aaa = fPsiSub[n][k][j] ; cout << aaa << " , " ; } 
3056    cout << endl ;
3057    cout << "                 " ; 
3058   }
3059  }
3060  cout << endl ; cout << "Delta_PsiSub[k][j] = " ; 
3061  for(int k=0;k<AliFlowConstants::kSels;k++) 
3062  {
3063   for(int j=0;j<AliFlowConstants::kHars;j++) { Float_t aaa = fPsiSub[0][k][j]-fPsiSub[1][k][j] ; cout << aaa << " , " ; } 
3064   cout << endl ;
3065   cout << "             " ; 
3066  }
3067  cout << endl ;
3068
3069  cout << " N. of selected tracks - SelPart(AliFlowTrack*) = " << fSelParts << endl ;                                    //! n. of tracks selected for correlation analysis
3070  cout << " N. of selected V0s    - SelPart(AliFlowV0*)    = " << fSelV0s << endl ;                                      //! n. of v0s selected for correlation analysis
3071
3072  cout << endl ;
3073 }
3074 //----------------------------------------------------------------------
3075 // ###
3076 //----------------------------------------------------------------------
3077 Bool_t AliFlowAnalyser::Resolution()
3078 {
3079  // Calculates the resolution from cos(dPsi) and corrects the observed flow.
3080  // New histograms are then created and filled with the corrected vn 
3081  // (see fVnResHistList), and saved in the otput file.
3082  
3083  cout << " AliFlowAnalyser::Resolution()" << endl ;
3084  
3085  // VnRes histogram collection
3086  fVnResHistList = new TOrdCollection(AliFlowConstants::kSels*AliFlowConstants::kHars);
3087  
3088  // Calculate resolution from sqrt(fHistCos)
3089  double cosPair[AliFlowConstants::kSels][AliFlowConstants::kHars];
3090  double cosPairErr[AliFlowConstants::kSels][AliFlowConstants::kHars];
3091  double content, contentV0;
3092  double error, errorV0;
3093  double totalError;
3094  TString* histTitle;
3095
3096  for(int k = 0; k < AliFlowConstants::kSels; k++)
3097  {
3098   // v for tracks (corrected for resolution)
3099   histTitle = new TString("Flow_v_Sel");
3100   *histTitle += k+1;
3101   fHistFull[k].fHistv = fHistFull[k].fHistvObs->ProjectionX(histTitle->Data());
3102   fHistFull[k].fHistv->SetTitle(histTitle->Data());
3103   fHistFull[k].fHistv->SetXTitle("Harmonic");
3104   fHistFull[k].fHistv->SetYTitle("v (%)");
3105   delete histTitle;
3106   fVnResHistList->AddLast(fHistFull[k].fHistv);
3107    
3108   // v for v0s (corrected for resolution)
3109   histTitle = new TString("FlowV0_v_Sel");
3110   *histTitle += k+1;
3111   fHistFull[k].fHistV0v = fHistFull[k].fHistV0vObs->ProjectionX(histTitle->Data());
3112   fHistFull[k].fHistV0v->SetTitle(histTitle->Data());
3113   fHistFull[k].fHistV0v->SetXTitle("Harmonic");
3114   fHistFull[k].fHistV0v->SetYTitle("v (%)");
3115   delete histTitle;
3116   fVnResHistList->AddLast(fHistFull[k].fHistV0v);
3117    
3118   for(int j = 0; j < AliFlowConstants::kHars; j++) 
3119   {
3120    double order = (double)(j+1);
3121    cosPair[k][j]    = fHistFull[k].fHistCos->GetBinContent(j+1);
3122    cosPairErr[k][j] = fHistFull[k].fHistCos->GetBinError(j+1);
3123    if(cosPair[k][j] > 0.) 
3124    {
3125     double resSub = 0. ;
3126     double resSubErr = 0. ;
3127     if(fV1Ep1Ep2 == kTRUE && order == 1) // calculate resolution of second order event plane first
3128     {
3129      double res2 = 0. ;
3130      double res2error = 0. ;
3131      if(fHistFull[k].fHistCos->GetBinContent(2) > 0.) 
3132      {
3133       if(fEtaSub)  // sub res only
3134       {
3135        res2 = TMath::Sqrt(fHistFull[k].fHistCos->GetBinContent(2));
3136        res2error = fHistFull[k].fHistCos->GetBinError(2) / (2. * res2);
3137       }
3138       else 
3139       {
3140        if (fHistFull[k].fHistCos->GetBinContent(2) > 0.92)  // resolution saturates
3141        {
3142         res2 = 0.99;
3143         res2error = 0.007;
3144        } 
3145        else 
3146        {
3147         double deltaRes2Sub = 0.005;  // differential for the error propagation
3148         double res2Sub = TMath::Sqrt(fHistFull[k].fHistCos->GetBinContent(2));
3149         double res2SubErr = fHistFull[k].fHistCos->GetBinError(2) / (2. * res2Sub);
3150         double chiSub2 = Chi(res2Sub);
3151         double chiSub2Delta = Chi(res2Sub + deltaRes2Sub);
3152         res2 = ResEventPlane(TMath::Sqrt(2.) * chiSub2); // full event plane res.
3153         double mRes2Delta = ResEventPlane(TMath::Sqrt(2.) * chiSub2Delta);
3154         res2error = res2SubErr * fabs((double)res2 - mRes2Delta) / deltaRes2Sub;
3155        }
3156       }
3157      }
3158      else 
3159      {
3160       res2 = 0.;
3161       res2error = 0.;
3162      }
3163      // now put everything together with first order event plane
3164      fRes[k][j]    = TMath::Sqrt(cosPair[k][0]*res2);
3165      fResErr[k][j] = 1./(2.*fRes[k][j]) * TMath::Sqrt(cosPairErr[k][0]*cosPairErr[k][0] + res2error*res2error); // Gaussian error propagation
3166     }
3167     else if(fEtaSub)  // sub res only
3168     {
3169      resSub = TMath::Sqrt(cosPair[k][j]);                
3170      resSubErr = cosPairErr[k][j] / (2. * resSub);
3171      fRes[k][j]    = resSub;
3172      fResErr[k][j] = resSubErr;
3173     } 
3174     else if(order==4. || order==6.|| order==8.)  // 2nd harmonic event plane
3175     {
3176      double deltaResSub = 0.005;  // differential for the error propagation
3177      double resSub = TMath::Sqrt(cosPair[k][1]);
3178      double resSubErr = cosPairErr[k][1] / (2. * resSub);
3179      double chiSub = Chi(resSub);
3180      double chiSubDelta = Chi(resSub + deltaResSub);
3181      double mResDelta;
3182      if (order==4.) 
3183      {
3184       fRes[k][j] = ResEventPlaneK2(TMath::Sqrt(2.) * chiSub); // full event plane res.
3185       mResDelta = ResEventPlaneK2(TMath::Sqrt(2.) * chiSubDelta);
3186      } 
3187      else if(order==6.) 
3188      {
3189       fRes[k][j] = ResEventPlaneK3(TMath::Sqrt(2.) * chiSub); // full event plane res.
3190       mResDelta = ResEventPlaneK3(TMath::Sqrt(2.) * chiSubDelta);
3191      } 
3192      else 
3193      {
3194       fRes[k][j] = ResEventPlaneK4(TMath::Sqrt(2.) * chiSub); // full event plane res.
3195       mResDelta = ResEventPlaneK4(TMath::Sqrt(2.) * chiSubDelta);
3196      }
3197      fResErr[k][j] = resSubErr * fabs((double)fRes[k][j] - mResDelta) / deltaResSub;
3198     }
3199     else 
3200     {
3201      if(cosPair[k][j] > 0.92) // resolution saturates
3202      {
3203       fRes[k][j]    = 0.99;
3204       fResErr[k][j] = 0.007;
3205      } 
3206      else 
3207      {
3208       double deltaResSub = 0.005;  // differential for the error propagation
3209       double resSub = TMath::Sqrt(cosPair[k][j]);
3210       double resSubErr = cosPairErr[k][j] / (2. * resSub);
3211       double chiSub = Chi(resSub);
3212       double chiSubDelta = Chi(resSub + deltaResSub);
3213       fRes[k][j] = ResEventPlane(TMath::Sqrt(2.) * chiSub); // full event plane res.
3214       double mResDelta = ResEventPlane(TMath::Sqrt(2.) * chiSubDelta);
3215       fResErr[k][j] = resSubErr * TMath::Abs((double)fRes[k][j] - mResDelta) / deltaResSub;
3216      }
3217     }
3218    }
3219    else  // subevent correlation must be positive
3220    {
3221     fRes[k][j]    = 0.;
3222     fResErr[k][j] = 0.;
3223    }
3224    fHistFull[k].fHistRes->SetBinContent(j+1, fRes[k][j]);
3225    fHistFull[k].fHistRes->SetBinError(j+1, fResErr[k][j]);
3226
3227    // Create the v 2D histogram (Flow corrected for res.) - v,Pt,Eta
3228    histTitle = new TString("Flow_v2D_Sel");
3229    *histTitle += k+1;
3230    histTitle->Append("_Har");
3231    *histTitle += j+1;
3232    fHistFull[k].fHistFullHar[j].fHistv2D = fHistFull[k].fHistFullHar[j].fHistvObs2D->ProjectionXY(histTitle->Data());
3233    fHistFull[k].fHistFullHar[j].fHistv2D->SetTitle(histTitle->Data());
3234    fHistFull[k].fHistFullHar[j].fHistv2D->SetXTitle((char*)fLabel.Data());
3235    fHistFull[k].fHistFullHar[j].fHistv2D->SetYTitle("Pt (GeV/c)");
3236    fHistFull[k].fHistFullHar[j].fHistv2D->SetZTitle("v (%)");
3237    delete histTitle;
3238    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistv2D);
3239
3240    // Create the 1D v histograms - v,Eta
3241    histTitle = new TString("Flow_vEta_Sel");
3242    *histTitle += k+1;
3243    histTitle->Append("_Har");
3244    *histTitle += j+1;
3245    fHistFull[k].fHistFullHar[j].fHistvEta = fHistFull[k].fHistFullHar[j].fHistvObsEta->ProjectionX(histTitle->Data());
3246    fHistFull[k].fHistFullHar[j].fHistvEta->SetTitle(histTitle->Data());
3247    fHistFull[k].fHistFullHar[j].fHistvEta->SetXTitle((char*)fLabel.Data());
3248    fHistFull[k].fHistFullHar[j].fHistvEta->SetYTitle("v (%)");
3249    delete histTitle;
3250    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistvEta);
3251    
3252    // v,Pt
3253    histTitle = new TString("Flow_vPt_Sel");
3254    *histTitle += k+1;
3255    histTitle->Append("_Har");
3256    *histTitle += j+1;
3257    fHistFull[k].fHistFullHar[j].fHistvPt = fHistFull[k].fHistFullHar[j].fHistvObsPt->ProjectionX(histTitle->Data());
3258    fHistFull[k].fHistFullHar[j].fHistvPt->SetTitle(histTitle->Data());
3259    fHistFull[k].fHistFullHar[j].fHistvPt->SetXTitle("Pt (GeV/c)");
3260    fHistFull[k].fHistFullHar[j].fHistvPt->SetYTitle("v (%)");
3261    delete histTitle;
3262    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistvPt);
3263
3264    // Create the v 2D histogram (V0s Flow corrected for res.) - v,Pt,Eta
3265    histTitle = new TString("FlowV0_v2D_Sel");
3266    *histTitle += k+1;
3267    histTitle->Append("_Har");
3268    *histTitle += j+1;
3269    fHistFull[k].fHistFullHar[j].fHistV0v2D = fHistFull[k].fHistFullHar[j].fHistV0vObs2D->ProjectionXY(histTitle->Data());
3270    fHistFull[k].fHistFullHar[j].fHistV0v2D->SetTitle(histTitle->Data());
3271    fHistFull[k].fHistFullHar[j].fHistV0v2D->SetXTitle((char*)fLabel.Data());
3272    fHistFull[k].fHistFullHar[j].fHistV0v2D->SetYTitle("Pt (GeV/c)");
3273    fHistFull[k].fHistFullHar[j].fHistV0v2D->SetZTitle("v (%)");
3274    delete histTitle;
3275    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0v2D);
3276
3277    // Create the 1D v histograms (V0s) - v,Eta
3278    histTitle = new TString("FlowV0_vEta_Sel");
3279    *histTitle += k+1;
3280    histTitle->Append("_Har");
3281    *histTitle += j+1;
3282    fHistFull[k].fHistFullHar[j].fHistV0vEta = fHistFull[k].fHistFullHar[j].fHistV0vObsEta->ProjectionX(histTitle->Data());
3283    fHistFull[k].fHistFullHar[j].fHistV0vEta->SetTitle(histTitle->Data());
3284    fHistFull[k].fHistFullHar[j].fHistV0vEta->SetXTitle((char*)fLabel.Data());
3285    fHistFull[k].fHistFullHar[j].fHistV0vEta->SetYTitle("v (%)");
3286    delete histTitle;
3287    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0vEta);
3288    
3289    // (V0s) v,Pt
3290    histTitle = new TString("FlowV0_vPt_Sel");
3291    *histTitle += k+1;
3292    histTitle->Append("_Har");
3293    *histTitle += j+1;
3294    fHistFull[k].fHistFullHar[j].fHistV0vPt = fHistFull[k].fHistFullHar[j].fHistV0vObsPt->ProjectionX(histTitle->Data());
3295    fHistFull[k].fHistFullHar[j].fHistV0vPt->SetTitle(histTitle->Data());
3296    fHistFull[k].fHistFullHar[j].fHistV0vPt->SetXTitle("Pt (GeV/c)");
3297    fHistFull[k].fHistFullHar[j].fHistV0vPt->SetYTitle("v (%)");
3298    delete histTitle;
3299    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0vPt);
3300
3301    // Create the v 2D histogram (V0s sidebands Flow corrected for res.) - v,Pt,Eta
3302    histTitle = new TString("FlowV0sb_v2D_Sel");
3303    *histTitle += k+1;
3304    histTitle->Append("_Har");
3305    *histTitle += j+1;
3306    fHistFull[k].fHistFullHar[j].fHistV0sbv2D = fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D->ProjectionXY(histTitle->Data());
3307    fHistFull[k].fHistFullHar[j].fHistV0sbv2D->SetTitle(histTitle->Data());
3308    fHistFull[k].fHistFullHar[j].fHistV0sbv2D->SetXTitle((char*)fLabel.Data());
3309    fHistFull[k].fHistFullHar[j].fHistV0sbv2D->SetYTitle("Pt (GeV/c)");
3310    fHistFull[k].fHistFullHar[j].fHistV0sbv2D->SetZTitle("v (%)");
3311    delete histTitle;
3312    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbv2D);
3313
3314    // Create the 1D v histograms (V0s sidebands) - v,Eta
3315    histTitle = new TString("FlowV0sb_vEta_Sel");
3316    *histTitle += k+1;
3317    histTitle->Append("_Har");
3318    *histTitle += j+1;
3319    fHistFull[k].fHistFullHar[j].fHistV0sbvEta = fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->ProjectionX(histTitle->Data());
3320    fHistFull[k].fHistFullHar[j].fHistV0sbvEta->SetTitle(histTitle->Data());
3321    fHistFull[k].fHistFullHar[j].fHistV0sbvEta->SetXTitle((char*)fLabel.Data());
3322    fHistFull[k].fHistFullHar[j].fHistV0sbvEta->SetYTitle("v (%)");
3323    delete histTitle;
3324    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbvEta);
3325    
3326    // (V0s sidebands) v,Pt
3327    histTitle = new TString("FlowV0sb_vPt_Sel");
3328    *histTitle += k+1;
3329    histTitle->Append("_Har");
3330    *histTitle += j+1;
3331    fHistFull[k].fHistFullHar[j].fHistV0sbvPt = fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->ProjectionX(histTitle->Data());
3332    fHistFull[k].fHistFullHar[j].fHistV0sbvPt->SetTitle(histTitle->Data());
3333    fHistFull[k].fHistFullHar[j].fHistV0sbvPt->SetXTitle("Pt (GeV/c)");
3334    fHistFull[k].fHistFullHar[j].fHistV0sbvPt->SetYTitle("v (%)");
3335    delete histTitle;
3336    fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbvPt);
3337
3338   // Calulate v = vObs / Resolution
3339    if(fRes[k][j]) 
3340    {
3341     cout << "***** Resolution of the " << j+1 << "th harmonic = " << fRes[k][j] << " +/- " << fResErr[k][j] << endl;
3342     // selected tracks
3343     fHistFull[k].fHistFullHar[j].fHistv2D->Scale(1. / fRes[k][j]);
3344     fHistFull[k].fHistFullHar[j].fHistvEta->Scale(1. / fRes[k][j]);
3345     fHistFull[k].fHistFullHar[j].fHistvPt->Scale(1. / fRes[k][j]);
3346     content = fHistFull[k].fHistv->GetBinContent(j+1);        
3347     content /=  fRes[k][j]; 
3348     fHistFull[k].fHistv->SetBinContent(j+1, content);            
3349     // selected V0s    
3350     fHistFull[k].fHistFullHar[j].fHistV0v2D->Scale(1. / fRes[k][j]);
3351     fHistFull[k].fHistFullHar[j].fHistV0vEta->Scale(1. / fRes[k][j]);
3352     fHistFull[k].fHistFullHar[j].fHistV0vPt->Scale(1. / fRes[k][j]);
3353     contentV0 = fHistFull[k].fHistV0v->GetBinContent(j+1);        
3354     contentV0 /=  fRes[k][j];
3355     fHistFull[k].fHistV0v->SetBinContent(j+1, contentV0);            
3356     // V0s sidebands    
3357     fHistFull[k].fHistFullHar[j].fHistV0sbv2D->Scale(1. / fRes[k][j]);
3358     fHistFull[k].fHistFullHar[j].fHistV0sbvEta->Scale(1. / fRes[k][j]);
3359     fHistFull[k].fHistFullHar[j].fHistV0sbvPt->Scale(1. / fRes[k][j]);
3360     // The systematic error of the resolution is folded in.
3361     // tracks
3362     error = fHistFull[k].fHistv->GetBinError(j+1) ; 
3363     error /= fRes[k][j];                                         
3364     totalError = fabs(content) * TMath::Sqrt((error/content)*(error/content) + (fResErr[k][j]/fRes[k][j])*(fResErr[k][j]/fRes[k][j]));
3365     fHistFull[k].fHistv->SetBinError(j+1, totalError);
3366     cout << "***** v" << j+1 << "= (" << content << " +/- " << error << " +/- " << totalError << "(with syst.)) %" << endl;
3367     // V0s
3368     errorV0 = fHistFull[k].fHistV0v->GetBinError(j+1) ; 
3369     errorV0 /= fRes[k][j]; 
3370     if (contentV0>0)
3371       totalError = fabs(contentV0) * TMath::Sqrt((errorV0/contentV0)*(errorV0/contentV0) + (fResErr[k][j]/fRes[k][j])*(fResErr[k][j]/fRes[k][j]));
3372     else
3373       totalError = 0;
3374     fHistFull[k].fHistV0v->SetBinError(j+1, totalError);
3375    } 
3376    else 
3377    {
3378     cout << "##### Resolution of the " << j+1 << "th harmonic was zero." << endl;
3379     // tracks hist
3380     fHistFull[k].fHistFullHar[j].fHistv2D->Reset();
3381     fHistFull[k].fHistFullHar[j].fHistvEta->Reset();
3382     fHistFull[k].fHistFullHar[j].fHistvPt->Reset();
3383     fHistFull[k].fHistv->SetBinContent(j+1, 0.);
3384     fHistFull[k].fHistv->SetBinError(j+1, 0.);
3385     // v0s hist
3386     fHistFull[k].fHistFullHar[j].fHistV0v2D->Reset();
3387     fHistFull[k].fHistFullHar[j].fHistV0vEta->Reset();
3388     fHistFull[k].fHistFullHar[j].fHistV0vPt->Reset();
3389     fHistFull[k].fHistFullHar[j].fHistV0sbv2D->Reset();
3390     fHistFull[k].fHistFullHar[j].fHistV0sbvEta->Reset();
3391     fHistFull[k].fHistFullHar[j].fHistV0sbvPt->Reset();
3392     fHistFull[k].fHistV0v->SetBinContent(j+1, 0.);
3393     fHistFull[k].fHistV0v->SetBinError(j+1, 0.);
3394    }
3395   }
3396  }
3397
3398  return kTRUE ;
3399 }
3400 //-----------------------------------------------------------------------
3401 Double_t AliFlowAnalyser::Chi(Double_t res)
3402 {          
3403  // chi from the event plane resolution
3404
3405  Double_t chi   = 2.0;
3406  Double_t delta = 1.0;
3407  for(int i = 0; i < 15; i++) 
3408  {
3409   if(ResEventPlane(chi) < res) { chi = chi + delta ; }
3410   else                         { chi = chi - delta ; }
3411   delta = delta / 2.;
3412  }
3413
3414  return chi ;
3415 }
3416 //-----------------------------------------------------------------------
3417 Double_t AliFlowAnalyser::ResEventPlane(Double_t chi)    
3418 {
3419  // plane resolution as function of chi
3420
3421  Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ;   // ~ 0.626657
3422  Double_t arg = chi * chi / 4.;
3423  Double_t res = con * chi * exp(-arg) * (TMath::BesselI0(arg) + TMath::BesselI1(arg)); 
3424
3425  return res ;
3426 }
3427 //-----------------------------------------------------------------------
3428 Double_t AliFlowAnalyser::ResEventPlaneK2(Double_t chi)  
3429 {
3430  // ... for the case k=2.
3431
3432  Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ;   // ~ 0.626657
3433  Double_t arg = chi * chi / 4.;
3434
3435  Double_t besselOneHalf = TMath::Sqrt(arg/(TMath::Pi()/2)) * sinh(arg)/arg;
3436  Double_t besselThreeHalfs = TMath::Sqrt(arg/(TMath::Pi()/2)) * (cosh(arg)/arg - sinh(arg)/(arg*arg));
3437  Double_t res = con * chi * exp(-arg) * (besselOneHalf + besselThreeHalfs); 
3438
3439  return res ;
3440 }
3441 //-----------------------------------------------------------------------
3442 Double_t AliFlowAnalyser::ResEventPlaneK3(Double_t chi) 
3443 {
3444  // ...for the case k=3.
3445
3446  Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ;   // ~ 0.626657
3447  Double_t arg = chi * chi / 4.;
3448  Double_t res = con * chi * exp(-arg) * (TMath::BesselI1(arg) + TMath::BesselI(2, arg)); 
3449
3450  return res ;
3451 }
3452 //-----------------------------------------------------------------------
3453 Double_t AliFlowAnalyser::ResEventPlaneK4(Double_t chi) 
3454 {
3455  // ...for the case k=4.
3456
3457  Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ;   // ~ 0.626657
3458  Double_t arg = chi * chi / 4.;
3459
3460  Double_t besselOneHalf = TMath::Sqrt(arg/(TMath::Pi()/2)) * sinh(arg)/arg;
3461  Double_t besselThreeHalfs = TMath::Sqrt(arg/(TMath::Pi()/2)) * (cosh(arg)/arg - sinh(arg)/(arg*arg));
3462  Double_t besselFiveHalfs = besselOneHalf - 3*besselThreeHalfs/arg;
3463  Double_t res = con * chi * exp(-arg) * (besselThreeHalfs + besselFiveHalfs); 
3464
3465  return res ;
3466 }
3467 //-----------------------------------------------------------------------
3468 // ###
3469 //-----------------------------------------------------------------------
3470
3471
3472 #endif