PWGCF/FLOW/Attic/AliFlowAnalyser.cxx

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