AliAnalysisTaskFlowModesMC.h

Go to the documentation of this file.

/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. */
/* See cxx source for full Copyright notice */
/* $Id$ */

#ifndef AliAnalysisTaskFlowModesMC_H
#define AliAnalysisTaskFlowModesMC_H

#include <vector>

#include "AliAnalysisTaskSE.h"
#include "AliAODVertex.h"
#include "AliAODTrack.h"
#include "AliMCParticle.h"
#include "AliPIDResponse.h"
#include "AliPIDCombined.h"
#include "AliESDpid.h"
#include "AliFlowBayesianPID.h"

#include "TComplex.h"
#include "TFile.h"
#include "TList.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TProfile.h"
#include "TProfile2D.h"



class AliAnalysisTaskFlowModesMC : public AliAnalysisTaskSE
{
public:
    enum    RunMode {kTest, kFillWeights, kFull}; // task running mode (NOT GRID MODE)
    enum    AnalType {kAOD, kESD, kMC, kMCAOD}; // tag for analysis type
    enum    PartSpecies {kUnknown, kCharged, kPion, kKaon, kProton}; // list of all particle species of interest
    enum    ColSystem {kPP, kPbPb};
    struct FlowPart // representation of selected particle (species independent) storing only basic properties for flow calculations
    {
        FlowPart(Double_t dPt = 0, Double_t dPhi = 0, Double_t dEta = 0, Short_t iCharge = 0, PartSpecies sSpecies = kUnknown, Double_t dMass = 0, Double_t dPx = 0, Double_t dPy = 0, Double_t dPz = 0) :
        pt(dPt), px(dPx), py(dPy), pz(dPz), phi(dPhi), eta(dEta), mass(dMass), charge(iCharge), species(sSpecies) {} // constructor
        
        void    PrintPart() const { printf("pt %g | px %g | py %g| pz %g | phi %g | eta %g | mass %g | charge %d | species %d \n",pt,px,py,pz,phi,eta,mass,charge,species); } // print struct members
        
        Double_t pt,px,py,pz,phi,eta,mass;
        Short_t charge;
        PartSpecies species;
    };
    AliAnalysisTaskFlowModesMC(); // constructor
    AliAnalysisTaskFlowModesMC(const char *name); // named (primary) constructor
    virtual                 ~AliAnalysisTaskFlowModesMC(); // destructor
    
    virtual void            UserCreateOutputObjects(); //
    virtual void            UserExec(Option_t* option); // main methond - called for each event
    virtual void            Terminate(Option_t* option); // called after all events are processed
    // analysis setters
    void                    SetColisionSystem(ColSystem colSystem = kPbPb) {fColSystem = colSystem; }
    void                    SetRunMode(RunMode mode = kFull) { fRunMode = mode; }
    void                    SetNumEventsAnalyse(Short_t num) { fNumEventsAnalyse = num; }
    void            SetAnalysisType(AnalType type = kAOD) { fAnalType = type; }
    void                    SetSampling(Bool_t sample = kTRUE) { fSampling = sample; }
    void                    SetImpactParameterRange(Double_t min=0., Double_t max=14.) { fImpactParameterMin=min;fImpactParameterMax=max;fCutImpactParameter=kTRUE; }
    void                    SetFillQAhistos(Bool_t fill = kTRUE) { fFillQA = fill; }
    void                    SetProcessCharged(Bool_t filter = kTRUE) { fProcessCharged = filter; }
    void                    SetProcessPID(Bool_t filter = kTRUE, Bool_t PIDbayesian = kFALSE) {
        fProcessPID = filter;
        if(PIDbayesian){fPIDbayesian = PIDbayesian; fPIDnsigma = kFALSE;}else{fPIDbayesian = kFALSE; fPIDnsigma = kTRUE;}
    }
    // flow related setters
    void                    SetFlowRFPsPtMin(Float_t pt) { fCutFlowRFPsPtMin = pt; }
    void                    SetFlowRFPsPtMax(Float_t pt) { fCutFlowRFPsPtMax = pt; }
    void                    SetFlowDoFourCorrelations(Bool_t four = kTRUE) { fCutFlowDoFourCorrelations = four; }
    void                    SetFlowDoOnlyMixedCorrelations(Bool_t b = kFALSE) { fDoOnlyMixedCorrelations = b; }
    void                    SetFlowFillWeights(Bool_t weights = kTRUE) { fFlowFillWeights = weights; }
    void                    SetUseNUAWeigthsFile(const char* file) { fFlowNUAWeightsPath = file; if(fFlowNUAWeightsPath){fFlowUseNUAWeights = kTRUE;} }
    void                    SetUseNUEWeigthsFile(const char* file) { fFlowNUEWeightsPath = file; if(fFlowNUEWeightsPath){fFlowUseNUEWeights = kTRUE;} }
    
