reduce memory usage
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 May 2012 10:25:20 +0000 (10:25 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 May 2012 10:25:20 +0000 (10:25 +0000)
PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx
PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.h
PWGCF/FLOW/macros/AddTaskPhiFlow.C

index 97567e7..d78bafa 100644 (file)
@@ -25,6 +25,7 @@
 #include "TH2F.h"
 #include "TCanvas.h"
 #include "TMath.h"
+#include "TObjArray.h"
 #include "AliAnalysisTask.h"
 #include "AliAnalysisTaskSE.h"
 #include "AliAnalysisManager.h"
 #include "AliESDVZERO.h"
 #include "AliAODVZERO.h"
 #include "AliPID.h"
-#include "AliESDpid.h"
-#include "AliAODpidUtil.h"
+#include "AliPIDResponse.h"
 
 class AliFlowTrackCuts;
 
 ClassImp(AliAnalysisTaskPhiFlow)
 
 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
-   fAODAnalysis(0), fEtaMinA(0), fEtaMaxA(0), fEtaMinB(0), fEtaMaxB(0), fCutsRP(NULL), fNullCuts(0), fBayesianResponse(0), fOldTrackParam(0), fRequireTPCStandAlone(0), fStrictKaonCuts(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fParticleProbability(0.5), fCentrality(0), fESD(0), fAOD(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDDeltaDip(0), fInvMNP03(0), fInvMPP03(0), fInvMNN03(0), fInvMNP36(0), fInvMPP36(0), fInvMNN36(0), fInvMNP69(0), fInvMPP69(0), fInvMNN69(0), fInvMNP912(0), fInvMPP912(0), fInvMNN912(0), fInvMNP1215(0), fInvMPP1215(0), fInvMNN1215(0), fInvMNP1518(0), fInvMPP1518(0), fInvMNN1518(0), fInvMNP1821(0), fInvMPP1821(0), fInvMNN1821(0), fInvMNP2124(0), fInvMPP2124(0), fInvMNN2124(0), fInvMNP2427(0), fInvMPP2427(0), fInvMNN2427(0), fInvMNP2730(0), fInvMPP2730(0), fInvMNN2730(0), fInvMNP3035(0), fInvMPP3035(0), fInvMNN3035(0), fInvMNP3540(0), fInvMPP3540(0), fInvMNN3540(0), fInvMNP4045(0), fInvMPP4045(0), fInvMNN4045(0), fInvMNP4550(0), fInvMPP4550(0), fInvMNN4550(0), fInvMNP5055(0), fInvMPP5055(0), fInvMNN5055(0), fInvMNP5560(0), fInvMPP5560(0), fInvMNN5560(0), fInvMNP6065(0), fInvMPP6065(0), fInvMNN6065(0), fInvMNP6570(0), fInvMPP6570(0), fInvMNN6570(0), fDeltaPhiPsiNP03(0), fDeltaPhiPsiNP36(0), fDeltaPhiPsiNP69(0), fDeltaPhiPsiNP912(0), fDeltaPhiPsiNP1215(0), fDeltaPhiPsiNP1518(0), fDeltaPhiPsiNP1821(0), fDeltaPhiPsiNP2124(0), fDeltaPhiPsiNP2427(0), fDeltaPhiPsiNP2730(0), fDeltaPhiPsiNP3035(0), fDeltaPhiPsiNP3540(0), fDeltaPhiPsiNP4045(0), fDeltaPhiPsiNP4550(0), fDeltaPhiPsiNP5055(0), fDeltaPhiPsiNP5560(0), fDeltaPhiPsiNP6065(0), fDeltaPhiPsiNP6570(0), fProfV2(0), fProfV2Sin(0), fProfV2InvM03(0), fProfV2InvM36(0), fProfV2InvM69(0), fProfV2InvM912(0), fProfV2InvM1215(0), fProfV2InvM1518(0), fProfV2InvM1821(0), fProfV2InvM2124(0), fProfV2InvM2427(0), fProfV2InvM2730(0), fProfV2InvM3035(0), fProfV2InvM3540(0), fProfV2InvM4045(0), fProfV2InvM4550(0), fProfV2InvM5055(0), fProfV2InvM5560(0), fProfV2InvM6065(0), fProfV2InvM6570(0), fProfV2SinInvM03(0), fProfV2SinInvM36(0), fProfV2SinInvM69(0), fProfV2SinInvM912(0), fProfV2SinInvM1215(0), fProfV2SinInvM1518(0), fProfV2SinInvM1821(0), fProfV2SinInvM2124(0), fProfV2SinInvM2427(0), fProfV2SinInvM2730(0), fProfV2SinInvM3035(0), fProfV2SinInvM3540(0), fProfV2SinInvM4045(0), fProfV2SinInvM4550(0), fProfV2SinInvM5055(0), fProfV2SinInvM5560(0), fProfV2SinInvM6065(0), fProfV2SinInvM6570(0), fPtSpectra03(0), fPtSpectra36(0), fPtSpectra69(0), fPtSpectra912(0), fPtSpectra1215(0), fPtSpectra1518(0), fPtSpectra1821(0), fPtSpectra2124(0), fPtSpectra2427(0), fPtSpectra2730(0), fPtSpectra3035(0), fPtSpectra3540(0), fPtSpectra4045(0), fPtSpectra4550(0), fPtSpectra5055(0), fPtSpectra5560(0), fPtSpectra6065(0), fPtSpectra6570(0), fEventPlaneSTAR(0), fEventPlaneResolutionRandom(0), fEventPlaneResolutionEta(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fRPCuts(0), fPOICuts(0), fVertexRange(0), fQx(0), fQy(0), fEventPlane(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fPIDtype(kCombined), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fPairLoss(0), fEventPlanePtCut(0), fSetEventPlanePtCut(0)
+   fDebug(0), fAODAnalysis(0), fMassBins(0), fMinMass(0), fMaxMass(0), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fOldTrackParam(0), fRequireTPCStandAlone(0), fStrictKaonCuts(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCentrality(0), fESD(0), fAOD(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fInvMNP03(0), fInvMPP03(0), fInvMNN03(0), fInvMNP36(0), fInvMPP36(0), fInvMNN36(0), fInvMNP69(0), fInvMPP69(0), fInvMNN69(0), fInvMNP912(0), fInvMPP912(0), fInvMNN912(0), fInvMNP1215(0), fInvMPP1215(0), fInvMNN1215(0), fInvMNP1518(0), fInvMPP1518(0), fInvMNN1518(0), fInvMNP1821(0), fInvMPP1821(0), fInvMNN1821(0), fInvMNP2124(0), fInvMPP2124(0), fInvMNN2124(0), fInvMNP2427(0), fInvMPP2427(0), fInvMNN2427(0), fInvMNP2730(0), fInvMPP2730(0), fInvMNN2730(0), fInvMNP3035(0), fInvMPP3035(0), fInvMNN3035(0), fInvMNP3540(0), fInvMPP3540(0), fInvMNN3540(0), fInvMNP4045(0), fInvMPP4045(0), fInvMNN4045(0), fInvMNP4550(0), fInvMPP4550(0), fInvMNN4550(0), fInvMNP5055(0), fInvMPP5055(0), fInvMNN5055(0), fInvMNP5560(0), fInvMPP5560(0), fInvMNN5560(0), fInvMNP6065(0), fInvMPP6065(0), fInvMNN6065(0), fInvMNP6570(0), fInvMPP6570(0), fInvMNN6570(0), fPtSpectra03(0), fPtSpectra36(0), fPtSpectra69(0), fPtSpectra912(0), fPtSpectra1215(0), fPtSpectra1518(0), fPtSpectra1821(0), fPtSpectra2124(0), fPtSpectra2427(0), fPtSpectra2730(0), fPtSpectra3035(0), fPtSpectra3540(0), fPtSpectra4045(0), fPtSpectra4550(0), fPtSpectra5055(0), fPtSpectra5560(0), fPtSpectra6065(0), fPtSpectra6570(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)
 {
-   // Dummy constructor
-   fNullCuts = new AliFlowTrackCuts("null_cuts");
-   for (Int_t i = 0; i < 2; i++)
-   {
-      for (Int_t j = 0; j < 30; j++)
-      {
-         fFlowBands[i][j] = 0;
-         if (i == 0)  fFlowEvent[j] = new AliFlowEvent(10000);
-      }
-   }
-   if (fPIDtype == kCombined) fBayesianResponse = new AliFlowBayesianPID();
+   // Default constructor
+   for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = 0.;
 }
 //_____________________________________________________________________________
 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow(const char *name) : AliAnalysisTaskSE(name),
-   fAODAnalysis(0), fEtaMinA(0), fEtaMaxA(0), fEtaMinB(0), fEtaMaxB(0), fCutsRP(NULL), fNullCuts(0), fBayesianResponse(0), fOldTrackParam(0), fRequireTPCStandAlone(0), fStrictKaonCuts(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fParticleProbability(0.5), fCentrality(0), fESD(0), fAOD(0),fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDDeltaDip(0), fInvMNP03(0), fInvMPP03(0), fInvMNN03(0), fInvMNP36(0), fInvMPP36(0), fInvMNN36(0), fInvMNP69(0), fInvMPP69(0), fInvMNN69(0), fInvMNP912(0), fInvMPP912(0), fInvMNN912(0), fInvMNP1215(0), fInvMPP1215(0), fInvMNN1215(0), fInvMNP1518(0), fInvMPP1518(0), fInvMNN1518(0), fInvMNP1821(0), fInvMPP1821(0), fInvMNN1821(0), fInvMNP2124(0), fInvMPP2124(0), fInvMNN2124(0), fInvMNP2427(0), fInvMPP2427(0), fInvMNN2427(0), fInvMNP2730(0), fInvMPP2730(0), fInvMNN2730(0), fInvMNP3035(0), fInvMPP3035(0), fInvMNN3035(0), fInvMNP3540(0), fInvMPP3540(0), fInvMNN3540(0), fInvMNP4045(0), fInvMPP4045(0), fInvMNN4045(0), fInvMNP4550(0), fInvMPP4550(0), fInvMNN4550(0), fInvMNP5055(0), fInvMPP5055(0), fInvMNN5055(0), fInvMNP5560(0), fInvMPP5560(0), fInvMNN5560(0), fInvMNP6065(0), fInvMPP6065(0), fInvMNN6065(0), fInvMNP6570(0), fInvMPP6570(0), fInvMNN6570(0), fDeltaPhiPsiNP03(0), fDeltaPhiPsiNP36(0), fDeltaPhiPsiNP69(0), fDeltaPhiPsiNP912(0), fDeltaPhiPsiNP1215(0), fDeltaPhiPsiNP1518(0), fDeltaPhiPsiNP1821(0), fDeltaPhiPsiNP2124(0), fDeltaPhiPsiNP2427(0), fDeltaPhiPsiNP2730(0), fDeltaPhiPsiNP3035(0), fDeltaPhiPsiNP3540(0), fDeltaPhiPsiNP4045(0), fDeltaPhiPsiNP4550(0), fDeltaPhiPsiNP5055(0), fDeltaPhiPsiNP5560(0), fDeltaPhiPsiNP6065(0), fDeltaPhiPsiNP6570(0), fProfV2(0), fProfV2Sin(0), fProfV2InvM03(0), fProfV2InvM36(0), fProfV2InvM69(0), fProfV2InvM912(0), fProfV2InvM1215(0), fProfV2InvM1518(0), fProfV2InvM1821(0), fProfV2InvM2124(0), fProfV2InvM2427(0), fProfV2InvM2730(0), fProfV2InvM3035(0), fProfV2InvM3540(0), fProfV2InvM4045(0), fProfV2InvM4550(0), fProfV2InvM5055(0), fProfV2InvM5560(0), fProfV2InvM6065(0), fProfV2InvM6570(0), fProfV2SinInvM03(0), fProfV2SinInvM36(0), fProfV2SinInvM69(0), fProfV2SinInvM912(0), fProfV2SinInvM1215(0), fProfV2SinInvM1518(0), fProfV2SinInvM1821(0), fProfV2SinInvM2124(0), fProfV2SinInvM2427(0), fProfV2SinInvM2730(0), fProfV2SinInvM3035(0), fProfV2SinInvM3540(0), fProfV2SinInvM4045(0), fProfV2SinInvM4550(0), fProfV2SinInvM5055(0), fProfV2SinInvM5560(0), fProfV2SinInvM6065(0), fProfV2SinInvM6570(0), fPtSpectra03(0), fPtSpectra36(0), fPtSpectra69(0), fPtSpectra912(0), fPtSpectra1215(0), fPtSpectra1518(0), fPtSpectra1821(0), fPtSpectra2124(0), fPtSpectra2427(0), fPtSpectra2730(0), fPtSpectra3035(0), fPtSpectra3540(0), fPtSpectra4045(0), fPtSpectra4550(0), fPtSpectra5055(0), fPtSpectra5560(0), fPtSpectra6065(0), fPtSpectra6570(0), fEventPlaneSTAR(0), fEventPlaneResolutionRandom(0), fEventPlaneResolutionEta(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fRPCuts(0), fPOICuts(0), fVertexRange(0), fQx(0), fQy(0), fEventPlane(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fPIDtype(kCombined), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fPairLoss(0), fEventPlanePtCut(0), fSetEventPlanePtCut(0)
+    fDebug(0), fAODAnalysis(0), fMassBins(0), fMinMass(0), fMaxMass(0), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fOldTrackParam(0), fRequireTPCStandAlone(0), fStrictKaonCuts(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCentrality(0), fESD(0), fAOD(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fInvMNP03(0), fInvMPP03(0), fInvMNN03(0), fInvMNP36(0), fInvMPP36(0), fInvMNN36(0), fInvMNP69(0), fInvMPP69(0), fInvMNN69(0), fInvMNP912(0), fInvMPP912(0), fInvMNN912(0), fInvMNP1215(0), fInvMPP1215(0), fInvMNN1215(0), fInvMNP1518(0), fInvMPP1518(0), fInvMNN1518(0), fInvMNP1821(0), fInvMPP1821(0), fInvMNN1821(0), fInvMNP2124(0), fInvMPP2124(0), fInvMNN2124(0), fInvMNP2427(0), fInvMPP2427(0), fInvMNN2427(0), fInvMNP2730(0), fInvMPP2730(0), fInvMNN2730(0), fInvMNP3035(0), fInvMPP3035(0), fInvMNN3035(0), fInvMNP3540(0), fInvMPP3540(0), fInvMNN3540(0), fInvMNP4045(0), fInvMPP4045(0), fInvMNN4045(0), fInvMNP4550(0), fInvMPP4550(0), fInvMNN4550(0), fInvMNP5055(0), fInvMPP5055(0), fInvMNN5055(0), fInvMNP5560(0), fInvMPP5560(0), fInvMNN5560(0), fInvMNP6065(0), fInvMPP6065(0), fInvMNN6065(0), fInvMNP6570(0), fInvMPP6570(0), fInvMNN6570(0), fPtSpectra03(0), fPtSpectra36(0), fPtSpectra69(0), fPtSpectra912(0), fPtSpectra1215(0), fPtSpectra1518(0), fPtSpectra1821(0), fPtSpectra2124(0), fPtSpectra2427(0), fPtSpectra2730(0), fPtSpectra3035(0), fPtSpectra3540(0), fPtSpectra4045(0), fPtSpectra4550(0), fPtSpectra5055(0), fPtSpectra5560(0), fPtSpectra6065(0), fPtSpectra6570(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)
 {
    // Constructor
-   fNullCuts = new AliFlowTrackCuts("null_cuts");
-   for (int i = 0; i < 2; i++)
-   {
-      for (int j = 0; j < 30; j++)
-      {
-         fFlowBands[i][j] = 0;
-         if (i == 0)  fFlowEvent[j] = new AliFlowEvent(10000);
-      }
-   }
-   if (fPIDtype == kCombined) fBayesianResponse = new AliFlowBayesianPID();
+   for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = 0.;
+
    DefineInput(0, TChain::Class());
    DefineOutput(1, TList::Class());
-   for (Int_t i = (0 + 2); i < (30 + 2); i++)
-   {
-      DefineOutput(i, AliFlowEventSimple::Class());
-   }
+   DefineOutput(2, AliFlowEventSimple::Class());
 }
 //_____________________________________________________________________________
 AliAnalysisTaskPhiFlow::~AliAnalysisTaskPhiFlow()
 {
    // Destructor
-   for (Int_t i = 0; i < 30; i++)
-   {
-      delete fFlowEvent[i];
-   }
-   delete fNullCuts;
-   if (fPIDtype == kCombined) delete fBayesianResponse;
+   if (fNullCuts) delete fNullCuts;
+   if (fOutputList) delete fOutputList;
+   if (fCandidates) delete fCandidates;
+   if (fFlowEvent) delete fFlowEvent;
+   if (fBayesianResponse) delete fBayesianResponse;
 }
 //_____________________________________________________________________________
 TH1F* AliAnalysisTaskPhiFlow::BookHistogram(const char* name)
@@ -118,18 +97,6 @@ TH1F* AliAnalysisTaskPhiFlow::BookHistogram(const char* name)
    return hist;
 }
 //_____________________________________________________________________________
-TH1F* AliAnalysisTaskPhiFlow::BookDPhiPsiHistogram(const char* name)
-{
-   // Return a pointer to a TH1 with predefined binning
-   TH1F *hist = new TH1F(Form("#phi - #Psi (%s)", name), Form("#phi - #Psi (%s)", name), 60, -7, 7);
-   hist->GetXaxis()->SetTitle("#phi - #Psi (rad)");
-   hist->GetYaxis()->SetTitle("No. of pairs");
-   hist->SetMarkerStyle(kFullCircle);
-   hist->Sumw2();
-   fOutputList->Add(hist);
-   return hist;
-}
-//_____________________________________________________________________________
 TH2F* AliAnalysisTaskPhiFlow::BookPIDHistogram(const char* name)
 {
    // Return a pointer to a TH2 with predefined binning
@@ -145,44 +112,14 @@ TH1F* AliAnalysisTaskPhiFlow::InitPtSpectraHistograms(Int_t i)
    // intialize p_t histograms for each p_t bin
    Double_t nmin(0);
    Double_t nmax(0);
-   (i < 11) ? nmin = 0.3 * (i-1) : nmin = 0.5 * (i-11) + 3;
-   (i < 11) ? nmax = 0.3 * (i-1) + 0.3 : nmax = 0.5 * (i-11) + 3.5;
+   (i < 11) ? nmin = 0.3 * (i - 1) : nmin = 0.5 * (i - 11) + 3;
+   (i < 11) ? nmax = 0.3 * (i - 1) + 0.3 : nmax = 0.5 * (i - 11) + 3.5;
    TH1F* hist = new TH1F(Form("%f p_{t} %f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 60, nmin, nmax);
    hist->GetXaxis()->SetTitle("p_{T} GeV / c");
    fOutputList->Add(hist);
    return hist;
 }
 //_____________________________________________________________________________
-TProfile* AliAnalysisTaskPhiFlow::BookV2Profile(const char* name, Bool_t pt, Bool_t cos)
-{
-   // Return a pointer to a v_2 TProfile with predefined binning
-   TProfile* profile = 0x0;
-   if (pt)
-   {
-      profile = new TProfile(name, Form("v_{c, 2}^{pair} (%s)", name), 100, 0., 4);
-      profile->GetXaxis()->SetTitle("p_{T} GeV / c");
-   }
-   if (!pt)
-   {
-      if (cos)
-      {
-         profile = new TProfile(Form("v_{c, 2}^{pair} (%s)", name), Form("v_{c, 2}^{pair} (%s)", name), 30, 0.99, 1.092);
-         profile->GetXaxis()->SetTitle("M_{INV} GeV / c^{2}");
-      }
-      if (!cos)
-      {
-         profile = new TProfile(Form("v_{s, 2}^{pair} (%s)", name), Form("v_{s, 2}^{pair} (%s)", name), 30, 0.99, 1.092);
-         profile->GetXaxis()->SetTitle("M_{INV} GeV / c^{2}");
-      }
-   }
-   profile->SetMarkerStyle(kOpenCross);
-   if (cos) profile->GetYaxis()->SetTitle("v_{2} <cos n[#phi_{pair} - #Psi_{RP}]>");
-   if (!cos) profile->GetYaxis()->SetTitle("v_{2} <sin n[#phi_{pair} - #Psi_{RP}]>");
-   profile->Sumw2();
-   fOutputList->Add(profile);
-   return profile;
-}
-//_____________________________________________________________________________
 TH1F* AliAnalysisTaskPhiFlow::BookPtHistogram(const char* name)
 {
    // Return a pointer to a p_T spectrum histogram
@@ -206,12 +143,6 @@ void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
    fOutputList->Add(fCentralityPass);
    fOutputList->Add(fCentralityNoPass);
 
-   fEventPlaneSTAR = new TH1F("fEventPlaneSTAR", "Event plane orientation", 200, -.2, 4);
-   fEventPlaneSTAR->GetXaxis()->SetTitle("#Psi_{2}");
-   fEventPlaneSTAR->GetYaxis()->SetTitle("Events");
-   fEventPlaneResolutionRandom = new TProfile("fEventPlaneResolutionRandom", "RP kaon pairs, Random", 1, 0.5, 1.5);
-   fEventPlaneResolutionEta = new TProfile("fEventPlaneResolutionEta", "RP kaon pairs, Eta", 1, 0.5, 1.5);
-
    fNOPID = BookPIDHistogram("TPC signal, all particles");
    fPIDk = BookPIDHistogram("TPC signal, kaons");
    fInvMNP03 = BookHistogram("NP, 0 < p_{T} < 0.3 GeV");
@@ -268,64 +199,6 @@ void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
    fInvMNP6570 = BookHistogram("NP, 6.5 < p_{T} < 7.0 GeV");
    fInvMPP6570 = BookHistogram("PP, 6.5 < p_{T} < 7.0 GeV");
    fInvMNN6570 = BookHistogram("NN, 6.5 < p_{T} < 7.0 GeV");
-   fProfV2 = BookV2Profile("v_{2} cos terms", kTRUE, kTRUE);
-   fProfV2Sin = BookV2Profile("v_{2} sin terms", kTRUE, kFALSE);
-   fProfV2InvM03 = BookV2Profile("0.0 < p_{T} < 0.3 GeV", kFALSE, kTRUE);
-   fProfV2InvM36 = BookV2Profile("0.3 < p_{T} < 0.6 GeV", kFALSE, kTRUE);
-   fProfV2InvM69 = BookV2Profile("0.6 < p_{T} < 0.9 GeV", kFALSE, kTRUE);
-   fProfV2InvM912 = BookV2Profile("0.9 < p_{T} < 1.2 GeV", kFALSE, kTRUE);
-   fProfV2InvM1215 = BookV2Profile("1.1 < p_{T} < 1.5 GeV", kFALSE, kTRUE);
-   fProfV2InvM1518 = BookV2Profile("1.5 < p_{T} < 1.8 GeV", kFALSE, kTRUE);
-   fProfV2InvM1821 = BookV2Profile("1.8 < p_{T} < 2.1 GeV", kFALSE, kTRUE);
-   fProfV2InvM2124 = BookV2Profile("2.1 < p_{T} < 2.4 GeV", kFALSE, kTRUE);
-   fProfV2InvM2427 = BookV2Profile("2.4 < p_{T} < 2.7 GeV", kFALSE, kTRUE);
-   fProfV2InvM2730 = BookV2Profile("2.7 < p_{T} < 3.0 GeV", kFALSE, kTRUE);
-   fProfV2InvM3035 = BookV2Profile("3.0 < p_{T} < 3.5 GeV", kFALSE, kTRUE);
-   fProfV2InvM3540 = BookV2Profile("3.5 < p_{T} < 4.0 GeV", kFALSE, kTRUE);
-   fProfV2InvM4045 = BookV2Profile("4.0 < p_{T} < 4.5 GeV", kFALSE, kTRUE);
-   fProfV2InvM4550 = BookV2Profile("4.5 < p_{T} < 5.0 GeV", kFALSE, kTRUE);
-   fProfV2InvM5055 = BookV2Profile("5.0 < p_{T} < 5.5 GeV", kFALSE, kTRUE);
-   fProfV2InvM5560 = BookV2Profile("5.5 < p_{T} < 6.0 GeV", kFALSE, kTRUE);
-   fProfV2InvM6065 = BookV2Profile("6.0 < p_{T} < 6.5 GeV", kFALSE, kTRUE);
-   fProfV2InvM6570 = BookV2Profile("6.5 < p_{T} < 7.0 GeV", kFALSE, kTRUE);
-   fProfV2SinInvM03 = BookV2Profile("0.0 < p_{T} < 0.3 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM36 = BookV2Profile("0.3 < p_{T} < 0.6 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM69 = BookV2Profile("0.6 < p_{T} < 0.9 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM912 = BookV2Profile("0.9 < p_{T} < 1.2 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM1215 = BookV2Profile("1.2 < p_{T} < 1.5 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM1518 = BookV2Profile("1.5 < p_{T} < 1.8 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM1821 = BookV2Profile("1.8 < p_{T} < 2.1 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM2124 = BookV2Profile("2.1 < p_{T} < 2.4 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM2427 = BookV2Profile("2.4 < p_{T} < 2.7 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM2730 = BookV2Profile("2.5 < p_{T} < 3.0 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM3035 = BookV2Profile("3.0 < p_{T} < 3.5 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM3540 = BookV2Profile("3.5 < p_{T} < 4.0 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM4045 = BookV2Profile("4.0 < p_{T} < 4.5 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM4550 = BookV2Profile("4.5 < p_{T} < 5.0 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM5055 = BookV2Profile("5.0 < p_{T} < 5.5 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM5560 = BookV2Profile("5.5 < p_{T} < 6.0 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM6065 = BookV2Profile("6.0 < p_{T} < 6.5 GeV", kFALSE, kFALSE);
-   fProfV2SinInvM6570 = BookV2Profile("6.5 < p_{T} < 7.0 GeV", kFALSE, kFALSE);
-
-   fDeltaPhiPsiNP03 = BookDPhiPsiHistogram("NP, 0 < p_{T} < 0.3 GeV");
-   fDeltaPhiPsiNP36 = BookDPhiPsiHistogram("NP, 0.3 < p_{T} < 0.6 GeV");
-   fDeltaPhiPsiNP69 = BookDPhiPsiHistogram("NP, 0.6 < p_{T} < 0.9 GeV");
-   fDeltaPhiPsiNP912 = BookDPhiPsiHistogram("NP, 0.9 < p_{T} < 1.2 GeV");
-   fDeltaPhiPsiNP1215 = BookDPhiPsiHistogram("NP, 1.2 < p_{T} < 1.5 GeV");
-   fDeltaPhiPsiNP1518 = BookDPhiPsiHistogram("NP, 1.5 < p_{T} < 1.8 GeV");
-   fDeltaPhiPsiNP1821 = BookDPhiPsiHistogram("NP, 1.8 < p_{T} < 2.1 GeV");
-   fDeltaPhiPsiNP2124 = BookDPhiPsiHistogram("NP, 2.1 < p_{T} < 2.4 GeV");
-   fDeltaPhiPsiNP2427 = BookDPhiPsiHistogram("NP, 2.4 < p_{T} < 2.7 GeV");
-   fDeltaPhiPsiNP2730 = BookDPhiPsiHistogram("NP, 2.7 < p_{T} < 3.0 GeV");
-   fDeltaPhiPsiNP3035 = BookDPhiPsiHistogram("NP, 3.0 < p_{T} < 3.5 GeV");
-   fDeltaPhiPsiNP3540 = BookDPhiPsiHistogram("NP, 3.5 < p_{T} < 4.0 GeV");
-   fDeltaPhiPsiNP4045 = BookDPhiPsiHistogram("NP, 4.0 < p_{T} < 4.5 GeV");
-   fDeltaPhiPsiNP4550 = BookDPhiPsiHistogram("NP, 4.5 < p_{T} < 5.0 GeV");
-   fDeltaPhiPsiNP5055 = BookDPhiPsiHistogram("NP, 5.0 < p_{T} < 5.5 GeV");
-   fDeltaPhiPsiNP5560 = BookDPhiPsiHistogram("NP, 5.5 < p_{T} < 6.0 GeV");
-   fDeltaPhiPsiNP6065 = BookDPhiPsiHistogram("NP, 6.0 < p_{T} < 6.5 GeV");
-   fDeltaPhiPsiNP6570 = BookDPhiPsiHistogram("NP, 6.5 < p_{T} < 7.0 GeV");
-
 
    fPtSpectra03 = InitPtSpectraHistograms(1);
    fPtSpectra36 = InitPtSpectraHistograms(2);
@@ -346,17 +219,11 @@ void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
    fPtSpectra6065 = InitPtSpectraHistograms(17);
    fPtSpectra6570 = InitPtSpectraHistograms(18);
 
-   fOutputList->Add(fEventPlaneSTAR);
    fPtP = BookPtHistogram("i^{+}");
    fPtN = BookPtHistogram("i^{-}");
    fPtKP = BookPtHistogram("K^{+}");
    fPtKN = BookPtHistogram("K^{-}");
 
-   fEventPlane = new TH1F("fEventPlane", "Event plane", 200, -0.2, 4.);
-   fEventPlane->GetXaxis()->SetTitle("#Psi_{2}");
-   fEventPlane->GetYaxis()->SetTitle("Events");
-   fOutputList->Add(fEventPlane);
-
    fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
    fOutputList->Add(fPhi);
 
@@ -374,20 +241,16 @@ void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
 
    fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
    fOutputList->Add(fTPCM);
-
-   fPairLoss = new TH2F("fPairLoss", "DphiDeta", 100, -0.2, 0.2, 100, -0.2, 0.2);
-   fPairLoss->GetXaxis()->SetTitle("#Delta #eta");
-   fPairLoss->GetYaxis()->SetTitle("#Delta #phi");
-
-   fOutputList->Add(fEventPlaneResolutionRandom);
-   fOutputList->Add(fEventPlaneResolutionEta);
-   fOutputList->Add(fPairLoss);
-   fPIDDeltaDip = BookPIDHistogram("TPC signal after delta dip cut");
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
 {
    // Create user defined output objects
+   fNullCuts = new AliFlowTrackCuts("null_cuts");
+   fBayesianResponse = new AliFlowBayesianPID();
+   Double_t t(0);
+   for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
+   if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis!!!");
 
    AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
    cc->SetNbinsMult(10000);
@@ -410,32 +273,41 @@ void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
    cc->SetQMin(0.0);
    cc->SetQMax(3.0);
 
+   cc->SetNbinsMass(fMassBins);
+   cc->SetMassMin(fMinMass);
+   cc->SetMassMax(fMaxMass);
+
    // setup initial state of PID response object
    if (!fOldTrackParam) fBayesianResponse->SetNewTrackParam();
 
+   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
+   if (man)
+   {
+      AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
+      if (inputHandler)   fPIDResponse = inputHandler->GetPIDResponse();
+   }
+
    // Create all output objects and store them to a list
    fOutputList = new TList();
-   fOutputList->SetOwner();
+   fOutputList->SetOwner(kTRUE);
 
    // Create and post the Phi identification output objects
    AddPhiIdentificationOutputObjects();
    PostData(1, fOutputList);
 
-   // post flow events (see ctor)
-   for (Int_t i = 0; i < 30; i++) PostData(2 + i, fFlowEvent[i]);
-}
-//_____________________________________________________________________________
-template <typename T> void AliAnalysisTaskPhiFlow::PairLoss(const T* track1, const T* track2) const
-{
-   if (track1->Pt() > track2->Pt()) fPairLoss->Fill((track1->Eta() - track2->Eta()), (track1->Phi() - track2->Phi() + TMath::ASin(0.075 * (1.2 / track2->Pt())) - TMath::ASin(0.075 * (1.2 / track1->Pt()))));
-   if (track2->Pt() > track1->Pt()) fPairLoss->Fill((track2->Eta() - track1->Eta()), (track2->Phi() - track1->Phi() + TMath::ASin(0.075 * (1.2 / track1->Pt())) - TMath::ASin(0.075 * (1.2 / track2->Pt()))));
+   // create candidate array
+   fCandidates = new TObjArray(1000);
+   fCandidates->SetOwner(kTRUE);
+
+   // create and post flowevent
+   fFlowEvent = new AliFlowEvent(10000);
+   PostData(2, fFlowEvent);
 }
 //_____________________________________________________________________________
 template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* track1, const T* track2) const
 {
    // Return the invariant mass of two tracks, assuming both tracks are kaons
    if ((!track2) || (!track1)) return 0.;
-   PairLoss(track1, track2);
    Double_t masss = TMath::Power(4.93676999999999977e-01, 2);
    Double_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
    Double_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
@@ -449,10 +321,7 @@ template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* tr
 //_____________________________________________________________________________
 template <typename T> Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const T* track1, const T* track2) const
 {
-   // Calculate the delta dip angle between two particles (the opening angle of a pair
-   // in the p_t p_z plane)
-   // A cut can be used to remove electron / positron contamination from photon conversion
-   // 999 (no cut) is returned when angle cannot be calculated
+   // Calculate the delta dip angle between two particles
    if (track1->P()*track2->P() == 0) return 999;
    return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
 }
@@ -461,20 +330,18 @@ template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const T*
 {
    // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
    if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PhiPt(track1, track2) < fDeltaDipPt)) return kFALSE;
-   fPIDDeltaDip->Fill(track1->P(), track1->GetTPCsignal());
-   fPIDDeltaDip->Fill(track2->P(), track2->GetTPCsignal());
    return kTRUE;
 }
 //_____________________________________________________________________________
 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
 {
-    // Check if pair passes eta and pt cut
-    if (fCandidateMinPt > PhiPt(track1, track2) || fCandidateMaxPt < PhiPt(track1, track2) ) return kFALSE;
-    TVector3 a(track1->Px(), track1->Py(), track1->Pz());
-    TVector3 b(track2->Px(), track2->Py(), track2->Pz());
-    TVector3 c = a + b;
-    if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
-    return kTRUE;
+   // Check if pair passes eta and pt cut
+   if (fCandidateMinPt > PhiPt(track1, track2) || fCandidateMaxPt < PhiPt(track1, track2)) return kFALSE;
+   TVector3 a(track1->Px(), track1->Py(), track1->Pz());
+   TVector3 b(track2->Px(), track2->Py(), track2->Pz());
+   TVector3 c = a + b;
+   if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
+   return kTRUE;
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskPhiFlow::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
@@ -528,12 +395,14 @@ void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliESDEvent* event)
 {
    // Initialize the Bayesian PID object for ESD analysis
    fBayesianResponse->SetDetResponse(event, fCentrality, AliESDpid::kTOF_T0, kTRUE);
+   if (fDebug) cout << " --> Initialized Bayesian ESD PID object " << endl;
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliAODEvent* event)
 {
    // Initialize the Bayesian PID object for AOD
    fBayesianResponse->SetDetResponse(event, fCentrality);
+   if (fDebug) cout << " --> Initialized Bayesian AOD PID object " << endl;
 }
 //_____________________________________________________________________________
 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* track) const
@@ -543,7 +412,7 @@ template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* tra
    fBayesianResponse->ComputeProb(track);
    if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
    Float_t *probabilities = fBayesianResponse->GetProb();
-   if (probabilities[3] > fParticleProbability)
+   if (probabilities[3] > fPIDConfig[6])
    {
       fPhi->Fill(track->Phi());
       fPt->Fill(track->Pt());
@@ -555,8 +424,7 @@ template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* tra
 //_____________________________________________________________________________
 Bool_t AliAnalysisTaskPhiFlow::PassesStrictKaonCuts(AliESDtrack* track) const
 {
-   // see if track passes additional kaon selection cuts
-   // with many thanks to francesco noferini
+   // propagate dca from TPC
    Double_t b[2] = { -99., -99.};
    Double_t bCov[3] = { -99., -99., -99.};
    if (!track->PropagateToDCA(fESD->GetPrimaryVertex(), fESD->GetMagneticField(), 100., b, bCov)) return kFALSE;
@@ -566,8 +434,7 @@ Bool_t AliAnalysisTaskPhiFlow::PassesStrictKaonCuts(AliESDtrack* track) const
 //_____________________________________________________________________________
 Bool_t AliAnalysisTaskPhiFlow::PassesStrictKaonCuts(AliAODTrack* track) const
 {
-   // see if track passes additional kaon selection cuts
-   // with many thanks to francesco noferini
+   // propagate dca from TPC
    Double_t b[2] = { -99., -99.};
    Double_t bCov[3] = { -99., -99., -99.};
    if (!track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
@@ -579,67 +446,64 @@ Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliESDtrack* track) const
 {
    // Check if particle is a kaon according to method set in steering macro
    if (fRequireTPCStandAlone && (track->GetStatus()&AliESDtrack::kTPCin) == 0) return kFALSE;
-   switch (fPIDtype)
+   if (fPIDConfig[1]!=0 || fPIDConfig[4]!=0) AliFatal("TPC || ITS PID not available in ESD anlaysis -- terminating analysis !!!");
+   if (PassesTPCbayesianCut(track))
    {
-      case kTPC:
-         // general tpc pid
-         fNOPID->Fill(track->P(), track->GetTPCsignal());
-         Double_t p[AliPID::kSPECIES];
-         track->GetTPCpid(p);
-         if(p[3] > fParticleProbability)
-         {
-            fPIDk->Fill(track->P(), track->GetTPCsignal());
-            return kTRUE;
-         }
-         break;
-      case kCombined:
-         // Use null_cuts dummy to check if track passes bayesian TPC TOF pid cut
-         fNOPID->Fill(track->P(), track->GetTPCsignal());
-         if (PassesTPCbayesianCut(track))
-         {
-            fPIDk->Fill(track->P(), track->GetTPCsignal());
-            return kTRUE;
-         }
-         break;
-      default:
-         AliFatal("No PID procedure set or available. Analysis will terminate!");
-         return kFALSE;
-         break;
+     fPIDk->Fill(track->P(), track->GetTPCsignal());
+     return kTRUE;
    }
    return kFALSE;
 }
 //_____________________________________________________________________________
 Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
 {
-   // Check if particle is a kaon according to method set in steering macro
-
+   // Kaon identification routine, based on multiple detectors and approaches
    if (fRequireTPCStandAlone && (!track->TestFilterBit(1))) return kFALSE;
-   switch (fPIDtype)
-   {
-      case kTPC:
-         fNOPID->Fill(track->P(), track->GetTPCsignal());
-//        Double_t p[AliPID::kSPECIES];
-//         AliAODpidUtil* pid;
-//         pid->MakeTPCPID(track, p);
-//         if(p[3] > fParticleProbability)
-//         {
-//            fPIDk->Fill(track->P(), track->GetTPCsignal());
-//            return kTRUE;
-//         }
-         return kFALSE;
-         break;
-      case kCombined:
-         fNOPID->Fill(track->P(), track->GetTPCsignal());
-         if (PassesTPCbayesianCut(track))
-         {
-            fPIDk->Fill(track->P(), track->GetTPCsignal());
-            return kTRUE;
-         }
-         break;
-      default:
-         AliFatal("No PID procedure set or available. Analysis will terminate!");
-         return kFALSE;
-         break;
+   fNOPID->Fill(track->P(), track->GetTPCsignal());
+   if(track->Pt() < fPIDConfig[1])
+   {
+       if(fDebug) cout << " ITS received track with p_t " << track->Pt() << endl;
+       // if tpc control is disabled, pure its pid
+       if(fPIDConfig[2] < 0.)
+       {
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) return kTRUE;
+           return kFALSE;
+       }
+       // else, switch to ITS pid with TPC rejection of protons and pions
+       if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) < fPIDConfig[3]) return kFALSE;
+       else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) < fPIDConfig[3]) return kFALSE;
+       else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0])
+       {
+           fPIDk->Fill(track->P(), track->GetTPCsignal());
+           return kTRUE;
+       }
+       return kFALSE;
+   }
+   if ((track->Pt() > fPIDConfig[1]) && (track->Pt() < fPIDConfig[4]))
+   {
+       if(fDebug) cout << " TPC received track with p_t " << track->Pt() << endl;
+       // if its control is disabled, pure tpc pid
+       if(fPIDConfig[5] < 0.)
+       {
+           if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) return kTRUE;
+           return kFALSE;
+       }
+       // else, switch to TPC pid with ITS rejection of protons and pions
+       if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kProton)) < fPIDConfig[0]) return kFALSE;
+       else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kPion)) < fPIDConfig[0]) return kFALSE;
+       else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3])
+       {
+           fPIDk->Fill(track->P(), track->GetTPCsignal());
+           return kTRUE;
+       }
+       return kFALSE;
+   }
+   if(fDebug) cout << " Bayesian method received track with p_t " << track->Pt() << endl;
+   // switch to bayesian PID
+   if (PassesTPCbayesianCut(track))
+   {
+       fPIDk->Fill(track->P(), track->GetTPCsignal());
+       return kTRUE;
    }
    return kFALSE;
 }
@@ -659,24 +523,96 @@ template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, c
    Double_t pt = PhiPt(track1, track2);
    if (tracktype == 0)
    {
-      if ((0.0 <= pt) && (0.3 > pt)) { fInvMNP03->Fill(InvariantMass(track1, track2)); fPtSpectra03->Fill(pt);}
-      if ((0.3 <= pt) && (0.6 > pt)) { fInvMNP36->Fill(InvariantMass(track1, track2)); fPtSpectra36->Fill(pt); }
-      if ((0.6 <= pt) && (0.9 > pt)) { fInvMNP69->Fill(InvariantMass(track1, track2)); fPtSpectra69->Fill(pt); }
-      if ((0.9 <= pt) && (1.2 > pt)) { fInvMNP912->Fill(InvariantMass(track1, track2)); fPtSpectra912->Fill(pt); }
-      if ((1.2 <= pt) && (1.5 > pt)) { fInvMNP1215->Fill(InvariantMass(track1, track2)); fPtSpectra1215->Fill(pt); }
-      if ((1.5 <= pt) && (1.8 > pt)) { fInvMNP1518->Fill(InvariantMass(track1, track2)); fPtSpectra1518->Fill(pt); }
-      if ((1.8 <= pt) && (2.1 > pt)) { fInvMNP1821->Fill(InvariantMass(track1, track2)); fPtSpectra1821->Fill(pt); }
-      if ((2.1 <= pt) && (2.4 > pt)) { fInvMNP2124->Fill(InvariantMass(track1, track2)); fPtSpectra2124->Fill(pt); }
-      if ((2.4 <= pt) && (2.7 > pt)) { fInvMNP2427->Fill(InvariantMass(track1, track2)); fPtSpectra2427->Fill(pt); }
-      if ((2.7 <= pt) && (3.0 > pt)) { fInvMNP2730->Fill(InvariantMass(track1, track2)); fPtSpectra2730->Fill(pt); }
-      if ((3.0 <= pt) && (3.5 > pt)) { fInvMNP3035->Fill(InvariantMass(track1, track2)); fPtSpectra3035->Fill(pt); }
-      if ((3.5 <= pt) && (4.0 > pt)) { fInvMNP3540->Fill(InvariantMass(track1, track2)); fPtSpectra3540->Fill(pt); }
-      if ((4.0 <= pt) && (4.5 > pt)) { fInvMNP4045->Fill(InvariantMass(track1, track2)); fPtSpectra4045->Fill(pt); }
-      if ((4.5 <= pt) && (5.0 > pt)) { fInvMNP4550->Fill(InvariantMass(track1, track2)); fPtSpectra4550->Fill(pt); }
-      if ((5.0 <= pt) && (5.5 > pt)) { fInvMNP5055->Fill(InvariantMass(track1, track2)); fPtSpectra5055->Fill(pt); }
-      if ((5.5 <= pt) && (6.0 > pt)) { fInvMNP5560->Fill(InvariantMass(track1, track2)); fPtSpectra5560->Fill(pt); }
-      if ((6.0 <= pt) && (6.5 > pt)) { fInvMNP6065->Fill(InvariantMass(track1, track2)); fPtSpectra6065->Fill(pt); }
-      if ((6.5 <= pt) && (7.0 > pt)) { fInvMNP6570->Fill(InvariantMass(track1, track2)); fPtSpectra6570->Fill(pt); }
+      if ((0.0 <= pt) && (0.3 > pt))
+      {
+         fInvMNP03->Fill(InvariantMass(track1, track2));
+         fPtSpectra03->Fill(pt);
+      }
+      if ((0.3 <= pt) && (0.6 > pt))
+      {
+         fInvMNP36->Fill(InvariantMass(track1, track2));
+         fPtSpectra36->Fill(pt);
+      }
+      if ((0.6 <= pt) && (0.9 > pt))
+      {
+         fInvMNP69->Fill(InvariantMass(track1, track2));
+         fPtSpectra69->Fill(pt);
+      }
+      if ((0.9 <= pt) && (1.2 > pt))
+      {
+         fInvMNP912->Fill(InvariantMass(track1, track2));
+         fPtSpectra912->Fill(pt);
+      }
+      if ((1.2 <= pt) && (1.5 > pt))
+      {
+         fInvMNP1215->Fill(InvariantMass(track1, track2));
+         fPtSpectra1215->Fill(pt);
+      }
+      if ((1.5 <= pt) && (1.8 > pt))
+      {
+         fInvMNP1518->Fill(InvariantMass(track1, track2));
+         fPtSpectra1518->Fill(pt);
+      }
+      if ((1.8 <= pt) && (2.1 > pt))
+      {
+         fInvMNP1821->Fill(InvariantMass(track1, track2));
+         fPtSpectra1821->Fill(pt);
+      }
+      if ((2.1 <= pt) && (2.4 > pt))
+      {
+         fInvMNP2124->Fill(InvariantMass(track1, track2));
+         fPtSpectra2124->Fill(pt);
+      }
+      if ((2.4 <= pt) && (2.7 > pt))
+      {
+         fInvMNP2427->Fill(InvariantMass(track1, track2));
+         fPtSpectra2427->Fill(pt);
+      }
+      if ((2.7 <= pt) && (3.0 > pt))
+      {
+         fInvMNP2730->Fill(InvariantMass(track1, track2));
+         fPtSpectra2730->Fill(pt);
+      }
+      if ((3.0 <= pt) && (3.5 > pt))
+      {
+         fInvMNP3035->Fill(InvariantMass(track1, track2));
+         fPtSpectra3035->Fill(pt);
+      }
+      if ((3.5 <= pt) && (4.0 > pt))
+      {
+         fInvMNP3540->Fill(InvariantMass(track1, track2));
+         fPtSpectra3540->Fill(pt);
+      }
+      if ((4.0 <= pt) && (4.5 > pt))
+      {
+         fInvMNP4045->Fill(InvariantMass(track1, track2));
+         fPtSpectra4045->Fill(pt);
+      }
+      if ((4.5 <= pt) && (5.0 > pt))
+      {
+         fInvMNP4550->Fill(InvariantMass(track1, track2));
+         fPtSpectra4550->Fill(pt);
+      }
+      if ((5.0 <= pt) && (5.5 > pt))
+      {
+         fInvMNP5055->Fill(InvariantMass(track1, track2));
+         fPtSpectra5055->Fill(pt);
+      }
+      if ((5.5 <= pt) && (6.0 > pt))
+      {
+         fInvMNP5560->Fill(InvariantMass(track1, track2));
+         fPtSpectra5560->Fill(pt);
+      }
+      if ((6.0 <= pt) && (6.5 > pt))
+      {
+         fInvMNP6065->Fill(InvariantMass(track1, track2));
+         fPtSpectra6065->Fill(pt);
+      }
+      if ((6.5 <= pt) && (7.0 > pt))
+      {
+         fInvMNP6570->Fill(InvariantMass(track1, track2));
+         fPtSpectra6570->Fill(pt);
+      }
    }
    if (tracktype == 1)
    {
@@ -722,501 +658,6 @@ template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, c
    }
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskPhiFlow::EventPlane(const AliESDEvent* event)
-{
-   // Calculate Q-vectors and event plane resolution (two different methods)
-   fQy = 0;
-   fQx = 0;
-   Double_t ceventsin(0);
-   Double_t ceventcos(0);
-   Double_t deventsin(0);
-   Double_t deventcos(0);
-   Double_t eeventsin(0);
-   Double_t eeventcos(0);
-   Double_t feventsin(0);
-   Double_t feventcos(0);
-   Int_t nTracks = event->GetNumberOfTracks();
-   Int_t tpcMultiplicity(0);
-   for (Int_t i = 0; i < nTracks; i++)
-   {
-      AliESDtrack* track = event->GetTrack(i);
-      if (!EventPlaneTrack(track)) continue;
-      if (fSetEventPlanePtCut && (track->Pt() > fEventPlanePtCut)) continue;
-      tpcMultiplicity++;
-      Double_t r = gRandom->Uniform(0, 1);
-      if (r < 0.5)
-      {
-         ceventsin += TMath::Sin(2 * track->Phi());
-         ceventcos += TMath::Cos(2 * track->Phi());
-      }
-      if (r > 0.5)
-      {
-         deventsin += TMath::Sin(2 * track->Phi());
-         deventcos += TMath::Cos(2 * track->Phi());
-      }
-      if (track->Eta() > 0)
-      {
-         eeventsin += TMath::Sin(2 * track->Phi());
-         eeventcos += TMath::Cos(2 * track->Phi());
-      }
-      if (track->Eta() < 0)
-      {
-         feventsin += TMath::Sin(2 * track->Phi());
-         feventcos += TMath::Cos(2 * track->Phi());
-      }
-      fQy += TMath::Sin(2 * track->Phi());
-      fQx += TMath::Cos(2 * track->Phi());
-   }
-   fTPCM->Fill(tpcMultiplicity);
-   EventPlaneResolution(kTRUE, ceventcos, ceventsin, deventcos, deventsin, fQx, fQy, eeventcos, eeventsin, feventcos, feventsin);
-   EventPlaneResolution(kFALSE, ceventcos, ceventsin, deventcos, deventsin, fQx, fQy, eeventcos, eeventsin, feventcos, feventsin);
-   fEventPlane->Fill((TMath::ATan2(fQy, fQx) / 2) + TMath::Pi() / 2.);
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskPhiFlow::EventPlane(const AliAODEvent* event)
-{
-   // Overloaded function - see comment at EventPlane(const AliESDEvent* event);
-   fQy = 0;
-   fQx = 0;
-   Double_t ceventsin(0);
-   Double_t ceventcos(0);
-   Double_t deventsin(0);
-   Double_t deventcos(0);
-   Double_t eeventsin(0);
-   Double_t eeventcos(0);
-   Double_t feventsin(0);
-   Double_t feventcos(0);
-   Int_t nTracks = event->GetNumberOfTracks();
-   Int_t tpcMultiplicity(0);
-   for (Int_t i = 0; i < nTracks; i++)
-   {
-      AliAODTrack* track = event->GetTrack(i);
-      if (!EventPlaneTrack(track)) continue;
-      if (fSetEventPlanePtCut && (track->Pt() > fEventPlanePtCut)) continue;
-      tpcMultiplicity++;
-      Double_t r = gRandom->Uniform(0, 1);
-      if (r < 0.5)
-      {
-         ceventsin += TMath::Sin(2 * track->Phi());
-         ceventcos += TMath::Cos(2 * track->Phi());
-      }
-      if (r > 0.5)
-      {
-         deventsin += TMath::Sin(2 * track->Phi());
-         deventcos += TMath::Cos(2 * track->Phi());
-      }
-      if (track->Eta() > 0)
-      {
-         eeventsin += TMath::Sin(2 * track->Phi());
-         eeventcos += TMath::Cos(2 * track->Phi());
-      }
-      if (track->Eta() < 0)
-      {
-         feventsin += TMath::Sin(2 * track->Phi());
-         feventcos += TMath::Cos(2 * track->Phi());
-      }
-      fQy += TMath::Sin(2 * track->Phi());
-      fQx += TMath::Cos(2 * track->Phi());
-   }
-   fTPCM->Fill(tpcMultiplicity);
-   EventPlaneResolution(kTRUE, ceventcos, ceventsin, deventcos, deventsin, fQx, fQy, eeventcos, eeventsin, feventcos, feventsin);
-   EventPlaneResolution(kFALSE, ceventcos, ceventsin, deventcos, deventsin, fQx, fQy, eeventcos, eeventsin, feventcos, feventsin);
-   fEventPlane->Fill((TMath::ATan2(fQy, fQx) / 2 + TMath::Pi() / 2.));
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskPhiFlow::EventPlaneResolution(Bool_t random,
-                                                     Double_t cossumQ1,
-                                                     Double_t sinsumQ1,
-                                                     Double_t cossumQ2,
-                                                     Double_t sinsumQ2,
-                                                     Double_t cossum,
-                                                     Double_t sinsum,
-                                                     Double_t cossumQwest,
-                                                     Double_t sinsumQwest,
-                                                     Double_t cossumQeast,
-                                                     Double_t sinsumQeast) const
-{
-   // Calculate the event plane resolution (called twice, for different methods)
-   if (random)
-   {
-      Double_t res2Sub = EventPlaneStar(kTRUE, cossumQ1, sinsumQ1, cossumQ2, sinsumQ2, cossum, sinsum, cossumQwest, sinsumQwest, cossumQeast, sinsumQeast);
-      if (res2Sub < -900) return;
-      Double_t chiSub2 = Chi(res2Sub);
-      Double_t res2 = ResolutionAsFunctionOfChi(TMath::Sqrt(2.) * chiSub2); // full event plane res.
-      fEventPlaneResolutionRandom->Fill(1, res2);
-   }
-   if (!random)
-   {
-      Double_t res2Sub = EventPlaneStar(kFALSE, cossumQ1, sinsumQ1, cossumQ2, sinsumQ2, cossum, sinsum, cossumQwest, sinsumQwest, cossumQeast, sinsumQeast);
-      if (res2Sub < -900) return;
-      Double_t chiSub2 = Chi(res2Sub);
-      Double_t res2 = ResolutionAsFunctionOfChi(TMath::Sqrt(2.) * chiSub2); // full event plane res.
-      fEventPlaneResolutionEta->Fill(1, res2);
-   }
-}
-//_____________________________________________________________________________
-Double_t AliAnalysisTaskPhiFlow::EventPlaneStar(Bool_t random,
-                                                   Double_t cossumQ1,
-                                                   Double_t sinsumQ1,
-                                                   Double_t cossumQ2,
-                                                   Double_t sinsumQ2,
-                                                   Double_t cossum,
-                                                   Double_t sinsum,
-                                                   Double_t cossumQwest,
-                                                   Double_t sinsumQwest,
-                                                   Double_t cossumQeast,
-                                                   Double_t sinsumQeast) const
-{
-   // Calculate subevent resolution
-   // This method also calculates the event plane
-   // This implementation is copied from STAR analysis code
-   TVector2 *q1 = new TVector2(cossumQ1, sinsumQ1);
-   TVector2 *q2 = new TVector2(cossumQ2, sinsumQ2);
-   TVector2 *qt = new TVector2(cossum, sinsum);
-   TVector2 *qwest  = new TVector2(cossumQwest, sinsumQwest);
-   TVector2 *qeast = new TVector2(cossumQeast, sinsumQeast);
-   if (qt->Mod() == 0)
-   {
-      delete qt;
-      delete q2;
-      delete q1;
-      delete qwest;
-      delete qeast;
-      return -999.; // make sure length > 0
-   }
-   Double_t twopi = TMath::TwoPi();
-
-   Double_t rxnPlane, rxnPlane1, rxnPlane2;
-   rxnPlane = qt->Phi() / 2.; //reaction plane angle ranges from 0 to pi
-   rxnPlane1 = q1->Phi() / 2.;
-   rxnPlane2 = q2->Phi() / 2.;
-   if ((qt->Mod()) != 0.0)
-   {
-      Double_t  eventPlane = 0.5 * qt->Phi();
-      if (random) fEventPlaneSTAR->Fill(eventPlane);
-   }
-   if (random)
-   {
-      if ((q1->Mod2()) != 0.0 && (q2->Mod2()) != 0.0)
-      {
-         Double_t eventDiff = 0.5 * q1->Phi() - 0.5 * q2->Phi();
-         if (eventDiff < 0.0) eventDiff += twopi / 2.0;
-         return TMath::Cos(2.*eventDiff);
-      }
-   }
-   if (!random)
-   {
-      if ((qwest->Mod2()) != 0.0 && (qeast->Mod2()) != 0.0)
-      {
-         Double_t eventDiff = 0.5 * qwest->Phi() - 0.5 * qeast->Phi();
-         if (eventDiff < 0.0) eventDiff += twopi / 2.0;
-         return TMath::Cos(2.*eventDiff);
-      }
-   }
-   return -999.;
-}
-//_____________________________________________________________________________
-Double_t AliAnalysisTaskPhiFlow::Chi(Double_t res) const
-{
-   // Iterative algorithm to solve R(chi) for chi
-   Double_t chi   = 2.0;
-   Double_t delta = 1.0;
-   for (Int_t i = 0; i < 15; i++)
-   {
-      chi   = (ResolutionAsFunctionOfChi(chi) < res) ? chi + delta : chi - delta;
-      delta = delta / 2.;
-   }
-   return chi;
-}
-//_____________________________________________________________________________
-Double_t AliAnalysisTaskPhiFlow::ResolutionAsFunctionOfChi(Double_t chi) const
-{
-   // Definition of R(chi)
-
-   Double_t con = 0.626657;                   // sqrt(pi/2)/2
-   Double_t arg = chi * chi / 4.;
-   Double_t res = con * chi * exp(-arg) * (TMath::BesselI0(arg) + TMath::BesselI1(arg));
-   return res;
-}
-//_____________________________________________________________________________
-template <typename T> void AliAnalysisTaskPhiFlow::EllipticFlow(Double_t* const v2, const T* track1, const T* track2) const
-{
-   // 1) Calculate elliptic flow coefficient v_2 - cos terms
-   // 2) Store v_2 in array for further analysis
-   // 3) Fill phi - Psi histograms per p_t bin
-   Double_t y = fQy - (TMath::Sin(2 * track1->Phi()) + TMath::Sin(2 * track2->Phi()));
-   Double_t x = fQx - (TMath::Cos(2 * track1->Phi()) + TMath::Cos(2 * track2->Phi()));
-   Double_t psi2 = (TMath::ATan2(y, x)) / 2.;
-   Double_t psi2_sub(0);
-   psi2 > 0 ? psi2_sub = psi2: psi2_sub = psi2 + TMath::Pi();
-   Double_t pt = PhiPt(track1, track2);
-   // to calculate the angle between the particles, firstly, construct the p_t vector of the phi
-   TVector3 a(track1->Px(), track1->Py(), track1->Pz());
-   TVector3 b(track2->Px(), track2->Py(), track2->Pz());
-   TVector3 c = a + b;
-   Double_t phi = c.Phi();
-   Double_t invm  = InvariantMass(track1, track2);
-   Double_t flow = TMath::Cos(2 * (phi - psi2));
-   if ((0 <= pt) && (0.3 > pt))
-   {
-      v2[0] += flow;
-      v2[18]++;
-      fProfV2InvM03->Fill(invm, flow);
-      fDeltaPhiPsiNP03->Fill(phi - psi2_sub);
-   }
-   if ((0.3 <= pt) && (0.6 > pt))
-   {
-      v2[1] += flow;
-      v2[19]++;
-      fProfV2InvM36->Fill(invm, flow);
-      fDeltaPhiPsiNP36->Fill(phi - psi2_sub);
-   }
-   if ((0.6 <= pt) && (0.9 > pt))
-   {
-      v2[2] += flow;
-      v2[20]++;
-      fProfV2InvM69->Fill(invm, flow);
-      fDeltaPhiPsiNP69->Fill(phi - psi2_sub);
-   }
-   if ((0.9 <= pt) && (1.2 > pt))
-   {
-      v2[3] += flow;
-      v2[21]++;
-      fProfV2InvM912->Fill(invm, flow);
-      fDeltaPhiPsiNP912->Fill(phi - psi2_sub);
-   }
-   if ((1.2 <= pt) && (1.5 > pt))
-   {
-      v2[4] += flow;
-      v2[22]++;
-      fProfV2InvM1215->Fill(invm, flow);
-      fDeltaPhiPsiNP1215->Fill(phi - psi2_sub);
-   }
-   if ((1.5 <= pt) && (1.8 > pt))
-   {
-      v2[5] += flow;
-      v2[23]++;
-      fProfV2InvM1518->Fill(invm, flow);
-      fDeltaPhiPsiNP1518->Fill(phi - psi2_sub);
-   }
-   if ((1.8 <= pt) && (2.1 > pt))
-   {
-      v2[6] += flow;
-      v2[24]++;
-      fProfV2InvM1821->Fill(invm, flow);
-      fDeltaPhiPsiNP1821->Fill(phi - psi2_sub);
-   }
-   if ((2.1 <= pt) && (2.4 > pt))
-   {
-      v2[7] += flow;
-      v2[25]++;
-      fProfV2InvM2124->Fill(invm, flow);
-      fDeltaPhiPsiNP2124->Fill(phi - psi2_sub);
-   }
-   if ((2.4 <= pt) && (2.7 > pt))
-   {
-      v2[8] += flow;
-      v2[26]++;
-      fProfV2InvM2427->Fill(invm, flow);
-      fDeltaPhiPsiNP2427->Fill(phi - psi2_sub);
-   }
-   if ((2.7 <= pt) && (3.0 > pt))
-   {
-      v2[9] += flow;
-      v2[27]++;
-      fProfV2InvM2730->Fill(invm, flow);
-      fDeltaPhiPsiNP2730->Fill(phi - psi2_sub);
-   }
-   if ((3.0 <= pt) && (3.5 > pt))
-   {
-      v2[10] += flow;
-      v2[28]++;
-      fProfV2InvM3035->Fill(invm, flow);
-      fDeltaPhiPsiNP3035->Fill(phi - psi2_sub);
-   }
-   if ((3.5 <= pt) && (4.0 > pt))
-   {
-      v2[11] += flow;
-      v2[29]++;
-      fProfV2InvM3540->Fill(invm, flow);
-      fDeltaPhiPsiNP3540->Fill(phi - psi2_sub);
-   }
-   if ((4.0 <= pt) && (4.5 > pt))
-   {
-      v2[12] += flow;
-      v2[30]++;
-      fProfV2InvM4045->Fill(invm, flow);
-      fDeltaPhiPsiNP4045->Fill(phi - psi2_sub);
-   }
-   if ((4.5 <= pt) && (5.0 > pt))
-   {
-      v2[13] += flow;
-      v2[31]++;
-      fProfV2InvM4550->Fill(invm, flow);
-      fDeltaPhiPsiNP4550->Fill(phi - psi2_sub);
-   }
-   if ((5.0 <= pt) && (5.5 > pt))
-   {
-      v2[14] += flow;
-      v2[32]++;
-      fProfV2InvM5055->Fill(invm, flow);
-      fDeltaPhiPsiNP5055->Fill(phi - psi2_sub);
-   }
-   if ((5.5 <= pt) && (6.0 > pt))
-   {
-      v2[15] += flow;
-      v2[33]++;
-      fProfV2InvM5560->Fill(invm, flow);
-      fDeltaPhiPsiNP5560->Fill(phi - psi2_sub);
-   }
-   if ((6.0 <= pt) && (6.5 > pt))
-   {
-      v2[16] += flow;
-      v2[34]++;
-      fProfV2InvM6065->Fill(invm, flow);
-      fDeltaPhiPsiNP6065->Fill(phi - psi2_sub);
-   }
-   if ((6.5 <= pt) && (7.0 > pt))
-   {
-      v2[17] += flow;
-      v2[35]++;
-      fProfV2InvM6570->Fill(invm, flow);
-      fDeltaPhiPsiNP6570->Fill(phi - psi2_sub);
-   }
-}
-//_____________________________________________________________________________
-template <typename T> void AliAnalysisTaskPhiFlow::EllipticFlowSin(Double_t* const v2Sin, const T* track1, const T* track2) const
-{
-   // Calculate elliptic flow coefficient v_2 - sin terms
-   // store v_2 in array for further analysis
-   Double_t y = fQy - (TMath::Sin(2 * track1->Phi()) + TMath::Sin(2 * track2->Phi()));
-   Double_t x = fQx - (TMath::Cos(2 * track1->Phi()) + TMath::Cos(2 * track2->Phi()));
-   Double_t psi2 = (TMath::ATan2(y, x)) / 2.;
-   Double_t pt = PhiPt(track1, track2);
-   TVector3 a(track1->Px(), track1->Py(), track1->Pz());
-   TVector3 b(track2->Px(), track2->Py(), track2->Pz());
-   TVector3 c = a + b;
-   Double_t phi = c.Phi();
-   Double_t invm = InvariantMass(track1, track2);
-   Double_t flow = TMath::Sin(2 * (phi - psi2));
-   if ((0 <= pt) && (0.3 > pt))
-   {
-      v2Sin[0] += flow;
-      v2Sin[18]++;
-      fProfV2SinInvM03->Fill(invm, flow);
-   }
-   if ((0.3 <= pt) && (0.6 > pt))
-   {
-      v2Sin[1] += flow;
-      v2Sin[19]++;
-      fProfV2SinInvM36->Fill(invm, flow);
-   }
-   if ((0.6 <= pt) && (0.9 > pt))
-   {
-      v2Sin[2] += flow;
-      v2Sin[20]++;
-      fProfV2SinInvM69->Fill(invm, flow);
-   }
-   if ((0.9 <= pt) && (1.2 > pt))
-   {
-      v2Sin[3] += flow;
-      v2Sin[21]++;
-      fProfV2SinInvM912->Fill(invm, flow);
-   }
-   if ((1.2 <= pt) && (1.5 > pt))
-   {
-      v2Sin[4] += flow;
-      v2Sin[22]++;
-      fProfV2SinInvM1215->Fill(invm, flow);
-   }
-   if ((1.5 <= pt) && (1.8 > pt))
-   {
-      v2Sin[5] += flow;
-      v2Sin[23]++;
-      fProfV2SinInvM1518->Fill(invm, flow);
-   }
-   if ((1.8 <= pt) && (2.1 > pt))
-   {
-      v2Sin[6] += flow;
-      v2Sin[24]++;
-      fProfV2SinInvM1821->Fill(invm, flow);
-   }
-   if ((2.1 <= pt) && (2.4 > pt))
-   {
-      v2Sin[7] += flow;
-      v2Sin[25]++;
-      fProfV2SinInvM2124->Fill(invm, flow);
-   }
-   if ((2.4 <= pt) && (2.7 > pt))
-   {
-      v2Sin[8] += flow;
-      v2Sin[26]++;
-      fProfV2SinInvM2427->Fill(invm, flow);
-   }
-   if ((2.7 <= pt) && (3.0 > pt))
-   {
-      v2Sin[9] += flow;
-      v2Sin[27]++;
-      fProfV2SinInvM2730->Fill(invm, flow);
-   }
-
-   if ((3.0 <= pt) && (3.5 > pt))
-   {
-      v2Sin[10] += flow;
-      v2Sin[28]++;
-      fProfV2SinInvM3035->Fill(invm, flow);
-   }
-   if ((3.5 <= pt) && (4.0 > pt))
-   {
-      v2Sin[11] += flow;
-      v2Sin[29]++;
-      fProfV2SinInvM3540->Fill(invm, flow);
-   }
-   if ((4.0 <= pt) && (4.5 > pt))
-   {
-      v2Sin[12] += flow;
-      v2Sin[30]++;
-      fProfV2SinInvM4045->Fill(invm, flow);
-   }
-   if ((4.5 <= pt) && (5.0 > pt))
-   {
-      v2Sin[13] += flow;
-      v2Sin[31]++;
-      fProfV2SinInvM4550->Fill(invm, flow);
-   }
-   if ((5.0 <= pt) && (5.5 > pt))
-   {
-      v2Sin[14] += flow;
-      v2Sin[32]++;
-      fProfV2SinInvM5055->Fill(invm, flow);
-   }
-   if ((5.5 <= pt) && (6.0 > pt))
-   {
-      v2Sin[15] += flow;
-      v2Sin[33]++;
-      fProfV2SinInvM5560->Fill(invm, flow);
-   }
-   if ((6.0 <= pt) && (6.5 > pt))
-   {
-      v2Sin[16] += flow;
-      v2Sin[34]++;
-      fProfV2SinInvM6065->Fill(invm, flow);
-   }
-   if ((6.5 <= pt) && (7.0 > pt))
-   {
-      v2Sin[17] += flow;
-      v2Sin[35]++;
-      fProfV2SinInvM6570->Fill(invm, flow);
-   }
-}
-//_____________________________________________________________________________
-template <typename T> Bool_t AliAnalysisTaskPhiFlow::EventPlaneTrack(T* track) const
-{
-   // Check if track is suitable for event plane estimation
-   if (!track) return kFALSE;
-   return fCutsRP->IsSelected(track);
-}
-//_____________________________________________________________________________
 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
 {
    // Check if track is suitable for phi flow analysis
@@ -1224,31 +665,10 @@ template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
    return fPOICuts->IsSelected(track);
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskPhiFlow::FlowFinish(Double_t* const v2) const
-{
-   // Push track averaged (event) flow into flow array - cos terms
-   Double_t pt[] = {150, 450, 750, 1050, 1350, 1650, 1950, 2250, 2550, 2850, 3250, 3750, 4250, 4750, 5250, 5750, 6250, 6750};
-   for (Int_t i = 0; i < 18; i++)
-   {
-      if (v2[i + 18] == 0) continue;
-      fProfV2->Fill(pt[i] / 1000., (v2[i] / v2[i + 18]));
-   }
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskPhiFlow::FlowFinishSin(Double_t* const v2Sin) const
-{
-   // Push track averaged (event) flow into flow array - sin terms
-   Double_t pt[] = {150, 450, 750, 1050, 1350, 1650, 1950, 2250, 2550, 2850, 3250, 3750, 4250, 4750, 5250, 5750, 6250, 6750};
-   for (Int_t i = 0; i < 18; i++)
-   {
-      if (v2Sin[i + 18] == 0) continue;
-      fProfV2Sin->Fill(pt[i] / 1000., (v2Sin[i] / v2Sin[i + 18]));
-   }
-}
-//_____________________________________________________________________________
 template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
 {
    // Set null cuts
+   if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
    fCutsRP->SetEvent(event, MCEvent());
    fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
    fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
@@ -1259,21 +679,21 @@ template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
 void AliAnalysisTaskPhiFlow::PrepareFlowEvent(Int_t iMulti)
 {
    // Prepare flow events
-   for (int r = 0; r != 30; ++r)
-   {
-      fFlowEvent[r]->ClearFast();
-      fFlowEvent[r]->Fill(fCutsRP, fNullCuts);
-      fFlowEvent[r]->SetReferenceMultiplicity(iMulti);
-      fFlowEvent[r]->DefineDeadZone(0, 0, 0, 0);
-      fFlowEvent[r]->TagSubeventsInEta(fEtaMinA, fEtaMaxA,
-                                       fEtaMinB, fEtaMaxB);
-   }
+   fFlowEvent->ClearFast();
+   fFlowEvent->Fill(fCutsRP, fNullCuts);
+   fFlowEvent->SetReferenceMultiplicity(iMulti);
+   fFlowEvent->DefineDeadZone(0, 0, 0, 0);
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
 {
-   // UserExec: execute for each event. Commented where necessary
+   // UserExec: called for each event. Commented where necessary
    // check for AOD data type
+   if (!fPIDResponse)
+   {
+      AliError("Cannot get pid response");
+      return;
+   }
 
    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
    if (fAOD)
@@ -1281,31 +701,24 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
       fAODAnalysis = kTRUE;
       // Check whether event passes event cuts
       if (!EventCut(fAOD)) return;
+      // initialize event objects
       InitializeBayesianPID(fAOD);
       SetNullCuts(fAOD);
       PrepareFlowEvent(fAOD->GetNumberOfTracks());
+      fCandidates->SetLast(-1);
       // Calculate event plane Q vectors and event plane resolution
-      EventPlane(fAOD);
       fEventStats->Fill(0);
       Int_t unTracks = fAOD->GetNumberOfTracks();
       AliAODTrack* un[unTracks];
       AliAODTrack* up[unTracks];
       Int_t unp(0);
       Int_t unn(0);
-      Double_t v2[36];
-      Double_t v2Sin[36];
-      // Flush flow arrays
-      for (Int_t i = 0; i < 36; i++)
-      {
-         v2[i] = 0.;
-         v2Sin[i] = 0.;
-      }
       // Loop through tracks, check for species (Kaons), fill arrays according to charge
       for (Int_t iTracks = 0; iTracks < unTracks; iTracks++)
       {
          AliAODTrack* track = fAOD->GetTrack(iTracks);
          if (!PhiTrack(track)) continue;
-         if (fStrictKaonCuts&&(!PassesStrictKaonCuts(track))) continue;
+         if (fStrictKaonCuts && (!PassesStrictKaonCuts(track))) continue;
          Bool_t charge = kFALSE;
          if (track->Charge() > 0)
          {
@@ -1344,8 +757,6 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
             if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
             PtSelector(0, up[pTracks], un[nTracks]);
-            EllipticFlow(v2, up[pTracks], un[nTracks]);
-            EllipticFlowSin(v2Sin, up[pTracks], un[nTracks]);
             Double_t pt = PhiPt(up[pTracks], un[nTracks]);
             Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
             TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
@@ -1356,26 +767,39 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
             Int_t nIDs[2];
             nIDs[0] = up[pTracks]->GetID();
             nIDs[1] = un[nTracks]->GetID();
-            for (Int_t r = 0; r != 30; ++r)
-               if ((mass >= fFlowBands[0][r]) && (mass < fFlowBands[1][r]))
+            MakeTrack(mass, pt, phi, eta, 2, nIDs);
+         }
+      }
+
+      if (fDebug)  printf("I received %d candidates\n", fCandidates->GetEntriesFast());
+      for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand)
+      {
+         AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
+         if (!cand) continue;
+         if (fDebug) printf(" >Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
+         for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau)
+         {
+            if (fDebug) printf("  >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
+            for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs)
+            {
+               AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
+               if (!iRP) continue;
+               if (!iRP->InRPSelection()) continue;
+               if (cand->GetIDDaughter(iDau) == iRP->GetID())
                {
-                  AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*)
-                                                  MakeTrack(mass, pt, phi, eta, 2, nIDs);
-                  for (Int_t iDau = 0; iDau != 2; ++iDau)
-                     for (Int_t iRPs = 0; iRPs != fFlowEvent[r]->NumberOfTracks(); ++iRPs)
-                     {
-                        AliFlowTrack *iRP = (AliFlowTrack*)(fFlowEvent[r]->GetTrack(iRPs));
-                        if (!iRP->InRPSelection()) continue;
-                        if (fabs(sTrack->GetIDDaughter(iDau)) == fabs(iRP->GetID()))
-                        {
-                           sTrack->SetDaughter(iDau, iRP);
-                           iRP->SetForRPSelection(kFALSE);
-                        }
-                     }
-                  fFlowEvent[r]->AddTrack(sTrack);
+                  if (fDebug) printf(" was in RP set");
+                  iRP->SetForRPSelection(kFALSE);
+                  fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
                }
+            }
+            if (fDebug) printf("\n");
          }
+         cand->SetForPOISelection(kTRUE);
+         fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
       }
+      if (fDebug) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
+
+
       for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
       {
          for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++)
@@ -1394,12 +818,8 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
             PtSelector(2, un[nTracks], un[pTracks]);
          }
       }
-      //Push event averaged track values into flow members
-      FlowFinish(v2);
-      FlowFinishSin(v2Sin);
-
       PostData(1, fOutputList);
-      for (int m = 0; m != 30; ++m) PostData(2 + m, fFlowEvent[m]);
+      PostData(2, fFlowEvent);
    }
 
    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
@@ -1412,21 +832,12 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
       SetNullCuts(fESD);
       PrepareFlowEvent(fESD->GetNumberOfTracks());
       // Calculate event plane Q vectors and event plane resolution
-      EventPlane(fESD);
       fEventStats->Fill(0);
       Int_t unTracks = fESD->GetNumberOfTracks();
       AliESDtrack* un[unTracks];
       AliESDtrack* up[unTracks];
       Int_t unp(0);
       Int_t unn(0);
-      Double_t v2[36];
-      Double_t v2Sin[36];
-      // Flush flow arrays
-      for (Int_t i = 0; i < 36; i++)
-      {
-         v2[i] = 0.;
-         v2Sin[i] = 0.;
-      }
       // Loop through tracks, check for species (Kaons), fill arrays according to charge
       for (Int_t iTracks = 0; iTracks < unTracks; iTracks++)
       {
@@ -1463,15 +874,13 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
          }
       }
       // Calculate invariant mass of like- and unlike sign pairs as a function of p_T, store data in histograms
-      for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
+        for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
       {
          for (Int_t nTracks = 0; nTracks < unn ; nTracks++)
          {
             if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
             PtSelector(0, up[pTracks], un[nTracks]);
-            EllipticFlow(v2, up[pTracks], un[nTracks]);
-            EllipticFlowSin(v2Sin, up[pTracks], un[nTracks]);
             Double_t pt = PhiPt(up[pTracks], un[nTracks]);
             Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
             TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
@@ -1479,31 +888,41 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
             TVector3 c = a + b;
             Double_t phi = c.Phi();
             Double_t eta = c.Eta();
-            int nIDs[2];
+            Int_t nIDs[2];
             nIDs[0] = up[pTracks]->GetID();
             nIDs[1] = un[nTracks]->GetID();
-            for (Int_t r = 0; r != 30; ++r)
-               if ((mass >= fFlowBands[0][r]) && (mass < fFlowBands[1][r]))
+            MakeTrack(mass, pt, phi, eta, 2, nIDs);
+         }
+      }
+
+      if (fDebug)  printf("I received %d candidates\n", fCandidates->GetEntriesFast());
+      for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand)
+      {
+         AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
+         if (!cand) continue;
+         if (fDebug) printf(" >Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
+         for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau)
+         {
+            if (fDebug) printf("  >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
+            for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs)
+            {
+               AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
+               if (!iRP) continue;
+               if (!iRP->InRPSelection()) continue;
+               if (cand->GetIDDaughter(iDau) == iRP->GetID())
                {
-                  AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*)
-                                                  MakeTrack(mass, pt, phi, eta, 2, nIDs);
-                  if (0) printf("   á¶«Injecting phi candidate on band %d \n", r);
-                  for (Int_t iDau = 0; iDau != 2; ++iDau)
-                     for (Int_t iRPs = 0; iRPs != fFlowEvent[r]->NumberOfTracks(); ++iRPs)
-                     {
-                        AliFlowTrack *iRP = (AliFlowTrack*)(fFlowEvent[r]->GetTrack(iRPs));
-                        if (!iRP->InRPSelection()) continue;
-                        if (fabs(sTrack->GetIDDaughter(iDau)) == fabs(iRP->GetID()))
-                        {
-                           sTrack->SetDaughter(iDau, iRP);
-                           iRP->SetForRPSelection(kFALSE);
-                           if (0) printf("    á¶«daughter%d with fID %d was removed from this RP set\n", iDau, sTrack->GetIDDaughter(iDau));
-                        }
-                     }
-                  fFlowEvent[r]->AddTrack(sTrack);
+                  if (fDebug) printf(" was in RP set");
+                  iRP->SetForRPSelection(kFALSE);
+                  fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
                }
+            }
+            if (fDebug) printf("\n");
          }
+         cand->SetForPOISelection(kTRUE);
+         fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
       }
+      if (fDebug) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
+
       for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
       {
          for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++)
@@ -1522,12 +941,8 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
             PtSelector(2, un[nTracks], un[pTracks]);
          }
       }
-      //Push event averaged track values into flow members
-      FlowFinish(v2);
-      FlowFinishSin(v2Sin);
-
       PostData(1, fOutputList);
-      for (int m = 0; m != 30; ++m) PostData(2 + m, fFlowEvent[m]);
+      PostData(2, fFlowEvent);
    }
 }
 //_____________________________________________________________________________
@@ -1536,21 +951,37 @@ void AliAnalysisTaskPhiFlow::Terminate(Option_t *)
    // Terminate
 }
 //______________________________________________________________________________
