PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // AliAnalysisTaskPhiFlow:
  17 // author: Redmer Alexander Bertens (rbertens@nikhef.nl)
  18 // analyis task for phi-meson reconstruction and determination of V2
  19 // handles aod's and esd's transparently
  20
  21 #include "TChain.h"
  22 #include "TTree.h"
  23 #include "TH1.h"
  24 #include "TH1F.h"
  25 #include "TH2F.h"
  26 #include "TCanvas.h"
  27 #include "TMath.h"
  28 #include "TObjArray.h"
  29 #include "AliAnalysisTask.h"
  30 #include "AliAnalysisTaskSE.h"
  31 #include "AliAnalysisManager.h"
  32 #include "AliESDEvent.h"
  33 #include "AliAODEvent.h"
  34 #include "AliESDInputHandler.h"
  35 #include "AliAODInputHandler.h"
  36 #include "AliCentrality.h"
  37 #include "AliVEvent.h"
  38 #include "AliAnalysisTaskPhiFlow.h"
  39 #include "AliFlowBayesianPID.h"
  40 #include "AliStack.h"
  41 #include "AliESDtrackCuts.h"
  42 #include "AliMCEvent.h"
  43 #include "TProfile.h"
  44 #include "AliFlowCandidateTrack.h"
  45 #include "AliFlowTrackCuts.h"
  46 #include "AliFlowEventSimple.h"
  47 #include "AliFlowCommonConstants.h"
  48 #include "AliFlowEvent.h"
  49 #include "TVector3.h"
  50 #include "TRandom2.h"
  51 #include "AliESDVZERO.h"
  52 #include "AliAODVZERO.h"
  53 #include "AliPID.h"
  54 #include "AliPIDResponse.h"
  55
  56 class AliFlowTrackCuts;
  57
  58 ClassImp(AliAnalysisTaskPhiFlow)
  59
  60 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
  61    fDebug(0), fAODAnalysis(0), fMassBins(0), fMinMass(0), fMaxMass(0), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fOldTrackParam(0), fRequireTPCStandAlone(0), fStrictKaonCuts(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCentrality(0), fESD(0), fAOD(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fInvMNP03(0), fInvMPP03(0), fInvMNN03(0), fInvMNP36(0), fInvMPP36(0), fInvMNN36(0), fInvMNP69(0), fInvMPP69(0), fInvMNN69(0), fInvMNP912(0), fInvMPP912(0), fInvMNN912(0), fInvMNP1215(0), fInvMPP1215(0), fInvMNN1215(0), fInvMNP1518(0), fInvMPP1518(0), fInvMNN1518(0), fInvMNP1821(0), fInvMPP1821(0), fInvMNN1821(0), fInvMNP2124(0), fInvMPP2124(0), fInvMNN2124(0), fInvMNP2427(0), fInvMPP2427(0), fInvMNN2427(0), fInvMNP2730(0), fInvMPP2730(0), fInvMNN2730(0), fInvMNP3035(0), fInvMPP3035(0), fInvMNN3035(0), fInvMNP3540(0), fInvMPP3540(0), fInvMNN3540(0), fInvMNP4045(0), fInvMPP4045(0), fInvMNN4045(0), fInvMNP4550(0), fInvMPP4550(0), fInvMNN4550(0), fInvMNP5055(0), fInvMPP5055(0), fInvMNN5055(0), fInvMNP5560(0), fInvMPP5560(0), fInvMNN5560(0), fInvMNP6065(0), fInvMPP6065(0), fInvMNN6065(0), fInvMNP6570(0), fInvMPP6570(0), fInvMNN6570(0), fPtSpectra03(0), fPtSpectra36(0), fPtSpectra69(0), fPtSpectra912(0), fPtSpectra1215(0), fPtSpectra1518(0), fPtSpectra1821(0), fPtSpectra2124(0), fPtSpectra2427(0), fPtSpectra2730(0), fPtSpectra3035(0), fPtSpectra3540(0), fPtSpectra4045(0), fPtSpectra4550(0), fPtSpectra5055(0), fPtSpectra5560(0), fPtSpectra6065(0), fPtSpectra6570(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)
  62 {
  63    // Default constructor
  64    for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = 0.;
  65 }
  66 //_____________________________________________________________________________
  67 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow(const char *name) : AliAnalysisTaskSE(name),
  68     fDebug(0), fAODAnalysis(0), fMassBins(0), fMinMass(0), fMaxMass(0), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fOldTrackParam(0), fRequireTPCStandAlone(0), fStrictKaonCuts(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCentrality(0), fESD(0), fAOD(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fInvMNP03(0), fInvMPP03(0), fInvMNN03(0), fInvMNP36(0), fInvMPP36(0), fInvMNN36(0), fInvMNP69(0), fInvMPP69(0), fInvMNN69(0), fInvMNP912(0), fInvMPP912(0), fInvMNN912(0), fInvMNP1215(0), fInvMPP1215(0), fInvMNN1215(0), fInvMNP1518(0), fInvMPP1518(0), fInvMNN1518(0), fInvMNP1821(0), fInvMPP1821(0), fInvMNN1821(0), fInvMNP2124(0), fInvMPP2124(0), fInvMNN2124(0), fInvMNP2427(0), fInvMPP2427(0), fInvMNN2427(0), fInvMNP2730(0), fInvMPP2730(0), fInvMNN2730(0), fInvMNP3035(0), fInvMPP3035(0), fInvMNN3035(0), fInvMNP3540(0), fInvMPP3540(0), fInvMNN3540(0), fInvMNP4045(0), fInvMPP4045(0), fInvMNN4045(0), fInvMNP4550(0), fInvMPP4550(0), fInvMNN4550(0), fInvMNP5055(0), fInvMPP5055(0), fInvMNN5055(0), fInvMNP5560(0), fInvMPP5560(0), fInvMNN5560(0), fInvMNP6065(0), fInvMPP6065(0), fInvMNN6065(0), fInvMNP6570(0), fInvMPP6570(0), fInvMNN6570(0), fPtSpectra03(0), fPtSpectra36(0), fPtSpectra69(0), fPtSpectra912(0), fPtSpectra1215(0), fPtSpectra1518(0), fPtSpectra1821(0), fPtSpectra2124(0), fPtSpectra2427(0), fPtSpectra2730(0), fPtSpectra3035(0), fPtSpectra3540(0), fPtSpectra4045(0), fPtSpectra4550(0), fPtSpectra5055(0), fPtSpectra5560(0), fPtSpectra6065(0), fPtSpectra6570(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)
  69 {
  70    // Constructor
  71    for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = 0.;
  72
  73    DefineInput(0, TChain::Class());
  74    DefineOutput(1, TList::Class());
  75    DefineOutput(2, AliFlowEventSimple::Class());
  76 }
  77 //_____________________________________________________________________________
  78 AliAnalysisTaskPhiFlow::~AliAnalysisTaskPhiFlow()
  79 {
  80    // Destructor
  81    if (fNullCuts) delete fNullCuts;
  82    if (fOutputList) delete fOutputList;
  83    if (fCandidates) delete fCandidates;
  84    if (fFlowEvent) delete fFlowEvent;
  85    if (fBayesianResponse) delete fBayesianResponse;
  86 }
  87 //_____________________________________________________________________________
  88 TH1F* AliAnalysisTaskPhiFlow::BookHistogram(const char* name)
  89 {
  90    // Return a pointer to a TH1 with predefined binning
  91    TH1F *hist = new TH1F(name, Form("M_{INV} (%s)", name), 60, .99, 1.092);
  92    hist->GetXaxis()->SetTitle("M_{INV} (GeV / c^{2})");
  93    hist->GetYaxis()->SetTitle("No. of pairs");
  94    hist->SetMarkerStyle(kFullCircle);
  95    hist->Sumw2();
  96    fOutputList->Add(hist);
  97    return hist;
  98 }
  99 //_____________________________________________________________________________
 100 TH2F* AliAnalysisTaskPhiFlow::BookPIDHistogram(const char* name)
 101 {
 102    // Return a pointer to a TH2 with predefined binning
 103    TH2F *hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
 104    hist->GetXaxis()->SetTitle("P (GeV / c)");
 105    hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
 106    fOutputList->Add(hist);
 107    return hist;
 108 }
 109 //_____________________________________________________________________________
 110 TH1F* AliAnalysisTaskPhiFlow::InitPtSpectraHistograms(Int_t i)
 111 {
 112    // intialize p_t histograms for each p_t bin
 113    Double_t nmin(0);
 114    Double_t nmax(0);
 115    (i < 11) ? nmin = 0.3 * (i - 1) : nmin = 0.5 * (i - 11) + 3;
 116    (i < 11) ? nmax = 0.3 * (i - 1) + 0.3 : nmax = 0.5 * (i - 11) + 3.5;
 117    TH1F* hist = new TH1F(Form("%f p_{t} %f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 60, nmin, nmax);
 118    hist->GetXaxis()->SetTitle("p_{T} GeV / c");
 119    fOutputList->Add(hist);
 120    return hist;
 121 }
 122 //_____________________________________________________________________________
 123 TH1F* AliAnalysisTaskPhiFlow::BookPtHistogram(const char* name)
 124 {
 125    // Return a pointer to a p_T spectrum histogram
 126    TH1F* ratio = new TH1F(name, name, 100, 0, 7);
 127    ratio->GetXaxis()->SetTitle("p_{T} ( GeV / c^{2} )");
 128    ratio->GetYaxis()->SetTitle("No. of events");
 129    ratio->Sumw2();
 130    fOutputList->Add(ratio);
 131    return ratio;
 132 }
 133 //_____________________________________________________________________________
 134 void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
 135 {
 136    // Add Phi Identification Output Objects
 137    fEventStats = new TH1F("fHistStats", "Event Statistics", 18, -.25, 4.25);
 138    fEventStats->GetXaxis()->SetTitle("No. of events");
 139    fEventStats->GetYaxis()->SetTitle("Statistics");
 140    fOutputList->Add(fEventStats);
 141    fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
 142    fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
 143    fOutputList->Add(fCentralityPass);
 144    fOutputList->Add(fCentralityNoPass);
 145
 146    fNOPID = BookPIDHistogram("TPC signal, all particles");
 147    fPIDk = BookPIDHistogram("TPC signal, kaons");
 148    fInvMNP03 = BookHistogram("NP, 0 < p_{T} < 0.3 GeV");
 149    fInvMPP03 = BookHistogram("PP, 0 < p_{T} < 0.3 GeV");
 150    fInvMNN03 = BookHistogram("NN, 0 < p_{T} < 0.3 GeV");
 151    fInvMNP36 = BookHistogram("NP, 0.3 < p_{T} < 0.6 GeV");
 152    fInvMPP36 = BookHistogram("PP, 0.3 < p_{T} < 0.6 GeV");
 153    fInvMNN36 = BookHistogram("NN, 0.3 < p_{T} < 0.6 GeV");
 154    fInvMNP69 = BookHistogram("NP, 0.6 < p_{T} < 0.9 GeV");
 155    fInvMPP69 = BookHistogram("PP, 0.6 < p_{T} < 0.9 GeV");
 156    fInvMNN69 = BookHistogram("NN, 0.6 < p_{T} < 0.9 GeV");
 157    fInvMNP912 = BookHistogram("NP, 0.9 < p_{T} < 1.2 GeV");
 158    fInvMPP912 = BookHistogram("PP, 0.9 < p_{T} < 1.2 GeV");
 159    fInvMNN912 = BookHistogram("NN, 0.9 < p_{T} < 1.2 GeV");
 160    fInvMNP1215 = BookHistogram("NP, 1.2 < p_{T} < 1.5 GeV");
 161    fInvMPP1215 = BookHistogram("PP, 1.2 < p_{T} < 1.5 GeV");
 162    fInvMNN1215 = BookHistogram("NN, 1.2 < p_{T} < 1.5 GeV");
 163    fInvMNP1518 = BookHistogram("NP, 1.5 < p_{T} < 1.8 GeV");
 164    fInvMPP1518 = BookHistogram("PP, 1.5 < p_{T} < 1.8 GeV");
 165    fInvMNN1518 = BookHistogram("NN, 1.5 < p_{T} < 1.8 GeV");
 166    fInvMNP1821 = BookHistogram("NP, 1.8 < p_{T} < 2.1 GeV");
 167    fInvMPP1821 = BookHistogram("PP, 1.8 < p_{T} < 2.1 GeV");
 168    fInvMNN1821 = BookHistogram("NN, 1.8 < p_{T} < 2.1 GeV");
 169    fInvMNP2124 = BookHistogram("NP, 2.1 < p_{T} < 2.4 GeV");
 170    fInvMPP2124 = BookHistogram("PP, 2.1 < p_{T} < 2.4 GeV");
 171    fInvMNN2124 = BookHistogram("NN, 2.1 < p_{T} < 2.4 GeV");
 172    fInvMNP2427 = BookHistogram("NP, 2.4 < p_{T} < 2.7 GeV");
 173    fInvMPP2427 = BookHistogram("PP, 2.4 < p_{T} < 2.7 GeV");
 174    fInvMNN2427 = BookHistogram("NN, 2.4 < p_{T} < 2.7 GeV");
 175    fInvMNP2730 = BookHistogram("NP, 2.7 < p_{T} < 3.0 GeV");
 176    fInvMPP2730 = BookHistogram("PP, 2.7 < p_{T} < 3.0 GeV");
 177    fInvMNN2730 = BookHistogram("NN, 2.7 < p_{T} < 3.0 GeV");
 178    fInvMNP3035 = BookHistogram("NP, 3.0 < p_{T} < 3.5 GeV");
 179    fInvMPP3035 = BookHistogram("PP, 3.0 < p_{T} < 3.5 GeV");
 180    fInvMNN3035 = BookHistogram("NN, 3.0 < p_{T} < 3.5 GeV");
 181    fInvMNP3540 = BookHistogram("NP, 3.5 < p_{T} < 4.0 GeV");
 182    fInvMPP3540 = BookHistogram("PP, 3.5 < p_{T} < 4.0 GeV");
 183    fInvMNN3540 = BookHistogram("NN, 3.5 < p_{T} < 4.0 GeV");
 184    fInvMNP4045 = BookHistogram("NP, 4.0 < p_{T} < 4.5 GeV");
 185    fInvMPP4045 = BookHistogram("PP, 4.0 < p_{T} < 4.5 GeV");
 186    fInvMNN4045 = BookHistogram("NN, 4.0 < p_{T} < 4.5 GeV");
 187    fInvMNP4550 = BookHistogram("NP, 4.5 < p_{T} < 5.0 GeV");
 188    fInvMPP4550 = BookHistogram("PP, 4.5 < p_{T} < 5.0 GeV");
 189    fInvMNN4550 = BookHistogram("NN, 4.5 < p_{T} < 5.0 GeV");
 190    fInvMNP5055 = BookHistogram("NP, 5.0 < p_{T} < 5.5 GeV");
 191    fInvMPP5055 = BookHistogram("PP, 5.0 < p_{T} < 5.5 GeV");
 192    fInvMNN5055 = BookHistogram("NN, 5.0 < p_{T} < 5.5 GeV");
 193    fInvMNP5560 = BookHistogram("NP, 5.5 < p_{T} < 6.0 GeV");
 194    fInvMPP5560 = BookHistogram("PP, 5.5 < p_{T} < 6.0 GeV");
 195    fInvMNN5560 = BookHistogram("NN, 5.5 < p_{T} < 6.0 GeV");
 196    fInvMNP6065 = BookHistogram("NP, 6.0 < p_{T} < 6.5 GeV");
 197    fInvMPP6065 = BookHistogram("PP, 6.0 < p_{T} < 6.5 GeV");
 198    fInvMNN6065 = BookHistogram("NN, 6.0 < p_{T} < 6.5 GeV");
 199    fInvMNP6570 = BookHistogram("NP, 6.5 < p_{T} < 7.0 GeV");
 200    fInvMPP6570 = BookHistogram("PP, 6.5 < p_{T} < 7.0 GeV");
 201    fInvMNN6570 = BookHistogram("NN, 6.5 < p_{T} < 7.0 GeV");
 202
 203    fPtSpectra03 = InitPtSpectraHistograms(1);
 204    fPtSpectra36 = InitPtSpectraHistograms(2);
 205    fPtSpectra69 = InitPtSpectraHistograms(3);
 206    fPtSpectra912 = InitPtSpectraHistograms(4);
 207    fPtSpectra1215 = InitPtSpectraHistograms(5);
 208    fPtSpectra1518 = InitPtSpectraHistograms(6);
 209    fPtSpectra1821 = InitPtSpectraHistograms(7);
 210    fPtSpectra2124 = InitPtSpectraHistograms(8);
 211    fPtSpectra2427 = InitPtSpectraHistograms(9);
 212    fPtSpectra2730 = InitPtSpectraHistograms(10);
 213    fPtSpectra3035 = InitPtSpectraHistograms(11);
 214    fPtSpectra3540 = InitPtSpectraHistograms(12);
 215    fPtSpectra4045 = InitPtSpectraHistograms(13);
 216    fPtSpectra4550 = InitPtSpectraHistograms(14);
 217    fPtSpectra5055 = InitPtSpectraHistograms(15);
 218    fPtSpectra5560 = InitPtSpectraHistograms(16);
 219    fPtSpectra6065 = InitPtSpectraHistograms(17);
 220    fPtSpectra6570 = InitPtSpectraHistograms(18);
 221
 222    fPtP = BookPtHistogram("i^{+}");
 223    fPtN = BookPtHistogram("i^{-}");
 224    fPtKP = BookPtHistogram("K^{+}");
 225    fPtKN = BookPtHistogram("K^{-}");
 226
 227    fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
 228    fOutputList->Add(fPhi);
 229
 230    fPt = new TH1F("fPt", "p_{T}", 100, 0, 5.5);
 231    fOutputList->Add(fPt);
 232
 233    fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
 234    fOutputList->Add(fEta);
 235
 236    fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
 237    fOutputList->Add(fVZEROA);
 238
 239    fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
 240    fOutputList->Add(fVZEROC);
 241
 242    fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
 243    fOutputList->Add(fTPCM);
 244 }
 245 //_____________________________________________________________________________
 246 void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
 247 {
 248    // Create user defined output objects
 249    fNullCuts = new AliFlowTrackCuts("null_cuts");
 250    fBayesianResponse = new AliFlowBayesianPID();
 251    Double_t t(0);
 252    for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
 253    if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis!!!");
 254
 255    AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 256    cc->SetNbinsMult(10000);
 257    cc->SetMultMin(0);
 258    cc->SetMultMax(10000);
 259
 260    cc->SetNbinsPt(100);
 261    cc->SetPtMin(0);
 262    cc->SetPtMax(10);
 263
 264    cc->SetNbinsPhi(180);
 265    cc->SetPhiMin(0.0);
 266    cc->SetPhiMax(TMath::TwoPi());
 267
 268    cc->SetNbinsEta(200);
 269    cc->SetEtaMin(-5.0);
 270    cc->SetEtaMax(+5.0);
 271
 272    cc->SetNbinsQ(500);
 273    cc->SetQMin(0.0);
 274    cc->SetQMax(3.0);
 275
 276    cc->SetNbinsMass(fMassBins);
 277    cc->SetMassMin(fMinMass);
 278    cc->SetMassMax(fMaxMass);
 279
 280    // setup initial state of PID response object
 281    if (!fOldTrackParam) fBayesianResponse->SetNewTrackParam();
 282
 283    AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 284    if (man)
 285    {
 286       AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 287       if (inputHandler)   fPIDResponse = inputHandler->GetPIDResponse();
 288    }
 289
 290    // Create all output objects and store them to a list
 291    fOutputList = new TList();
 292    fOutputList->SetOwner(kTRUE);
 293
 294    // Create and post the Phi identification output objects
 295    AddPhiIdentificationOutputObjects();
 296    PostData(1, fOutputList);
 297
 298    // create candidate array
 299    fCandidates = new TObjArray(1000);
 300    fCandidates->SetOwner(kTRUE);
 301
 302    // create and post flowevent
 303    fFlowEvent = new AliFlowEvent(10000);
 304    PostData(2, fFlowEvent);
 305 }
 306 //_____________________________________________________________________________
 307 template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* track1, const T* track2) const
 308 {
 309    // Return the invariant mass of two tracks, assuming both tracks are kaons
 310    if ((!track2) || (!track1)) return 0.;
 311    Double_t masss = TMath::Power(4.93676999999999977e-01, 2);
 312    Double_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
 313    Double_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
 314    Double_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
 315    Double_t e1 = TMath::Sqrt(track1->P() * track1->P() + masss);
 316    Double_t e2 = TMath::Sqrt(track2->P() * track2->P() + masss);
 317    Double_t es = TMath::Power((e1 + e2), 2);
 318    if ((es - (pxs + pys + pzs)) < 0) return 0.;
 319    return TMath::Sqrt((es - (pxs + pys + pzs)));
 320 }
 321 //_____________________________________________________________________________
 322 template <typename T> Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const T* track1, const T* track2) const
 323 {
 324    // Calculate the delta dip angle between two particles
 325    if (track1->P()*track2->P() == 0) return 999;
 326    return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
 327 }
 328 //_____________________________________________________________________________
 329 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const T* track1, const T* track2) const
 330 {
 331    // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
 332    if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PhiPt(track1, track2) < fDeltaDipPt)) return kFALSE;
 333    return kTRUE;
 334 }
 335 //_____________________________________________________________________________
 336 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
 337 {
 338    // Check if pair passes eta and pt cut
 339    if (fCandidateMinPt > PhiPt(track1, track2) || fCandidateMaxPt < PhiPt(track1, track2)) return kFALSE;
 340    TVector3 a(track1->Px(), track1->Py(), track1->Pz());
 341    TVector3 b(track2->Px(), track2->Py(), track2->Pz());
 342    TVector3 c = a + b;
 343    if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
 344    return kTRUE;
 345 }
 346 //_____________________________________________________________________________
 347 void AliAnalysisTaskPhiFlow::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
 348 {
 349    // Set a centrality range ]min, max] and define the method to use for centrality selection
 350    fCentralityMin = CentralityMin;
 351    fCentralityMax = CentralityMax;
 352    fkCentralityMethod = CentralityMethod;
 353 }
 354 //_____________________________________________________________________________
 355 template <typename T> Bool_t AliAnalysisTaskPhiFlow::EventCut(T* event)
 356 {
 357    // Impose event cuts
 358    if (!event) return kFALSE;
 359    if (!CheckVertex(event)) return kFALSE;
 360    if (!CheckCentrality(event)) return kFALSE;
 361    PlotVZeroMultiplcities(event);
 362    return kTRUE;
 363 }
 364 //_____________________________________________________________________________
 365 template <typename T> void AliAnalysisTaskPhiFlow::PlotVZeroMultiplcities(const T* event) const
 366 {
 367    // QA multiplicity plots
 368    fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
 369    fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
 370 }
 371 //_____________________________________________________________________________
 372 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckVertex(const T* event) const
 373 {
 374    // Check if event vertex is within given range
 375    if (!event->GetPrimaryVertex()) return 0x0;
 376    if (TMath::Abs((event->GetPrimaryVertex())->GetZ()) > fVertexRange) return 0x0;
 377    return kTRUE;
 378 }
 379 //_____________________________________________________________________________
 380 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCentrality(T* event)
 381 {
 382    // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
 383    if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
 384    fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
 385    if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
 386    {
 387       fCentralityNoPass->Fill(fCentrality) ;
 388       return kFALSE;
 389    }
 390    fCentralityPass->Fill(fCentrality);
 391    return kTRUE;
 392 }
 393 //_____________________________________________________________________________
 394 void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliESDEvent* event)
 395 {
 396    // Initialize the Bayesian PID object for ESD analysis
 397    fBayesianResponse->SetDetResponse(event, fCentrality, AliESDpid::kTOF_T0, kTRUE);
 398    if (fDebug) cout << " --> Initialized Bayesian ESD PID object " << endl;
 399 }
 400 //_____________________________________________________________________________
 401 void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliAODEvent* event)
 402 {
 403    // Initialize the Bayesian PID object for AOD
 404    fBayesianResponse->SetDetResponse(event, fCentrality);
 405    if (fDebug) cout << " --> Initialized Bayesian AOD PID object " << endl;
 406 }
 407 //_____________________________________________________________________________
 408 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* track) const
 409 {
 410    // Check if the particle passes the TPC TOF bayesian cut.
 411    if ((!fStrictKaonCuts) && (!PassesStrictKaonCuts(track))) return kFALSE;
 412    fBayesianResponse->ComputeProb(track);
 413    if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
 414    Float_t *probabilities = fBayesianResponse->GetProb();
 415    if (probabilities[3] > fPIDConfig[6])
 416    {
 417       fPhi->Fill(track->Phi());
 418       fPt->Fill(track->Pt());
 419       fEta->Fill(track->Eta());
 420       return kTRUE;
 421    }
 422    return kFALSE;
 423 }
 424 //_____________________________________________________________________________
 425 Bool_t AliAnalysisTaskPhiFlow::PassesStrictKaonCuts(AliESDtrack* track) const
 426 {
 427    // propagate dca from TPC
 428    Double_t b[2] = { -99., -99.};
 429    Double_t bCov[3] = { -99., -99., -99.};
 430    if (!track->PropagateToDCA(fESD->GetPrimaryVertex(), fESD->GetMagneticField(), 100., b, bCov)) return kFALSE;
 431    if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)) return kFALSE;
 432    return kTRUE;
 433 }
 434 //_____________________________________________________________________________
 435 Bool_t AliAnalysisTaskPhiFlow::PassesStrictKaonCuts(AliAODTrack* track) const
 436 {
 437    // propagate dca from TPC
 438    Double_t b[2] = { -99., -99.};
 439    Double_t bCov[3] = { -99., -99., -99.};
 440    if (!track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
 441    if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)) return kFALSE;
 442    return kTRUE;
 443 }
 444 //_____________________________________________________________________________
 445 Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliESDtrack* track) const
 446 {
 447    // Check if particle is a kaon according to method set in steering macro
 448    if (fRequireTPCStandAlone && (track->GetStatus()&AliESDtrack::kTPCin) == 0) return kFALSE;
 449    if (fPIDConfig[1]!=0 || fPIDConfig[4]!=0) AliFatal("TPC || ITS PID not available in ESD anlaysis -- terminating analysis !!!");
 450    if (PassesTPCbayesianCut(track))
 451    {
 452      fPIDk->Fill(track->P(), track->GetTPCsignal());
 453      return kTRUE;
 454    }
 455    return kFALSE;
 456 }
 457 //_____________________________________________________________________________
 458 Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
 459 {
 460    // Kaon identification routine, based on multiple detectors and approaches
 461    if (fRequireTPCStandAlone && (!track->TestFilterBit(1))) return kFALSE;
 462    fNOPID->Fill(track->P(), track->GetTPCsignal());
 463    if(track->Pt() < fPIDConfig[1])
 464    {
 465        if(fDebug) cout << " ITS received track with p_t " << track->Pt() << endl;
 466        // if tpc control is disabled, pure its pid
 467        if(fPIDConfig[2] < 0.)
 468        {
 469            if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) return kTRUE;
 470            return kFALSE;
 471        }
 472        // else, switch to ITS pid with TPC rejection of protons and pions
 473        if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) < fPIDConfig[3]) return kFALSE;
 474        else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) < fPIDConfig[3]) return kFALSE;
 475        else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0])
 476        {
 477            fPIDk->Fill(track->P(), track->GetTPCsignal());
 478            return kTRUE;
 479        }
 480        return kFALSE;
 481    }
 482    if ((track->Pt() > fPIDConfig[1]) && (track->Pt() < fPIDConfig[4]))
 483    {
 484        if(fDebug) cout << " TPC received track with p_t " << track->Pt() << endl;
 485        // if its control is disabled, pure tpc pid
 486        if(fPIDConfig[5] < 0.)
 487        {
 488            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) return kTRUE;
 489            return kFALSE;
 490        }
 491        // else, switch to TPC pid with ITS rejection of protons and pions
 492        if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kProton)) < fPIDConfig[0]) return kFALSE;
 493        else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kPion)) < fPIDConfig[0]) return kFALSE;
 494        else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3])
 495        {
 496            fPIDk->Fill(track->P(), track->GetTPCsignal());
 497            return kTRUE;
 498        }
 499        return kFALSE;
 500    }
 501    if(fDebug) cout << " Bayesian method received track with p_t " << track->Pt() << endl;
 502    // switch to bayesian PID
 503    if (PassesTPCbayesianCut(track))
 504    {
 505        fPIDk->Fill(track->P(), track->GetTPCsignal());
 506        return kTRUE;
 507    }
 508    return kFALSE;
 509 }
 510 //_____________________________________________________________________________
 511 template <typename T> Double_t AliAnalysisTaskPhiFlow::PhiPt(const T* track1, const T* track2) const
 512 {
 513    // Calculate transverse momentum (p_T) of two tracks
 514    TVector3 a(track1->Px(), track1->Py(), track1->Pz());
 515    TVector3 b(track2->Px(), track2->Py(), track2->Pz());
 516    TVector3 c = a + b;
 517    return c.Pt();
 518 }
 519 //_____________________________________________________________________________
 520 template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, const T* track1, const T* track2) const
 521 {
 522    // Request transverse momentum (p_T), then fill invariant mass histograms as a function of p_T
 523    Double_t pt = PhiPt(track1, track2);
 524    if (tracktype == 0)
 525    {
 526       if ((0.0 <= pt) && (0.3 > pt))
 527       {
 528          fInvMNP03->Fill(InvariantMass(track1, track2));
 529          fPtSpectra03->Fill(pt);
 530       }
 531       if ((0.3 <= pt) && (0.6 > pt))
 532       {
 533          fInvMNP36->Fill(InvariantMass(track1, track2));
 534          fPtSpectra36->Fill(pt);
 535       }
 536       if ((0.6 <= pt) && (0.9 > pt))
 537       {
 538          fInvMNP69->Fill(InvariantMass(track1, track2));
 539          fPtSpectra69->Fill(pt);
 540       }
 541       if ((0.9 <= pt) && (1.2 > pt))
 542       {
 543          fInvMNP912->Fill(InvariantMass(track1, track2));
 544          fPtSpectra912->Fill(pt);
 545       }
 546       if ((1.2 <= pt) && (1.5 > pt))
 547       {
 548          fInvMNP1215->Fill(InvariantMass(track1, track2));
 549          fPtSpectra1215->Fill(pt);
 550       }
 551       if ((1.5 <= pt) && (1.8 > pt))
 552       {
 553          fInvMNP1518->Fill(InvariantMass(track1, track2));
 554          fPtSpectra1518->Fill(pt);
 555       }
 556       if ((1.8 <= pt) && (2.1 > pt))
 557       {
 558          fInvMNP1821->Fill(InvariantMass(track1, track2));
 559          fPtSpectra1821->Fill(pt);
 560       }
 561       if ((2.1 <= pt) && (2.4 > pt))
 562       {
 563          fInvMNP2124->Fill(InvariantMass(track1, track2));
 564          fPtSpectra2124->Fill(pt);
 565       }
 566       if ((2.4 <= pt) && (2.7 > pt))
 567       {
 568          fInvMNP2427->Fill(InvariantMass(track1, track2));
 569          fPtSpectra2427->Fill(pt);
 570       }
 571       if ((2.7 <= pt) && (3.0 > pt))
 572       {
 573          fInvMNP2730->Fill(InvariantMass(track1, track2));
 574          fPtSpectra2730->Fill(pt);
 575       }
 576       if ((3.0 <= pt) && (3.5 > pt))
 577       {
 578          fInvMNP3035->Fill(InvariantMass(track1, track2));
 579          fPtSpectra3035->Fill(pt);
 580       }
 581       if ((3.5 <= pt) && (4.0 > pt))
 582       {
 583          fInvMNP3540->Fill(InvariantMass(track1, track2));
 584          fPtSpectra3540->Fill(pt);
 585       }
 586       if ((4.0 <= pt) && (4.5 > pt))
 587       {
 588          fInvMNP4045->Fill(InvariantMass(track1, track2));
 589          fPtSpectra4045->Fill(pt);
 590       }
 591       if ((4.5 <= pt) && (5.0 > pt))
 592       {
 593          fInvMNP4550->Fill(InvariantMass(track1, track2));
 594          fPtSpectra4550->Fill(pt);
 595       }
 596       if ((5.0 <= pt) && (5.5 > pt))
 597       {
 598          fInvMNP5055->Fill(InvariantMass(track1, track2));
 599          fPtSpectra5055->Fill(pt);
 600       }
 601       if ((5.5 <= pt) && (6.0 > pt))
 602       {
 603          fInvMNP5560->Fill(InvariantMass(track1, track2));
 604          fPtSpectra5560->Fill(pt);
 605       }
 606       if ((6.0 <= pt) && (6.5 > pt))
 607       {
 608          fInvMNP6065->Fill(InvariantMass(track1, track2));
 609          fPtSpectra6065->Fill(pt);
 610       }
 611       if ((6.5 <= pt) && (7.0 > pt))
 612       {
 613          fInvMNP6570->Fill(InvariantMass(track1, track2));
 614          fPtSpectra6570->Fill(pt);
 615       }
 616    }
 617    if (tracktype == 1)
 618    {
 619       if ((0.0 <= pt) && (0.3 > pt)) fInvMPP03->Fill(InvariantMass(track1, track2));
 620       if ((0.3 <= pt) && (0.6 > pt)) fInvMPP36->Fill(InvariantMass(track1, track2));
 621       if ((0.6 <= pt) && (0.9 > pt)) fInvMPP69->Fill(InvariantMass(track1, track2));
 622       if ((0.9 <= pt) && (1.2 > pt)) fInvMPP912->Fill(InvariantMass(track1, track2));
 623       if ((1.2 <= pt) && (1.5 > pt)) fInvMPP1215->Fill(InvariantMass(track1, track2));
 624       if ((1.5 <= pt) && (1.8 > pt)) fInvMPP1518->Fill(InvariantMass(track1, track2));
 625       if ((1.8 <= pt) && (2.1 > pt)) fInvMPP1821->Fill(InvariantMass(track1, track2));
 626       if ((2.1 <= pt) && (2.4 > pt)) fInvMPP2124->Fill(InvariantMass(track1, track2));
 627       if ((2.4 <= pt) && (2.7 > pt)) fInvMPP2427->Fill(InvariantMass(track1, track2));
 628       if ((2.7 <= pt) && (3.0 > pt)) fInvMPP2730->Fill(InvariantMass(track1, track2));
 629       if ((3.0 <= pt) && (3.5 > pt)) fInvMPP3035->Fill(InvariantMass(track1, track2));
 630       if ((3.5 <= pt) && (4.0 > pt)) fInvMPP3540->Fill(InvariantMass(track1, track2));
 631       if ((4.0 <= pt) && (4.5 > pt)) fInvMPP4045->Fill(InvariantMass(track1, track2));
 632       if ((4.5 <= pt) && (5.0 > pt)) fInvMPP4550->Fill(InvariantMass(track1, track2));
 633       if ((5.0 <= pt) && (5.5 > pt)) fInvMPP5055->Fill(InvariantMass(track1, track2));
 634       if ((5.5 <= pt) && (6.0 > pt)) fInvMPP5560->Fill(InvariantMass(track1, track2));
 635       if ((6.0 <= pt) && (6.5 > pt)) fInvMPP6065->Fill(InvariantMass(track1, track2));
 636       if ((6.5 <= pt) && (7.0 > pt)) fInvMPP6570->Fill(InvariantMass(track1, track2));
 637    }
 638    if (tracktype == 2)
 639    {
 640       if ((0.0 <= pt) && (0.3 > pt)) fInvMNN03->Fill(InvariantMass(track1, track2));
 641       if ((0.3 <= pt) && (0.6 > pt)) fInvMNN36->Fill(InvariantMass(track1, track2));
 642       if ((0.6 <= pt) && (0.9 > pt)) fInvMNN69->Fill(InvariantMass(track1, track2));
 643       if ((0.9 <= pt) && (1.2 > pt)) fInvMNN912->Fill(InvariantMass(track1, track2));
 644       if ((1.2 <= pt) && (1.5 > pt)) fInvMNN1215->Fill(InvariantMass(track1, track2));
 645       if ((1.5 <= pt) && (1.8 > pt)) fInvMNN1518->Fill(InvariantMass(track1, track2));
 646       if ((1.8 <= pt) && (2.1 > pt)) fInvMNN1821->Fill(InvariantMass(track1, track2));
 647       if ((2.1 <= pt) && (2.4 > pt)) fInvMNN2124->Fill(InvariantMass(track1, track2));
 648       if ((2.4 <= pt) && (2.7 > pt)) fInvMNN2427->Fill(InvariantMass(track1, track2));
 649       if ((2.7 <= pt) && (3.0 > pt)) fInvMNN2730->Fill(InvariantMass(track1, track2));
 650       if ((3.0 <= pt) && (3.5 > pt)) fInvMNN3035->Fill(InvariantMass(track1, track2));
 651       if ((3.5 <= pt) && (4.0 > pt)) fInvMNN3540->Fill(InvariantMass(track1, track2));
 652       if ((4.0 <= pt) && (4.5 > pt)) fInvMNN4045->Fill(InvariantMass(track1, track2));
 653       if ((4.5 <= pt) && (5.0 > pt)) fInvMNN4550->Fill(InvariantMass(track1, track2));
 654       if ((5.0 <= pt) && (5.5 > pt)) fInvMNN5055->Fill(InvariantMass(track1, track2));
 655       if ((5.5 <= pt) && (6.0 > pt)) fInvMNN5560->Fill(InvariantMass(track1, track2));
 656       if ((6.0 <= pt) && (6.5 > pt)) fInvMNN6065->Fill(InvariantMass(track1, track2));
 657       if ((6.5 <= pt) && (7.0 > pt)) fInvMNN6570->Fill(InvariantMass(track1, track2));
 658    }
 659 }
 660 //_____________________________________________________________________________
 661 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
 662 {
 663    // Check if track is suitable for phi flow analysis
 664    if (!track) return kFALSE;
 665    return fPOICuts->IsSelected(track);
 666 }
 667 //_____________________________________________________________________________
 668 template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
 669 {
 670    // Set null cuts
 671    if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
 672    fCutsRP->SetEvent(event, MCEvent());
 673    fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
 674    fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
 675    fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
 676    fNullCuts->SetEvent(event, MCEvent());
 677 }
 678 //_____________________________________________________________________________
 679 void AliAnalysisTaskPhiFlow::PrepareFlowEvent(Int_t iMulti)
 680 {
 681    // Prepare flow events
 682    fFlowEvent->ClearFast();
 683    fFlowEvent->Fill(fCutsRP, fNullCuts);
 684    fFlowEvent->SetReferenceMultiplicity(iMulti);
 685    fFlowEvent->DefineDeadZone(0, 0, 0, 0);
 686 }
 687 //_____________________________________________________________________________
 688 void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
 689 {
 690    // UserExec: called for each event. Commented where necessary
 691    // check for AOD data type
 692    if (!fPIDResponse)
 693    {
 694       AliError("Cannot get pid response");
 695       return;
 696    }
 697
 698    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 699    if (fAOD)
 700    {
 701       fAODAnalysis = kTRUE;
 702       // Check whether event passes event cuts
 703       if (!EventCut(fAOD)) return;
 704       // initialize event objects
 705       InitializeBayesianPID(fAOD);
 706       SetNullCuts(fAOD);
 707       PrepareFlowEvent(fAOD->GetNumberOfTracks());
 708       fCandidates->SetLast(-1);
 709       // Calculate event plane Q vectors and event plane resolution
 710       fEventStats->Fill(0);
 711       Int_t unTracks = fAOD->GetNumberOfTracks();
 712       AliAODTrack* un[unTracks];
 713       AliAODTrack* up[unTracks];
 714       Int_t unp(0);
 715       Int_t unn(0);
 716       // Loop through tracks, check for species (Kaons), fill arrays according to charge
 717       for (Int_t iTracks = 0; iTracks < unTracks; iTracks++)
 718       {
 719          AliAODTrack* track = fAOD->GetTrack(iTracks);
 720          if (!PhiTrack(track)) continue;
 721          if (fStrictKaonCuts && (!PassesStrictKaonCuts(track))) continue;
 722          Bool_t charge = kFALSE;
 723          if (track->Charge() > 0)
 724          {
 725             charge = kTRUE;
 726             fEventStats->Fill(1);
 727             fPtP->Fill(track->Pt());
 728          }
 729          if (track->Charge() < 0)
 730          {
 731             fEventStats->Fill(2);
 732             fPtN->Fill(track->Pt());
 733          }
 734          if (IsKaon(track))
 735          {
 736             if (charge)
 737             {
 738                up[unp] = track;
 739                unp++;
 740                fEventStats->Fill(3);
 741                fPtKP->Fill(track->Pt());
 742             }
 743             if (!charge)
 744             {
 745                un[unn] = track;
 746                unn++;
 747                fEventStats->Fill(4);
 748                fPtKN->Fill(track->Pt());
 749             }
 750          }
 751       }
 752       // Calculate invariant mass of like- and unlike sign pairs as a function of p_T, store data in histograms
 753       for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
 754       {
 755          for (Int_t nTracks = 0; nTracks < unn ; nTracks++)
 756          {
 757             if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
 758             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
 759             PtSelector(0, up[pTracks], un[nTracks]);
 760             Double_t pt = PhiPt(up[pTracks], un[nTracks]);
 761             Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
 762             TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
 763             TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
 764             TVector3 c = a + b;
 765             Double_t phi = c.Phi();
 766             Double_t eta = c.Eta();
 767             Int_t nIDs[2];
 768             nIDs[0] = up[pTracks]->GetID();
 769             nIDs[1] = un[nTracks]->GetID();
 770             MakeTrack(mass, pt, phi, eta, 2, nIDs);
 771          }
 772       }
 773
 774       if (fDebug)  printf("I received %d candidates\n", fCandidates->GetEntriesFast());
 775       for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand)
 776       {
 777          AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
 778          if (!cand) continue;
 779          if (fDebug) printf(" >Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
 780          for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau)
 781          {
 782             if (fDebug) printf("  >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
 783             for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs)
 784             {
 785                AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
 786                if (!iRP) continue;
 787                if (!iRP->InRPSelection()) continue;
 788                if (cand->GetIDDaughter(iDau) == iRP->GetID())
 789                {
 790                   if (fDebug) printf(" was in RP set");
 791                   iRP->SetForRPSelection(kFALSE);
 792                   fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
 793                }
 794             }
 795             if (fDebug) printf("\n");
 796          }
 797          cand->SetForPOISelection(kTRUE);
 798          fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
 799       }
 800       if (fDebug) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
 801
 802
 803       for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
 804       {
 805          for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++)
 806          {
 807             if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
 808             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
 809             PtSelector(1, up[pTracks], up[nTracks]);
 810          }
 811       }
 812       for (Int_t nTracks = 0; nTracks < unn ; nTracks++)
 813       {
 814          for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++)
 815          {
 816             if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
 817             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
 818             PtSelector(2, un[nTracks], un[pTracks]);
 819          }
 820       }
 821       PostData(1, fOutputList);
 822       PostData(2, fFlowEvent);
 823    }
 824
 825    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 826    if (fESD)
 827    {
 828       fAODAnalysis = kFALSE;
 829       // Check whether event passes event cuts
 830       if (!EventCut(fESD)) return;
 831       InitializeBayesianPID(fESD);
 832       SetNullCuts(fESD);
 833       PrepareFlowEvent(fESD->GetNumberOfTracks());
 834       // Calculate event plane Q vectors and event plane resolution
 835       fEventStats->Fill(0);
 836       Int_t unTracks = fESD->GetNumberOfTracks();
 837       AliESDtrack* un[unTracks];
 838       AliESDtrack* up[unTracks];
 839       Int_t unp(0);
 840       Int_t unn(0);
 841       // Loop through tracks, check for species (Kaons), fill arrays according to charge
 842       for (Int_t iTracks = 0; iTracks < unTracks; iTracks++)
 843       {
 844          AliESDtrack* track = fESD->GetTrack(iTracks);
 845          if (!PhiTrack(track)) continue;
 846          Bool_t charge = kFALSE;
 847          if (track->Charge() > 0)
 848          {
 849             charge = kTRUE;
 850             fEventStats->Fill(1);
 851             fPtP->Fill(track->Pt());
 852          }
 853          if (track->Charge() < 0)
 854          {
 855             fEventStats->Fill(2);
 856             fPtN->Fill(track->Pt());
 857          }
 858          if (IsKaon(track))
 859          {
 860             if (charge)
 861             {
 862                up[unp] = track;
 863                unp++;
 864                fEventStats->Fill(3);
 865                fPtKP->Fill(track->Pt());
 866             }
 867             if (!charge)
 868             {
 869                un[unn] = track;
 870                unn++;
 871                fEventStats->Fill(4);
 872                fPtKN->Fill(track->Pt());
 873             }
 874          }
 875       }
 876       // Calculate invariant mass of like- and unlike sign pairs as a function of p_T, store data in histograms
 877         for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
 878       {
 879          for (Int_t nTracks = 0; nTracks < unn ; nTracks++)
 880          {
 881             if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
 882             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
 883             PtSelector(0, up[pTracks], un[nTracks]);
 884             Double_t pt = PhiPt(up[pTracks], un[nTracks]);
 885             Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
 886             TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
 887             TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
 888             TVector3 c = a + b;
 889             Double_t phi = c.Phi();
 890             Double_t eta = c.Eta();
 891             Int_t nIDs[2];
 892             nIDs[0] = up[pTracks]->GetID();
 893             nIDs[1] = un[nTracks]->GetID();
 894             MakeTrack(mass, pt, phi, eta, 2, nIDs);
 895          }
 896       }
 897
 898       if (fDebug)  printf("I received %d candidates\n", fCandidates->GetEntriesFast());
 899       for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand)
 900       {
 901          AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
 902          if (!cand) continue;
 903          if (fDebug) printf(" >Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
 904          for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau)
 905          {
 906             if (fDebug) printf("  >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
 907             for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs)
 908             {
 909                AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
 910                if (!iRP) continue;
 911                if (!iRP->InRPSelection()) continue;
 912                if (cand->GetIDDaughter(iDau) == iRP->GetID())
 913                {
 914                   if (fDebug) printf(" was in RP set");
 915                   iRP->SetForRPSelection(kFALSE);
 916                   fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
 917                }
 918             }
 919             if (fDebug) printf("\n");
 920          }
 921          cand->SetForPOISelection(kTRUE);
 922          fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
 923       }
 924       if (fDebug) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
 925
 926       for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
 927       {
 928          for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++)
 929          {
 930             if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
 931             if (fCandidateMinEta && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
 932             PtSelector(1, up[pTracks], up[nTracks]);
 933          }
 934       }
 935       for (Int_t nTracks = 0; nTracks < unn ; nTracks++)
 936       {
 937          for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++)
 938          {
 939             if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
 940             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
 941             PtSelector(2, un[nTracks], un[pTracks]);
 942          }
 943       }
 944       PostData(1, fOutputList);
 945       PostData(2, fFlowEvent);
 946    }
 947 }
 948 //_____________________________________________________________________________
 949 void AliAnalysisTaskPhiFlow::Terminate(Option_t *)
 950 {
 951    // Terminate
 952 }
 953 //______________________________________________________________________________
 954 void  AliAnalysisTaskPhiFlow::MakeTrack(Double_t mass,
 955                                               Double_t pt, Double_t phi, Double_t eta,
 956                                               Int_t nDau, Int_t iID[]) const
 957 {
 958    // Consruct Flow Candidate Track from two selected candidates
 959    Bool_t overwrite = kTRUE;
 960    AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
 961    if (!sTrack)
 962    {
 963       sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
 964       overwrite = kFALSE;
 965    }
 966    else sTrack->ClearMe();
 967    sTrack->SetMass(mass);
 968    sTrack->SetPt(pt);
 969    sTrack->SetPhi(phi);
 970    sTrack->SetEta(eta);
 971    for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
 972    sTrack->SetForPOISelection(kTRUE);
 973    sTrack->SetForRPSelection(kFALSE);
 974    if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
 975    else fCandidates->AddLast(sTrack);
 976    return;
 977 }
 978 //_____________________________________________________________________________
 979 void AliAnalysisTaskPhiFlow::SetCommonConstants(Int_t massBins, Double_t minMass, Double_t maxMass)
 980 {
 981   // set common constants
 982   fMassBins = massBins;
 983   fMinMass = minMass;
 984   fMaxMass = maxMass;
 985 }
 986 //_____________________________________________________________________________
 987