package for phi meson flow analysis
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Feb 2012 13:20:05 +0000 (13:20 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Feb 2012 13:20:05 +0000 (13:20 +0000)
PWG/CMakelibPWGflowTasks.pkg
PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx [new file with mode: 0644]
PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.h [new file with mode: 0644]
PWG/PWGflowTasksLinkDef.h
PWGCF/FLOW/macros/AddTaskPhiFlow.C [new file with mode: 0644]

index c85bed5..2edf5f6 100644 (file)
@@ -49,6 +49,7 @@ set ( SRCS
     FLOW/Tasks/AliAnalysisTaskPIDflowQA.cxx 
     FLOW/Tasks/AliAnalysisTaskQAPmdflow.cxx 
     FLOW/Tasks/AliFlowBayesianPID.cxx
+    FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
diff --git a/PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx b/PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx
new file mode 100644 (file)
index 0000000..14a5161
--- /dev/null
@@ -0,0 +1,1371 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+// AliAnalysisTaskPhiFlow:
+// author: Redmer Alexander Bertens (rbertens@nikhef.nl)
+// analyis task for phi-meson reconstruction and determination of V2
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliCentrality.h"
+#include "AliVEvent.h"
+#include "AliAnalysisTaskPhiFlow.h"
+#include "AliPID.h"
+#include "AliVEvent.h"
+#include "AliPIDResponse.h"
+#include "AliStack.h"
+#include "AliESDtrackCuts.h"
+#include "AliMCEvent.h"
+#include "TProfile.h"
+#include "AliFlowCandidateTrack.h"
+#include "AliFlowTrackCuts.h"
+#include "AliFlowEventSimple.h"
+#include "AliFlowCommonConstants.h"
+#include "AliFlowEvent.h"
+#include "TVector3.h"
+#include "TRandom2.h"
+#include "AliESDVZERO.h"
+
+class AliFlowTrackCuts;
+
+ClassImp(AliAnalysisTaskPhiFlow)
+
+AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
+   fEtaMinA(0), fEtaMaxA(0), fEtaMinB(0), fEtaMaxB(0), fCutsRP(NULL), fNullCuts(0), fESD(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), 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), 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), fKaonCuts(NULL), fVZEROA(0), fVZEROC(0), fTPCM(0), fPIDtype(kCombined), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fPairLoss(0), fEventPlanePtCut(0), fSetEventPlanePtCut(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);
+      }
+   }
+}
+//_____________________________________________________________________________
+AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow(const char *name) : AliAnalysisTaskSE(name),
+   fEtaMinA(0), fEtaMaxA(0), fEtaMinB(0), fEtaMaxB(0), fCutsRP(NULL), fNullCuts(0), fESD(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), 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), 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), fKaonCuts(NULL), fVZEROA(0), fVZEROC(0), fTPCM(0), fPIDtype(kCombined), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fPairLoss(0), fEventPlanePtCut(0), fSetEventPlanePtCut(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);
+      }
+   }
+   DefineInput(0, TChain::Class());
+   DefineOutput(1, TList::Class());
+   for (Int_t i = (0 + 2); i < (30 + 2); i++)
+   {
+      DefineOutput(i, AliFlowEventSimple::Class());
+   }
+}
+//_____________________________________________________________________________
+AliAnalysisTaskPhiFlow::~AliAnalysisTaskPhiFlow()
+{
+   // Destructor
+   for (Int_t i = 0; i < 30; i++)
+   {
+      delete fFlowEvent[i];
+   }
+   delete fNullCuts;
+}
+//_____________________________________________________________________________
+TH1F* AliAnalysisTaskPhiFlow::BookHistogram(const char* name)
+{
+   // Return a pointer to a TH1 with predefined binning
+   TH1F *hist = new TH1F(name, Form("M_{INV} (%s)", name), 60, .99, 1.092);
+   hist->GetXaxis()->SetTitle("M_{INV} (GeV / c^{2})");
+   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
+   TH2F *hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
+   hist->GetXaxis()->SetTitle("P (GeV / c)");
+   hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
+   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
+   TH1F* ratio = new TH1F(name, name, 100, 0, 7);
+   ratio->GetXaxis()->SetTitle("p_{T} ( GeV / c^{2} )");
+   ratio->GetYaxis()->SetTitle("No. of events");
+   ratio->Sumw2();
+   fOutputList->Add(ratio);
+   return ratio;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
+{
+   // Add Phi Identification Output Objects
+   fEventStats = new TH1F("fHistStats", "Event Statistics", 18, -.25, 4.25);
+   fEventStats->GetXaxis()->SetTitle("No. of events");
+   fEventStats->GetYaxis()->SetTitle("Statistics");
+   fOutputList->Add(fEventStats);
+   fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
+   fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
+   fOutputList->Add(fCentralityPass);
+   fOutputList->Add(fCentralityNoPass);
+
+   fEventPlaneSTAR = new TH1F("fEventPlaneSTAR", "Event plane orientation", 200, 0, 4);
+   fEventPlaneSTAR->GetXaxis()->SetTitle("#Psi_{2} (degrees)");
+   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");
+   fInvMPP03 = BookHistogram("PP, 0 < p_{T} < 0.3 GeV");
+   fInvMNN03 = BookHistogram("NN, 0 < p_{T} < 0.3 GeV");
+   fInvMNP36 = BookHistogram("NP, 0.3 < p_{T} < 0.6 GeV");
+   fInvMPP36 = BookHistogram("PP, 0.3 < p_{T} < 0.6 GeV");
+   fInvMNN36 = BookHistogram("NN, 0.3 < p_{T} < 0.6 GeV");
+   fInvMNP69 = BookHistogram("NP, 0.6 < p_{T} < 0.9 GeV");
+   fInvMPP69 = BookHistogram("PP, 0.6 < p_{T} < 0.9 GeV");
+   fInvMNN69 = BookHistogram("NN, 0.6 < p_{T} < 0.9 GeV");
+   fInvMNP912 = BookHistogram("NP, 0.9 < p_{T} < 1.2 GeV");
+   fInvMPP912 = BookHistogram("PP, 0.9 < p_{T} < 1.2 GeV");
+   fInvMNN912 = BookHistogram("NN, 0.9 < p_{T} < 1.2 GeV");
+   fInvMNP1215 = BookHistogram("NP, 1.2 < p_{T} < 1.5 GeV");
+   fInvMPP1215 = BookHistogram("PP, 1.2 < p_{T} < 1.5 GeV");
+   fInvMNN1215 = BookHistogram("NN, 1.2 < p_{T} < 1.5 GeV");
+   fInvMNP1518 = BookHistogram("NP, 1.5 < p_{T} < 1.8 GeV");
+   fInvMPP1518 = BookHistogram("PP, 1.5 < p_{T} < 1.8 GeV");
+   fInvMNN1518 = BookHistogram("NN, 1.5 < p_{T} < 1.8 GeV");
+   fInvMNP1821 = BookHistogram("NP, 1.8 < p_{T} < 2.1 GeV");
+   fInvMPP1821 = BookHistogram("PP, 1.8 < p_{T} < 2.1 GeV");
+   fInvMNN1821 = BookHistogram("NN, 1.8 < p_{T} < 2.1 GeV");
+   fInvMNP2124 = BookHistogram("NP, 2.1 < p_{T} < 2.4 GeV");
+   fInvMPP2124 = BookHistogram("PP, 2.1 < p_{T} < 2.4 GeV");
+   fInvMNN2124 = BookHistogram("NN, 2.1 < p_{T} < 2.4 GeV");
+   fInvMNP2427 = BookHistogram("NP, 2.4 < p_{T} < 2.7 GeV");
+   fInvMPP2427 = BookHistogram("PP, 2.4 < p_{T} < 2.7 GeV");
+   fInvMNN2427 = BookHistogram("NN, 2.4 < p_{T} < 2.7 GeV");
+   fInvMNP2730 = BookHistogram("NP, 2.7 < p_{T} < 3.0 GeV");
+   fInvMPP2730 = BookHistogram("PP, 2.7 < p_{T} < 3.0 GeV");
+   fInvMNN2730 = BookHistogram("NN, 2.7 < p_{T} < 3.0 GeV");
+   fInvMNP3035 = BookHistogram("NP, 3.0 < p_{T} < 3.5 GeV");
+   fInvMPP3035 = BookHistogram("PP, 3.0 < p_{T} < 3.5 GeV");
+   fInvMNN3035 = BookHistogram("NN, 3.0 < p_{T} < 3.5 GeV");
+   fInvMNP3540 = BookHistogram("NP, 3.5 < p_{T} < 4.0 GeV");
+   fInvMPP3540 = BookHistogram("PP, 3.5 < p_{T} < 4.0 GeV");
+   fInvMNN3540 = BookHistogram("NN, 3.5 < p_{T} < 4.0 GeV");
+   fInvMNP4045 = BookHistogram("NP, 4.0 < p_{T} < 4.5 GeV");
+   fInvMPP4045 = BookHistogram("PP, 4.0 < p_{T} < 4.5 GeV");
+   fInvMNN4045 = BookHistogram("NN, 4.0 < p_{T} < 4.5 GeV");
+   fInvMNP4550 = BookHistogram("NP, 4.5 < p_{T} < 5.0 GeV");
+   fInvMPP4550 = BookHistogram("PP, 4.5 < p_{T} < 5.0 GeV");
+   fInvMNN4550 = BookHistogram("NN, 4.5 < p_{T} < 5.0 GeV");
+   fInvMNP5055 = BookHistogram("NP, 5.0 < p_{T} < 5.5 GeV");
+   fInvMPP5055 = BookHistogram("PP, 5.0 < p_{T} < 5.5 GeV");
+   fInvMNN5055 = BookHistogram("NN, 5.0 < p_{T} < 5.5 GeV");
+   fInvMNP5560 = BookHistogram("NP, 5.5 < p_{T} < 6.0 GeV");
+   fInvMPP5560 = BookHistogram("PP, 5.5 < p_{T} < 6.0 GeV");
+   fInvMNN5560 = BookHistogram("NN, 5.5 < p_{T} < 6.0 GeV");
+   fInvMNP6065 = BookHistogram("NP, 6.0 < p_{T} < 6.5 GeV");
+   fInvMPP6065 = BookHistogram("PP, 6.0 < p_{T} < 6.5 GeV");
+   fInvMNN6065 = BookHistogram("NN, 6.0 < p_{T} < 6.5 GeV");
+   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);
+
+   fOutputList->Add(fEventPlaneSTAR);
+   fPtP = BookPtHistogram("i^{+}");
+   fPtN = BookPtHistogram("i^{-}");
+   fPtKP = BookPtHistogram("K^{+}");
+   fPtKN = BookPtHistogram("K^{-}");
+
+   fEventPlane = new TH1F("fEventPlane", "Event plane", 100, -2., 2.);
+   fEventPlane->GetXaxis()->SetTitle("#Psi_{2}");
+   fEventPlane->GetYaxis()->SetTitle("Events");
+   fOutputList->Add(fEventPlane);
+
+   fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
+   fOutputList->Add(fPhi);
+
+   fPt = new TH1F("fPt", "p_{T}", 100, 0, 5.5);
+   fOutputList->Add(fPt);
+
+   fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
+   fOutputList->Add(fEta);
+
+   fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
+   fOutputList->Add(fVZEROA);
+
+   fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
+   fOutputList->Add(fVZEROC);
+
+   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);
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
+{
+   // Create user defined output objects
+
+   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
+   cc->SetNbinsMult(10000);
+   cc->SetMultMin(0);
+   cc->SetMultMax(10000);
+
+   cc->SetNbinsPt(100);
+   cc->SetPtMin(0);
+   cc->SetPtMax(10);
+
+   cc->SetNbinsPhi(180);
+   cc->SetPhiMin(0.0);
+   cc->SetPhiMax(TMath::TwoPi());
+
+   cc->SetNbinsEta(200);
+   cc->SetEtaMin(-5.0);
+   cc->SetEtaMax(+5.0);
+
+   cc->SetNbinsQ(500);
+   cc->SetQMin(0.0);
+   cc->SetQMax(3.0);
+
+   // Create all output objects and store them to a list
+   fOutputList = new TList();
+   fOutputList->SetOwner();
+
+   // 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]); // TODO put flow events in a TList
+   }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::PairLoss(const AliESDtrack* track1, const AliESDtrack* 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()))));
+   }
+}
+//_____________________________________________________________________________
+Double_t AliAnalysisTaskPhiFlow::InvariantMass(const AliESDtrack* track1, const AliESDtrack* 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);
+   Double_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
+   Double_t e1 = TMath::Sqrt(track1->P() * track1->P() + masss);
+   Double_t e2 = TMath::Sqrt(track2->P() * track2->P() + masss);
+   Double_t es = TMath::Power((e1 + e2), 2);
+   if ((es - (pxs + pys + pzs)) < 0)
+   {
+      return 0.;
+   }
+   return TMath::Sqrt((es - (pxs + pys + pzs)));
+}
+//_____________________________________________________________________________
+Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const AliESDtrack* track1, const AliESDtrack* 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
+   if (track1->P()*track2->P() == 0) return 999;
+   return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const AliESDtrack* track1, const AliESDtrack* track2) const
+{
+   // 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;
+   return kTRUE;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
+{
+   // Set a centrality range ]min, max] and define the method to use for centrality selection
+   fCentralityMin = CentralityMin;
+   fCentralityMax = CentralityMax;
+   fkCentralityMethod = CentralityMethod;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPhiFlow::EventCut(AliESDEvent* esd) const
+{
+   // Impose event cuts
+   if (!esd) return kFALSE;
+   if (!CheckVertex(esd)) return kFALSE;
+   if (!CheckCentrality(esd)) return kFALSE;
+   PlotVZeroMultiplcities(esd);
+   return kTRUE;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::PlotVZeroMultiplcities(const AliESDEvent* esd) const
+{
+   // QA multiplicity plots
+   AliESDVZERO* vzero = esd->GetVZEROData();
+   fVZEROA->Fill(vzero->GetMTotV0A());
+   fVZEROC->Fill(vzero->GetMTotV0C());
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPhiFlow::CheckVertex(const AliESDEvent* ESD) const
+{
+   // Check if event vertex is within given range
+   const AliESDVertex * vertex = ESD->GetPrimaryVertex();
+   if (!vertex) return 0x0;
+   if (TMath::Abs(vertex->GetZ()) > fVertexRange) return 0x0;
+   return kTRUE;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPhiFlow::CheckCentrality(AliESDEvent* esd) const
+{
+   // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
+   if (!fkCentralityMethod)
+   {
+      AliFatal("No centrality method set! FATAL ERROR!");
+   }
+   AliCentrality* centr = esd->GetCentrality();
+   Double_t centrality = centr->GetCentralityPercentile(fkCentralityMethod);
+   Bool_t pass = kTRUE;
+   if ((centrality <= fCentralityMin) || (centrality > fCentralityMax))
+   {
+      pass = kFALSE;
+      fCentralityNoPass->Fill(centrality) ;
+   }
+   if (pass)
+   {
+      fCentralityPass->Fill(centrality);
+   }
+   return pass;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPhiFlow::IsKaon(const AliESDtrack* track) const
+{
+   // Check if particle is a kaon according to method set in steering macro
+   switch (fPIDtype)
+   {
+      case kTPC:
+         // Check if particle is flagged as Kaon in the ESD (TPC pid)
+         fNOPID->Fill(track->P(), track->GetTPCsignal());
+         if (track->GetPID() == AliPID::kKaon)
+         {
+            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 (fKaonCuts->PassesESDpidCut(track))
+         {
+            fPIDk->Fill(track->P(), track->GetTPCsignal());
+            return kTRUE;
+         }
+         break;
+      default:
+         AliFatal("No PID procedure set or available. Analysis will terminate!");
+         return kFALSE;
+         break;
+   }
+   return kFALSE;
+}
+//_____________________________________________________________________________
+Double_t AliAnalysisTaskPhiFlow::PhiPt(const AliESDtrack* track1, const AliESDtrack* track2) const
+{
+   // Calculate transverse momentum (p_T) of two tracks
+   TVector3 a(track1->Px(), track1->Py(), track1->Pz());
+   TVector3 b(track2->Px(), track2->Py(), track2->Pz());
+   TVector3 c = a + b;
+   return c.Pt();
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, const AliESDtrack* track1, const AliESDtrack* track2) const
+{
+   // Request transverse momentum (p_T), then fill invariant mass histograms as a function of p_T
+   Double_t pt = PhiPt(track1, track2);
+   if (tracktype == 0)
+   {
+      if ((0 <= pt) && (0.3 > pt))
+      {
+         fInvMNP03->Fill(InvariantMass(track1, track2));
+      }
+      if ((0.3 <= pt) && (0.6 > pt))
+      {
+         fInvMNP36->Fill(InvariantMass(track1, track2));
+      }
+      if ((0.6 <= pt) && (0.9 > pt))
+      {
+         fInvMNP69->Fill(InvariantMass(track1, track2));
+      }
+      if ((0.9 <= pt) && (1.2 > pt))
+      {
+         fInvMNP912->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.2 <= pt) && (1.5 > pt))
+      {
+         fInvMNP1215->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.5 <= pt) && (1.8 > pt))
+      {
+         fInvMNP1518->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.8 <= pt) && (2.1 > pt))
+      {
+         fInvMNP1821->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.1 <= pt) && (2.4 > pt))
+      {
+         fInvMNP2124->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.4 <= pt) && (2.7 > pt))
+      {
+         fInvMNP2427->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.7 <= pt) && (3.0 > pt))
+      {
+         fInvMNP2730->Fill(InvariantMass(track1, track2));
+      }
+      if ((3.0 <= pt) && (3.5 > pt))
+      {
+         fInvMNP3035->Fill(InvariantMass(track1, track2));
+      }
+      if ((3.5 <= pt) && (4.0 > pt))
+      {
+         fInvMNP3540->Fill(InvariantMass(track1, track2));
+      }
+      if ((4.0 <= pt) && (4.5 > pt))
+      {
+         fInvMNP4045->Fill(InvariantMass(track1, track2));
+      }
+      if ((4.5 <= pt) && (5.0 > pt))
+      {
+         fInvMNP4550->Fill(InvariantMass(track1, track2));
+      }
+      if ((5.0 <= pt) && (5.5 > pt))
+      {
+         fInvMNP5055->Fill(InvariantMass(track1, track2));
+      }
+      if ((5.5 <= pt) && (6.0 > pt))
+      {
+         fInvMNP5560->Fill(InvariantMass(track1, track2));
+      }
+      if ((6.0 <= pt) && (6.5 > pt))
+      {
+         fInvMNP6065->Fill(InvariantMass(track1, track2));
+      }
+      if ((6.5 <= pt) && (7.0 > pt))
+      {
+         fInvMNP6570->Fill(InvariantMass(track1, track2));
+      }
+   }
+   if (tracktype == 1)
+   {
+      if ((0 <= pt) && (0.3 > pt))
+      {
+         fInvMPP03->Fill(InvariantMass(track1, track2));
+      }
+      if ((0.3 <= pt) && (0.6 > pt))
+      {
+         fInvMPP36->Fill(InvariantMass(track1, track2));
+      }
+      if ((.6 <= pt) && (0.9 > pt))
+      {
+         fInvMPP69->Fill(InvariantMass(track1, track2));
+      }
+      if ((0.9 <= pt) && (1.2 > pt))
+      {
+         fInvMPP912->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.2 <= pt) && (1.5 > pt))
+      {
+         fInvMPP1215->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.5 <= pt) && (1.8 > pt))
+      {
+         fInvMPP1518->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.8 <= pt) && (2.1 > pt))
+      {
+         fInvMPP1821->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.1 <= pt) && (2.4 > pt))
+      {
+         fInvMPP2124->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.4 <= pt) && (2.7 > pt))
+      {
+         fInvMPP2427->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.7 <= pt) && (3.0 > pt))
+      {
+         fInvMPP2730->Fill(InvariantMass(track1, track2));
+      }
+      if ((3.0 <= pt) && (3.5 > pt))
+      {
+         fInvMPP3035->Fill(InvariantMass(track1, track2));
+      }
+      if ((3.5 <= pt) && (4.0 > pt))
+      {
+         fInvMPP3540->Fill(InvariantMass(track1, track2));
+      }
+      if ((4.0 <= pt) && (4.5 > pt))
+      {
+         fInvMPP4045->Fill(InvariantMass(track1, track2));
+      }
+      if ((4.5 <= pt) && (5.0 > pt))
+      {
+         fInvMPP4550->Fill(InvariantMass(track1, track2));
+      }
+      if ((5.0 <= pt) && (5.5 > pt))
+      {
+         fInvMPP5055->Fill(InvariantMass(track1, track2));
+      }
+      if ((5.5 <= pt) && (6.0 > pt))
+      {
+         fInvMPP5560->Fill(InvariantMass(track1, track2));
+      }
+      if ((6.0 <= pt) && (6.5 > pt))
+      {
+         fInvMPP6065->Fill(InvariantMass(track1, track2));
+      }
+      if ((6.5 <= pt) && (7.0 > pt))
+      {
+         fInvMPP6570->Fill(InvariantMass(track1, track2));
+      }
+   }
+   if (tracktype == 2)
+   {
+      if ((0 <= pt) && (0.3 > pt))
+      {
+         fInvMNN03->Fill(InvariantMass(track1, track2));
+      }
+      if ((0.3 <= pt) && (0.6 > pt))
+      {
+         fInvMNN36->Fill(InvariantMass(track1, track2));
+      }
+      if ((.6 <= pt) && (0.9 > pt))
+      {
+         fInvMNN69->Fill(InvariantMass(track1, track2));
+      }
+      if ((0.9 <= pt) && (1.2 > pt))
+      {
+         fInvMNN912->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.2 <= pt) && (1.5 > pt))
+      {
+         fInvMNN1215->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.5 <= pt) && (1.8 > pt))
+      {
+         fInvMNN1518->Fill(InvariantMass(track1, track2));
+      }
+      if ((1.8 <= pt) && (2.1 > pt))
+      {
+         fInvMNN1821->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.1 <= pt) && (2.4 > pt))
+      {
+         fInvMNN2124->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.4 <= pt) && (2.7 > pt))
+      {
+         fInvMNN2427->Fill(InvariantMass(track1, track2));
+      }
+      if ((2.7 <= pt) && (3.0 > pt))
+      {
+         fInvMNN2730->Fill(InvariantMass(track1, track2));
+      }
+      if ((3.0 <= pt) && (3.5 > pt))
+      {
+         fInvMNN3035->Fill(InvariantMass(track1, track2));
+      }
+      if ((3.5 <= pt) && (4.0 > pt))
+      {
+         fInvMNN3540->Fill(InvariantMass(track1, track2));
+      }
+      if ((4.0 <= pt) && (4.5 > pt))
+      {
+         fInvMNN4045->Fill(InvariantMass(track1, track2));
+      }
+      if ((4.5 <= pt) && (5.0 > pt))
+      {
+         fInvMNN4550->Fill(InvariantMass(track1, track2));
+      }
+      if ((5.0 <= pt) && (5.5 > pt))
+      {
+         fInvMNN5055->Fill(InvariantMass(track1, track2));
+      }
+      if ((5.5 <= pt) && (6.0 > pt))
+      {
+         fInvMNN5560->Fill(InvariantMass(track1, track2));
+      }
+      if ((6.0 <= pt) && (6.5 > pt))
+      {
+         fInvMNN6065->Fill(InvariantMass(track1, track2));
+      }
+      if ((6.5 <= pt) && (7.0 > pt))
+      {
+         fInvMNN6570->Fill(InvariantMass(track1, track2));
+      }
+   }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::EventPlane(const AliESDEvent* esd)
+{
+   // 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 = esd->GetNumberOfTracks();
+   Int_t tpcMultiplicity(0);
+   for (Int_t i = 0; i < nTracks; i++)
+   {
+      AliESDtrack* track = fESD->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));
+}
+//_____________________________________________________________________________
+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();
+      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;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::EllipticFlow(Double_t* const v2, const AliESDtrack* track1, const AliESDtrack* track2) const
+{
+   // Calculate elliptic flow coefficient v_2 - cos 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);
+   // 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);
+   }
+   if ((0.3 <= pt) && (0.6 > pt))
+   {
+      v2[1] += flow;
+      v2[19]++;
+      fProfV2InvM36->Fill(invm, flow);
+   }
+   if ((0.6 <= pt) && (0.9 > pt))
+   {
+      v2[2] += flow;
+      v2[20]++;
+      fProfV2InvM69->Fill(invm, flow);
+   }
+   if ((0.9 <= pt) && (1.2 > pt))
+   {
+      v2[3] += flow;
+      v2[21]++;
+      fProfV2InvM912->Fill(invm, flow);
+   }
+   if ((1.2 <= pt) && (1.5 > pt))
+   {
+      v2[4] += flow;
+      v2[22]++;
+      fProfV2InvM1215->Fill(invm, flow);
+   }
+   if ((1.5 <= pt) && (1.8 > pt))
+   {
+      v2[5] += flow;
+      v2[23]++;
+      fProfV2InvM1518->Fill(invm, flow);
+   }
+   if ((1.8 <= pt) && (2.1 > pt))
+   {
+      v2[6] += flow;
+      v2[24]++;
+      fProfV2InvM1821->Fill(invm, flow);
+   }
+   if ((2.1 <= pt) && (2.4 > pt))
+   {
+      v2[7] += flow;
+      v2[25]++;
+      fProfV2InvM2124->Fill(invm, flow);
+   }
+   if ((2.4 <= pt) && (2.7 > pt))
+   {
+      v2[8] += flow;
+      v2[26]++;
+      fProfV2InvM2427->Fill(invm, flow);
+   }
+   if ((2.7 <= pt) && (3.0 > pt))
+   {
+      v2[9] += flow;
+      v2[27]++;
+      fProfV2InvM2730->Fill(invm, flow);
+   }
+
+   if ((3.0 <= pt) && (3.5 > pt))
+   {
+      v2[10] += flow;
+      v2[28]++;
+      fProfV2InvM3035->Fill(invm, flow);
+   }
+   if ((3.5 <= pt) && (4.0 > pt))
+   {
+      v2[11] += flow;
+      v2[29]++;
+      fProfV2InvM3540->Fill(invm, flow);
+   }
+   if ((4.0 <= pt) && (4.5 > pt))
+   {
+      v2[12] += flow;
+      v2[30]++;
+      fProfV2InvM4045->Fill(invm, flow);
+   }
+   if ((4.5 <= pt) && (5.0 > pt))
+   {
+      v2[13] += flow;
+      v2[31]++;
+      fProfV2InvM4550->Fill(invm, flow);
+   }
+   if ((5.0 <= pt) && (5.5 > pt))
+   {
+      v2[14] += flow;
+      v2[32]++;
+      fProfV2InvM5055->Fill(invm, flow);
+   }
+   if ((5.5 <= pt) && (6.0 > pt))
+   {
+      v2[15] += flow;
+      v2[33]++;
+      fProfV2InvM5560->Fill(invm, flow);
+   }
+   if ((6.0 <= pt) && (6.5 > pt))
+   {
+      v2[16] += flow;
+      v2[34]++;
+      fProfV2InvM6065->Fill(invm, flow);
+   }
+   if ((6.5 <= pt) && (7.0 > pt))
+   {
+      v2[17] += flow;
+      v2[35]++;
+      fProfV2InvM6570->Fill(invm, flow);
+   }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::EllipticFlowSin(Double_t* const v2Sin, const AliESDtrack* track1, const AliESDtrack* 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);
+   }
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPhiFlow::EventPlaneTrack(AliESDtrack* track) const
+{
+   // Check if track is suitable for event plane estimation
+   if (!track) return kFALSE;
+   Bool_t sel = fCutsRP->IsSelected(track);
+   if (sel)
+   {
+      fPhi->Fill(track->Phi());
+      fPt->Fill(track->Pt());
+      fEta->Fill(track->Eta());
+      return sel;
+   }
+   return kFALSE;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskPhiFlow::PhiTrack(AliESDtrack* track) const
+{
+   // Check if track is suitable for phi flow analysis
+   if (!track) return kFALSE;
+   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]));
+   }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::SetNullCuts(AliESDEvent* esd)
+{
+   // Set null cuts
+   fCutsRP->SetEvent(esd, MCEvent());
+   fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
+   fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
+   fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
+   fNullCuts->SetEvent(esd, MCEvent());
+}
+//_____________________________________________________________________________
+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);
+   }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
+{
+   // UserExec: execute for each event. Commented where necessary
+   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+   // Check whether event passes event cuts
+   if (!EventCut(fESD))
+   {
+      return;
+   }
+   SetNullCuts(fESD);
+   fKaonCuts->SetEvent(fESD, MCEvent());
+   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++)
+   {
+      AliESDtrack* track = fESD->GetTrack(iTracks);
+      if (!PhiTrack(track))
+      {
+         continue;
+      }
+      Bool_t charge = kFALSE;
+      if (track->Charge() > 0)
+      {
+         charge = kTRUE;
+         fEventStats->Fill(1);
+         fPtP->Fill(track->Pt());
+      }
+      if (track->Charge() < 0)
+      {
+         fEventStats->Fill(2);
+         fPtN->Fill(track->Pt());
+      }
+      if (IsKaon(track))
+      {
+         if (charge)
+         {
+            up[unp] = track;
+            unp++;
+            fEventStats->Fill(3);
+            fPtKP->Fill(track->Pt());
+         }
+         if (!charge)
+         {
+            un[unn] = track;
+            unn++;
+            fEventStats->Fill(4);
+            fPtKN->Fill(track->Pt());
+         }
+      }
+   }
+   // 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 nTracks = 0; nTracks < unn ; nTracks++)
+      {
+         if(fApplyDeltaDipCut&&(!CheckDeltaDipAngle(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());
+         TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
+         TVector3 c = a + b;
+         Double_t phi = c.Phi();
+         Double_t eta = c.Eta();
+         int 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]))
+            {
+               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);
+            }
+      }
+   }
+   for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
+   {
+      for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++)
+      {
+         if(fApplyDeltaDipCut&&(!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
+         PtSelector(1, up[pTracks], up[nTracks]);
+      }
+   }
+   for (Int_t nTracks = 0; nTracks < unn ; nTracks++)
+   {
+      for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++)
+      {
+         if(fApplyDeltaDipCut&&(!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
+         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]);
+}
+//_____________________________________________________________________________
+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
+{
+   // Consruct Flow Candidate Track from two selected candidates
+   AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
+   sTrack->SetMass(mass);
+   sTrack->SetPt(pt);
+   sTrack->SetPhi(phi);
+   sTrack->SetEta(eta);
+   for (Int_t iDau = 0; iDau != nDau; ++iDau)
+      sTrack->AddDaughter(iID[iDau]);
+   sTrack->SetForPOISelection(kTRUE);
+   sTrack->SetForRPSelection(kFALSE);
+   return sTrack;
+}
+//_____________________________________________________________________________
+
+
diff --git a/PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.h b/PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.h
new file mode 100644 (file)
index 0000000..41141e3
--- /dev/null
@@ -0,0 +1,262 @@
+#ifndef ALIANALYSISTASKPHIFLOW_H
+#define ALIANALYSISTASKPHIFLOW_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. */
+/* See cxx source for full Copyright notice */
+/* $Id$ */
+
+// AliAnalysisTaskPhiFlow:
+// origin: Redmer Alexander Bertens (rbertens@nikhef.nl)
+// analyis task for phi-meson reconstruction and determination of V2
+
+class TH1F;
+class TH1F;
+class TH2F;
+class TProfile;
+class AliESDEvent;
+class AliESDtrackCuts;
+class AliFlowTrackCuts;
+class AliFlowEvent;
+class AliFlowCandidateTrack;
+
+#include "AliAnalysisTaskSE.h"
+
+enum ESDEventStats_t
+{
+   kNUnlikePairs = 0,
+   kNLikeNPairs,
+   kNLikePPairs,
+   kNUnlikeKPairs,
+   kNLikeNKPairs,
+   kNLikePKPairs,
+   kNStats = kNLikePKPairs,
+};
+
+class AliAnalysisTaskPhiFlow : public AliAnalysisTaskSE
+{
+public:
+   AliAnalysisTaskPhiFlow();
+   AliAnalysisTaskPhiFlow(const char *name);
+   virtual ~AliAnalysisTaskPhiFlow();
+
+   enum PIDtype
+   {
+       kTPC = 0,
+       kCombined,
+   };
+
+
+   TH1F*          BookHistogram(const char * name);
+   TH2F*          BookPIDHistogram(const char * name);
+   TProfile*      BookV2Profile(const char * name, Bool_t pt, Bool_t cos);
+   TH1F*          BookPtHistogram(const char* name);
+   void           AddPhiIdentificationOutputObjects();
+   virtual void   UserCreateOutputObjects();
+   void           PairLoss(const AliESDtrack* track1, const AliESDtrack* track2) const;
+   Double_t       InvariantMass(const AliESDtrack* track_1, const AliESDtrack* track_2) const;
+   Double_t       DeltaDipAngle(const AliESDtrack* track1, const AliESDtrack* track2) const;
+   Bool_t         CheckDeltaDipAngle(const AliESDtrack* track1, const AliESDtrack* track2) const;
+   void           SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod);
+   void           SetVertexZ(Float_t z)  { fVertexRange = z; };
+   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; } ;
+   Bool_t         EventCut( AliESDEvent* esd) const;
+   void           PlotVZeroMultiplcities(const AliESDEvent* esd) const;
+   Bool_t         CheckVertex(const AliESDEvent* esd) const ;
+   Bool_t         CheckCentrality(AliESDEvent* esd) const;
+   Bool_t         IsKaon(const AliESDtrack* track) const;
+   Double_t       PhiPt(const AliESDtrack* track_1, const AliESDtrack* track_2) const;
+   void           PtSelector(Int_t _track_type, const AliESDtrack* track_1, const AliESDtrack* track_2) const;
+   void           EventPlane(const AliESDEvent* esd);
+   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;
+   void           EllipticFlow(Double_t* const v2, const AliESDtrack* track_1, const AliESDtrack* track_2) const;
+   Bool_t         EventPlaneTrack(AliESDtrack* track) const;
+   Bool_t         PhiTrack(AliESDtrack* track) const;
+   void           FlowFinish(Double_t* const v2) const;
+   void           EllipticFlowSin(Double_t* const v2Sin, const AliESDtrack* track_1, const AliESDtrack* track_2) const;
+   void           FlowFinishSin(Double_t* const v2Sin) const;
+   void           SetNullCuts(AliESDEvent* 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           SetKaonCuts(AliFlowTrackCuts *cutsKaon) { fKaonCuts = cutsKaon; }
+   void           SetIdentificationType(PIDtype type) {fPIDtype = type;}
+   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; }
+
+private:
+
+   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)
+   AliFlowTrackCuts     *fNullCuts; // dummy cuts for flow event tracks
+   AliFlowEvent         *fFlowEvent[30]; //! flow events (one for each inv mass band)
+   AliESDEvent          *fESD;    //! ESD object
+   TList                *fOutputList; // ! Output list
+   TH1F                 *fEventStats; // ! Histogram for event statistics
+   TH1F                 *fCentralityPass; // ! QA histogram of events that pass centrality cut
+   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
+   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
+   TH1F                 *fInvMNP36; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP36; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN36; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP69; //!  Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP69; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN69; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP912; //!  Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP912; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN912; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP1215; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP1215; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN1215; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP1518; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP1518; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN1518; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP1821; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP1821; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN1821; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP2124; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP2124; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN2124; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP2427; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP2427; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN2427; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP2730; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP2730; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN2730; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP3035; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP3035; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN3035; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP3540; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP3540; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN3540; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP4045; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP4045; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN4045; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP4550; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP4550; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN4550; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP5055; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP5055; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN5055; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP5560; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP5560; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN5560; //! Invariant mass of like sign (--) kaon pairs
+   TH1F                 *fInvMNP6065; //! Invariant mass of unlike sign kaon pairs
+   TH1F                 *fInvMPP6065; //! Invariant mass of like sign (++) kaon pairs
+   TH1F                 *fInvMNN6065; //! Invariant mass of like sign (--) kaon pairs
+   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
+   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                 *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                 *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
+   TH1F                 *fPtKN; //! QA histogram of p_t distribution of negative kaons
+   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)
+   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
+   AliFlowTrackCuts     *fKaonCuts; // dummy flowtrackcuts object for bayesian pid method
+   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&);
+   AliAnalysisTaskPhiFlow& operator=(const AliAnalysisTaskPhiFlow&);
+
+   AliFlowCandidateTrack*  MakeTrack(Double_t, Double_t, Double_t, Double_t, Int_t , Int_t[]) const;
+
+   ClassDef(AliAnalysisTaskPhiFlow, 1);
+
+};
+
+#endif
+
+
index 954158b..f816ff3 100644 (file)
@@ -27,5 +27,7 @@
 #pragma link C++ class AliAnalysisTaskPIDflowQA+;
 #pragma link C++ class AliAnalysisTaskQAPmdflow+;
 #pragma link C++ class AliFlowBayesianPID+;
