]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
reduce memory usage
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskPhiFlow.cxx
CommitLineData
19d77c98 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
48077518 19// handles aod's and esd's transparently
19d77c98 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"
b4f29a18 28#include "TObjArray.h"
19d77c98 29#include "AliAnalysisTask.h"
30#include "AliAnalysisTaskSE.h"
31#include "AliAnalysisManager.h"
32#include "AliESDEvent.h"
48077518 33#include "AliAODEvent.h"
19d77c98 34#include "AliESDInputHandler.h"
48077518 35#include "AliAODInputHandler.h"
19d77c98 36#include "AliCentrality.h"
37#include "AliVEvent.h"
38#include "AliAnalysisTaskPhiFlow.h"
48077518 39#include "AliFlowBayesianPID.h"
19d77c98 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"
48077518 52#include "AliAODVZERO.h"
7c29ebab 53#include "AliPID.h"
b4f29a18 54#include "AliPIDResponse.h"
19d77c98 55
56class AliFlowTrackCuts;
57
58ClassImp(AliAnalysisTaskPhiFlow)
59
60AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
b4f29a18 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)
19d77c98 62{
b4f29a18 63 // Default constructor
64 for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = 0.;
19d77c98 65}
66//_____________________________________________________________________________
67AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow(const char *name) : AliAnalysisTaskSE(name),
b4f29a18 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)
19d77c98 69{
70 // Constructor
b4f29a18 71 for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = 0.;
72
19d77c98 73 DefineInput(0, TChain::Class());
74 DefineOutput(1, TList::Class());
b4f29a18 75 DefineOutput(2, AliFlowEventSimple::Class());
19d77c98 76}
77//_____________________________________________________________________________
78AliAnalysisTaskPhiFlow::~AliAnalysisTaskPhiFlow()
79{
80 // Destructor
b4f29a18 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;
19d77c98 86}
87//_____________________________________________________________________________
88TH1F* 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//_____________________________________________________________________________
100TH2F* 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//_____________________________________________________________________________
48077518 110TH1F* 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);
b4f29a18 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;
48077518 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//_____________________________________________________________________________
19d77c98 123TH1F* 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//_____________________________________________________________________________
134void 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
19d77c98 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");
48077518 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
19d77c98 222 fPtP = BookPtHistogram("i^{+}");
223 fPtN = BookPtHistogram("i^{-}");
224 fPtKP = BookPtHistogram("K^{+}");
225 fPtKN = BookPtHistogram("K^{-}");
226
19d77c98 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);
19d77c98 244}
245//_____________________________________________________________________________
246void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
247{
248 // Create user defined output objects
b4f29a18 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!!!");
19d77c98 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
b4f29a18 276 cc->SetNbinsMass(fMassBins);
277 cc->SetMassMin(fMinMass);
278 cc->SetMassMax(fMaxMass);
279
48077518 280 // setup initial state of PID response object
281 if (!fOldTrackParam) fBayesianResponse->SetNewTrackParam();
282
b4f29a18 283 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
284 if (man)
285 {
286 AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
287 if (inputHandler) fPIDResponse = inputHandler->GetPIDResponse();
288 }
289
19d77c98 290 // Create all output objects and store them to a list
291 fOutputList = new TList();
b4f29a18 292 fOutputList->SetOwner(kTRUE);
19d77c98 293
294 // Create and post the Phi identification output objects
295 AddPhiIdentificationOutputObjects();
296 PostData(1, fOutputList);
297
b4f29a18 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);
19d77c98 305}
306//_____________________________________________________________________________
48077518 307template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* track1, const T* track2) const
19d77c98 308{
309 // Return the invariant mass of two tracks, assuming both tracks are kaons
310 if ((!track2) || (!track1)) return 0.;
19d77c98 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);
48077518 318 if ((es - (pxs + pys + pzs)) < 0) return 0.;
19d77c98 319 return TMath::Sqrt((es - (pxs + pys + pzs)));
320}
321//_____________________________________________________________________________
48077518 322template <typename T> Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const T* track1, const T* track2) const
19d77c98 323{
b4f29a18 324 // Calculate the delta dip angle between two particles
19d77c98 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//_____________________________________________________________________________
48077518 329template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const T* track1, const T* track2) const
19d77c98 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//_____________________________________________________________________________
48077518 336template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
337{
b4f29a18 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;
48077518 345}
346//_____________________________________________________________________________
19d77c98 347void 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//_____________________________________________________________________________
48077518 355template <typename T> Bool_t AliAnalysisTaskPhiFlow::EventCut(T* event)
19d77c98 356{
357 // Impose event cuts
48077518 358 if (!event) return kFALSE;
48077518 359 if (!CheckVertex(event)) return kFALSE;
360 if (!CheckCentrality(event)) return kFALSE;
361 PlotVZeroMultiplcities(event);
19d77c98 362 return kTRUE;
363}
364//_____________________________________________________________________________
48077518 365template <typename T> void AliAnalysisTaskPhiFlow::PlotVZeroMultiplcities(const T* event) const
19d77c98 366{
367 // QA multiplicity plots
48077518 368 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
369 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
19d77c98 370}
371//_____________________________________________________________________________
48077518 372template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckVertex(const T* event) const
19d77c98 373{
374 // Check if event vertex is within given range
48077518 375 if (!event->GetPrimaryVertex()) return 0x0;
376 if (TMath::Abs((event->GetPrimaryVertex())->GetZ()) > fVertexRange) return 0x0;
19d77c98 377 return kTRUE;
378}
379//_____________________________________________________________________________
48077518 380template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCentrality(T* event)
19d77c98 381{
382 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
48077518 383 if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
384 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
385 if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
19d77c98 386 {
48077518 387 fCentralityNoPass->Fill(fCentrality) ;
388 return kFALSE;
19d77c98 389 }
48077518 390 fCentralityPass->Fill(fCentrality);
391 return kTRUE;
392}
393//_____________________________________________________________________________
394void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliESDEvent* event)
395{
396 // Initialize the Bayesian PID object for ESD analysis
397 fBayesianResponse->SetDetResponse(event, fCentrality, AliESDpid::kTOF_T0, kTRUE);
b4f29a18 398 if (fDebug) cout << " --> Initialized Bayesian ESD PID object " << endl;
48077518 399}
400//_____________________________________________________________________________
401void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliAODEvent* event)
402{
403 // Initialize the Bayesian PID object for AOD
404 fBayesianResponse->SetDetResponse(event, fCentrality);
b4f29a18 405 if (fDebug) cout << " --> Initialized Bayesian AOD PID object " << endl;
48077518 406}
407//_____________________________________________________________________________
408template <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();
b4f29a18 415 if (probabilities[3] > fPIDConfig[6])
19d77c98 416 {
48077518 417 fPhi->Fill(track->Phi());
418 fPt->Fill(track->Pt());
419 fEta->Fill(track->Eta());
420 return kTRUE;
19d77c98 421 }
48077518 422 return kFALSE;
423}
424//_____________________________________________________________________________
425Bool_t AliAnalysisTaskPhiFlow::PassesStrictKaonCuts(AliESDtrack* track) const
426{
b4f29a18 427 // propagate dca from TPC
48077518 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//_____________________________________________________________________________
435Bool_t AliAnalysisTaskPhiFlow::PassesStrictKaonCuts(AliAODTrack* track) const
436{
b4f29a18 437 // propagate dca from TPC
48077518 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;
19d77c98 443}
444//_____________________________________________________________________________
48077518 445Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliESDtrack* track) const
19d77c98 446{
447 // Check if particle is a kaon according to method set in steering macro
48077518 448 if (fRequireTPCStandAlone && (track->GetStatus()&AliESDtrack::kTPCin) == 0) return kFALSE;
b4f29a18 449 if (fPIDConfig[1]!=0 || fPIDConfig[4]!=0) AliFatal("TPC || ITS PID not available in ESD anlaysis -- terminating analysis !!!");
450 if (PassesTPCbayesianCut(track))
19d77c98 451 {
b4f29a18 452 fPIDk->Fill(track->P(), track->GetTPCsignal());
453 return kTRUE;
48077518 454 }
455 return kFALSE;
456}
457//_____________________________________________________________________________
458Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
459{
b4f29a18 460 // Kaon identification routine, based on multiple detectors and approaches
48077518 461 if (fRequireTPCStandAlone && (!track->TestFilterBit(1))) return kFALSE;
b4f29a18 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;
19d77c98 507 }
508 return kFALSE;
509}
510//_____________________________________________________________________________
48077518 511template <typename T> Double_t AliAnalysisTaskPhiFlow::PhiPt(const T* track1, const T* track2) const
19d77c98 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//_____________________________________________________________________________
48077518 520template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, const T* track1, const T* track2) const
19d77c98 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 {
b4f29a18 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 }
19d77c98 616 }
617 if (tracktype == 1)
618 {
48077518 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));
19d77c98 637 }
638 if (tracktype == 2)
639 {
48077518 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//_____________________________________________________________________________
48077518 661template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
19d77c98 662{
663 // Check if track is suitable for phi flow analysis
664 if (!track) return kFALSE;
665 return fPOICuts->IsSelected(track);
666}
667//_____________________________________________________________________________
48077518 668template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
19d77c98 669{
670 // Set null cuts
b4f29a18 671 if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
48077518 672 fCutsRP->SetEvent(event, MCEvent());
19d77c98 673 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
674 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
675 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
48077518 676 fNullCuts->SetEvent(event, MCEvent());
19d77c98 677}
678//_____________________________________________________________________________
679void AliAnalysisTaskPhiFlow::PrepareFlowEvent(Int_t iMulti)
680{
681 // Prepare flow events
b4f29a18 682 fFlowEvent->ClearFast();
683 fFlowEvent->Fill(fCutsRP, fNullCuts);
684 fFlowEvent->SetReferenceMultiplicity(iMulti);
685 fFlowEvent->DefineDeadZone(0, 0, 0, 0);
19d77c98 686}
687//_____________________________________________________________________________
688void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
689{
b4f29a18 690 // UserExec: called for each event. Commented where necessary
48077518 691 // check for AOD data type
b4f29a18 692 if (!fPIDResponse)
693 {
694 AliError("Cannot get pid response");
695 return;
696 }
48077518 697
698 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
699 if (fAOD)
700 {
48077518 701 fAODAnalysis = kTRUE;
702 // Check whether event passes event cuts
703 if (!EventCut(fAOD)) return;
b4f29a18 704 // initialize event objects
48077518 705 InitializeBayesianPID(fAOD);
706 SetNullCuts(fAOD);
707 PrepareFlowEvent(fAOD->GetNumberOfTracks());
b4f29a18 708 fCandidates->SetLast(-1);
48077518 709 // Calculate event plane Q vectors and event plane resolution
48077518 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);
48077518 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;
b4f29a18 721 if (fStrictKaonCuts && (!PassesStrictKaonCuts(track))) continue;
48077518 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 }
19d77c98 751 }
48077518 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++)
19d77c98 754 {
48077518 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]);
48077518 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();
b4f29a18 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())
48077518 789 {
b4f29a18 790 if (fDebug) printf(" was in RP set");
791 iRP->SetForRPSelection(kFALSE);
792 fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
48077518 793 }
b4f29a18 794 }
795 if (fDebug) printf("\n");
48077518 796 }
b4f29a18 797 cand->SetForPOISelection(kTRUE);
798 fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
19d77c98 799 }
b4f29a18 800 if (fDebug) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
801
802
48077518 803 for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
19d77c98 804 {
48077518 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 }
19d77c98 811 }
48077518 812 for (Int_t nTracks = 0; nTracks < unn ; nTracks++)
19d77c98 813 {
48077518 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 }
48077518 821 PostData(1, fOutputList);
b4f29a18 822 PostData(2, fFlowEvent);
48077518 823 }
824
825 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
826 if (fESD)
827 {
48077518 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
48077518 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);
48077518 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)
19d77c98 848 {
48077518 849 charge = kTRUE;
850 fEventStats->Fill(1);
851 fPtP->Fill(track->Pt());
19d77c98 852 }
48077518 853 if (track->Charge() < 0)
19d77c98 854 {
48077518 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 }
19d77c98 874 }
875 }
48077518 876 // Calculate invariant mass of like- and unlike sign pairs as a function of p_T, store data in histograms
b4f29a18 877 for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
19d77c98 878 {
48077518 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]);
48077518 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();
b4f29a18 891 Int_t nIDs[2];
48077518 892 nIDs[0] = up[pTracks]->GetID();
893 nIDs[1] = un[nTracks]->GetID();
b4f29a18 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())
48077518 913 {
b4f29a18 914 if (fDebug) printf(" was in RP set");
915 iRP->SetForRPSelection(kFALSE);
916 fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
48077518 917 }
b4f29a18 918 }
919 if (fDebug) printf("\n");
48077518 920 }
b4f29a18 921 cand->SetForPOISelection(kTRUE);
922 fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
19d77c98 923 }
b4f29a18 924 if (fDebug) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
925
48077518 926 for (Int_t pTracks = 0; pTracks < unp ; pTracks++)
19d77c98 927 {
48077518 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 }
19d77c98 934 }
48077518 935 for (Int_t nTracks = 0; nTracks < unn ; nTracks++)
19d77c98 936 {
48077518 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 }
19d77c98 943 }
48077518 944 PostData(1, fOutputList);
b4f29a18 945 PostData(2, fFlowEvent);
48077518 946 }
19d77c98 947}
948//_____________________________________________________________________________
949void AliAnalysisTaskPhiFlow::Terminate(Option_t *)
950{
951 // Terminate
952}
953//______________________________________________________________________________
b4f29a18 954void AliAnalysisTaskPhiFlow::MakeTrack(Double_t mass,
955 Double_t pt, Double_t phi, Double_t eta,
956 Int_t nDau, Int_t iID[]) const
19d77c98 957{
958 // Consruct Flow Candidate Track from two selected candidates
b4f29a18 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();
19d77c98 967 sTrack->SetMass(mass);
968 sTrack->SetPt(pt);
969 sTrack->SetPhi(phi);
970 sTrack->SetEta(eta);
b4f29a18 971 for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
19d77c98 972 sTrack->SetForPOISelection(kTRUE);
973 sTrack->SetForRPSelection(kFALSE);
b4f29a18 974 if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
975 else fCandidates->AddLast(sTrack);
976 return;
977}
978//_____________________________________________________________________________
979void 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;
19d77c98 985}
986//_____________________________________________________________________________
987