PWG2/FLOW/AliFlowAnalyser.cxx

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