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