+#pragma link C++ class AliAnalysisTaskPhiFlow+;
 
 #endif
+
diff --git a/PWGCF/FLOW/macros/AddTaskPhiFlow.C b/PWGCF/FLOW/macros/AddTaskPhiFlow.C
new file mode 100644 (file)
index 0000000..d0fb559
--- /dev/null
@@ -0,0 +1,168 @@
+/////////////////////////////////////////////////////////////////////////////////////////////
+//
+// AddTaskPhiFlow macro
+// Author: Redmer A. Bertens, Utrecht University, 2012
+/////////////////////////////////////////////////////////////////////////////////////////////
+AliAnalysisTaskPhiFlow* AddTaskPhiFlow( Float_t centrMin = 30.,
+                                        Float_t centrMax = 40.,
+                                        Float_t vertexZ = 10.,
+                                        Float_t bayesThresholdP = 0.5,
+                                        Float_t minEtaA = -0.8,
+                                        Float_t maxEtaA = 0.,
+                                        Float_t minEtaB = 0.,
+                                        Float_t maxEtaB = 0.8,
+                                        Float_t deltaDip = 0.0,
+                                        Float_t deltaDipMaxPt = 0.0,
+                                        Float_t eventPlanePtCut = 10,
+                                        AliAnalysisTaskPhiFlow::PIDtype pidtype = AliAnalysisTaskPhiFlow::kCombined)
+{
+   TString centralityName("");
+   centralityName += Form("%.0f", centrMin);
+   centralityName += "-";
+   centralityName += Form("%.0f", centrMax);
+   centralityName += "_";
+   centralityName += Form("vZ%.1f", vertexZ);
+   centralityName += "_";
+   centralityName += Form("bayP%.1f", bayesThresholdP);
+   centralityName += "_";
+   centralityName += "_EtaA";
+   centralityName += Form("%.1f", minEtaA);
+   centralityName += "-";
+   centralityName += Form("%.1f", maxEtaA);
+   centralityName += "_EtaB";
+   centralityName += Form("%.1f", minEtaB);
+   centralityName += "-";
+   centralityName += Form("%.1f", maxEtaB);
+   centralityName += "-";
+   centralityName += Form("dDip%.f", deltaDip);
+   centralityName += "-";
+   centralityName += Form("dDipPt%.f", deltaDipMaxPt);
+   centralityName += "-";
+   centralityName += Form("epPtCut%.f", eventPlanePtCut);
+
+   TString fileName = AliAnalysisManager::GetCommonFileName();
+   fileName += ":PhiV2FlowEvents";
+
+   // get the manager
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   if (!mgr)
+   {
+      ::Error("AddTaskPhiV2", "No analysis manager to connect to.");
+      return NULL;
+   }
+   if (!mgr->GetInputEventHandler())
+   {
+      ::Error("AddTaskPhiV2", "This task requires an input event handler");
+      return NULL;
+   }
+
+   // create the main task
+   AliAnalysisTaskPhiFlow *task = new AliAnalysisTaskPhiFlow("TaskPhiFlow");
+   task->SelectCollisionCandidates(AliVEvent::kMB);
+
+   // set RP cuts for EP method
+   AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts("GlobalRP");
+   cutsRP->GetStandardTPCStandaloneTrackCuts();
+
+   // set POI cuts for EP method
+   AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts("GlobalPOI");
+   cutsPOI->GetStandardGlobalTrackCuts2010();
+
+   // set and configure kaon track cuts (for the bayesian PID method)
+   // select this method by setting PIDtype to kCombined
+   AliFlowTrackCuts* cutsKaon = new AliFlowTrackCuts("kaon_cuts");
+   AliFlowTrackCuts::PIDsource pid = AliFlowTrackCuts::kTPCbayesian;
+   cutsKaon->SetPID(AliPID::kKaon, pid, bayesThresholdP);
+   cutsKaon->SetAllowTOFmismatchFlag(kTRUE);
+   cutsKaon->SetRequireStrictTOFTPCagreement(kFALSE);
+
+   // set some more cuts ...
+   task->EventPlanePtCut(kTRUE); // set to true to enforce pt cut on event plane tracks
+   task->SetEventPlanePtCut(eventPlanePtCut); // set the upper value for Pt
+   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);
+   task->SetPOICuts(cutsPOI);
+
+   // set the kaon cuts, and specify the PID procedure which will be used
+   task->SetKaonCuts(cutsKaon);
+   task->SetIdentificationType(pidtype);
+
+   //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.);
+
+   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;
+   }
+   task->SetMassRanges(flowBands);
+   task->SetEtaRanges(minEtaA, maxEtaA, minEtaB, maxEtaB);
+   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]);
+   }
+
+   // Scalar Product
+   AliAnalysisTaskScalarProduct *taskSP[30];
+   AliAnalysisDataContainer *coutputSP[30];
+   for (int r = 0; r != 30; ++r)
+   {
+      taskSP[r] = new AliAnalysisTaskScalarProduct(Form("SP_MassBand_%d", r), kFALSE);
+      taskSP[r]->SetRelDiffMsub(1.0);
+      taskSP[r]->SetApplyCorrectionForNUA(kTRUE);
+      mgr->AddTask(taskSP[r]);
+      coutputSP[r] = mgr->CreateContainer(Form("cobjSP_MassBand_%d", r), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
+      mgr->ConnectInput(taskSP[r], 0, coutputCandidates[r]);
+      mgr->ConnectOutput(taskSP[r], 1, coutputSP[r]);
+   }
+   // Q-Cumulants
+   AliAnalysisTaskQCumulants *taskQC[30];
+   AliAnalysisDataContainer *coutputQC[30];
+   for (int r = 0; r != 30; ++r)
+   {
+      taskQC[r] = new AliAnalysisTaskQCumulants(Form("QC_MassBand_%d", r), kFALSE);
+      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", r), TList::Class(), AliAnalysisManager::kOutputContainer, fileName);
+      mgr->ConnectInput(taskQC[r], 0, coutputCandidates[r]);
+      mgr->ConnectOutput(taskQC[r], 1, coutputQC[r]);
+   }
+
+   mgr->ConnectInput(task, 0, cinput);
+   mgr->ConnectOutput(task, 1, coutput);
+
+   // return the task
+   return task;
+};
+