    // events setters
    void                    SetMultEstimator(const char* mult = "CHARGED") { fMultEstimator = mult; }
    void                    SetTrigger(Short_t trigger = 0) { fTrigger = trigger; }
    void                      SetPVtxZMax(Double_t z) { fPVtxCutZ = z; }
    void                    SetCentralityRange(Bool_t kCentRange){fFullCentralityRange = kCentRange;}
    // track setters
    void                    SetChargedEtaMax(Double_t eta) { fCutChargedEtaMax = eta; }
    void                    SetChargedPtMax(Double_t pt) { fCutChargedPtMax = pt; }
    void                    SetChargedPtMin(Double_t pt) { fCutChargedPtMin = pt; }
    void                    SetChargedDCAzMax(Double_t dcaz) {  fCutChargedDCAzMax = dcaz; }
    void                    SetChargedDCAxyMax(Double_t dcaxy) {  fCutChargedDCAxyMax = dcaxy; }
    void                    SetChargedNumTPCclsMin(UShort_t tpcCls) { fCutChargedNumTPCclsMin = tpcCls; }
    void                      SetMaxChi2perTPCcls(Double_t MaxChi2 = 4){ fMaxChi2perTPCcls = MaxChi2;}
    void                    SetChargedTrackFilterBit(UInt_t filter) { fCutChargedTrackFilterBit = filter; }
    // PID (pi,K,p) setters
    
    void                    SetPIDUseAntiProtonOnly(Bool_t use = kTRUE) { fCutPIDUseAntiProtonOnly = use; }
    void                    SetPIDNumSigmasPionMax(Double_t numSigmas) { fCutPIDnSigmaPionMax = numSigmas; }
    void                    SetPIDNumSigmasKaonMax(Double_t numSigmas) { fCutPIDnSigmaKaonMax = numSigmas; }
    void                    SetPIDNumSigmasProtonMax(Double_t numSigmas) { fCutPIDnSigmaProtonMax = numSigmas; }
    void                    SetPIDNumSigmasCombinedNoTOFrejection(Bool_t reject = kTRUE) { fCutPIDnSigmaCombinedNoTOFrejection = reject; }
    void                    SetPIDnsigmaCombination(Int_t Comb =2){fPIDnsigmaCombination = Comb;}
    void            SetExtraPileUpCut(){fExtraPileUp = kTRUE;}
    void                    SetPositivelyChargedRef(Bool_t Pos=kFALSE){fPositivelyChargedRef = Pos;}
    void                    SetNegativelyChargedRef(Bool_t Neg=kFALSE){fNegativelyChargedRef = Neg;}
    void                    SetPositivelyChargedPOI(Bool_t Pos=kFALSE){fPositivelyChargedPOI = Pos;}
    void                    SetNegativelyChargedPOI(Bool_t Neg=kFALSE){fNegativelyChargedPOI = Neg;}
    void                    SetBayesianProbability(Double_t prob=0.9){fParticleProbability = prob;}
    void                    SetPriors(Float_t centr = 0); // set Noferini's favourite priors for Bayesian PID (requested if Bayesian PID is used)
    AliESDpid&              GetESDpid() {return fESDpid;}
    Bool_t                  TPCTOFagree(const AliVTrack *track);
    
private:
    
    // properties
    AliAODEvent*            fEventAOD; 
    AliMCEvent*             fEventMC; 
    AliPIDResponse*         fPIDResponse; 
    AliPIDCombined*         fPIDCombined; 
    TFile*                  fFlowNUAWeightsFile; 
    TFile*                  fFlowNUEWeightsFile; 
    Bool_t                  fInit; // initialization check
    Short_t                 fIndexSampling; // sampling index (randomly generated)
    Short_t                 fIndexCentrality; // centrality bin index (based on centrality est. or number of selected tracks)
    Double_t                fIndexImpactParameter; // impact parameter of the collision in MC studies
    Short_t                 fEventCounter; // event counter (used for local test runmode purpose)
    Short_t                 fNumEventsAnalyse; // [50] number of events to be analysed / after passing selection (only in test mode)
    Int_t                   fRunNumber; // [-1] run number obtained from AliVHeader
    Bool_t                  fExtraPileUp; // extra pile-up cuts
    // array lenghts & constants
    const Double_t          fPDGMassPion; // [DPGMass] DPG mass of charged pion
    const Double_t          fPDGMassKaon; // [DPGMass] DPG mass of charged kaon
    const Double_t          fPDGMassProton; // [DPGMass] DPG mass of proton
    