-AliFlowCandidateTrack*  AliAnalysisTaskPhiFlow::MakeTrack(Double_t mass,
-      Double_t pt, Double_t phi, Double_t eta,
-      Int_t nDau, Int_t iID[]) const
+void  AliAnalysisTaskPhiFlow::MakeTrack(Double_t mass,
+                                              Double_t pt, Double_t phi, Double_t eta,
+                                              Int_t nDau, Int_t iID[]) const
 {
    // Consruct Flow Candidate Track from two selected candidates
-   AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
+   Bool_t overwrite = kTRUE;
+   AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
+   if (!sTrack)
+   {
+      sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
+      overwrite = kFALSE;
+   }
+   else sTrack->ClearMe();
    sTrack->SetMass(mass);
    sTrack->SetPt(pt);
    sTrack->SetPhi(phi);
    sTrack->SetEta(eta);
-   for (Int_t iDau = 0; iDau != nDau; ++iDau)
-      sTrack->AddDaughter(iID[iDau]);
+   for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
    sTrack->SetForPOISelection(kTRUE);
    sTrack->SetForRPSelection(kFALSE);
-   return sTrack;
+   if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
+   else fCandidates->AddLast(sTrack);
+   return;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::SetCommonConstants(Int_t massBins, Double_t minMass, Double_t maxMass) 
+{
+  // set common constants
+  fMassBins = massBins;
+  fMinMass = minMass;
+  fMaxMass = maxMass;
 }
 //_____________________________________________________________________________
 
index dd6488f..483816b 100644 (file)
@@ -10,7 +10,6 @@
 #ifndef ALIANALYSISTASKPHIFLOW_H
 #define ALIANALYSISTASKPHIFLOW_H
 
-
 class TH1F;
 class TH1F;
 class TH2F;
@@ -42,21 +41,13 @@ public:
    AliAnalysisTaskPhiFlow(const char *name);
    virtual ~AliAnalysisTaskPhiFlow();
 
-   enum PIDtype
-   {
-      kTPC = 0,
-      kCombined,
-   };
-
+   void                                 SetEnableDebugMode() {fDebug = kTRUE; };
    TH1F*                                BookHistogram(const char * name);
-   TH1F*                                BookDPhiPsiHistogram(const char * name);
    TH2F*                                BookPIDHistogram(const char * name);
    TH1F*                                InitPtSpectraHistograms(Int_t i);
-   TProfile*                            BookV2Profile(const char * name, Bool_t pt, Bool_t cos);
    TH1F*                                BookPtHistogram(const char* name);
    void                                 AddPhiIdentificationOutputObjects();
    virtual void                         UserCreateOutputObjects();
-   template <typename T> void           PairLoss(const T* track1, const T* track2) const;
    template <typename T> Double_t       InvariantMass(const T* track1, const T* track2) const;
    template <typename T> Double_t       DeltaDipAngle(const T* track1, const T* track2) const;
    template <typename T> Bool_t         CheckDeltaDipAngle(const T* track1, const T* track2) const;
@@ -66,8 +57,6 @@ public:
    void                                 SetMaxDeltaDipAngleAndPt(Float_t a, Float_t pt) { fDeltaDipAngle = a;
                                                                                           fDeltaDipPt = pt;
                                                                                           fApplyDeltaDipCut = kTRUE; };
-   void                                 EventPlanePtCut(Bool_t a = kFALSE) { fSetEventPlanePtCut = a; };
-   void                                 SetEventPlanePtCut(Float_t c) { fEventPlanePtCut = c; };
    template <typename T> Bool_t         EventCut(T* event);
    template <typename T> void           PlotVZeroMultiplcities(const T* event) const;
    template <typename T> Bool_t         CheckVertex(const T* event) const ;
@@ -81,71 +70,37 @@ public:
    Bool_t                               IsKaon(AliAODTrack* track) const;
    template <typename T> Double_t       PhiPt(const T* track_1, const T* track_2) const;
    template <typename T> void           PtSelector(Int_t _track_type, const T* track_1, const T* track_2) const;
-   void                                 EventPlane(const AliESDEvent* event);
-   void                                 EventPlane(const AliAODEvent* event);
-   void                                 EventPlaneResolution(Bool_t random,
-                                                        Double_t cossumQ1,
-                                                        Double_t sinsumQ1,
-                                                        Double_t cossumQ2,
-                                                        Double_t sinsumQ2,
-                                                        Double_t cossum,
-                                                        Double_t sinsum,
-                                                        Double_t cossumQwest,
-                                                        Double_t sinsumQwest,
-                                                        Double_t cossumQeast,
-                                                        Double_t sinsumQeast) const;
-   Double_t                             Chi(Double_t res) const;
-   Double_t                             ResolutionAsFunctionOfChi(Double_t chi) const;
-   Double_t                             EventPlaneStar(Bool_t random,
-                                                        Double_t cossumQ1,
-                                                        Double_t sinsumQ1,
-                                                        Double_t cossumQ2,
-                                                        Double_t sinsumQ2,
-                                                        Double_t cossum,
-                                                        Double_t sinsum,
-                                                        Double_t cossumQwest,
-                                                        Double_t sinsumQwest,
-                                                        Double_t cossumQeast,
-                                                        Double_t sinsumQeast) const;
-   template <typename T> void           EllipticFlow(Double_t* const v2, const T* track_1, const T* track_2) const;
-   template <typename T> Bool_t         EventPlaneTrack(T* track) const;
    template <typename T> Bool_t         PhiTrack(T* track) const;
-   void                                 FlowFinish(Double_t* const v2) const;
-   template <typename T> void           EllipticFlowSin(Double_t* const v2Sin, const T* track_1, const T* track_2) const;
-   void                                 FlowFinishSin(Double_t* const v2Sin) const;
    template <typename T> void           SetNullCuts(T* esd);
    void                                 PrepareFlowEvent(Int_t iMulti);
    virtual void                         UserExec(Option_t *option);
    virtual void                         Terminate(Option_t *);
-   void                                 SetRPCuts(AliFlowTrackCuts *cutsRP) { fRPCuts = cutsRP; }
    void                                 SetPOICuts(AliFlowTrackCuts *cutsPOI) { fPOICuts = cutsPOI; }
-   void                                 SetIdentificationType(PIDtype type) { fPIDtype = type; }
+   void                                 SetRPCuts(AliFlowTrackCuts *cutsRP) { fCutsRP = cutsRP; }
    void                                 SetRequireStrictKaonCuts() { fStrictKaonCuts= kTRUE; }
    void                                 SetRequireTPCStandAloneKaons() { fRequireTPCStandAlone = kTRUE; }
    void                                 SetOlderThanPbPbPass2() { fOldTrackParam = kTRUE; }
-   void                                 SetBayesianProbability(Double_t prob) { fParticleProbability = prob; }
-   void                                 SetMassRanges(Double_t flowBands[2][30]) { for (int i = 0; i != 2; ++i) for (int j = 0; j != 30; ++j) fFlowBands[i][j] = flowBands[i][j]; }
-   void                                 SetEtaRanges(Double_t mina, Double_t maxa, Double_t minb, Double_t maxb) { 
-                                                                                        fEtaMinA = mina;
-                                                                                        fEtaMaxA = maxa; 
-                                                                                        fEtaMinB = minb; 
-                                                                                        fEtaMaxB = maxb; }
-   void                                 SetFlowRPCuts(AliFlowTrackCuts *cutsRP) { fCutsRP = cutsRP; }
+   void                                 SetPIDConfiguration(Double_t prob[7]) { for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = prob[i]; }
    void                                 SetCandidateEtaAndPt(Double_t mineta, Double_t maxeta, Double_t minpt, Double_t maxpt) { fCandidateMinEta = mineta, 
                                                                                                                                fCandidateMaxEta = maxeta,
                                                                                                                                fCandidateMinPt = minpt, 
                                                                                                                                fCandidateMaxPt = maxpt,
                                                                                                                                fCandidateEtaPtCut = kTRUE;}
+   void                                 SetCommonConstants(Int_t massBins, Double_t minMass, Double_t maxMass);
 
 private:
 
+   Bool_t               fDebug; //! enable debug mode
    Bool_t               fAODAnalysis; // set aod analysis
-   Double_t             fFlowBands[2][30]; // Array containing the boundaries (upper , lower) of invariant mass bands, gev / c^2
-   Double_t             fEtaMinA, fEtaMaxA, fEtaMinB, fEtaMaxB; // upper and lower bounds of included eta regions
-   AliFlowTrackCuts     *fCutsRP; // track cuts for reference particles (event plane method)
+   Int_t                fMassBins; // mass bins
+   Double_t             fMinMass; // mass range
+   Double_t             fMaxMass; // mass range
+   AliFlowTrackCuts     *fCutsRP; // track cuts for reference particles
    AliFlowTrackCuts     *fNullCuts; // dummy cuts for flow event tracks
-   AliFlowEvent         *fFlowEvent[30]; //! flow events (one for each inv mass band)
+   AliPIDResponse       *fPIDResponse; //! pid response object
+   AliFlowEvent         *fFlowEvent; //! flow events (one for each inv mass band)
    AliFlowBayesianPID   *fBayesianResponse; //!PID response object
+   TObjArray            *fCandidates; // candidate array
    Bool_t               fOldTrackParam; // set to true if data is older than pbpb pass 2 production
    Bool_t               fRequireTPCStandAlone; // set TPC standalone cut for kaon selection
    Bool_t               fStrictKaonCuts; // require strict kaon cuts
@@ -154,7 +109,7 @@ private:
    Double_t             fCandidateMaxEta; // maximum eta for candidates
    Double_t             fCandidateMinPt; // minimum pt for candidates
    Double_t             fCandidateMaxPt; // maximum pt for candidates
-   Double_t             fParticleProbability; // set cutoff for bayesian probability
+   Double_t             fPIDConfig[7]; // set cutoff for bayesian probability
    Double_t             fCentrality; // event centrality
    AliESDEvent          *fESD;    //! ESD object
    AliAODEvent          *fAOD;    //! AOD oject
@@ -164,7 +119,6 @@ private:
    TH1F                 *fCentralityNoPass; //! QA histogram of events that do not pass centrality cut
    TH2F                 *fNOPID;//! QA histogram of TPC response of all charged particles
    TH2F                 *fPIDk;//! QA histogram of TPC response of kaons
-   TH2F                 *fPIDDeltaDip; //! QA histogram of TPC signal after delta dip cut
    TH1F                 *fInvMNP03; //! Invariant mass of unlike sign kaon pairs
    TH1F                 *fInvMPP03; //! Invariant mass of like sign (++) kaon pairs
    TH1F                 *fInvMNN03; //! Invariant mass of like sign (--) kaon pairs
@@ -219,83 +173,24 @@ private:
    TH1F                 *fInvMNP6570; //! Invariant mass of unlike sign kaon pairs
    TH1F                 *fInvMPP6570; //! Invariant mass of like sign (++) kaon pairs
    TH1F                 *fInvMNN6570; //! Invariant mass of like sign (--) kaon pairs
-   TH1F                 *fDeltaPhiPsiNP03; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP36; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP69; //!  Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP912; //!  Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP1215; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP1518; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP1821; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP2124; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP2427; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP2730; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP3035; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP3540; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP4045; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP4550; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP5055; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP5560; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP6065; //! Delta Phi Psi of unlike sign kaon pairs
-   TH1F                 *fDeltaPhiPsiNP6570; //! Delta Phi Psi of unlike sign kaon pairs
-   TProfile             *fProfV2; //! charged particle v2 (cos terms, event plane method)
-   TProfile             *fProfV2Sin; //! charged particle v2 (sin terms, event plane method)
-   TProfile             *fProfV2InvM03; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM36; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM69; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM912; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM1215; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM1518; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM1821; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM2124; //! cos component of v2 (eevent plane method)
-   TProfile             *fProfV2InvM2427; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM2730; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM3035; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM3540; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM4045; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM4550; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM5055; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM5560; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM6065; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2InvM6570; //! cos component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM03; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM36; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM69; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM912; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM1215; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM1518; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM1821; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM2124; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM2427; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM2730; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM3035; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM3540; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM4045; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM4550; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM5055; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM5560; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM6065; //! sin component of v2 (event plane method)
-   TProfile             *fProfV2SinInvM6570; //! sin component of v2 (event plane method)
-   TH1F                 *fPtSpectra03; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra36; //! Kaon Pt spectrum 
-   TH1F                 *fPtSpectra69; //!  Kaon Pt spectrum
-   TH1F                 *fPtSpectra912; //!  Kaon Pt spectrum
-   TH1F                 *fPtSpectra1215; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra1518; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra1821; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra2124; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra2427; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra2730; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra3035; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra3540; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra4045; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra4550; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra5055; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra5560; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra6065; //! Kaon Pt spectrum
-   TH1F                 *fPtSpectra6570; //! Kaon Pt spectrum
-   TH1F                 *fEventPlaneSTAR; //! Distribution of the orientation of the event plane, taken from STAR
-   TProfile             *fEventPlaneResolutionRandom; //! event plane resolution from randomized subevents (event plane method)
-   TProfile             *fEventPlaneResolutionEta; //! event plane resolution from eta separated subevents (event plane method)
+   TH1F                 *fPtSpectra03; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra36; //! kaon pair Pt spectrum 
+   TH1F                 *fPtSpectra69; //!  kaon pair Pt spectrum
+   TH1F                 *fPtSpectra912; //!  kaon pair Pt spectrum
+   TH1F                 *fPtSpectra1215; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra1518; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra1821; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra2124; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra2427; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra2730; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra3035; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra3540; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra4045; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra4550; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra5055; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra5560; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra6065; //! kaon pair Pt spectrum
+   TH1F                 *fPtSpectra6570; //! kaon pair Pt spectrum
    TH1F                 *fPtP; //! QA histogram of p_t distribution of positive particles
    TH1F                 *fPtN; //! QA histogram of p_t distribution of negative particles
    TH1F                 *fPtKP; //! QA histogram of p_t distribution of positive kaons
@@ -303,35 +198,27 @@ private:
    Double_t             fCentralityMin; // lower bound of cenrality bin
    Double_t             fCentralityMax; // upper bound of centrality bin
    const char           *fkCentralityMethod; // method used to determine centrality (V0 by default)
-   AliFlowTrackCuts     *fRPCuts; // cuts for reference paricles (flow package)
-   AliFlowTrackCuts     *fPOICuts; // cus for particles of interest (flow package)
+   AliFlowTrackCuts     *fPOICuts; // cuts for particles of interest (flow package)
    Float_t              fVertexRange; // absolute value of maximum distance of vertex along the z-axis
-   Double_t             fQx; // dummy variable for the Q_x vector (event plane method)
-   Double_t             fQy; // dummy variable for the Q_y vector (event plane method)
-   TH1F                 *fEventPlane; //! Distribution of the orientation of the event palen (event plane method)
    TH1F                 *fPhi; //! QA plot of azimuthal distribution of tracks used for event plane estimation
    TH1F                 *fPt; //! QA plot of p_t sectrum of tracks used for event plane estimation
    TH1F                 *fEta; //! QA plot of eta distribution of tracks used for event plane estimation
    TH1F                 *fVZEROA; //! QA plot vzeroa multiplicity (all tracks in event)
    TH1F                 *fVZEROC; //! QA plot vzeroc multiplicity (all tracks in event)
    TH1F                 *fTPCM; //! QA plot TPC multiplicity (tracks used for event plane estimation)
-   PIDtype              fPIDtype; // selector to indicate which PID process will be used
    Float_t              fDeltaDipAngle; // absolute value of delta dip angle to be excluded
    Float_t              fDeltaDipPt; // upper value of pt range in which delta dip angle must be applied
    Bool_t               fApplyDeltaDipCut; // enforce delta dip cut
-   TH2F                 *fPairLoss; //! pair loss histo
-   Float_t              fEventPlanePtCut; // Pt cut on tracks used for event plane estimation
-   Bool_t               fSetEventPlanePtCut; // kTRUE for Pt cut on event plane estimation tracks
 
    AliAnalysisTaskPhiFlow(const AliAnalysisTaskPhiFlow&); // Not implemented
    AliAnalysisTaskPhiFlow& operator=(const AliAnalysisTaskPhiFlow&); // Not implemented
 
-   AliFlowCandidateTrack*  MakeTrack(Double_t, Double_t, Double_t, Double_t, Int_t , Int_t[]) const;
-
-   ClassDef(AliAnalysisTaskPhiFlow, 2);
+   void                 MakeTrack(Double_t, Double_t, Double_t, Double_t, Int_t , Int_t[]) const;
+     ClassDef(AliAnalysisTaskPhiFlow, 2);
 
 };
 
 #endif
 
 
