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