    // selected POIs containers
    std::vector<FlowPart>*  fVectorCharged; 
    std::vector<FlowPart>*  fVectorPion; 
    std::vector<FlowPart>*  fVectorKaon; 
    std::vector<FlowPart>*  fVectorProton; 
    
    //cuts & selection: analysis
    RunMode                 fRunMode; // running mode (not grid related)
    AnalType                fAnalType; // analysis type: AOD / ESD
    Bool_t                  fSampling;      // Do random sampling ? (estimation of vn stat. uncertanity)
    Bool_t                  fCutImpactParameter; // Cut on the impact parameter
    Bool_t                  fFillQA; //[kTRUE] flag for filling the QA plots
    Bool_t                  fProcessCharged; // flag for processing charged tracks (both RPF and POIs)
    Bool_t                  fProcessPID; // flag for processing PID tracks (pi,K,p)
    AliFlowBayesianPID*     fBayesianResponse; 
    
    // cuts & selection: flow related
    Float_t                 fCutFlowRFPsPtMin; // [0] (GeV/c) min pT threshold for RFPs particle for reference flow
    Float_t                 fCutFlowRFPsPtMax; // [0] (GeV/c) max pT threshold for RFPs particle for reference flow
    const Float_t           fFlowPOIsPtMin; // [0] (GeV/c) min pT treshold for POIs for differential flow
    const Float_t           fFlowPOIsPtMax; // [20] (GeV/c) max pT treshold for POIs for differential flow
    const Int_t             fFlowCentMin; // [0] min range for centrality/multiplicity histos
    const Int_t             fFlowCentMax; // [150] min range for centrality/multiplicity histos
    const Int_t             fFlowCentNumBins; // [150] min range for centrality/multiplicity histos
    Double_t          fImpactParameterMin; // [0] min range for impact parameter
    Double_t          fImpactParameterMax; // [50] max range for impact parameter
    Bool_t                  fCutFlowDoFourCorrelations; // [kTRUE] flag for processing <4>
    Bool_t                  fDoOnlyMixedCorrelations; // [kTRUE] flag if I only want to analyse mixed harmonics
    Bool_t                  fFlowFillWeights; //[kFALSE] flag for filling weights
    Bool_t                  fFlowUseNUAWeights; //[kFALSE] flag for using the previously filled NUA weights (NOTE: this is turned on only when path to file is applied via fFlowWeightsPath)
    Bool_t                  fFlowUseNUEWeights; //[kFALSE] flag for using the previously filled NUE weights (NOTE: this is turned on only when path to file is applied via fFlowWeightsPath)
    TString                 fFlowNUAWeightsPath;//[] path to source root file with weigthts (if empty unit weights are applied) "alice/cern.ch/user/n/nmohamma/CorrectionMaps/fb96/NUACorrectionMap.root"
    
    TString                 fFlowNUEWeightsPath; //[] path to source root file with weigthts (if empty unit weights are applied) e.g. "alice/cern.ch/user/k/kgajdoso/EfficienciesWeights/2016/PhiWeight_LHC16kl.root"
    Bool_t                  fPositivelyChargedRef; //for same charged reference particle studies
    Bool_t                  fNegativelyChargedRef; //for same charged reference particle studies
    Bool_t                  fPositivelyChargedPOI; //for like sign reference particles and POIs studies
    Bool_t                  fNegativelyChargedPOI; //for unlike sign reference particle and POIs studies
    
    
    //cuts & selection: events
    Float_t                    fPVtxCutZ; // (cm) PV z cut
    ColSystem                  fColSystem; // collision system: pp or PbPb both with Run2 data (2016 data for pp and 2015 for PbPb)
    TString                 fMultEstimator; // [''] multiplicity estimator (suported: ''/Charged,VOA,V0C,V0M,CL0,CL1,ZNA,ZNC)
    Short_t                 fTrigger; // physics selection trigger
    Bool_t                  fFullCentralityRange; // flag for running over the full centrality range (0-100%). Otherwise runs from 50-100%
    
    //cuts & selection: tracks
    Float_t                 fCutChargedEtaMax; // (-) Maximum pseudorapidity range
    Float_t                 fCutChargedPtMax; // (GeV/c) Maximal track pT
    Float_t                 fCutChargedPtMin; // (GeV/c) Minimal track pT
    Float_t                 fCutChargedDCAzMax; // (cm) Maximal DCA-z cuts for tracks (pile-up rejection suggested for LHC16)
    Float_t                 fCutChargedDCAxyMax; // (cm) Maximal DCA-xy cuts for tracks (pile-up rejection suggested for LHC16)
    UInt_t                  fCutChargedTrackFilterBit; // (-) tracks filter bit
    UShort_t                fCutChargedNumTPCclsMin;  // (-) Minimal number of TPC clusters used for track reconstruction
    Double_t                fMaxChi2perTPCcls; // max chi2 per TPC clusters
    
