]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
test task to replace phi and K* analysis
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Oct 2012 13:59:20 +0000 (13:59 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Oct 2012 13:59:20 +0000 (13:59 +0000)
PWG/CMakelibPWGflowTasks.pkg
PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.cxx [new file with mode: 0644]
PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h [new file with mode: 0644]
PWG/PWGflowTasksLinkDef.h
PWGCF/FLOW/macros/AddTwoParticleResonanceFlowTask.C [new file with mode: 0644]

index ecb8a7636fb926ca399efa24b78c74783611d704..d3601ddc2406dced490a72aad72d1743cdb659b6 100644 (file)
@@ -50,6 +50,7 @@ set ( SRCS
     FLOW/Tasks/AliAnalysisTaskQAPmdflow.cxx 
     FLOW/Tasks/AliFlowBayesianPID.cxx
     FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx
+    FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.cxx
     FLOW/Tasks/AliAnalysisTaskFilterFE.cxx
     FLOW/Tasks/AliAnalysisTaskVnV0.cxx
     FLOW/Tasks/AliFlowVZEROResults.cxx
diff --git a/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.cxx b/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.cxx
new file mode 100644 (file)
index 0000000..6494a5f
--- /dev/null
@@ -0,0 +1,972 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+// AliAnalysisTwoParticleResonanceFlowTask:
+// author: Redmer Alexander Bertens (rbertens@cern.ch)
+//         You Zhou                 (you.zhou@cern.ch) 
+// analyis task for 
+// AliResonanceFlowHelperTrack provides a lightweight helper track for reconstruction
+// This analysis DOES NOT support ESDs
+// Version: v2012.9.7
+
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TObjArray.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliCentrality.h"
+#include "AliVEvent.h"
+#include "AliAnalysisTwoParticleResonanceFlowTask.h"
+#include "AliFlowBayesianPID.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "TProfile.h"
+#include "AliFlowCandidateTrack.h"
+#include "AliFlowTrackCuts.h"
+#include "AliFlowEventSimple.h"
+#include "AliFlowTrackSimple.h"
+#include "AliFlowCommonConstants.h"
+#include "AliFlowEvent.h"
+#include "TVector3.h"
+#include "AliAODVZERO.h"
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliAODMCParticle.h"
+#include "AliAnalysisTaskVnV0.h"
+#include "AliEventPoolManager.h"
+
+class AliFlowTrackCuts;
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAnalysisTwoParticleResonanceFlowTask)
+
+AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask() : AliAnalysisTaskSE(),
+  fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fDebug(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0),  fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fNOPIDTOF(0), fPIDTOF(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), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fAnalysisSummary(0)
+{
+   // Default constructor
+   for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
+   for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
+   for(Int_t i(0); i < 20; i++) {
+       fVertexMixingBins[i] = 0;
+       fCentralityMixingBins[i] = 0;
+   }
+   fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
+   for(Int_t i(0); i < 18; i++) {
+       for(Int_t j(0); j < 2; j++) {
+           fV0Data[i][j] = 0;
+           for(Int_t k(0); k < 15; k++) fPhiMinusPsiDataContainer[k][i][j] = 0x0;
+           for(Int_t k(0); k < 15; k++) fPhiMinusPsiBackgroundContainer[k][i][j] = 0x0;
+       }
+       fResonanceSignal[i] = 0; fResonanceBackground[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.; fdPhiBins[i] = 0.;
+   }
+}
+//_____________________________________________________________________________
+AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask(const char *name) : AliAnalysisTaskSE(name), fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fDebug(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0),  fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fNOPIDTOF(0), fPIDTOF(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), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fAnalysisSummary(0)
+{
+   // Constructor
+  for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
+  for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
+  for(Int_t i(0); i < 20; i++) {
+      fVertexMixingBins[i] = 0;
+      fCentralityMixingBins[i] = 0;
+  }
+  fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
+  for(Int_t i(0); i < 18; i++) {
+      for(Int_t j(0); j < 2; j++) {
+          fV0Data[i][j] = 0;
+          for(Int_t k(0); k < 15; k++) fPhiMinusPsiDataContainer[k][i][j] = 0x0;
+          for(Int_t k(0); k < 15; k++) fPhiMinusPsiBackgroundContainer[k][i][j] = 0x0;
+      }
+      fResonanceSignal[i] = 0; fResonanceBackground[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.; fdPhiBins[i] = 0.;
+  }
+  DefineInput(0, TChain::Class());
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, AliFlowEventSimple::Class());
+  if(fDebug > 0) cout << " === Created instance of AliAnaysisTaskGenericResonanceFlowAnalysis === " << endl;
+}
+//_____________________________________________________________________________
+AliAnalysisTwoParticleResonanceFlowTask::~AliAnalysisTwoParticleResonanceFlowTask()
+{
+   // Destructor
+   if (fNullCuts) delete fNullCuts;
+   if (fOutputList) delete fOutputList;
+   if (fCandidates) delete fCandidates;
+   if (fFlowEvent) delete fFlowEvent;
+   if (fBayesianResponse) delete fBayesianResponse;
+   if (fEventMixing) delete fPoolManager;
+   if (fDebug > 0) cout << " === Deleted instance of AliAnalysisTwoParticleResonanceFlowTask === " << endl;
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::ForceExit(Int_t type, const char* message)
+{
+    // force exit
+    if(type==0) cout << " --> Illegal options in AddTask <-- " << endl;
+    if(type==1) cout << " --> Unknown error <-- " << endl;
+    AliFatal(message);
+}
+//_____________________________________________________________________________
+TH1F* AliAnalysisTwoParticleResonanceFlowTask::BookHistogram(const char* name)
+{
+   // Return a pointer to a TH1 with predefined binning
+   if(fDebug > 0) cout << " *** BookHistogram() *** " << name << endl;
+   TH1F *hist = new TH1F(name, Form("M_{INV} (%s)", name), 2*fMassBins, fMinMass, fMaxMass);
+   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* AliAnalysisTwoParticleResonanceFlowTask::BookPIDHistogram(const char* name, Bool_t TPC)
+{
+   // Return a pointer to a TH2 with predefined binning
+   if(fDebug > 0) cout << " *** BookPIDHisotgram() *** " << endl;
+   TH2F *hist = 0x0;
+   if(TPC) {
+       hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
+       hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
+   }
+   if(!TPC) {
+       hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1.5);
+       hist->GetYaxis()->SetTitle("#beta");
+   }
+   hist->GetXaxis()->SetTitle("P (GeV / c)");
+   fOutputList->Add(hist);
+   return hist;
+}
+//_____________________________________________________________________________
+TH1F* AliAnalysisTwoParticleResonanceFlowTask::InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
+{
+   // intialize p_t histograms for each p_t bin
+   if(fDebug > 0) cout << " *** InitPtSpectraHistograms() *** " << endl;
+   TH1F* hist = new TH1F(Form("%4.2f p_{t} %4.2f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 2*fMassBins, nmin, nmax);
+   hist->GetXaxis()->SetTitle("p_{T} GeV / c");
+   fOutputList->Add(hist);
+   return hist;
+}
+//_____________________________________________________________________________
+TH1F* AliAnalysisTwoParticleResonanceFlowTask::BookPtHistogram(const char* name)
+{
+   // Return a pointer to a p_T spectrum histogram
+   if(fDebug > 0) cout << " *** BookPtHistogram() *** " << endl;
+   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;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTwoParticleResonanceFlowTask::InitializeAnalysis()
+{
+    // set up the analysis
+   fNullCuts = new AliFlowTrackCuts("null_cuts");
+   fBayesianResponse = new AliFlowBayesianPID();
+   Float_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 !!!");
+   if(fNPtBins > 18) AliFatal("Invalid number of pt bins initialied ( > 18 ) -- terminating analysis !!!");
+   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
+   cc->SetNbinsQ(500);           cc->SetNbinsPhi(180);           cc->SetNbinsMult(10000);
+   cc->SetQMin(0.0);             cc->SetPhiMin(0.0);             cc->SetMultMin(0);
+   cc->SetQMax(3.0);             cc->SetPhiMax(TMath::TwoPi());  cc->SetMultMax(10000);
+   cc->SetNbinsMass(fMassBins);  cc->SetNbinsEta(200);           (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
+   cc->SetMassMin(fMinMass);     cc->SetEtaMin(-5.0);            cc->SetPtMin(0);
+   cc->SetMassMax(fMaxMass);     cc->SetEtaMax(+5.0);            (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
+   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
+   if (man) {
+      AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
+      if (inputHandler)   fPIDResponse = inputHandler->GetPIDResponse();
+   }
+   fMassA = fMassA*fMassA;  //we always use mass squared in calculations, so doing it once in the init reduces cpu load
+   fMassB = fMassB*fMassB;
+   return kTRUE;
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::UserCreateOutputObjects()
+{
+   // Create user defined output objects
+   if(!InitializeAnalysis()) ForceExit(0, " > Couldn't initialize the analysis !!! <");
+   if(fDebug > 0) cout << " *** UserCreateOutputObjects() *** " << endl;
+     // Create all output objects and store them to a list
+   fOutputList = new TList();
+   fOutputList->SetOwner(kTRUE);
+   if(fQA) {
+       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);
+   fOutputList->Add(fCentralityPass);
+   if(fQA) {
+       fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
+       fOutputList->Add(fCentralityNoPass);
+       fNOPID = BookPIDHistogram("TPC signal, all particles", kTRUE);
+       fPIDk = BookPIDHistogram("TPC signal, kaons", kTRUE);
+       fPIDp = BookPIDHistogram("TPC signal, pions", kTRUE);
+       fNOPIDTOF = BookPIDHistogram("TOF signal, all particles", kFALSE);
+       fPIDTOF = BookPIDHistogram("TOF signal, kaons", kFALSE);
+   }
+   for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
+       if(!fPhiMinusPsiMethod) {
+           fResonanceSignal[ptbin] = BookHistogram(Form("NP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
+           fResonanceBackground[ptbin] = BookHistogram(Form("PP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
+       }
+       else if(fPhiMinusPsiBackgroundContainer) {
+           for(Int_t i(0); i<fNdPhiBins; i++) {// for all delta phi bins
+                   fPhiMinusPsiDataContainer[i][ptbin][0] = BookHistogram(Form("%4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{A}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
+                   fPhiMinusPsiDataContainer[i][ptbin][1] = BookHistogram(Form("%4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{C}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
+                   fPhiMinusPsiBackgroundContainer[i][ptbin][0] = BookHistogram(Form("BG, %4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{A}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
+                   fPhiMinusPsiBackgroundContainer[i][ptbin][1] = BookHistogram(Form("BG, %4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{C}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
+           }
+       }
+   }
+   if(fQA) {
+       for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) fPtSpectra[ptbin] = InitPtSpectraHistograms(fPtBins[ptbin], fPtBins[ptbin+1]);
+       fPtP = BookPtHistogram("i^{+}");
+       fPtN = BookPtHistogram("i^{-}");
+       fPtKP = BookPtHistogram("K^{+}");
+       fPtKN = BookPtHistogram("K^{-}");
+       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);
+       fDCAXYQA = new TH1F("fDCAXYQA", "fDCAXYQA", 1000, -5, 5);
+       fOutputList->Add(fDCAXYQA);
+       fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
+       fOutputList->Add(fDCAZQA);
+   }
+   if(fIsMC || fQA) {
+       fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
+       fOutputList->Add(fDCAAll);
+       fDCAPrim = new TH2F("fDCAprim","fDCAprim", 1000, 0, 10, 1000, -5, 5);
+       fOutputList->Add(fDCAPrim);
+       fDCASecondaryWeak = new TH2F("fDCASecondaryWeak","fDCASecondaryWeak", 1000, 0, 10, 1000, -5, 5);
+       fOutputList->Add(fDCASecondaryWeak);
+       fDCAMaterial = new TH2F("fDCAMaterial","fDCAMaterial", 1000, 0, 10, 1000, -5, 5);
+       fOutputList->Add(fDCAMaterial);
+   }
+   if(fV0||fPhiMinusPsiMethod) {
+       fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
+       fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
+       fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
+       fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
+       fOutputList->Add(fSubEventDPhiv2);
+       const char* V0[] = {"V0A", "V0C"};
+       for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++)
+           for(Int_t iV0(0); iV0 < 2; iV0++) {
+                   fV0Data[ptbin][iV0] = new TProfile(Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), fMassBins, fMinMass, fMaxMass);
+                   fOutputList->Add(fV0Data[ptbin][iV0]);
+           }
+   }
+   fAnalysisSummary = new TH1F("Analysis summary","Analysis summary", 26, 0, 26);
+   fAnalysisSummary->GetXaxis()->SetBinLabel(1, "fIsMC");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(2, "fEventMixing");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(3, "fAODAnalysis");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(4, "bayesian purity");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(5, "dca mode");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(6, "fVertex");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(7, "fCentralityMin");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(8, "fCentralityMax");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(9, "V0 centrality");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(10, "TPC centrality");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(11, "fApplyDeltaDipCut");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(12, "POI filterbit");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(13, "RP filterbit");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(14, "PhiMinusPsiMethod");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(15, "SP");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(16, "SPSUB");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(17, "QC");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(18, "EP");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(19, "EP3sub");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(20, "V0_SP");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(21, "eta gap");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(22, "harm");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(23, "highPtMode");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(24, "deltaDipAngle");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(25, "POImaxPt");
+   fAnalysisSummary->GetXaxis()->SetBinLabel(26, "iterator");
+   fOutputList->Add(fAnalysisSummary);
+   fAnalysisSummary->SetBinContent(1, (Float_t)fIsMC);
+   fAnalysisSummary->SetBinContent(2, (Float_t)fEventMixing);
+   fAnalysisSummary->SetBinContent(3, (Float_t)1);
+   fAnalysisSummary->SetBinContent(4, (Float_t)fPIDConfig[6]);
+   fAnalysisSummary->SetBinContent(5, (Float_t)fDCAConfig[0]);
+   fAnalysisSummary->SetBinContent(6, (Float_t)fVertexRange);
+   fAnalysisSummary->SetBinContent(7, (Float_t)fCentralityMin);
+   fAnalysisSummary->SetBinContent(8, (Float_t)fCentralityMax);
+   const char* v = "V0M";   const char* t = "TRK"; // dont compare string literals
+   if(fkCentralityMethod==v) fAnalysisSummary->SetBinContent(9, (Float_t)1);
+   if(fkCentralityMethod==t) fAnalysisSummary->SetBinContent(10, (Float_t)1);
+   fAnalysisSummary->SetBinContent(11, (Float_t)fApplyDeltaDipCut);
+   fAnalysisSummary->SetBinContent(12, (Float_t)fPOICuts->GetAODFilterBit());
+   fAnalysisSummary->SetBinContent(13, (Float_t)fCutsRP->GetAODFilterBit());
+   for(Int_t i(0); i<12;i++)  fAnalysisSummary->SetBinContent(14+i, fAddTaskMacroSummary[i]);
+   fAnalysisSummary->SetBinContent(26, (Float_t)1);
+   PostData(1, fOutputList);
+   // create candidate array
+   fCandidates = new TObjArray(1000);
+   fCandidates->SetOwner(kTRUE);
+   // create and post flowevent
+   fFlowEvent = new AliFlowEvent(10000);
+   PostData(2, fFlowEvent);
+   if(fEventMixing) fPoolManager = InitializeEventMixing();
+}
+//_____________________________________________________________________________
+AliEventPoolManager* AliAnalysisTwoParticleResonanceFlowTask::InitializeEventMixing()
+{
+   // initialize event mixing
+  if(fDebug > 0) cout << " *** InitializeEventMixing() *** " << endl;
+  Int_t _c(0), _v(0);
+  for(Int_t i(0); i < 19; i++) {
+      if (fCentralityMixingBins[i+1] < fCentralityMixingBins[i]) { _c = i; break; }
+      else _c = 19;
+  }
+  for(Int_t i(0); i < 19; i++) {
+      if (fVertexMixingBins[i+1] < fVertexMixingBins[i]) { _v = i; break; }
+      else _v = 19;
+  }
+  if(fDebug > 0 ) cout << Form("   --> found %d centrality bins for mixing, %d vertex bins for mixing", _c, _v) << endl;
+  Float_t centralityBins[_c];
+  Float_t vertexBins[_v];
+  for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
+  for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
+  return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
+}
+//_____________________________________________________________________________
+template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::InvariantMass(const T* track1, const T* track2) const
+{
+   // Return the invariant mass of two tracks
+   if(fDebug > 1) cout << " *** InvariantMass() *** " << endl;
+   if ((!track2) || (!track1)) return 0.;
+   Float_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
+   Float_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
+   Float_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
+   Float_t e1 = TMath::Sqrt(track1->P() * track1->P() + track1->Mass());
+   Float_t e2 = TMath::Sqrt(track2->P() * track2->P() + track2->Mass());
+   Float_t es = TMath::Power((e1 + e2), 2);
+   if ((es - (pxs + pys + pzs)) < 0) return 0.;
+   return TMath::Sqrt((es - (pxs + pys + pzs)));
+}
+//_____________________________________________________________________________
+template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::DeltaDipAngle(const T* track1, const T* track2) const
+{
+   // Calculate the delta dip angle between two particles
+   if(fDebug > 1) cout << " *** DeltaDipAngle() *** " << endl;
+   if (track1->P()*track2->P() == 0) return 999;
+   return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckDeltaDipAngle(const T* track1, const T* track2) const
+{
+   // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
+   if(fDebug > 1) cout << " *** CheckDeltaDipAngle() *** " << endl;
+   if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PairPt(track1, track2) < fDeltaDipPt)) return kFALSE;
+   return kTRUE;
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
+{
+   // Check if pair passes eta and pt cut
+   if(fDebug > 1) cout << " *** CheckCandidateEtaPtCut() *** " << endl;
+   if (fCandidateMinPt > PairPt(track1, track2) || fCandidateMaxPt < PairPt(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;
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::EventCut(T* event)
+{
+   // Impose event cuts
+   if(fDebug > 0) cout << " *** EventCut() *** " << endl;
+   if (!event) return kFALSE;
+   if (!CheckVertex(event)) return kFALSE;
+   if (!CheckCentrality(event)) return kFALSE;
+   if(fQA) PlotMultiplcities(event);
+   return kTRUE;
+}
+//_____________________________________________________________________________
+template <typename T> void AliAnalysisTwoParticleResonanceFlowTask::PlotMultiplcities(const T* event) const
+{
+   // QA multiplicity plots
+   if(fDebug > 1) cout << " *** PlotMultiplcities() *** " << endl;
+   fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
+   fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
+   fTPCM->Fill(event->GetNumberOfTracks());
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckVertex(const T* event)
+{
+   // Check if event vertex is within given range
+   if(fDebug > 0) cout << " *** CheckVertex() *** " << endl;
+   if (!event->GetPrimaryVertex()) return 0x0;
+   fVertex = event->GetPrimaryVertex()->GetZ();
+   if (TMath::Abs(fVertex) > fVertexRange) return 0x0;
+   return kTRUE;
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckCentrality(T* event)
+{
+   // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
+   if(fDebug > 0) cout << " *** CheckCentrality() *** " << endl;
+   if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
+   fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
+   if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax)) {
+      if(fQA) fCentralityNoPass->Fill(fCentrality) ;
+      return kFALSE;
+   }
+   fCentralityPass->Fill(fCentrality);
+   return kTRUE;
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::InitializeBayesianPID(AliAODEvent* event)
+{
+   // Initialize the Bayesian PID object for AOD
+   if(fDebug > 0) cout << " *** InitializeBayesianPID() *** " << endl;
+   fBayesianResponse->SetDetResponse(event, fCentrality);
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::PassesTPCbayesianCut(T* track, Int_t species) const
+{
+   // Check if the particle passes the TPC TOF bayesian cut.
+   if(fDebug > 1) cout << " *** PassesTPCbayesianCut() *** " << endl;
+   fBayesianResponse->ComputeProb(track);
+   if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
+   Int_t contamination(0); //type of particle we don't want in the sample
+   if(species==3)  contamination = 2;  //we want to reject pions when selection kaons
+   if(species==2)  contamination = 3; // we want to reject kaons when selecting pions
+   Float_t *probabilities = fBayesianResponse->GetProb();
+   if (probabilities[species] > fPIDConfig[6]  && (probabilities[contamination] < probabilities[species])) { // 3 is kaon, 2 is pion
+     if(fQA) {fPhi->Fill(track->Phi()); fPt->Fill(track->Pt()); fEta->Fill(track->Eta());}
+     return kTRUE;
+   }
+   return kFALSE;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTwoParticleResonanceFlowTask::PassesDCACut(AliAODTrack* track) const
+{
+    // check if track passes dca cut according to dca routine
+    // setup the routine as follows:
+    // fDCAConfig[0] < -1 no pt dependence
+    // fDCAConfig[0] =  0 do nothing
+    // fDCAConfig[0] >  1 pt dependent dca cut
+    if(fDebug > 1) cout << " *** PassesDCACut() *** " << endl;
+    if(fIsMC) return kTRUE;
+    if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
+        if(fQA) {
+            Double_t b[2] = { -99., -99.};
+            Double_t bCov[3] = { -99., -99., -99.};
+            track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov);
+            fDCAXYQA->Fill(b[0]);
+            fDCAZQA->Fill(b[1]);
+            fDCAPrim->Fill(track->Pt(), b[0]);
+        }
+        return kTRUE;
+    }
+    Double_t b[2] = { -99., -99.};
+    Double_t bCov[3] = { -99., -99., -99.};
+    track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov);
+    if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
+    if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
+    if(fDCAConfig[0] > .9) {
+        if(fDCAConfig[4] < TMath::Abs(b[1])) return kFALSE;
+        Float_t denom = TMath::Power(track->Pt(), fDCAConfig[3]);
+        if( denom < 0.0000001 ) return kFALSE; // avoid division by zero
+        if( (fDCAConfig[1] + fDCAConfig[2] / denom) < TMath::Abs(b[0]) ) return kFALSE;
+    }
+    if(fQA) {
+        fDCAXYQA->Fill(b[0]);
+        fDCAZQA->Fill(b[1]);
+        fDCAPrim->Fill(track->Pt(), b[0]);
+    }
+    return kTRUE;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTwoParticleResonanceFlowTask::AcceptTrack(AliAODTrack* track, Int_t species) const
+{
+   if(fDebug > 1) cout << " *** IsTrack() *** " << endl;
+   if (!PassesDCACut(track)) return kFALSE;
+   if (PassesTPCbayesianCut(track, species)) {
+       if(fQA) if(species==3) fPIDk->Fill(track->P(), track->GetTPCsignal());
+       if(fQA) if(species==2) fPIDp->Fill(track->P(), track->GetTPCsignal());
+       return kTRUE;
+   }
+   return kFALSE;
+}
+//_____________________________________________________________________________
+template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::PairPt(const T* track1, const T* track2, Bool_t phi) const
+{
+   // return p_t of track pair, select phi to return the azimuthal angle instead
+   TVector3 a(track1->Px(), track1->Py(), track1->Pz());
+   TVector3 b(track2->Px(), track2->Py(), track2->Pz());
+   TVector3 c = a + b;
+   if(phi) return c.Phi();
+   return c.Pt();
+}
+//_____________________________________________________________________________
+template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::PtSelector(Int_t tracktype, const T* track1, const T* track2, Float_t mass) const
+{
+   // plot m_inv spectra of like- and unlike sign pairs
+   Float_t pt = PairPt(track1, track2);
+   if (tracktype == 0) { // signal histograms
+       for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
+           if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
+               fResonanceSignal[ptbin]->Fill(mass);
+               if(fQA) fPtSpectra[ptbin]->Fill(pt);
+           }
+       }
+   }
+   if (tracktype == 1) for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) fResonanceBackground[ptbin]->Fill(mass);
+   return pt;
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::QualityCheck(T* track) const
+{
+   // check if track meets quality cuts
+   if(!track) return kFALSE;
+   return fPOICuts->IsSelected(track);
+}
+//_____________________________________________________________________________
+template <typename T> void AliAnalysisTwoParticleResonanceFlowTask::SetNullCuts(T* event)
+{
+   // Set null cuts
+   if (fDebug > 0) cout << " *** SetNullCuts() *** " << fCutsRP << endl;
+   fCutsRP->SetEvent(event, MCEvent());
+   fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
+   fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
+   fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
+   fNullCuts->SetEvent(event, MCEvent());
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::PrepareFlowEvent(Int_t iMulti)
+{
+   // Prepare flow events
+   if (fDebug > 0 ) cout << " *** PrepareFlowEvent() *** " << endl;
+   fFlowEvent->ClearFast();
+   fFlowEvent->Fill(fCutsRP, fNullCuts);
+   fFlowEvent->SetReferenceMultiplicity(iMulti);
+   fFlowEvent->DefineDeadZone(0, 0, 0, 0);
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::PhiMinusPsiMethod(TObjArray* MixingCandidates)
+{
+    // extract v2 using the phi minus psi method
+    if(fDebug > 0) cout << " ** PhiMinusPsiMethod() *** " << endl;
+    if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
+        if(fDebug > 0 ) cout << "Fatal error: unable to retrieve VZERO task output !!! Exiting ..." << endl;
+        return;
+    }
+    // retrieve data
+    Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
+    fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1]))); // vzeroa - tpc
+    fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2]))); // vzeroa - vzeroc
+    fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2]))); // tpc - vzeroc
+//    if(abcPsi2[0] < 0) abcPsi2[0] += TMath::Pi();
+//    if(abcPsi2[2] < 0) abcPsi2[2] += TMath::Pi();
+    // if event plane stuff went ok, fill the histograms necessary for flow analysis with phi - psi method
+    if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
+    AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
+    if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
+    if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
+        Int_t nEvents = pool->GetCurrentNEvents();
+        if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
+        for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
+            TObjArray* mixed_candidates = pool->GetEvent(iEvent);
+            if(!mixed_candidates) continue; // this should NEVER happen
+            Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
+            Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
+            if(fDebug > 0) cout << Form("   - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
+            TObjArray* SpeciesA = new TObjArray();
+            SpeciesA->SetOwner(kTRUE);
+            TObjArray* SpeciesB = new TObjArray();
+            SpeciesB->SetOwner(kTRUE);
+            TObjArray* SpeciesAFromBuffer = new TObjArray();
+            SpeciesAFromBuffer->SetOwner(kTRUE);
+            TObjArray* SpeciesBFromBuffer = new TObjArray();
+            SpeciesBFromBuffer->SetOwner(kTRUE);
+            for (Int_t iTracks = 0; iTracks < candidates; iTracks++) { 
+                AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)MixingCandidates->At(iTracks);
+                if (!track) continue;
+                if (track->Charge() == fChargeA && track->Species() == fSpeciesA) SpeciesA->Add(track);
+                else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesB->Add(track);
+            }
+            for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) { 
+                AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)mixed_candidates->At(iTracks);
+                if (!track) continue;
+                if (track->Charge() == fChargeA && track->Species() == fSpeciesA ) SpeciesAFromBuffer->Add(track);
+                else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesBFromBuffer->Add(track);
+            }
+            PhiMinusPsiMethodWriteData(kTRUE, SpeciesA, SpeciesB, abcPsi2); // write signal information to histograms
+            PhiMinusPsiMethodWriteData(kFALSE, SpeciesA, SpeciesBFromBuffer, abcPsi2);
+            PhiMinusPsiMethodWriteData(kFALSE, SpeciesAFromBuffer, SpeciesB, abcPsi2);
+        } // end of mixed events loop
+    } // end of checking to see whether pool is filled correctly
+    pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
+    if(fDebug > 0 ) pool->PrintInfo();
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::PhiMinusPsiMethodWriteData(Bool_t signal, TObjArray* SpeciesA, TObjArray* SpeciesB, Float_t* abcPsi2)
+{
+    // write data for phi minus psi method into correct  histograms
+    if(fDebug > 0 ) cout << " *** PhiMinusPsiMethodWriteData *** " << endl;
+    for (Int_t i = 0; i < SpeciesA->GetEntriesFast(); i++) { // create pairs
+        for(Int_t j = 0; j < SpeciesB->GetEntriesFast(); j++) { 
+            AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)(SpeciesA->At(i));
+            AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)(SpeciesB->At(j));
+            if(!(trackA||trackB)) continue; // shouldn't happen
+            if(trackA->ID() == trackB->ID()) continue; // remove autocorrelations 
+            if(fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
+           if(fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
+            Float_t minv = InvariantMass(trackA, trackB);
+            Float_t pt = PtSelector(2, trackA, trackB, minv);
+            TVector3 a(trackA->Px(), trackA->Py(), trackA->Pz());
+           TVector3 b(trackB->Px(), trackB->Py(), trackB->Pz());
+           TVector3 c = a + b;
+           Float_t phi = c.Phi();
+//            if(phi < 0) phi += TMath::Pi();
+            Float_t dPhi_a = phi - abcPsi2[0];
+            Float_t dPhi_c = phi - abcPsi2[2];
+//            if (dPhi_a < 0) dPhi_a += TMath::Pi();
+//            if (dPhi_c < 0) dPhi_c += TMath::Pi();
+            // now all necessary numbers are here, so we can put them into histograms
+            for(Int_t k(0); k < fNdPhiBins; k++) { // dPhi bin loop
+                if(dPhi_a >= fdPhiBins[k] && dPhi_a < fdPhiBins[k+1]) // see if track is in desired dPhi_a range
+                for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
+                    if(pt <= fPtBins[ptbin] || pt > fPtBins[ptbin+1]) continue;
+                    signal ? fPhiMinusPsiDataContainer[k][ptbin][0]->Fill(minv) : fPhiMinusPsiBackgroundContainer[k][ptbin][0]->Fill(minv);
+                }
+                if(dPhi_c >= fdPhiBins[k] && dPhi_c < fdPhiBins[k+1]) // see if track is in desired dPhi_c range
+                for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
+                    if(pt <= fPtBins[ptbin] || pt > fPtBins[ptbin+1]) continue;
+                    signal ? fPhiMinusPsiDataContainer[k][ptbin][1]->Fill(minv) : fPhiMinusPsiBackgroundContainer[k][ptbin][1]->Fill(minv);
+                }
+            }
+        } // end of loop on i- tracks
+    } // end of loop on i+ tracks
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::VZEROSubEventAnalysis()
+{
+    // vzero event plane analysis using three subevents
+    if(fDebug > 0) cout << " ** VZEROSubEventAnalysis() *** " << endl;
+    if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
+        if(fDebug > 0 ) cout << "Fatal error: unable to retrieve VZERO task output !!! Exiting ..." << endl;
+        return;
+    }
+    // retrieve data
+    Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
+    for(Int_t i(0); i < 3; i++) if(abcPsi2[i] == 0)  {
+        if(fDebug > 0) cout << " Warning: I've encountered 0 in one of the symmetry panes (TPC, VZERO) - skipping event !" << endl;
+            return;
+    }
+    // save info for resolution calculations
+    fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1]))); // vzeroa - tpc
+    fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2]))); // vzeroa - vzeroc
+    fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2]))); // tpc - vzeroc
+    Float_t minv[fMassBins+1];
+    Float_t _dphi[fMassBins][fNPtBins][2]; //minv, pt, v0a-c
+    Int_t occurence[fMassBins][fNPtBins]; //minv, pt
+    Float_t _inc = (fMaxMass - fMinMass) / (Float_t)fMassBins;
+    if(_inc==0) return; // avoid division by zero later on
+    for(Int_t mb(0); mb < fMassBins+1; mb++) minv[mb] = fMinMass + mb * _inc;
+    for(Int_t i(0); i < fMassBins; i++) for (Int_t j(0); j < fNPtBins; j++) for(Int_t k(0); k < 2; k ++) {
+        _dphi[i][j][k] = 0;
+        if(k==0) occurence[i][j] = 0;
+    }
+    for(Int_t iCand(0); iCand < fCandidates->GetEntriesFast(); iCand++) { // loop over unlike sign tracks
+        AliFlowTrackSimple *track = dynamic_cast<AliFlowTrackSimple*>(fCandidates->At(iCand));
+        if(!track) {
+            if(fDebug > 1) cout << " --> dynamic_cast returned null-pointer, skipping candidate <-- " << endl;
+            continue;
+        }
+        for(Int_t mb(0); mb < fMassBins; mb++) { // loop over massbands
+            if(track->Mass() <= minv[mb] || track->Mass() > minv[mb+1]) continue;
+            for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
+                if(track->Pt() <= fPtBins[ptbin] || track->Pt() > fPtBins[ptbin+1]) continue;
+                _dphi[mb][ptbin][0]+=TMath::Cos(2.*(track->Phi() - abcPsi2[0]));
+                _dphi[mb][ptbin][1]+=TMath::Cos(2.*(track->Phi() - abcPsi2[2]));
+                occurence[mb][ptbin]+=1;
+            }
+        }
+        for(Int_t mb(0); mb < fMassBins; mb++) // write vn values to tprofiles
+            for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
+                if(occurence[mb][ptbin]==0) continue;
+                fV0Data[ptbin][0]->Fill(mb*_inc+fMinMass+0.5*_inc, _dphi[mb][ptbin][0]/(Float_t)occurence[mb][ptbin]);
+                fV0Data[ptbin][1]->Fill(mb*_inc+fMinMass+0.5*_inc, _dphi[mb][ptbin][1]/(Float_t)occurence[mb][ptbin]);
+            }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::UserExec(Option_t *)
+{
+   // UserExec: called for each event. Commented where necessary
+   if(fDebug > 0 ) cout << " *** UserExec() *** " << endl;
+   TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
+   if(fEventMixing) {
+       MixingCandidates = new TObjArray();
+       MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
+   }
+
+   TObjArray* SpeciesA = new TObjArray(); // create arrays for the helper tracks
+   SpeciesA->SetOwner(kTRUE);
+   TObjArray* ocSpeciesA = new TObjArray(); // opposite charge particles
+   ocSpeciesA->SetOwner(kTRUE);
+   TObjArray* SpeciesB = new TObjArray();
+   SpeciesB->SetOwner(kTRUE);
+   TObjArray* ocSpeciesB = new TObjArray();
+   ocSpeciesB->SetOwner(kTRUE);
+
+   if (!fPIDResponse) { // kill the event if there isn't a pid response
+      if(fDebug > 0 ) cout << " Could not get PID response " << endl;
+      return;
+   }
+
+   fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
+   if (fAOD) {
+      if (!EventCut(fAOD)) return; // check for event cuts
+      InitializeBayesianPID(fAOD); // init event objects
+      SetNullCuts(fAOD);
+      Int_t iTracks = fAOD->GetNumberOfTracks();
+      PrepareFlowEvent(iTracks); // does number of tracks correspond to the set filterbit ??? FIXME !!!
+      fCandidates->SetLast(-1);
+      if(fIsMC) IsMC(); // launch mc stuff FIXME
+      if(fQA) fEventStats->Fill(0);
+      for (Int_t i = 0; i < iTracks; i++) { // select analysis candidates
+         AliAODTrack* track = fAOD->GetTrack(i);
+         if (!QualityCheck(track)) continue; // reject poor quality tracks
+         if (fQA) { // fill qa histo's
+            if(track->Charge() > 0) {fEventStats->Fill(1); fPtP->Fill(track->Pt());}
+            if(track->Charge() < 0) {fEventStats->Fill(2); fPtN->Fill(track->Pt());}
+         }
+         if (track->Charge() == fChargeA && AcceptTrack(track, fSpeciesA)) { // store species a
+             SpeciesA->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(), 
+                                                           track->Pt(), track->Charge(), fMassA, track->GetID(), fSpeciesA));
+            if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(),  track->Px(), 
+                                                                                    track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassA, 
+                                                                                    track->GetID(), fSpeciesA));
+        }
+         if (track->Charge() == -1*fChargeA && AcceptTrack(track, fSpeciesA)) { // store opposite charge species a
+             ocSpeciesA->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(), 
+                                                             track->Pt(), track->Charge(), fMassA, track->GetID(), fSpeciesA));
+            if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), 
+                                                                                    track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassA, 
+                                                                                    track->GetID(), fSpeciesA));
+        }
+        if (track->Charge() == fChargeB && AcceptTrack(track, fSpeciesB)) { // store species b
+             SpeciesB->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(), 
+                                                           track->Pt(), track->Charge(), fMassB, track->GetID(), fSpeciesB));
+            if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), 
+                                                                                    track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassB, 
+                                                                                    track->GetID(), fSpeciesB));
+        }
+         if (track->Charge() == -1*fChargeB && AcceptTrack(track, fSpeciesB)) { // store opposite charge species b
+             ocSpeciesB->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(), 
+                                                             track->Pt(), track->Charge(), fMassB, track->GetID(), fSpeciesB));
+            if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), 
+                                                                                    track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassB, 
+                                                                                    track->GetID(), fSpeciesB));
+        }
+      }
+//      cout << SpeciesA->GetEntriesFast() << endl;
+//      cout << SpeciesB->GetEntriesFast() << endl;
+//      cout << ocSpeciesA->GetEntriesFast() << endl;
+//      cout << ocSpeciesB->GetEntriesFast() << endl;
+      // do the phi minus psi method if that's specified FIXME currenlty only supports event mixing
+      if (fPhiMinusPsiMethod) { // for the phi minus psi method, no need for flow package etc 
+          PhiMinusPsiMethod(MixingCandidates); // call the method
+          PostData(1, fOutputList); // save the output data
+          return; // return to retrieve next event
+      }
+      // start the resonance reconstruction, this fills the fCandidates array with flow tracks 
+      ResonanceSignal(SpeciesA, SpeciesB);
+      // 
+      if(fV0) VZEROSubEventAnalysis();
+      if (fDebug > 0)  printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates into flow events
+      for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
+         AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
+         if (!cand) continue;
+         if (fDebug > 1) 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>1) 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()) {
+                  if (fDebug > 1) printf("      was in RP set");
+                  iRP->SetForRPSelection(kFALSE);
+                  fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
+               }
+            }
+            if (fDebug > 1) printf("\n");
+         }
+         cand->SetForPOISelection(kTRUE);
+         fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
+      }
+      if(!fEventMixing) { // do the combinatorial background
+          ResonanceBackground(ocSpeciesA, SpeciesB); // mix opposite charge species A with species B
+          ResonanceBackground(SpeciesA, ocSpeciesB); // mix species A with opposite charge species B
+      }
+      delete SpeciesA;
+      delete SpeciesB;
+      delete ocSpeciesA;
+      delete ocSpeciesB; // clean heap memory
+      if(fEventMixing) ReconstructionWithEventMixing(MixingCandidates);
+      PostData(1, fOutputList);
+      PostData(2, fFlowEvent);
+   }
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::ResonanceSignal(TObjArray* SpeciesA, TObjArray* SpeciesB) const
+{
+    // fill signal histograms
+      for (Int_t i(0); i < SpeciesA->GetEntriesFast(); i++) { //track loop over species A
+       for (Int_t j(0); j < SpeciesB->GetEntriesFast(); j++) { //track loop over species B
+          AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)SpeciesA->At(i);
+          AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)SpeciesB->At(j);
+          if(!(trackA||trackB)) continue; // shouldn't happen
+          if(trackA->ID()==trackB->ID()) continue; // beware of autocorrelations
+         if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
+         if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
+         Float_t mass = InvariantMass(trackA, trackB);
+         Float_t pt = PtSelector(0, trackA, trackB, mass);
+         TVector3 a(trackA->Px(), trackA->Py(), trackA->Pz());
+         TVector3 b(trackB->Px(), trackB->Py(), trackB->Pz());
+         TVector3 c = a + b;
+         Float_t phi = c.Phi();
+         Float_t eta = c.Eta();
+         Int_t nIDs[2];
+         nIDs[0] = trackA->ID();
+         nIDs[1] = trackB->ID();
+         MakeTrack(mass, pt, phi, eta, 2, nIDs);
+        }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::ResonanceBackground(TObjArray* SpeciesA, TObjArray* SpeciesB, Bool_t checkAutoCorrelations) const
+{
+    // fill background histograms
+     for (Int_t i(0); i < SpeciesA->GetEntriesFast(); i++) { //track loop over species A
+       for (Int_t j(0); j < SpeciesB->GetEntriesFast(); j++) { //track loop over species B
+          AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)SpeciesA->At(i);
+          AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)SpeciesB->At(j);
+          if(!(trackA||trackB)) continue; // shouldn't happen
+          if((trackA->ID() == trackB->ID()) && checkAutoCorrelations) continue; // remove autocorrelations 
+         if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
+         if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
+          PtSelector(1, trackA, trackB, InvariantMass(trackA, trackB));
+        }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::ReconstructionWithEventMixing(TObjArray* MixingCandidates) const
+{
+    // perform reconstruction with event mixing
+    if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
+    AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
+    if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
+    if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
+        Int_t nEvents = pool->GetCurrentNEvents();
+        if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
+        for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
+            TObjArray* mixed_candidates = pool->GetEvent(iEvent);
+            if(!mixed_candidates) continue; // this should NEVER happen
+            Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
+            Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
+            if(fDebug > 0) cout << Form("   - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
+            TObjArray* SpeciesA = new TObjArray();
+            SpeciesA->SetOwner(kTRUE);
+            TObjArray* SpeciesB = new TObjArray();
+            SpeciesB->SetOwner(kTRUE);
+            TObjArray* SpeciesAFromBuffer = new TObjArray();
+            SpeciesAFromBuffer->SetOwner(kTRUE);
+            TObjArray* SpeciesBFromBuffer = new TObjArray();
+            SpeciesBFromBuffer->SetOwner(kTRUE);
+            for (Int_t iTracks = 0; iTracks < candidates; iTracks++) { 
+                AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)MixingCandidates->At(iTracks);
+                if (!track) continue;
+                if (track->Charge() == fChargeA && track->Species() == fSpeciesA) SpeciesA->Add(track);
+                else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesB->Add(track);
+            }
+            for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) { 
+                AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)mixed_candidates->At(iTracks);
+                if (!track) continue;
+                if (track->Charge() == fChargeA && track->Species() == fSpeciesA ) SpeciesAFromBuffer->Add(track);
+                else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesBFromBuffer->Add(track);
+            }
+            ResonanceBackground(SpeciesA, SpeciesBFromBuffer, kFALSE);
+            ResonanceBackground(SpeciesB, SpeciesAFromBuffer, kFALSE); 
+        } // end of mixed events loop
+    } // end of checking to see whether pool is filled correctly
+    pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
+    if(fDebug > 0 ) pool->PrintInfo();
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::Terminate(Option_t *)
+{
+    // terminate
+    if(fDebug > 0) cout << " *** Terminate() *** " << endl;
+}
+//______________________________________________________________________________
+void  AliAnalysisTwoParticleResonanceFlowTask::MakeTrack(Float_t mass, Float_t pt, Float_t phi, Float_t eta, Int_t nDau, Int_t iID[]) const
+{
+   // Construct Flow Candidate Track from two selected candidates
+   if(fDebug > 1 ) cout << " *** MakeTrack() *** " << endl;
+   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]);
+   sTrack->SetForPOISelection(kTRUE);
+   sTrack->SetForRPSelection(kFALSE);
+   if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
+   else fCandidates->AddLast(sTrack);
+   return;
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::IsMC()
+{
+    ForceExit(1,"No MC mode available at this stage of developement ... ");
+}
+//_____________________________________________________________________________
diff --git a/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h b/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h
new file mode 100644 (file)
index 0000000..c4f78da
--- /dev/null
@@ -0,0 +1,243 @@
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. */
+/* See cxx source for full Copyright notice */
+/* $Id$ */
+
+// AliAnalysisTwoParticleResonanceFlowTask:
+// origin: Redmer Alexander Bertens (rbertens@nikhef.nl)
+// analyis task for Resonance-meson reconstruction and estimation of v_n
+
+#ifndef ALIANALYSISTWOPARTICLERESONANCEFLOWTASK_H
+#define ALIANALYSISTWOPARTICLERESONANCEFLOWTASK_H
+
+class TH1F;
+class TH2F;
+class TProfile;
+class AliFlowTrackCuts;
+class AliFlowEvent;
+class AliFlowCandidateTrack;
+class AliFlowBayesianPID;
+class AliEventPoolManager;
+class AliResonanceFlowHelperTrack;
+
+#include "AliAnalysisTaskSE.h"
+class AliResonanceFlowHelperTrack : public TObject
+{
+public:
+        AliResonanceFlowHelperTrack(Double_t eta, Double_t phi, Double_t p, Double_t px, Double_t py, Double_t pz, Double_t pt, Int_t charge, Double_t mass, Int_t id, Int_t species) : fEta(eta), fPhi(phi), fp(p), fpX(px), fpY(py), fpZ(pz), fpT(pt), fCharge(charge), fMass(mass), fID(id), fSpecies(species) {  }
+    ~AliResonanceFlowHelperTrack() {}
+    virtual Double_t P()                const { return fp; }
+    virtual Double_t Px()               const { return fpX; }
+    virtual Double_t Py()               const { return fpY; }
+    virtual Double_t Pz()               const { return fpZ; }
+    virtual Double_t Pt()               const { return fpT; }
+    virtual Double_t Phi()              const { return fPhi; }
+    virtual Double_t Eta()              const { return fEta; }
+    virtual Int_t Charge()              const { return fCharge; }
+    virtual Double_t Mass()             const { return fMass; }
+    virtual Short_t ID()                const { return fID; }
+    virtual Int_t Species()             const { return fSpecies; }
+    void    InitializeHelperTrack(Double_t eta, Double_t phi, Double_t p, Double_t px, Double_t py, Double_t pz, Double_t pt, Int_t charge, Double_t mass, Int_t id, Int_t species) { fEta = eta; fPhi = phi; fp = p; fpX = px; fpY = py; fpZ = pz; fpT = pt; fCharge = charge; fMass = mass; fID = id; fSpecies = species; }
+private:
+    Double_t                             fEta;      // eta
+    Double_t                             fPhi;      // phi
+    Double_t                             fp;        // p
+    Double_t                             fpX;       // pX
+    Double_t                             fpY;       // pY
+    Double_t                             fpZ;       // pZ
+    Double_t                             fpT;       // pT
+    Int_t                                fCharge;   // charge
+    Double_t                             fMass;     // mass
+    Short_t                              fID;       // id
+    Int_t                                fSpecies;  // species
+    ClassDef(AliResonanceFlowHelperTrack, 1); // lightweight helper track for phi reconstruction
+};
+
+class AliAnalysisTwoParticleResonanceFlowTask : public AliAnalysisTaskSE
+{
+public:
+   AliAnalysisTwoParticleResonanceFlowTask();
+   AliAnalysisTwoParticleResonanceFlowTask(const char *name);
+   virtual ~AliAnalysisTwoParticleResonanceFlowTask();
+   // technical aspects of the analysis
+   void                                 ForceExit(Int_t type, const char* message);
+   Int_t                                SetDebugLevelResonanceTask(Int_t debug) {fDebug = debug; return fDebug; }
+   Bool_t                               SetIsMC(Bool_t ismc) {fIsMC = ismc; return fIsMC; }
+   void                                 IsMC();
+   Bool_t                               UseEventMixing(Bool_t mix) { fEventMixing = mix; return mix; }
+   Bool_t                               UsePhiMinusPsiMethod(Bool_t p) {fPhiMinusPsiMethod = p; return p;}
+   Bool_t                               SetVZEROSubEvents(Bool_t v0) { fV0 = v0; return v0; }
+   // configure the output of the analysis
+   TH1F*                                BookHistogram(const char * name);
+   TH2F*                                BookPIDHistogram(const char * name, Bool_t TPC);
+   TH1F*                                InitPtSpectraHistograms(Float_t nmin, Float_t nmax);
+   TH1F*                                BookPtHistogram(const char* name);
+   Bool_t                               InitializeAnalysis();
+   void                                 AddResonanceIdentificationOutputObjects();
+   virtual void                         UserCreateOutputObjects();
+   // setters
+   void                                 SetPtBins(Float_t bin[19], Int_t n) { for(Int_t i = 0; i < n+1; i++) fPtBins[i] = bin[i]; fNPtBins = n; }
+   void                                 SetdPhiBins(Float_t bin[19], Int_t n) { for(Int_t i = 0; i < n+1; i++) fdPhiBins[i] = bin[i]; fNdPhiBins = n;}
+   void                                 SetCentralityParameters(Float_t CentralityMin, Float_t CentralityMax, const char* CentralityMethod) { 
+                                                                                          fCentralityMin = CentralityMin; 
+                                                                                          fCentralityMax = CentralityMax; 
+                                                                                          fkCentralityMethod = CentralityMethod; }
+   void                                 SetPOICuts(AliFlowTrackCuts *cutsPOI) { fPOICuts = cutsPOI; }
+   void                                 SetRPCuts(AliFlowTrackCuts *cutsRP) { fCutsRP = cutsRP; }
+   void                                 SetPIDConfiguration(Float_t prob[7]) { for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = prob[i]; }
+   Bool_t                               SetQA(Bool_t qa) {fQA = qa; return fQA;}
+   void                                 SetAddTaskMacroSummary(Float_t m[12]) {for(Int_t i(0); i < 12; i++) fAddTaskMacroSummary[i] = m[i];}
+   void                                 SetPOIDCAXYZ(Float_t dca[5]) { for(Int_t i = 0; i < 5; i++) fDCAConfig[i] = dca[i]; }
+   void                                 SetMixingParameters(Int_t p[3]) { for(Int_t i = 0; i < 3; i++) fMixingParameters[i] = p[i]; }
+   void                                 SetupSpeciesA(Int_t species, Int_t charge, Float_t mass) {fSpeciesA = species; fChargeA = charge; fMassA = mass;}
+   void                                 SetupSpeciesB(Int_t species, Int_t charge, Float_t mass) {fSpeciesB = species; fChargeB = charge; fMassB = mass;}   
+   void                                 SetCandidateEtaAndPt(Float_t mineta, Float_t maxeta, Float_t minpt, Float_t maxpt) { fCandidateMinEta = mineta;
+                                                                                                                                 fCandidateMaxEta = maxeta;
+                                                                                                                                 fCandidateMinPt = minpt;
+                                                                                                                                 fCandidateMaxPt = maxpt;
+                                                                                                                                 fCandidateEtaPtCut = kTRUE;}
+   void                                 SetCommonConstants(Int_t massBins, Float_t minMass, Float_t maxMass) { fMassBins = massBins;
+                                                                                                                 fMinMass = minMass;
+                                                                                                                 fMaxMass = maxMass; }
+   void                                 SetVertexZ(Float_t z) { fVertexRange = z; }
+   void                                 SetMaxDeltaDipAngleAndPt(Float_t a, Float_t pt) { fDeltaDipAngle = a;
+                                                                                          fDeltaDipPt = pt;
+                                                                                          fApplyDeltaDipCut = kTRUE; };
+   void                                 SetMixingBins(Int_t c[20], Int_t v[20]) {for(Int_t i = 0; i < 20; i++) { fCentralityMixingBins[i] = c[i];
+                                                                                                                 fVertexMixingBins[i] = v[i]; } }
+   //getters
+   void                                 GetMixingParameters(Int_t p[3]) const { for(Int_t i = 0; i < 3; i++) p[i] = fMixingParameters[i]; } 
+   Float_t                              GetCenMin() const {return fCentralityMin; }
+   Float_t                              GetCenMax() const {return fCentralityMax; }
+   const char*                          GetCentralityMethod() const {return fkCentralityMethod; }
+   Float_t                              GetVertexZ() const { return fVertexRange; }
+   Float_t                              GetDeltaDipAngle() const {return fDeltaDipAngle; }
+   Float_t                              GetDeltaDipPt() const {return fDeltaDipPt; }
+   void                                 GetPIDConfiguration(Float_t prob[7]) const {for(Int_t i = 0; i < 7; i++) prob[i] = fPIDConfig[i]; }
+   void                                 GetPOIDCZXYZ(Float_t dca[5]) const { for(Int_t i = 0; i < 5; i++) dca[i] = fDCAConfig[i]; }
+   void                                 GetCandidateEtaAndPt(Float_t etapt[4]) const { etapt[0] = fCandidateMinEta;
+                                                                                        etapt[1] = fCandidateMaxEta;
+                                                                                        etapt[2] = fCandidateMinPt;
+                                                                                        etapt[3] = fCandidateMaxPt; }
+   // the analysis itself
+   AliEventPoolManager*                 InitializeEventMixing();
+   template <typename T> Float_t        InvariantMass(const T* track1, const T* track2) const;
+   template <typename T> Float_t        DeltaDipAngle(const T* track1, const T* track2) const;
+   template <typename T> Bool_t         CheckDeltaDipAngle(const T* track1, const T* track2) const;
+   template <typename T> Bool_t         CheckCandidateEtaPtCut(const T* track1, const T* track2) const;
+   template <typename T> Bool_t         EventCut(T* event);
+   template <typename T> void           PlotMultiplcities(const T* event) const;
+   template <typename T> Bool_t         CheckVertex(const T* event);
+   template <typename T> Bool_t         CheckCentrality(T* event);
+   void                                 InitializeBayesianPID(AliAODEvent* event);
+   template <typename T> Bool_t         PassesTPCbayesianCut(T* track, Int_t species) const;
+   Bool_t                               PassesDCACut(AliAODTrack* track) const;
+   Bool_t                               AcceptTrack(AliAODTrack* track, Int_t species) const;
+   template <typename T> Float_t        PairPt(const T* track_1, const T* track_2, Bool_t phi = kFALSE) const;
+   template <typename T> Float_t        PtSelector(Int_t _track_type, const T* track_1, const T* track_2, Float_t mass) const;
+   template <typename T> Bool_t         QualityCheck(T* track) const;
+   template <typename T> void           SetNullCuts(T* esd);
+   void                                 PrepareFlowEvent(Int_t iMulti);
+   void                                 PhiMinusPsiMethod(TObjArray* MixingCandidates);
+   void                                 PhiMinusPsiMethodWriteData(Bool_t signal, TObjArray* SpeciesA, TObjArray* SpeciesB, Float_t* abcPsi2);
+   void                                 VZEROSubEventAnalysis();
+   virtual void                         UserExec(Option_t *option);
+   void                                 ResonanceSignal(TObjArray* SpeciesA, TObjArray* SpeciesB) const;
+   void                                 ResonanceBackground(TObjArray* SpeciesA, TObjArray* SpeciesB, Bool_t checkAutoCorrelations = kTRUE) const;
+   void                                 ReconstructionWithEventMixing(TObjArray* MixingCandidates) const;
+   virtual void                         Terminate(Option_t *);
+
+private:
+
+   Int_t                fSpeciesA; // particle species a
+   Int_t                fSpeciesB; // species b
+   Int_t                fChargeA; // charge for species a
+   Int_t                fChargeB; // charge for species b
+   Float_t              fMassA; // mass for species a
+   Float_t              fMassB; // mass  for species b
+   Int_t                fDebug; // debug level (0 none, 1 fcn calls, 2 verbose)
+   Bool_t               fIsMC; // use mc mode
+   Bool_t               fEventMixing; // use event mixing
+   Bool_t               fPhiMinusPsiMethod; //  use phi minus psi method (default is invariant mass fit method)
+   Bool_t               fQA; // make qa plots
+   Bool_t               fV0; // use three subevents including vzero
+   Int_t                fMassBins; // mass bins
+   Float_t              fMinMass; // mass range
+   Float_t              fMaxMass; // mass range
+   AliFlowTrackCuts     *fCutsRP; // track cuts for reference particles
+   AliFlowTrackCuts     *fNullCuts; // dummy cuts for flow event tracks
+   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               fCandidateEtaPtCut; // set eta and pt cut for candidate tracks and combinatorial background
+   Float_t              fCandidateMinEta; // minimum eta for candidates
+   Float_t              fCandidateMaxEta; // maximum eta for candidates
+   Float_t              fCandidateMinPt; // minimum pt for candidates
+   Float_t              fCandidateMaxPt; // maximum pt for candidates
+   Float_t              fPIDConfig[7]; // configure pid routine
+   Float_t              fDCAConfig[5]; // configure dca routine
+   Int_t                fMixingParameters[3]; // mixing: poolsize, mixing tracks, pool buffer
+   Int_t                fCentralityMixingBins[20]; // configure centrality bins for event mixing
+   Int_t                fVertexMixingBins[20]; // configure vertex bins for event mixing
+   Float_t              fPtBins[19]; // pt bin borders
+   Float_t              fdPhiBins[19]; // dPhi bin borders
+   Int_t                fNPtBins; // no of pt bins + 1
+   Int_t                fNdPhiBins; // no of dphi bins + 1
+   Float_t              fCentrality; // event centrality
+   Float_t              fVertex; // event vertex z 
+   AliAODEvent          *fAOD;    //! AOD oject
+   AliEventPoolManager  *fPoolManager; //! event pool manager
+   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
+   TH2F                 *fPIDp; //! QA histogram of TPC response of pions
+   TH2F                 *fNOPIDTOF; //! QA histo of TOF repsonse charged particles
+   TH2F                 *fPIDTOF; //! QA histo of TOF response kaons
+   TH1F                 *fResonanceSignal[18]; //! signal histograms
+   TH1F                 *fResonanceBackground[18]; //! like-sign kaon pairs
+   TH1F                 *fPtSpectra[18]; //! pt spectra
+   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
+   Float_t              fCentralityMin; // lower bound of cenrality bin
+   Float_t              fCentralityMax; // upper bound of centrality bin
+   const char           *fkCentralityMethod; // method used to determine centrality (V0 by default)
+   AliFlowTrackCuts     *fPOICuts; // cuts for particles of interest (flow package)
+   Float_t              fVertexRange; // absolute value of maximum distance of vertex along the z-axis
+   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)
+   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                 *fDCAAll;//! qa dca of all charged particles
+   TH1F                 *fDCAXYQA; //! qa plot of dca xz
+   TH1F                 *fDCAZQA; //!qa plot of dca z
+   TH2F                 *fDCAPrim; //!dca of primaries (mc) or kaons (data)
+   TH2F                 *fDCASecondaryWeak; //! dca of weak (mc)
+   TH2F                 *fDCAMaterial; //!dca material (mc) all (data)
+   TProfile             *fSubEventDPhiv2; //! subevent resolution info for v2
+   TProfile             *fV0Data[18][2]; //! profiles for vzero vn(minv)
+   TH1F                 *fPhiMinusPsiDataContainer[18][18][2]; //! histograms for phi minus psi results
+   TH1F                 *fPhiMinusPsiBackgroundContainer[18][18][2]; //! histograms for phi minus psi background
+   TH1F                 *fAnalysisSummary; //! plot analysis flags
+   Float_t              fAddTaskMacroSummary[12]; // add task macro summary
+
+   AliAnalysisTwoParticleResonanceFlowTask(const AliAnalysisTwoParticleResonanceFlowTask&); // Not implemented
+   AliAnalysisTwoParticleResonanceFlowTask& operator=(const AliAnalysisTwoParticleResonanceFlowTask&); // Not implemented
+   void                 MakeTrack(Float_t, Float_t, Float_t, Float_t, Int_t , Int_t[]) const;
+
+   ClassDef(AliAnalysisTwoParticleResonanceFlowTask, 1);
+};
+
+#endif
+
+
+
index bd38edb07470ff706112eaa633d11b6ca4d927a8..e5c7de23ca0b277736d14344587fb53a531c6381 100644 (file)
@@ -29,6 +29,8 @@
 #pragma link C++ class AliFlowBayesianPID+;
 #pragma link C++ class AliAnalysisTaskPhiFlow+;
 #pragma link C++ class AliPhiMesonHelperTrack+;
