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