    // cuts & selection: PID selection
    AliESDpid               fESDpid; 
    Bool_t                  fCutPIDUseAntiProtonOnly; // [kFALSE] check proton PID charge to select AntiProtons only
    Bool_t                  fPIDnsigma; // [kTRUE] default pid method
    Int_t                   fPIDnsigmaCombination; // 1,2,3 for PID nsigma combinations 1,2,3
    Bool_t                  fPIDbayesian; // [kFALSE] bayesian pid method
    Double_t                fCutPIDnSigmaPionMax; // [3] maximum of nSigmas (TPC or TPC & TOF combined) for pion candidates
    Double_t                fCutPIDnSigmaKaonMax; // [3] maximum of nSigmas (TPC or TPC & TOF combined) for kaon candidates
    Double_t                fCutPIDnSigmaProtonMax; // [3] maximum of nSigmas (TPC or TPC & TOF combined) for proton candidates
    Double_t                fCutPIDnSigmaTPCRejectElectron; // [3] number of TPC nSigma for electron rejection
    Bool_t                  fCutPIDnSigmaCombinedNoTOFrejection; // [kFALSE] flag for rejection candidates in TPC+TOF pt region if TOF is not available (if true and no TOF, only TPC is used)
    Float_t                 fCurrCentr; // current centrality used for set the priors
    Double_t                fParticleProbability; // Minimum Bayesian probability
    
    static const Int_t      fNumCentralities = 5; // number of centrality ranges for NUE maps (for now...)
    static const Int_t      fgkPIDptBin = 32; // pT bins for priors
    Float_t                 fC[fgkPIDptBin][5],fBinLimitPID[fgkPIDptBin]; // pt bin limit and priors
    static const Short_t    fFlowNumHarmonicsMax = 10; // maximum harmonics length of flow vector array
    static const Short_t    fFlowNumWeightPowersMax = 10; // maximum weight power length of flow vector array
    static const Short_t    fFlowPOIsPtNumBins = 200; // number of pT bins for POIs
    
    static const Short_t    fiNumIndexQA = 2; // QA indexes: 0: before cuts // 1: after cuts
    
    const static Short_t    fNumSamples = 6; // overall number of samples (from random sampling) used
    const static Int_t      fNumHarmonics = 5; // number of harmonics
    const static Int_t      fNumMixedHarmonics = 4; //4 number of mixed harmonics: 4{psi2}, 6{psi3} and 5{psi2,3}
    static Int_t            fHarmonics[fNumHarmonics]; // values of used harmonics
    static Int_t            fMixedHarmonics[fNumMixedHarmonics]; // values of used harmonics
    const static Int_t      fNumEtaGap = 1; // number of harmonics
    static Double_t         fEtaGap[fNumEtaGap]; // values of used harmonics
    
    
    // output lists
    TList*      fQAEvents; 
    TList*      fQACharged; 
    TList*      fQAPID; 
    TList*      fFlowWeights; 
    TList*      fFlowRefs; 
    TList*      fFlowCharged; 
    TList*      fFlowPID; 
    
    // histograms & profiles
    
    // Flow
    TH3D*           fh3BeforeNUAWeightsRefs; 
    TH3D*           fh3BeforeNUAWeightsCharged; 
    TH3D*           fh3BeforeNUAWeightsPion; 
    TH3D*           fh3BeforeNUAWeightsKaon; 
    TH3D*           fh3BeforeNUAWeightsProton; 
    
    TH3D*           fh3AfterNUAWeightsRefs; 
    TH3D*           fh3AfterNUAWeightsCharged; 
    TH3D*           fh3AfterNUAWeightsPion; 
    TH3D*           fh3AfterNUAWeightsKaon; 
    TH3D*           fh3AfterNUAWeightsProton; 
    
    TH1D*           fhBeforeNUEWeightsRefs; 
    TH1D*           fhBeforeNUEWeightsCharged; 
    TH1D*           fhBeforeNUEWeightsPion; 
    TH1D*           fhBeforeNUEWeightsKaon; 
    TH1D*           fhBeforeNUEWeightsProton; 
    