+#pragma link C++ class AliAnalysisTwoParticleResonanceFlowTask+;
+#pragma link C++ class AliResonanceFlowHelperTrack+;
 #pragma link C++ class AliAnalysisTaskFilterFE+;
 #pragma link C++ class AliAnalysisTaskVnV0+;
 #pragma link C++ class AliFlowVZEROResults+;
diff --git a/PWGCF/FLOW/macros/AddTwoParticleResonanceFlowTask.C b/PWGCF/FLOW/macros/AddTwoParticleResonanceFlowTask.C
new file mode 100644 (file)
index 0000000..ab0f12d
--- /dev/null
@@ -0,0 +1,561 @@
+/////////////////////////////////////////////////////////////////////////////////////////////
+//
+// AddTaskPhiFlow macro
+// Author: Redmer A. Bertens, Utrecht University, 2012
+//         rbertens@cern.ch, r.a.bertens@uu.nl
+//         Commented where necessary
+/////////////////////////////////////////////////////////////////////////////////////////////
+
+class AliAnalysisDataContainer;
+class AliFlowTrackCuts;
+class AliFlowTrackSimpleCuts;
+class AliFlowEventCuts;
+class AliFlowEventSimpleCuts;
+class AliAnalysisDataContainer;
+
+AliAnalysisTwoParticleResonanceFlowTask* AddTwoParticleResonanceFlowTask(Bool_t SP = 1, // select flow analysis methods
+                                       Bool_t SPSUB = 0,
+                                       Bool_t QC = 0,
+                                       Bool_t EP = 0,
+                                       Bool_t EP3sub = 0,
+                                       Bool_t VZERO_SP = 0, // use vzero sp method
+                                       Float_t centrMin = 20., // centrality selection
+                                       Float_t centrMax = 30.,
+                                       Float_t ITSsigma = 0., // pid mode (see task implementation)
+                                       Float_t ITSrange = 0.,
+                                       Float_t TPCcontrol = 1.,
+                                       Float_t TPCsigma = 3.,
+                                       Float_t TPCrange = 0.,
+                                       Float_t ITScontrol = -1.,
+                                       Float_t Bpurity = 0.3,
+                                       TString suffixName = "UniqueID", // unique id for output objects
+                                       Bool_t bCentralTrigger = kTRUE, // trigger selection
+                                       Float_t EtaGap = 0., // eta gap for SPSUB
+                                       TString DCA = "pt", // dca mode (see task implementation)
+                                       Int_t harm = 2, // harmonic vn
+                                       UInt_t poi_filter = 1, // aod filterbits
+                                       UInt_t rp_filter = 1,
+                                       Float_t minMass = .71, // min mass for flow analysis
+                                       Float_t maxMaxx = 1.13, // max mass for flow analysis
+                                       Int_t speciesA = 3, // species a
+                                       Int_t speciesB = 2, // species b
+                                       Int_t chargeA = -1, // charge a
+                                       Int_t chargeB = 1, // charge b
+                                       Float_t massA = 4.93676999999999977e-01, // mass species a
+                                       Float_t massB = 1.39570e-01, // mass species b
+                                       Bool_t PhiMinusPsiMethod = 0, // use phi minus psi method
+                                       Bool_t event_mixing = kTRUE, // use event mixing        
+                                       Bool_t highPtMode = kFALSE, // use with caution !!! disables invariant mass fit method
+                                       Float_t deltaPhiMass = 0.0003, // dM in which to look for phi 
+                                       Float_t POIPtMax = 5., // max pt of daughterp particles
+                                       Bool_t shrinkSP = kTRUE, // shrink output
+                                       Bool_t debug = kTRUE) // macro debug mode, for task's debug mode see header
+{
+   /* some defaults that are used frequently:
+    * for phi:
+    *  - minMass = 0.99
+    *  - maxMass = 1.092
+    * for kstar
+    *  - minMass = 0.71
+    *  - maxMass = 1.13
+    * kaon mass = 4.93676999999999977e-01
+    * pion mass = 1.39570e-01
+    */
+   // some defaults that have been removed as function arguments (august 30 2012)
+   Float_t POIPtMin = 0.2;  // pt of daughters
+   Float_t deltaDip = 0.;
+   Float_t deltaDipMaxPt = 0.;
+   Bool_t TPCStandAloneTracks = kFALSE; // deprecated
+   Bool_t useGlobalRPCuts = kTRUE; // deprecated
+   Float_t vertexZ = 10.;
+   Float_t POIEtaMin = -0.8;
+   Float_t POIEtaMax = 0.8;
+   // start of macro
+   Float_t PIDconfig[] = {ITSsigma, ITSrange, TPCcontrol, TPCsigma, TPCrange, ITScontrol, Bpurity};
+   // main function, create and add tasks
+   if(debug) cout << " === Adding Task PhiFlow === " << endl;
+   // set up main output container's name
+   TString fileName = AliAnalysisManager::GetCommonFileName();
+   fileName += ":PhiReconstruction";
+   suffixName += Form("%.0f", centrMin);
+   suffixName += Form("%.0f", centrMax);
+   fileName+=suffixName;
+   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||EP3sub)&&(EtaGap > 0.)) {
+       if(debug) cout << " --> Fatal error: one or more of the flow analyses is not compatible with the Eta Gap !!! <--" << endl;
+       return 0x0;
+   }
+   else if(EP3sub) {
+       if(highPtMode) {
+           if(debug) cout << " --> Can't launch 3 subevent method for high pt analysis, exiting ... <--" << endl;
+           return 0x0;
+       }
+       if(debug) cout << " --> Starting 3 subevent plane method - try at your own risk !!! <-- " << endl;
+       if(!(harm!=2||harm!=3)) {
+           if(debug) cout << " --> Fatal error: can only return v2 and v3 with 3 subevent method! " << endl;
+           return 0x0;
+       }
+   }
+   if((EP3sub||PhiMinusPsiMethod)&&debug) {
+       gROOT->LoadMacro("$ALICE_ROOT/PWGCF/FLOW/macros/AddTaskVZERO.C");
+       AddTaskVZERO(0,0,0,0);
+   }
+   // get the manager and input event handler
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   if (!mgr) {
+      if(debug) cout << " Fatal error: no analysis manager found! " << endl;
+      return 0x0;
+   }
+   if (!mgr->GetInputEventHandler()) {
+      if(debug) cout << " Fatal error: no imput event handler found!" << endl;
+      return 0x0;
+   }
+   if(EP3sub&&debug) { // don't use this on train ! this is why it's only enabled in macro debug mode
+         gROOT->LoadMacro("$ALICE_ROOT/PWGCF/FLOW/macros/AddTaskVZERO.C");
+         AddTaskVZERO(0,0,0,0);
+   }
+   // create the main task
+   AliAnalysisTwoParticleResonanceFlowTask *task = new AliAnalysisTwoParticleResonanceFlowTask("TaskPhiFlow");
+   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;
+   }
+   if(task->UsePhiMinusPsiMethod(PhiMinusPsiMethod)) {
+       cout << " --> Using phi - Psi method <-- ... " << endl;
+       Float_t dphibins[] = {0., 0.63, 1.26, 1.89, 2.52, 3.15};
+       task->SetdPhiBins(dphibins, (Int_t)(sizeof(dphibins)/sizeof(dphibins[1])-1));
+       event_mixing = kTRUE;
+   }
+   task->SetupSpeciesA(speciesA, chargeA, massA);
+   task->SetupSpeciesB(speciesB, chargeB, massB);
+
+   if(task->SetVZEROSubEvents(EP3sub)) cout << " --> Setting up VZERO subevents method ... " << endl;
+   if(event_mixing) {
+      if(debug) cout << " --> Enabeling event mixing for phi reconstruction - try at your own risk !!! <-- " << endl;
+      // set vertex and mixing bins - arrays MUST have length 20!
+      Int_t c[] = {0, 2, 4, 6, 8, 10, 20, 30, 40, 50, 60, 70, 80, 90, 101, 0, 0, 0, 0, 0};
+      Int_t v[] = {-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+      if(((Int_t)(sizeof(c)/sizeof(c[1]))!=20)||((Int_t)(sizeof(v)/sizeof(v[1]))!=20)) {
+          cout << " --> Fatal error: check mixing parameters ! <-- " << endl;
+          return 0x0;
+      }
+      else task->SetMixingBins(c, v);
+   }
+   // 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
+   AliFlowTrackCuts* 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&&(!VZERO_SP)) {
+       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(rp_filter < 9999 ) {
+           if(debug) cout << "  > set RP filterbit " << rp_filter << endl;     
+           cutsRP->SetAODfilterBit(rp_filter);
+       }
+       if(debug) cout << "    --> kGlobal RP's " << cutsRP << endl;
+   }
+   if(VZERO_SP) { // use vzero sub analysis
+       cutsRP = cutsRP->GetStandardVZEROOnlyTrackCuts(); // select vzero tracks
+       SP = kFALSE; // disable other methods
+       SPSUB = kTRUE; // calculate sp_qa and sp_qb
+       QC = kFALSE;
+       EP = kFALSE;
+       EP3sub = kFALSE;
+       EtaGap = 0.; // no eta gap, full tpc poi's
+       if(debug) cout << "    --> VZERO 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;
+   }
+   cutsPOI = cutsPOI->GetStandardGlobalTrackCuts2010();
+   cutsPOI->SetPtRange(POIPtMin, POIPtMax); // pt range of DAUGHTERS !!!
+   cutsPOI->SetMaxDCAToVertexXY(0.3); // FIXME not implemented in aod086 aod095 see PassesDCACut() in implementation
+   cutsPOI->SetMaxDCAToVertexZ(0.3);
+   if(poi_filter < 9999 ) {
+       if(debug) cout << "  > set POI filterbit " << poi_filter << endl;     
+       cutsPOI->SetAODfilterBit(poi_filter);
+   }
+   if(debug) cout << "    --> cutsPOI " << cutsPOI << endl;
+   task->SetPOICuts(cutsPOI);
+   //set POI cuts for aods XY Z - 3 distinct cases
+   Float_t POIDCA[] = {0., 0., 0., 0., 0.};
+   if(DCA == "none" ) { // 1 --- do nothing
+       if (debug) cout << " --> No DCA cut on POI's <-- " << endl;
+       for (Int_t i = 0; i < 5; i++) POIDCA[i] = 0.;
+   }
+   if(DCA == "fix" ) { // 2 --- use fixed values for xy z
+       if (debug) cout << " --> Fixed DCA cut on POI's <-- " << endl;
+       POIDCA[0] = -1.; POIDCA[1] = 0.3; POIDCA[2] = 0.3; POIDCA[3] = 0.; POIDCA[4] = 0.;
+   }
+   if(DCA == "pt" ) { // 3 --- use pt dependent cut
+       if (debug) cout << " --> Pt dependent DCA cut on POI's <-- " << endl;
+       POIDCA[0] = 1.; POIDCA[1] = 0.0105; POIDCA[2] = 0.0350; POIDCA[3] = 1.1; POIDCA[4] = 2.;
+   }
+   task->SetPOIDCAXYZ(POIDCA);
+   Float_t AddTaskMacroSummary[] = {(Float_t)0, (Float_t)SP, (Float_t)SPSUB, (Float_t)QC, (Float_t)EP, (Float_t)EP3sub, (Float_t)VZERO_SP, (Float_t)EtaGap, (Float_t)harm, (Float_t)highPtMode, (Float_t)deltaDip, (Float_t)POIPtMax}; 
+   task->SetAddTaskMacroSummary(AddTaskMacroSummary);
+
+   if(highPtMode) { // high pt loop - loop will end macro
+       Float_t _pt[] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 15.};
+       task->SetPtBins(_pt, (Int_t)(sizeof(_pt)/sizeof(_pt[1]))-1);
+       // general approach: use kinematic filters which select kaon pairs with a certain mass window
+       AliFlowTrackSimpleCuts* HighPtSubEventFilterA = new AliFlowTrackSimpleCuts("HighPtSubEventFilterA"); 
+       HighPtSubEventFilterA->SetEtaMin(-0.8);
+       HighPtSubEventFilterA->SetEtaMax(0.0);
+       HighPtSubEventFilterA->SetMassMin(1.019445 - deltaPhiMass);
+       HighPtSubEventFilterA->SetMassMax(1.019445 + deltaPhiMass);
+       AliFlowTrackSimpleCuts* HighPtSubEventFilterB = new AliFlowTrackSimpleCuts("HighPtSubEventFilterB"); 
+       HighPtSubEventFilterB->SetEtaMin(0.0);
+       HighPtSubEventFilterB->SetEtaMax(+0.8);
+       HighPtSubEventFilterB->SetMassMin(1.019445 - deltaPhiMass);
+       HighPtSubEventFilterB->SetMassMax(1.019445 + deltaPhiMass);
+       AliFlowTrackSimpleCuts* HighPtGenericFilter = new AliFlowTrackSimpleCuts("HighPtGenericFilter");
+       HighPtGenericFilter->SetEtaMin(-0.8);
+       HighPtGenericFilter->SetEtaMax(+0.8);
+       HighPtGenericFilter->SetMassMin(1.019445 - deltaPhiMass);
+       HighPtGenericFilter->SetMassMax(1.019445 + deltaPhiMass);
+       if(debug) cout << "    --> Created poi filters " << 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, 0., 15.);
+       task->SetCentralityParameters(centrMin, centrMax, "TRK");
+       task->SetVertexZ(vertexZ);
+       if(debug) cout << "    --> Set pair cuts and event cuts" << endl;
+       // specify the PID procedure which will be used
+       task->SetPIDConfiguration(PIDconfig);
+       AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+       AliAnalysisDataContainer *coutHist = mgr->CreateContainer(Form("PhiV2_%s", OutputName(centrMin, centrMax,PIDconfig, suffixName, bCentralTrigger, EtaGap, POIEtaMin, POIEtaMax, 0., 15., 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
+       {
+          Int_t mb(999);
+          if(debug) cout << " === RECEIVED REQUEST FOR FLOW ANALYSIS === " << endl;
+          AliAnalysisDataContainer *flowEvent = mgr->CreateContainer(Form("FC%s", suffixName.Data()), AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer);
+          mgr->ConnectOutput(task, 2, flowEvent);
+          if(debug) cout << "    --> Created IO containers " << flowEvent << endl;
+          if(debug) cout << "    --> suffixName " << suffixName << endl;
+          if (QC) {  // add qc tasks
+             AddQCmethod(Form("QCTPCMB_%d_%s", mb, suffixName.Data()), harm, flowEvent, HighPtGenericFilter, debug, 0x0, suffixName.Data());
+             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, HighPtSubEventFilterB, NULL, false, shrinkSP, debug, true, suffixName.Data());
+             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, HighPtSubEventFilterA, NULL, false, shrinkSP, debug, true, suffixName.Data());
+             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, HighPtGenericFilter, NULL, false, shrinkSP, debug, 0x0, suffixName.Data());
+             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, HighPtGenericFilter, NULL, true, shrinkSP, debug, 0x0, suffixName.Data());
+             if(debug) cout << "    --> Hanging EP task ... " << mb << " succes!" << endl;
+          }
+       }
+       // print summary to screen
+       cout << endl << endl << "       ==== AddTaskPhiFlow launched  ===== " << 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 << " ************ Configured DCA setup ************** " << endl;
+       cout << "  DCA type: " << DCA;
+       if (DCA == "") cout << "default";
+       cout << endl << " ************* Task statisti:q!cs ****************** " << endl;
+       cout << "   -> Launched PHI reconstruction " << endl;
+       if(SP) cout << "   --> Launched QaQb SP filters and corresponding SP task " << endl;
+       if(EP) cout << "   --> Launched QaQb QC filters and corresponding EP task " << endl;
+       if(SPSUB) cout << "   --> Launched Qa&&Qb SP filters and corresponding SP tasks " << endl;
+       if(QC) cout << "   --> Launched QaQb QC filters and corresponding QC task " << endl;
+       if(EP3sub) cout << " --> Launched VZERO subevent analysis alongside reconstruction - USE WITH CAUTION!" << endl;
+       cout << " ************************************************ " << endl;
+       TString condit = "";
+       (task->SetQA(kFALSE)) ? condit+= " --> Enabled QA plots <-- " : condit+= " --> Disabled QA plots <-- ";
+       (task->SetIsMC(kFALSE)) ? condit+= " --> MC mode <-- " : condit+= " --> DATA mode <-- ";
+       (task->UseEventMixing(event_mixing)) ? condit+= " --> Using EVENT MIXING <--" : condit+= "--> Combinatorial background <--";
+       cout << condit << endl;
+       cout << "           --> Now go for a coffee! <-- " << endl;
+       cout << " ************************************************ " << endl;
+       return task;
+   } //  end of high pt loop - high-pt task is ready to go at this point
+   Float_t _pt[] = {0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0};
+   task->SetPtBins(_pt, (Int_t)(sizeof(_pt)/sizeof(_pt[1]))-1);
+   // POI filter cuts, will filter invm mass bands and subevents
+   AliFlowTrackSimpleCuts* POIfilterQC[30];
+   AliFlowTrackSimpleCuts* POIfilterSP[30][2];
+   Double_t flowBands[2][30];
+   Double_t _inc = (maxMaxx-minMass)/30.;
+   for (Int_t mb = 0; mb < 30; mb++) {
+      flowBands[0][mb] = minMass + mb * _inc;
+      flowBands[1][mb] = minMass + (mb + 1) * _inc;
+      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.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.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]);
+   }
+   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, 0., 10.);
+   task->SetCentralityParameters(centrMin, centrMax, "TRK");
+   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;
+   }
+   if (TPCStandAloneTracks)
+   { // switch to tpc standalone tracks for POI selection
+      if(debug) cout << "    --> Switching to TPC standalone analsis " << endl;
+      task->SetRequireStrictKaonCuts();
+      task->SetRequireTPCStandAloneKaons();
+   }
+   // 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, 0., 10., 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
+   {
+      if(debug) cout << " === RECEIVED REQUEST FOR FLOW ANALYSIS === " << endl;
+      AliAnalysisDataContainer *flowEvent = mgr->CreateContainer(Form("FC%s", suffixName.Data()), AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer);
+      mgr->ConnectOutput(task, 2, flowEvent);
+      if(debug) cout << "    --> Created IO containers " << flowEvent << endl;
+      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, 0x0, suffixName.Data());
+            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, suffixName.Data(), VZERO_SP);
+            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, suffixName.Data(), VZERO_SP);
+            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, 0x0, suffixName.Data());
+            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, 0x0, suffixName.Data());
+            if(debug) cout << "    --> Hanging EP task ... " << mb << " succes!" << endl;
+         }
+      }
+   }
+   // print summary to screen
+   cout << endl << endl << "       ==== AddTaskPhiFlow launched  ===== " << 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 << " ************ Configured DCA setup ************** " << endl;
+   cout << "  DCA type: " << DCA;
+   if (DCA == "") cout << "default";
+   cout << endl << " ************* Task statisti:q!cs ****************** " << endl;
+   cout << "   -> Launched PHI reconstruction " << 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;
+   if(EP3sub) cout << " --> Launched VZERO subevent analysis alongside reconstruction - USE WITH CAUTION!" << endl;
+   cout << " ************************************************ " << endl;
+   TString condit = "";
+   (task->SetQA(kFALSE)) ? condit+= " --> Enabled QA plots <-- " : condit+= " --> Disabled QA plots <-- ";
+   (task->SetIsMC(kFALSE)) ? condit+= " --> MC mode <-- " : condit+= " --> DATA mode <-- ";
+   (task->UseEventMixing(event_mixing)) ? condit+= " --> Using EVENT MIXING <--" : condit+= "--> Combinatorial background <--";
+   cout << condit << endl;
+   cout << "           --> Now go for a coffee! <-- " << endl;
+   cout << " ************************************************ " << 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, char *suffixName, Bool_t VZERO_SP = kFALSE)
+{
+   // 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";
+   fileName+=suffixName;
+   if(etagap) {
+       fileName+="_SUBEVENTS";
+       if(debug) cout << "    --> Setting up subevent analysis <-- " << endl;
+   }
+   if(debug) cout << "    --> fileName " << fileName << endl;
+   TString myNameSP;
+   (bEP) ? myNameSP = Form("%sEPv%d%s", name, harmonic, Qvector): myNameSP = Form("%sSPv%d%s", name, harmonic, Qvector);
+   if(debug) cout << " Task and filter name: " << 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);
+   if(VZERO_SP) tskFilter->SetSubeventEtaRange(-10, 0, 0, 10);
+   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, char *suffixName)
+{
+   // 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";
+   fileName+=suffixName;
+   if(debug) cout << "    --> Common filename: " << fileName << 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,
+                    Float_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,
+                    TString DCA,
+                    Int_t harm,
+                    Bool_t TPCStandAloneTracks,
+                    Float_t vertexZ,
+                    Bool_t debug,
+                    Bool_t useGlobalRPCuts)
+{
+   // generate output name
+   TString centralityName = "";
+   centralityName += suffixName;
+   centralityName += "_DCA";
+   centralityName += DCA;
+   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;
+}
+//_____________________________________________________________________________
+