+
index 9d50f78..da805f5 100644 (file)
 //
 // AddTaskPhiFlow macro
 // Author: Redmer A. Bertens, Utrecht University, 2012
+//         Many thanks to Carlos Perez Lara
+//         Commented where necessary
 /////////////////////////////////////////////////////////////////////////////////////////////
+
+class AliAnalysisDataContainer;
+class AliFlowTrackCuts;
+class AliFlowTrackSimpleCuts;
+class AliFlowEventCuts;
+class AliFlowEventSimpleCuts;
+class AliAnalysisDataContainer;
+
 AliAnalysisTaskPhiFlow* AddTaskPhiFlow(Bool_t SP = kTRUE,
+                                       Bool_t SPSUB = kTRUE,
                                        Bool_t QC = kTRUE,
                                        Bool_t EP = kTRUE,
-                                       Int_t harmonic = 2,
-                                       Float_t centrMin = 30.,
-                                       Float_t centrMax = 40.,
-                                       Float_t bayesThresholdP = 0.5,
+                                       Float_t centrMin = 20.,
+                                       Float_t centrMax = 30.,
+                                       Double_t PIDconfig[7],
                                        TString suffixName = "",
                                        Bool_t bCentralTrigger = kFALSE,
-                                       Float_t RPEtaMinA = -0.8,
-                                       Float_t RPEtaMaxA = 0.0,
-                                       Float_t RPEtaMinB = 0.0,
-                                       Float_t RPEtaMaxB = 0.8,
+                                       Float_t EtaGap = 0.,
                                        Float_t POIEtaMin = -0.8,
                                        Float_t POIEtaMax = 0.8,
                                        Float_t POIPtMin = 0.15,
                                        Float_t POIPtMax = 10.,
-                                       Float_t deltaDip = 0.04,
-                                       Float_t deltaDipMaxPt = 1.2,
-                                       AliAnalysisTaskPhiFlow::PIDtype pidtype = AliAnalysisTaskPhiFlow::kCombined,
-                                       Bool_t strictKaonCuts = kFALSE,
+                                       Float_t deltaDip = 0.,
+                                       Float_t deltaDipMaxPt = 0.,
+                                       Float_t DCA = 0.3,
+                                       Int_t harm = 2,
                                        Bool_t TPCStandAloneTracks = kFALSE,
-                                       Float_t vertexZ = 10.)
+                                       Bool_t useGlobalRPCuts = kTRUE,
+                                       Float_t vertexZ = 10.,
+                                       Bool_t shrinkSP = kTRUE,
+                                       Bool_t debug = kTRUE)
 {
+   // main function, create and add tasks
+   if(debug) cout << " === Adding Task PhiFlow === " << endl;
    // set up main output container's name
-   TString centralityName("");
-   centralityName += Form("%.0f", centrMin);
-   centralityName += "-";
-   centralityName += Form("%.0f", centrMax);
-   centralityName += "_";
-   centralityName += Form("vZ%.f", vertexZ);
-   centralityName += "_";
-   centralityName += Form("bayP%.1f", bayesThresholdP);
-   centralityName += "_";
-   centralityName += "_POIEta";
-   centralityName += Form("%.1f", POIEtaMin);
-   centralityName += "-";
-   centralityName += Form("%.1f", POIEtaMax);
-   centralityName += "_RPEta";
-   centralityName += Form("%.1f", RPEtaMinA);
-   centralityName += "-";
-   centralityName += Form("%.1f", RPEtaMaxA);
-   centralityName += "-";
-   centralityName += Form("%.1f", RPEtaMinB);
-   centralityName += "-";
-   centralityName += Form("%.1f", RPEtaMaxB);
-   centralityName += "-";
-   centralityName += Form("dDip%.2f", deltaDip);
-   centralityName += "-";
-   centralityName += Form("dDipPt%.2f", deltaDipMaxPt);
-   if (strictKaonCuts)
-   {
-      centralityName += "-";
-      centralityName += "StrictKaonCuts";
-   }
-   if (bCentralTrigger)
-   {
-       centralityName += "-";
-       centralityName += "kMBkCkSC";
-   }
-
    TString fileName = AliAnalysisManager::GetCommonFileName();
    fileName += ":PhiV2FlowEvents";
-
-   // get the manager
+   if(debug) cout << "    --> Reconstruction data container: " << fileName << endl;
+   // check validity of arguments
+   if((!SPSUB)&&(EtaGap > 0.)) {
+       if(debug) cout << " --> Fatal error: Eta gap is introduced but method SPSUB is not enabled !!! <-- " << endl;
+       return 0x0;
+   }
+   else if ((QC||SP||EP)&&(EtaGap > 0.)) {
+       if(debug) cout << " --> Fatal error: one or more of the flow analyses is not compatible with the Eta Gap !!! <--" << endl;
+       return 0x0;
+   }
+   // get the manager and input event handler
    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
-   if (!mgr)
-   {
-      ::Error("AddTaskPhiV2", "No analysis manager to connect to.");
-      return NULL;
+   if (!mgr) {
+      if(debug) cout << " Fatal error: no analysis manager found! " << endl;
+      return 0x0;
    }
-   if (!mgr->GetInputEventHandler())
-   {
-      ::Error("AddTaskPhiV2", "This task requires an input event handler");
-      return NULL;
+   if (!mgr->GetInputEventHandler()) {
+      if(debug) cout << " Fatal error: no imput event handler found!" << endl;
+      return 0x0;
    }
-
    // create the main task
    AliAnalysisTaskPhiFlow *task = new AliAnalysisTaskPhiFlow("TaskPhiFlow");
-
-    // set triggers
-   if(bCentralTrigger) task->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);
-   else                task->SelectCollisionCandidates(AliVEvent::kMB);
-
-   // set RP cuts for EP method
-   AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts("GlobalRP");
-   AliFlowTrackCuts* cutsRP = cutsRP->GetStandardTPCStandaloneTrackCuts();
-
+   if(debug) cout << " === AliAnalysisTaskPhiFlow === " << task << endl;
+   if(!task) {
+       if(debug) cout << " --> Unexpected error occurred: NO TASK WAS CREATED! (could be a library problem!) " << endl;
+       return 0x0;
+   }
+   // set triggers
+   if (bCentralTrigger) task->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);
+   else                 task->SelectCollisionCandidates(AliVEvent::kMB);
+   if(debug) cout << "    --> Set trigger selection to ";
+   if(debug&&bCentralTrigger) cout << " kMB, kCentral, kSemiCentral " << endl;
+   if(debug&&(!bCentralTrigger)) cout << " kMB " << endl;
+   //set RP cuts for flow package analysis
+   cutsRP = new AliFlowTrackCuts("RFPcuts");
+   if(!cutsRP) {
+       if(debug) cout << " Fatal error: no RP cuts found, could be a library problem! " << endl;
+       return 0x0;
+   }
+   if(useGlobalRPCuts) {
+       AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;
+       cutsRP->SetParamType(rptype);
+       cutsRP->SetPtRange(0.2, 5.0);
+       cutsRP->SetEtaRange(-0.8, 0.8);
+       cutsRP->SetMinNClustersTPC(70);
+       cutsRP->SetMinChi2PerClusterTPC(0.1);
+       cutsRP->SetMaxChi2PerClusterTPC(4.0);
+       cutsRP->SetRequireTPCRefit(kTRUE);
+       cutsRP->SetMaxDCAToVertexXY(0.3);
+       cutsRP->SetMaxDCAToVertexZ(0.3);
+       cutsRP->SetAcceptKinkDaughters(kFALSE);
+       cutsRP->SetMinimalTPCdedx(10.);
+       if(debug) cout << "    --> kGlobal RP's " << cutsRP << endl;
+   }
+   if(!useGlobalRPCuts) {
+       AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kTPCstandalone;
+       cutsRP->SetParamType(rptype);
+       cutsRP->SetPtRange(0.2,5.);
+       cutsRP->SetEtaRange(-0.8,0.8);
+       cutsRP->SetMinNClustersTPC(70);
+       cutsRP->SetMinChi2PerClusterTPC(0.2);
+       cutsRP->SetMaxChi2PerClusterTPC(4.0);
+       cutsRP->SetMaxDCAToVertexXY(3.0);
+       cutsRP->SetMaxDCAToVertexZ(3.0);
+       cutsRP->SetDCAToVertex2D(kTRUE);
+       cutsRP->SetAcceptKinkDaughters(kFALSE);
+       cutsRP->SetMinimalTPCdedx(10.);
+       if(debug) cout << "    --> kTPCStandAlone RP's " << cutsRP << endl;
+   }
+   task->SetRPCuts(cutsRP);
    // set POI cuts for kaon selection
    AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts("GlobalPOI");
+   if(!cutsPOI) {
+       if(debug) cout << " Fatal error: no POI cuts (could be a library problem)!" << endl;
+       return 0x0;
+   }
    AliFlowTrackCuts* cutsPOI = cutsPOI->GetStandardGlobalTrackCuts2010();
-
-   // set some more cuts ...
-   task->EventPlanePtCut(kFALSE); // set to true to enforce pt cut on event plane tracks
-   task->SetMaxDeltaDipAngleAndPt(deltaDip, deltaDipMaxPt); // set an upper value for d-dip angle and pt of tracks on which this should be applied
-   task->SetCentralityParameters(centrMin, centrMax, "V0M");
-   task->SetVertexZ(vertexZ);
-   task->SetRPCuts(cutsRP);
+   cutsPOI->SetPtRange(0.2, 5); //?
+   cutsPOI->SetMaxDCAToVertexXY(DCA);
+   cutsPOI->SetMaxDCAToVertexZ(DCA);
+   if(debug) cout << "    --> cutsPOI " << cutsPOI << endl;
    task->SetPOICuts(cutsPOI);
-
-   // set the kaon cuts, and specify the PID procedure which will be used
-   task->SetIdentificationType(pidtype);
-   task->SetBayesianProbability(bayesThresholdP);
-   if (strictKaonCuts) task->SetRequireStrictKaonCuts();
-   if (TPCStandAloneTracks) task->SetRequireTPCStandAloneKaons();
-   task->SetCandidateEtaAndPt(POIEtaMin, POIEtaMax, POIPtMin, POIPtMax);
-
-   //set RP cuts for flow package analysis
-   cutsFlowRP = new AliFlowTrackCuts("rp_cuts_111");
-   AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;
-   cutsFlowRP->SetParamType(rptype);
-   cutsFlowRP->SetPtRange(0.2, 5.0);
-   cutsFlowRP->SetEtaRange(-0.8, 0.8);
-   cutsFlowRP->SetMinNClustersTPC(70);
-   cutsFlowRP->SetMinChi2PerClusterTPC(0.1);
-   cutsFlowRP->SetMaxChi2PerClusterTPC(4.0);
-   cutsFlowRP->SetRequireTPCRefit(kTRUE);
-   cutsFlowRP->SetMaxDCAToVertexXY(0.3);
-   cutsFlowRP->SetMaxDCAToVertexZ(0.3);
-   cutsFlowRP->SetAcceptKinkDaughters(kFALSE);
-   cutsFlowRP->SetMinimalTPCdedx(10.);
-
+   // POI filter cuts, will filter invm mass bands and subevents
+   AliFlowTrackSimpleCuts* POIfilterQC[30];
+   AliFlowTrackSimpleCuts* POIfilterSP[30][2];
    Double_t flowBands[2][30];
-   for (Int_t i = 0; i < 30; i++)
-   {
-      flowBands[0][i] = 0.99 + i * 0.0034;
-      flowBands[1][i] = 0.99 + (i + 1) * 0.0034;
+   for (Int_t mb = 0; mb < 30; mb++) {
+      flowBands[0][mb] = 0.99 + mb * 0.0034;
+      flowBands[1][mb] = 0.99 + (mb + 1) * 0.0034;
+      POIfilterSP[mb][0] = new AliFlowTrackSimpleCuts(Form("FilterPOISP_MB%d_ETANEG", mb));
+      if(!POIfilterSP[mb][0]) {
+          if(debug) cout << " FAILED to initialize POI filters, could be a library problem!" << endl;
+          return 0x0;
+      }
+      POIfilterSP[mb][0]->SetEtaMin(-0.8);
+      POIfilterSP[mb][0]->SetEtaMax(0.);
+      POIfilterSP[mb][0]->SetMassMin(flowBands[0][mb]);
+      POIfilterSP[mb][0]->SetMassMax(flowBands[1][mb]);
+      POIfilterSP[mb][1] = new AliFlowTrackSimpleCuts(Form("FilterPOISP_MB%d_ETAPOS", mb));
+      POIfilterSP[mb][1]->SetEtaMin(0.);
+      POIfilterSP[mb][1]->SetEtaMax(+0.8);
+      POIfilterSP[mb][1]->SetMassMin(flowBands[0][mb]);
+      POIfilterSP[mb][1]->SetMassMax(flowBands[1][mb]);
+      POIfilterQC[mb] = new AliFlowTrackSimpleCuts(Form("FilterPOIQC_MB%d", mb));
+      POIfilterQC[mb]->SetEtaMin(-0.8);
+      POIfilterQC[mb]->SetEtaMax(+0.8);
+      POIfilterQC[mb]->SetMassMin(flowBands[0][mb]);
+      POIfilterQC[mb]->SetMassMax(flowBands[1][mb]);
    }
-   task->SetMassRanges(flowBands);
-   task->SetEtaRanges(RPEtaMinA, RPEtaMaxA, RPEtaMinB, RPEtaMaxB);
-   task->SetFlowRPCuts(cutsFlowRP);
-   mgr->AddTask(task);
-
-   AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
-   AliAnalysisDataContainer *coutput = mgr->CreateContainer(Form("PhiV2_%s", centralityName.Data()), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
-
-   AliAnalysisDataContainer *coutputCandidates[30];
-
-   for (int r = 0; r != 30; ++r)
-   {
-      coutputCandidates[r] = mgr->CreateContainer(
-                                Form("Flow_MassBand_%d", r),
-                                AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer);
-      mgr->ConnectOutput(task, 2 + r, coutputCandidates[r]);
+   if(debug) cout << "    --> Created poi filters " << endl;
+   task->SetCommonConstants(30, flowBands[0][0], flowBands[1][29]);
+   if(debug) cout << "    --> Set common constants " << endl;
+   // set pair and event cuts
+   if((deltaDip>0.005)&&(deltaDipMaxPt>0.005)) task->SetMaxDeltaDipAngleAndPt(deltaDip, deltaDipMaxPt);
+   else cout << " --> Disabled Delta-Dip exclusion. <-- " << endl;
+   task->SetCandidateEtaAndPt(POIEtaMin, POIEtaMax, POIPtMin, POIPtMax);
+   task->SetCentralityParameters(centrMin, centrMax, "V0M");
+   task->SetVertexZ(vertexZ);
+   if(debug) cout << "    --> Set pair cuts and event cuts" << endl;
+   // set the kaon cuts, and specify the PID procedure which will be used
+   task->SetPIDConfiguration(PIDconfig);
+   if(debug) {
+       cout << " ************ Configured PID routine ************ " << endl;
+       cout << "    0 < " << PIDconfig[1] << " p_t, ITS || TPC with s < " << PIDconfig[0] << endl;
+       if(PIDconfig[2] < 0.) cout << "  --> TPC control disabled " << endl;
+       if(PIDconfig[2] > 0.) cout << "  --> TPC control enabled " << endl;
+       cout << "    " << PIDconfig[1] << " < " << PIDconfig[4] << " p_t, TPC || ITS with s < " << PIDconfig[3] << endl;
+       if(PIDconfig[5] < 0.) cout << "  --> ITS control disabled " << endl;
+       if(PIDconfig[5] > 0.) cout << "  --> ITS control enabled " << endl;
+       cout << "    " << PIDconfig[4] << " < 7 p_t, TPC / TOF Bayesian with p < " << PIDconfig[6] << endl;
+       cout << " ************************************************ " << endl;
    }
-
-   // set, if necessary, a suffixname as a unique identifier for each wagon
-   if(suffixName == "") {
-       suffixName += Form("%.0f", centrMin);
-       suffixName += Form("%.0f", centrMax);
+   if (TPCStandAloneTracks)
+   { // switch to tpc standalone tracks for POI selection
+      if(debug) cout << "    --> Switching to TPC standalone analsis " << endl;
+      task->SetRequireStrictKaonCuts();
+      task->SetRequireTPCStandAloneKaons();
    }
-
-   if (SP) // set up scalar product (SP) tasks
-   {
-      AliAnalysisTaskScalarProduct *taskSP[30];
-      AliAnalysisDataContainer *coutputSP[30];
-      for (int r = 0; r != 30; ++r)
-      {
-         taskSP[r] = new AliAnalysisTaskScalarProduct(Form("SP_MassBand_%d_%s", r, suffixName.Data()), kFALSE);
-         taskSP[r]->SetHarmonic(harmonic);
-         taskSP[r]->SetRelDiffMsub(1.0);
-         taskSP[r]->SetApplyCorrectionForNUA(kTRUE);
-         mgr->AddTask(taskSP[r]);
-         coutputSP[r] = mgr->CreateContainer(Form("cobjSP_MassBand_%d_%s", r, suffixName.Data()), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
-         mgr->ConnectInput(taskSP[r], 0, coutputCandidates[r]);
-         mgr->ConnectOutput(taskSP[r], 1, coutputSP[r]);
-      }
-   }
-
-   if (EP) // set up event plane (ep) tasks
+   // create and configure IO containers
+   AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+   AliAnalysisDataContainer *coutHist = mgr->CreateContainer(Form("PhiV2_%s", OutputName(centrMin, centrMax,PIDconfig, suffixName, bCentralTrigger, EtaGap, POIEtaMin, POIEtaMax, POIPtMin, POIPtMax, deltaDip, deltaDipMaxPt, DCA, harm, TPCStandAloneTracks, vertexZ, debug, useGlobalRPCuts).Data()), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
+   if(debug) cout << "    --> Created IO containers " << cinput << ", " << coutHist << endl;
+   mgr->AddTask(task);
+   if(debug) cout << " === ADDING MAIN TASK == " << endl;
+   mgr->ConnectInput(task, 0, cinput);
+   mgr->ConnectOutput(task, 1, coutHist);
+   if(debug) cout << "    --> Connected IO containers " << endl;
+   if (SP || EP || QC || SPSUB) // if flow analysis should be done after reconstruction
    {
-      AliAnalysisTaskScalarProduct *taskEP[30];
-      AliAnalysisDataContainer *coutputEP[30];
-      for (int r = 0; r != 30; ++r)
+      if(debug) cout << " === RECEIVED REQUEST FOR FLOW ANALYSIS === " << endl;
+      AliAnalysisDataContainer *flowEvent = mgr->CreateContainer("FlowContainer", AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer);
+      mgr->ConnectOutput(task, 2, flowEvent);
+      if(debug) cout << "    --> Created IO containers " << flowEvent << endl;
+      // set a suffixname as a unique identifier for each wagon
+      if (suffixName == "")
       {
-         taskEP[r] = new AliAnalysisTaskScalarProduct(Form("EP_MassBand_%d_%s", r, suffixName.Data()), kFALSE);
-         taskEP[r]->SetBehaveAsEP();
-         taskEP[r]->SetHarmonic(harmonic);
-         taskEP[r]->SetRelDiffMsub(1.0);
-         taskEP[r]->SetApplyCorrectionForNUA(kTRUE);
-         mgr->AddTask(taskEP[r]);
-         coutputEP[r] = mgr->CreateContainer(Form("cobjEP_MassBand_%d_%s", r, suffixName.Data()), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
-         mgr->ConnectInput(taskEP[r], 0, coutputCandidates[r]);
-         mgr->ConnectOutput(taskEP[r], 1, coutputEP[r]);
+         suffixName += Form("%.0f", centrMin);
+         suffixName += Form("%.0f", centrMax);
       }
-   }
-
-   if (QC) // set up q-cumulants task
-   {
-      AliAnalysisTaskQCumulants *taskQC[30];
-      AliAnalysisDataContainer *coutputQC[30];
-      for (int r = 0; r != 30; ++r)
-      {
-         taskQC[r] = new AliAnalysisTaskQCumulants(Form("QC_MassBand_%d_%s", r, suffixName.Data()), kFALSE);
-         taskQC[r]->SetHarmonic(harmonic);
-         taskQC[r]->SetCalculateCumulantsVsM(kFALSE);
-         taskQC[r]->SetnBinsMult(10000);
-         taskQC[r]->SetMinMult(0.);
-         taskQC[r]->SetMaxMult(10000.);
-         taskQC[r]->SetApplyCorrectionForNUA(kTRUE);
-         taskQC[r]->SetFillMultipleControlHistograms(kFALSE);
-         mgr->AddTask(taskQC[r]);
-         coutputQC[r] = mgr->CreateContainer(Form("cobjQC_MassBand_%d_%s", r, suffixName.Data()), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
-         mgr->ConnectInput(taskQC[r], 0, coutputCandidates[r]);
-         mgr->ConnectOutput(taskQC[r], 1, coutputQC[r]);
+      if(debug) cout << "    --> suffixName " << suffixName << endl;
+      for (int mb = 0; mb != 30; ++mb) {
+         if (QC) {  // add qc tasks
+            AddQCmethod(Form("QCTPCMB_%d_%s", mb, suffixName.Data()), harm, flowEvent, POIfilterQC[mb], debug);
+            if(debug) cout << "    --> Hanging QC task ... " << mb << " succes! "<< endl;
+         }
+         if (SPSUB) {  // add sp subevent tasks
+            AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -.5*EtaGap, .5*EtaGap, +0.8, "Qa", harm, flowEvent, POIfilterSP[mb][1], NULL, false, shrinkSP, debug, true);
+            if(debug) cout << "    --> Hanging SP Qa task " << mb << " succes!" << endl;
+            AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -.5*EtaGap, .5*EtaGap, +0.8, "Qb", harm, flowEvent, POIfilterSP[mb][0], NULL, false, shrinkSP, debug, true);
+            if(debug) cout << "    --> Hanging SP Qb task ..." << mb << " succes!"<< endl;
+         }
+         if (SP) { // add sp tasks
+            AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -0.0, +0.0, +0.8, "QaQb", harm, flowEvent, POIfilterQC[mb], NULL, false, shrinkSP, debug);
+            if(debug) cout << "    --> Hanging SP task ... " << mb << " succes!" << endl;
+         }
+         if (EP) { // add ep tasks
+            AddSPmethod(Form("EPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -0.0, +0.0, +0.8, "QaQb", harm, flowEvent, POIfilterQC[mb], NULL, true, shrinkSP, debug);
+            if(debug) cout << "    --> Hanging EP task ... " << mb << " succes!" << endl;
+         }
       }
    }
-   mgr->ConnectInput(task, 0, cinput);
-   mgr->ConnectOutput(task, 1, coutput);
-
-   // return the task
+   // print summary to screen
+   cout << endl << endl << " ========= AddTaskPhiFlow launched succesfully - SUMMARY: ========== " << endl;
+   cout << " ************ Configured PID routine ************ " << endl;
+   cout << "      0 < " << PIDconfig[1] << " p_t, ITS || TPC with s < " << PIDconfig[0] << endl;
+   if(PIDconfig[2] < 0.) cout << "    --> TPC control disabled " << endl;
+   if(PIDconfig[2] > 0.) cout << "    --> TPC control enabled " << endl;
+   cout << "    " << PIDconfig[1] << " < " << PIDconfig[4] << " p_t, TPC || ITS with s < " << PIDconfig[3] << endl;
+   if(PIDconfig[5] < 0.) cout << "    --> ITS control disabled " << endl;
+   if(PIDconfig[5] > 0.) cout << "    --> ITS control enabled " << endl;
+   cout << "      " << PIDconfig[4] << " < 7 p_t, TPC / TOF Bayesian with p < " << PIDconfig[6] << endl;
+   cout <<   " ************************************************ " << endl << endl;
+   cout << " ************* Task statistics ****************** " << endl;
+   if(SP) cout << "   --> Launched 30 QaQb SP filters and corresponding 30 SP tasks " << endl;
+   if(EP) cout << "   --> Launched 30 QaQb QC filters and corresponding 30 EP tasks " << endl;
+   if(SPSUB) cout << "   --> Launched 30+30 Qa&&Qb SP filters and corresponding 30+30 SP tasks " << endl;
+   if(QC) cout << "   --> Launched 30 QaQb QC filters and corresponding 30 QC tasks " << endl;
+   cout << " ************************************************** " << endl;
+   cout << "           --> Now go for a coffee! <-- " << endl;
    return task;
-};
+}
+//_____________________________________________________________________________
+void AddSPmethod(char *name, double minEtaA, double maxEtaA, double minEtaB, double maxEtaB, char *Qvector, int harmonic, AliAnalysisDataContainer *flowEvent, AliFlowTrackSimpleCuts *cutsPOI = NULL, AliFlowTrackSimpleCuts *cutsRFP = NULL, bool bEP, bool shrink = false, bool debug, bool etagap = false)
+{
+   // add sp task and invm filter tasks
+   if(debug) (bEP) ? cout << " ****** Reveived request for EP task ****** " << endl : cout << " ******* Switching to SP task ******* " << endl;
+   TString fileName = AliAnalysisManager::GetCommonFileName();
+   (bEP) ? fileName+=":EP" : fileName+=":SP";
+   if(etagap) {
+       fileName+="_SUBEVENTS";
+       if(debug) cout << "    --> Setting up subevent analysis <-- " << endl;
+   }
+   if(debug) cout << "    --> fileName " << fileName << endl;
+   TString myFolder = fileName;
+   if(debug) cout << "    --> myFolder " << myFolder << endl;
+   TString myNameSP;
+   (bEP) ? myNameSP = Form("%sEPv%d%s", name, harmonic, Qvector): myNameSP = Form("%sSPv%d%s", name, harmonic, Qvector);
+   if(debug) cout << " myNameSP " << myNameSP << endl;
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   AliAnalysisDataContainer *flowEvent2 = mgr->CreateContainer(Form("Filter_%s", myNameSP.Data()), AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer);
+   AliAnalysisTaskFilterFE *tskFilter = new AliAnalysisTaskFilterFE(Form("TaskFilter_%s", myNameSP.Data()), cutsRFP, cutsPOI);
+   tskFilter->SetSubeventEtaRange(minEtaA, maxEtaA, minEtaB, maxEtaB);
+   mgr->AddTask(tskFilter);
+   mgr->ConnectInput(tskFilter, 0, flowEvent);
+   mgr->ConnectOutput(tskFilter, 1, flowEvent2);
+   AliAnalysisDataContainer *outSP = mgr->CreateContainer(myNameSP.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
+   AliAnalysisTaskScalarProduct *tskSP = new AliAnalysisTaskScalarProduct(Form("TaskScalarProduct_%s", myNameSP.Data()), kFALSE);
+   tskSP->SetApplyCorrectionForNUA(kTRUE);
+   tskSP->SetHarmonic(harmonic);
+   tskSP->SetTotalQvector(Qvector);
+   if (bEP) tskSP->SetBehaveAsEP();
+  if (shrink) tskSP->SetBookOnlyBasicCCH(kTRUE);
+   mgr->AddTask(tskSP);
+   mgr->ConnectInput(tskSP, 0, flowEvent2);
+   mgr->ConnectOutput(tskSP, 1, outSP);
+}
+//_____________________________________________________________________________
+void AddQCmethod(char *name, int harmonic, AliAnalysisDataContainer *flowEvent, AliFlowTrackSimpleCuts *cutsPOI = NULL, Bool_t debug, AliFlowTrackSimpleCuts *cutsRFP = NULL)
+{
+   // add qc task and invm filter tasks
+   if(debug) cout << " ****** Received request for QC v" << harmonic << " task " << name << ", POI " << cutsPOI << ", IO ****** " << flowEvent << endl;
+   TString fileName = AliAnalysisManager::GetCommonFileName();
+   fileName+=":QC";
+   if(debug) cout << "    --> Common filename: " << fileName << endl;
+   TString myFolder = Form("v%d", harmonic);
+   if(debug) cout << "    --> myFolder: " << myFolder << endl;
+   TString myName = Form("%s", name);
+   if(debug) cout << "    --> myName: " << myName << endl;
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   AliAnalysisDataContainer *flowEvent2 = mgr->CreateContainer(Form("Filter_%s", myName.Data()), AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
+   AliAnalysisTaskFilterFE *tskFilter = new AliAnalysisTaskFilterFE(Form("TaskFilter_%s", myName.Data()), cutsRFP, cutsPOI);
+   mgr->AddTask(tskFilter);
+   mgr->ConnectInput(tskFilter, 0, flowEvent);
+   mgr->ConnectOutput(tskFilter, 1, flowEvent2);
+   AliAnalysisDataContainer *outQC = mgr->CreateContainer(myName.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
+   AliAnalysisTaskQCumulants *tskQC = new AliAnalysisTaskQCumulants(Form("TaskQCumulants_%s", myName.Data()), kFALSE);
+   tskQC->SetApplyCorrectionForNUA(kTRUE);
+   tskQC->SetHarmonic(harmonic);
+   tskQC->SetBookOnlyBasicCCH(kTRUE);
+   mgr->AddTask(tskQC);
+   mgr->ConnectInput(tskQC, 0, flowEvent2);
+   mgr->ConnectOutput(tskQC, 1, outQC);
+}
+//_____________________________________________________________________________
+TString OutputName( Float_t centrMin,
+                    Float_t centrMax,
+                    Double_t PIDconfig[7],
+                    TString suffixName,
+                    Bool_t bCentralTrigger,
+                    Float_t EtaGap,
+                    Float_t POIEtaMin,
+                    Float_t POIEtaMax,
+                    Float_t POIPtMin,
+                    Float_t POIPtMax,
+                    Float_t deltaDip,
+                    Float_t deltaDipMaxPt,
+                    Float_t DCA,
+                    Int_t harm,
+                    Bool_t TPCStandAloneTracks,
+                    Float_t vertexZ,
+                    Bool_t debug,
+                    Bool_t useGlobalRPCuts)
+{
+   // generate output name
+   TString centralityName = ("");
+   centralityName += Form("%.0f", centrMin);
+   centralityName += Form("%.0f", centrMax);
+   centralityName += "_";
+   centralityName += Form("vZ%.f", vertexZ);
+   centralityName += "_";
+   for(Int_t i = 0; i < 7; i++) centralityName += Form("%.1f_", PIDconfig[i]);
+   centralityName += "_POIEta";
+   centralityName += Form("%.1f", POIEtaMin);
+   centralityName += Form("%.1f", POIEtaMax);
+   centralityName += "_gap";
+   centralityName += Form("%.1f", -0.5*EtaGap);
+   centralityName += Form("%.1f", 0.5*EtaGap);
+   centralityName += "_";
+   centralityName += Form("dDip%.2f", deltaDip);
+   centralityName += "-";
+   centralityName += Form("dDipPt%.2f", deltaDipMaxPt);
+   if (TPCStandAloneTracks) {
+      centralityName += "-";
+      centralityName += "TPCStandAloneTracks";
+   }
+   if (bCentralTrigger) {
+      centralityName += "-";
+      centralityName += "kMBkCkSC";
+   }
+   if (!useGlobalRPCuts) {
+       centralityName += "-";
+       centralityName += "TPCRP";
+   }
+   if(debug) cout << "    --> centralityName " << centralityName << endl;
+   return centralityName;
+}
+//_____________________________________________________________________________