    TH1D*           fhAfterNUEWeightsRefs; 
    TH1D*           fhAfterNUEWeightsCharged; 
    TH1D*           fhAfterNUEWeightsPion; 
    TH1D*           fhAfterNUEWeightsKaon; 
    TH1D*           fhAfterNUEWeightsProton; 
    
    TH3D*           fh3NUAWeightRefsPlus; 
    TH3D*           fh3NUAWeightChargedPlus; 
    TH3D*           fh3NUAWeightPionPlus; 
    TH3D*           fh3NUAWeightKaonPlus; 
    TH3D*           fh3NUAWeightProtonPlus; 
    
    TH3D*           fh3NUAWeightRefsMinus; 
    TH3D*           fh3NUAWeightChargedMinus; 
    TH3D*           fh3NUAWeightPionMinus; 
    TH3D*           fh3NUAWeightKaonMinus; 
    TH3D*           fh3NUAWeightProtonMinus; 
    
    // Events
    TH2D*           fhEventSampling; 
    TH1D*           fhEventCentrality; 
    TH2D*           fh2EventCentralityNumSelCharged; 
    TH1D*           fhEventCounter; 
    // Charged
    TH2D*           fh2RefsMult; 
    TH2D*           fh2RefsPt; 
    TH2D*           fh2RefsEta; 
    TH2D*           fh2RefsPhi; 
    TH1D*           fhChargedCounter; 
    // PID
    TH2D*           fh2PIDPionMult; 
    TH2D*           fh2PIDPionPt; 
    TH2D*           fh2PIDPionPhi; 
    TH2D*           fh2PIDPionEta; 
    TH1D*           fhPIDPionCharge; 
    
    TH2D*           fh2PIDKaonMult; 
    TH2D*           fh2PIDKaonPt; 
    TH2D*           fh2PIDKaonPhi; 
    TH2D*           fh2PIDKaonEta; 
    TH1D*           fhPIDKaonCharge; 
    
    TH2D*           fh2PIDProtonMult; 
    TH2D*           fh2PIDProtonPt; 
    TH2D*           fh2PIDProtonPhi; 
    TH2D*           fh2PIDProtonEta; 
    TH1D*           fhPIDProtonCharge; 
    
    TH2D*           fh2PIDPionTPCdEdx; 
    TH2D*           fh2PIDPionTOFbeta; 
    
    TH2D*           fh2PIDKaonTPCdEdx; 
    TH2D*           fh2PIDKaonTOFbeta; 
    
    TH2D*           fh2PIDProtonTPCdEdx; 
    TH2D*           fh2PIDProtonTOFbeta; 
    
    TH2D*           fh2PIDPionTPCnSigmaPion; 
    TH2D*           fh2PIDPionTOFnSigmaPion; 
    TH2D*           fh2PIDPionTPCnSigmaKaon; 
    TH2D*           fh2PIDPionTOFnSigmaKaon; 
    TH2D*           fh2PIDPionTPCnSigmaProton; 
    TH2D*           fh2PIDPionTOFnSigmaProton; 
    TH2D*           fh2PIDKaonTPCnSigmaPion; 
    TH2D*           fh2PIDKaonTOFnSigmaPion; 
    TH2D*           fh2PIDKaonTPCnSigmaKaon; 
    TH2D*           fh2PIDKaonTOFnSigmaKaon; 
    TH2D*           fh2PIDKaonTPCnSigmaProton; 
    TH2D*           fh2PIDKaonTOFnSigmaProton; 
    TH2D*           fh2PIDProtonTPCnSigmaPion; 
    TH2D*           fh2PIDProtonTOFnSigmaPion; 
    TH2D*           fh2PIDProtonTPCnSigmaKaon; 
    TH2D*           fh2PIDProtonTOFnSigmaKaon; 
    TH2D*           fh2PIDProtonTPCnSigmaProton; 
    TH2D*           fh2PIDProtonTOFnSigmaProton; 
    
    
    TH1F*           fhNUEWeightRefsPlus[fNumCentralities]; 
    TH1F*           fhNUEWeightChargedPlus[fNumCentralities]; 
    TH1F*           fhNUEWeightPionPlus[fNumCentralities]; 
    TH1F*           fhNUEWeightKaonPlus[fNumCentralities]; 
    TH1F*           fhNUEWeightProtonPlus[fNumCentralities]; 
    
    TH1F*           fhNUEWeightRefsMinus[fNumCentralities]; 
    TH1F*           fhNUEWeightChargedMinus[fNumCentralities]; 
    TH1F*           fhNUEWeightPionMinus[fNumCentralities]; 
    TH1F*           fhNUEWeightKaonMinus[fNumCentralities]; 
    TH1F*           fhNUEWeightProtonMinus[fNumCentralities]; 
    
    
    TComplex                fFlowVecQpos[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax]; // flow vector array for flow calculation
    TComplex                fFlowVecQneg[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax]; // flow vector array for flow calculation
    TComplex                fFlowVecPpos[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax][fFlowPOIsPtNumBins]; // flow vector array for flow calculation
    TComplex                fFlowVecPneg[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax][fFlowPOIsPtNumBins]; // flow vector array for flow calculation
    TComplex                fFlowVecS[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax][fFlowPOIsPtNumBins]; // flow vector array for flow calculation
    
    
    // histograms & profiles
    
    TProfile*       fpMeanQxRefsPos[fNumEtaGap][fNumHarmonics]; 
    TProfile*       fpMeanQxRefsNeg[fNumEtaGap][fNumHarmonics]; 
    TProfile*       fpMeanQyRefsPos[fNumEtaGap][fNumHarmonics]; 
    TProfile*       fpMeanQyRefsNeg[fNumEtaGap][fNumHarmonics]; 
    
    TProfile*       fpRefsCor2[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    TProfile2D*     fp2ChargedCor2Pos[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    TProfile2D*     fp2ChargedCor2Neg[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    TProfile2D*     fp2PionCor2Pos[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    TProfile2D*     fp2PionCor2Neg[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    TProfile2D*     fp2KaonCor2Pos[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    TProfile2D*     fp2KaonCor2Neg[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    TProfile2D*     fp2ProtonCor2Pos[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    TProfile2D*     fp2ProtonCor2Neg[fNumSamples][fNumEtaGap][fNumHarmonics]; 
    
    //TProfile*       fpRefsCor4[fNumSamples][fNumHarmonics]; //! <4> correlations for RFPs
    //TProfile2D*     fp2ChargedCor4[fNumSamples][fNumHarmonics]; //! <4'> correlations for Charged tracks POIs
    //TProfile2D*     fp2PionCor4[fNumSamples][fNumHarmonics]; //! <4'> correlations for pion POIs
    //TProfile2D*     fp2KaonCor4[fNumSamples][fNumHarmonics]; //! <4'> correlations for kaon POIs
    //TProfile2D*     fp2ProtonCor4[fNumSamples][fNumHarmonics]; //! <4'> correlations for proton POIs
    
    //Mixed harmonics:
    TProfile*       fpMixedRefsCor4[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile*       fpMixedRefsCor6[fNumSamples][fNumEtaGap]; 
    TProfile2D*     fpMixedChargedCor3Pos[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile2D*     fpMixedChargedCor3Neg[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile2D*     fpMixedPionCor3Pos[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile2D*     fpMixedPionCor3Neg[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile2D*     fpMixedKaonCor3Pos[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile2D*     fpMixedKaonCor3Neg[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile2D*     fpMixedProtonCor3Pos[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile2D*     fpMixedProtonCor3Neg[fNumSamples][fNumEtaGap][fNumMixedHarmonics]; 
    TProfile2D*     fpMixedChargedCor4Pos[fNumSamples][fNumEtaGap]; 
    TProfile2D*     fpMixedChargedCor4Neg[fNumSamples][fNumEtaGap]; 
    TProfile2D*     fpMixedPionCor4Pos[fNumSamples][fNumEtaGap]; 
    TProfile2D*     fpMixedPionCor4Neg[fNumSamples][fNumEtaGap]; 
    TProfile2D*     fpMixedKaonCor4Pos[fNumSamples][fNumEtaGap]; 
    TProfile2D*     fpMixedKaonCor4Neg[fNumSamples][fNumEtaGap]; 
    TProfile2D*     fpMixedProtonCor4Pos[fNumSamples][fNumEtaGap]; 
    TProfile2D*     fpMixedProtonCor4Neg[fNumSamples][fNumEtaGap]; 
    
    // QA: events
    TH1D*           fhQAEventsPVz[fiNumIndexQA]; 
    TH1D*           fhQAEventsNumContrPV[fiNumIndexQA]; 
    TH1D*           fhQAEventsNumSPDContrPV[fiNumIndexQA]; 
    TH1D*           fhQAEventsDistPVSPD[fiNumIndexQA]; 
    TH1D*           fhQAEventsSPDresol[fiNumIndexQA]; 
    TH2D*           fhQAEventsCentralityOutliers[fiNumIndexQA]; 
    TH2D*           fhQAEventsPileUp[fiNumIndexQA]; 
    TH2D*           fhEventsMultTOFFilterbit32[fiNumIndexQA]; 
    // QA: charged tracks
    TH1D*           fhQAChargedMult[fiNumIndexQA];       
    TH1D*           fhQAChargedPt[fiNumIndexQA];         
    TH1D*           fhQAChargedEta[fiNumIndexQA];        
    TH1D*           fhQAChargedPhi[fiNumIndexQA];        
    TH1D*           fhQAChargedCharge[fiNumIndexQA];     
    TH1D*           fhQAChargedFilterBit[fiNumIndexQA];  
    TH1D*           fhQAChargedNumTPCcls[fiNumIndexQA];  
    TH1D*           fhQAChargedDCAxy[fiNumIndexQA];      
    TH1D*           fhQAChargedDCAz[fiNumIndexQA];       
    // QA: PID tracks
    TH1D*           fhQAPIDTPCstatus[fiNumIndexQA];  
    TH1D*           fhQAPIDTOFstatus[fiNumIndexQA];  
    TH2D*           fhQAPIDTPCdEdx[fiNumIndexQA];    
    TH2D*           fhQAPIDTOFbeta[fiNumIndexQA];    
    
    TH3D*           fh3PIDPionTPCTOFnSigmaPion[fiNumIndexQA]; 
    TH3D*           fh3PIDPionTPCTOFnSigmaKaon[fiNumIndexQA]; 
    TH3D*           fh3PIDPionTPCTOFnSigmaProton[fiNumIndexQA]; 
    
    TH3D*           fh3PIDKaonTPCTOFnSigmaPion[fiNumIndexQA]; 
    TH3D*           fh3PIDKaonTPCTOFnSigmaKaon[fiNumIndexQA]; 
    TH3D*           fh3PIDKaonTPCTOFnSigmaProton[fiNumIndexQA]; 
    
    TH3D*           fh3PIDProtonTPCTOFnSigmaPion[fiNumIndexQA]; 
    TH3D*           fh3PIDProtonTPCTOFnSigmaKaon[fiNumIndexQA]; 
    TH3D*           fh3PIDProtonTPCTOFnSigmaProton[fiNumIndexQA]; 
    
    
    // functions:
    Bool_t                  InitializeTask(); // called once on beginning of task (within CreateUserObjects method)
    void                    ListParameters(); // list all task parameters
    
    Bool_t                  EventSelection(); // main method for event selection (specific event selection is applied within)
    Bool_t                  IsEventSelected_PbPb(); // event selection for LHC2015 PbPb data
    Bool_t                  IsEventSelected_pp(); // event selection for LHC2016 MB pp data
    void                    FillEventsQA(const Short_t iQAindex); // filling QA plots related to event selection
    Short_t                 GetSamplingIndex(); // returns sampling index based on sampling selection (number of samples)
    Short_t                 GetCentralityIndex(); // returns centrality index based centrality estimator or number of selected tracks
    Double_t                GetImpactParameter(); // returns the impact parameter or the number of selected tracks
    Double_t                GetWDist(const AliAODVertex* v0, const AliAODVertex* v1); // gets the distance between the two vertices
    Bool_t                  ProcessEvent(); // main (envelope) method for processing events passing selection
    
    Bool_t                  Filtering(); // main (envelope) method for filtering all POIs in event
    void                    FilterCharged(); // charged tracks filtering
    Bool_t                  IsChargedSelected(const AliAODTrack* track = 0x0); // charged track selection
    void                    FillQARefs(const Short_t iQAindex, const AliAODTrack* track = 0x0); // filling QA plots for RFPs selection
    void                    FillQARefs(const Short_t iQAindex, const AliMCParticle* particle = 0x0); // filling QA plots for RFPs selection

    void                    FillQACharged(const Short_t iQAindex, const AliAODTrack* track = 0x0); // filling QA plots for charged track selection
    void                    FillQACharged(const Short_t iQAindex, const AliMCParticle* particle = 0x0); // filling QA plots for charged track selection

    void                    FilterPID(); // pi,K,p filtering
    PartSpecies             IsPIDSelected(const AliAODTrack* track); // PID tracks selections
    void                    FillPIDQA(const Short_t iQAindex, const AliAODTrack* track = 0x0, const PartSpecies species = kUnknown); // filling pi,K,p QA histograms
    void                    FillPIDQA(const Short_t iQAindex, const AliMCParticle* particle = 0x0, const PartSpecies species = kUnknown); // filling pi,K,p QA histograms

    // Flow related methods
    void                    DoFlowRefs(const Short_t iEtaGapIndex = 0); // Estimate <2> for reference flow
    void                    DoFlowCharged(const Short_t iEtaGapIndex = 0); // Estimate <2'> for pt diff. flow of charged hadrons
    void                    DoFlowPID(const Short_t iEtaGapIndex = 0, const PartSpecies species = kUnknown); // Estimate <2'> for pt diff. flow of PID (pi,K,p) hadrons
    void                    FillRefsVectors(const Short_t iEtaGapIndex = 0); // fill flow vector Q with RFPs for reference flow
    void                    FillPOIsVectors(const Short_t iEtaGapIndex = 0, const PartSpecies species = kUnknown, const Short_t iMassIndex = 0); // fill flow vectors p,q and s with POIs (for given species) for differential flow calculations
    Short_t                 GetPOIsPtBinIndex(const Double_t pt); // return pT bin index based on momenta value
    void                    ResetRFPsVector(TComplex (&array)[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax]); // set values to TComplex(0,0,0) for given array
    void                    ResetPOIsVector(TComplex (&array)[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax][fFlowPOIsPtNumBins]); // set values to TComplex(0,0,0) for given array
    void                    ListFlowVector(TComplex (&array)[fFlowNumHarmonicsMax][fFlowNumWeightPowersMax]); // printf all values of given Flow vector array
    
    TComplex                Q(const Short_t n, const Short_t p);
    TComplex                QGapPos(const Short_t n, const Short_t p);
    TComplex                QGapNeg(const Short_t n, const Short_t p);
    TComplex                P(const Short_t n, const Short_t p, const Short_t pt);
    TComplex                PGapPos(const Short_t n, const Short_t p, const Short_t pt);
    TComplex                PGapNeg(const Short_t n, const Short_t p, const Short_t pt);
    TComplex                S(const Short_t n, const Short_t p, const Short_t pt);
    
    TComplex                Two(const Short_t n1, const Short_t n2); // Two particle reference correlation calculations (no eta gap)
    TComplex                TwoGap(const Short_t n1, const Short_t n2); // Two particle reference correlation calculations (with eta gap)
    TComplex                TwoDiff(const Short_t n1, const Short_t n2, const Short_t pt); // Two particle diff. correlation calculations (no eta gap)
    TComplex                TwoDiffGapPos(const Short_t n1, const Short_t n2, const Short_t pt); // Two particle diff. correlation calculations (with eta gap)
    TComplex                TwoDiffGapNeg(const Short_t n1, const Short_t n2, const Short_t pt); // Two particle diff. correlation calculations (with eta gap)
    
    TComplex                ThreeDiffGapPos(const Short_t n1, const Short_t n2, const Short_t n3,const Short_t pt); // Three particle diff. correlation calculations (with eta gap);
    TComplex                ThreeDiffGapNeg(const Short_t n1, const Short_t n2, const Short_t n3,const Short_t pt); // Three particle diff. correlation calculations (with eta gap);
    
    TComplex                Four(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4); // Four particle reference correlation calculations (no eta gap)
    TComplex                FourDiff(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t pt); // Four particle reference correlation calculations (no eta gap)
    TComplex                FourGapPos(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4); // Four particle reference correlation calculations (with eta gap); n1+n2 = n3+n4;  n1, n2 from P & n3, n4 from M
    TComplex                FourGapNeg(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4); // Four particle reference correlation calculations (with eta gap); n1+n2 = n3+n4;  n1, n2 from M & n3, n4 from P
    TComplex                Four13DiffGapPos(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t pt); // Four particle diff. correlation calculations (with eta gap); n1 = n2+n3+n4;  n1 from P & n2, n3, n4 from M
    TComplex                Four13DiffGapNeg(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t pt); // Four particle diff. correlation calculations (with eta gap); n1 = n2+n3+n4;  n1 from M & n2, n3, n4 from P
    TComplex                SixGapPos(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t n5, const Short_t n6); // Six particle reference correlation calculations (with eta gap); n1 + n2 + n3 = n4 + n5 + n6 ;n1, n2, n3 from P and n4, n5 and n6 from M
    TComplex                SixGapNeg(const Short_t n1, const Short_t n2, const Short_t n3, const Short_t n4, const Short_t n5, const Short_t n6); // Six particle reference correlation calculations (with eta gap); n1 + n2 + n3 = n4 + n5 + n6; n1, n2, n3 from M and n4, n5 and n6 from P
    
    
    AliAnalysisTaskFlowModesMC(const AliAnalysisTaskFlowModesMC&); // not implemented
    AliAnalysisTaskFlowModesMC& operator=(const AliAnalysisTaskFlowModesMC&); // not implemented
    
    ClassDef(AliAnalysisTaskFlowModesMC, 4);
};

#endif