]> 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
Completed changes needed because of previous commit
[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:
5725c281 17// author: Redmer Alexander Bertens (rbertens@cern.ch)
18// analyis task for phi-meson reconstruction and determination of v_n
19// AliPhiMesonHelperTrack provides a lightweight helper track for reconstruction
20// new in this version (use with caution): vzero event plane, event mixing
19d77c98 21
22#include "TChain.h"
19d77c98 23#include "TH1.h"
24#include "TH1F.h"
25#include "TH2F.h"
19d77c98 26#include "TMath.h"
b4f29a18 27#include "TObjArray.h"
19d77c98 28#include "AliAnalysisTaskSE.h"
29#include "AliAnalysisManager.h"
48077518 30#include "AliAODEvent.h"
48077518 31#include "AliAODInputHandler.h"
19d77c98 32#include "AliCentrality.h"
19d77c98 33#include "AliAnalysisTaskPhiFlow.h"
48077518 34#include "AliFlowBayesianPID.h"
9fb90465 35#include "AliPIDCombined.h"
19d77c98 36#include "AliMCEvent.h"
37#include "TProfile.h"
38#include "AliFlowCandidateTrack.h"
39#include "AliFlowTrackCuts.h"
40#include "AliFlowEventSimple.h"
5725c281 41#include "AliFlowTrackSimple.h"
19d77c98 42#include "AliFlowCommonConstants.h"
43#include "AliFlowEvent.h"
44#include "TVector3.h"
48077518 45#include "AliAODVZERO.h"
b4f29a18 46#include "AliPIDResponse.h"
863bc87f 47#include "AliAODMCParticle.h"
5725c281 48#include "AliAnalysisTaskVnV0.h"
49#include "AliEventPoolManager.h"
19d77c98 50
51class AliFlowTrackCuts;
52
3a7af7bd 53using std::cout;
54using std::endl;
5725c281 55
19d77c98 56ClassImp(AliAnalysisTaskPhiFlow)
5725c281 57ClassImp(AliPhiMesonHelperTrack)
19d77c98 58
59AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
0d91ada7 60 fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCandidateYCut(kFALSE), fCandidateMinY(-.5), fCandidateMaxY(.5), fNPtBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0),fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut2010(0), fCentralityCut2011(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0)
19d77c98 61{
b4f29a18 62 // Default constructor
5725c281 63 for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
64 for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
65 for(Int_t i(0); i < 20; i++) {
66 fVertexMixingBins[i] = 0;
67 fCentralityMixingBins[i] = 0;
68 }
69 fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
70 for(Int_t i(0); i < 18; i++) {
ffef9156 71 for(Int_t j(0); j < 2; j++) fV0Data[i][j] = 0;
72 fInvMNP[i] = 0; fInvMNN[i] = 0; fInvMPP[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.;
5725c281 73 }
19d77c98 74}
75//_____________________________________________________________________________
76AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow(const char *name) : AliAnalysisTaskSE(name),
0d91ada7 77 fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCandidateYCut(kFALSE), fCandidateMinY(-.5), fCandidateMaxY(.5), fNPtBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut2010(0), fCentralityCut2011(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0)
19d77c98 78{
79 // Constructor
5725c281 80 for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
81 for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
82 for(Int_t i(0); i < 20; i++) {
83 fVertexMixingBins[i] = 0;
84 fCentralityMixingBins[i] = 0;
85 }
86 fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
87 for(Int_t i(0); i < 18; i++) {
ffef9156 88 for(Int_t j(0); j < 2; j++) fV0Data[i][j] = 0;
89 fInvMNP[i] = 0; fInvMNN[i] = 0; fInvMPP[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.;
5725c281 90 }
19d77c98 91 DefineInput(0, TChain::Class());
92 DefineOutput(1, TList::Class());
b4f29a18 93 DefineOutput(2, AliFlowEventSimple::Class());
5725c281 94 if(fDebug > 0) cout << " === Created instance of AliAnaysisTaskPhiFlow === " << endl;
19d77c98 95}
96//_____________________________________________________________________________
97AliAnalysisTaskPhiFlow::~AliAnalysisTaskPhiFlow()
98{
99 // Destructor
b4f29a18 100 if (fNullCuts) delete fNullCuts;
101 if (fOutputList) delete fOutputList;
102 if (fCandidates) delete fCandidates;
103 if (fFlowEvent) delete fFlowEvent;
104 if (fBayesianResponse) delete fBayesianResponse;
5725c281 105 if (fEventMixing) delete fPoolManager;
9fb90465 106 if (fPIDCombined) delete fPIDCombined;
5725c281 107 if (fDebug > 0) cout << " === Deleted instance of AliAnalysisTaskPhiFlow === " << endl;
19d77c98 108}
109//_____________________________________________________________________________
110TH1F* AliAnalysisTaskPhiFlow::BookHistogram(const char* name)
111{
112 // Return a pointer to a TH1 with predefined binning
ffef9156 113 if(fDebug > 0) cout << " *** BookHistogram() *** " << name << endl;
19d77c98 114 TH1F *hist = new TH1F(name, Form("M_{INV} (%s)", name), 60, .99, 1.092);
115 hist->GetXaxis()->SetTitle("M_{INV} (GeV / c^{2})");
116 hist->GetYaxis()->SetTitle("No. of pairs");
117 hist->SetMarkerStyle(kFullCircle);
118 hist->Sumw2();
119 fOutputList->Add(hist);
120 return hist;
121}
122//_____________________________________________________________________________
863bc87f 123TH2F* AliAnalysisTaskPhiFlow::BookPIDHistogram(const char* name, Bool_t TPC)
19d77c98 124{
125 // Return a pointer to a TH2 with predefined binning
5725c281 126 if(fDebug > 0) cout << " *** BookPIDHisotgram() *** " << endl;
863bc87f 127 TH2F *hist = 0x0;
128 if(TPC) {
129 hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
130 hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
131 }
132 if(!TPC) {
133 hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1.5);
134 hist->GetYaxis()->SetTitle("#beta");
135 }
19d77c98 136 hist->GetXaxis()->SetTitle("P (GeV / c)");
19d77c98 137 fOutputList->Add(hist);
138 return hist;
139}
140//_____________________________________________________________________________
ffef9156 141TH1F* AliAnalysisTaskPhiFlow::InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
48077518 142{
143 // intialize p_t histograms for each p_t bin
5725c281 144 if(fDebug > 0) cout << " *** InitPtSpectraHistograms() *** " << endl;
ffef9156 145 TH1F* hist = new TH1F(Form("%4.2f p_{t} %4.2f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 60, nmin, nmax);
48077518 146 hist->GetXaxis()->SetTitle("p_{T} GeV / c");
147 fOutputList->Add(hist);
148 return hist;
149}
150//_____________________________________________________________________________
19d77c98 151TH1F* AliAnalysisTaskPhiFlow::BookPtHistogram(const char* name)
152{
153 // Return a pointer to a p_T spectrum histogram
5725c281 154 if(fDebug > 0) cout << " *** BookPtHistogram() *** " << endl;
19d77c98 155 TH1F* ratio = new TH1F(name, name, 100, 0, 7);
156 ratio->GetXaxis()->SetTitle("p_{T} ( GeV / c^{2} )");
157 ratio->GetYaxis()->SetTitle("No. of events");
158 ratio->Sumw2();
159 fOutputList->Add(ratio);
160 return ratio;
161}
162//_____________________________________________________________________________
163void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
164{
165 // Add Phi Identification Output Objects
5725c281 166 if(fDebug > 0) cout << " ** AddPhiIdentificationOutputObjects() *** " << endl;
863bc87f 167 if(fQA) {
168 fEventStats = new TH1F("fHistStats", "Event Statistics", 18, -.25, 4.25);
169 fEventStats->GetXaxis()->SetTitle("No. of events");
170 fEventStats->GetYaxis()->SetTitle("Statistics");
171 fOutputList->Add(fEventStats);
172 }
19d77c98 173 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
19d77c98 174 fOutputList->Add(fCentralityPass);
863bc87f 175 if(fQA) {
176 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
177 fOutputList->Add(fCentralityNoPass);
178 fNOPID = BookPIDHistogram("TPC signal, all particles", kTRUE);
179 fPIDk = BookPIDHistogram("TPC signal, kaons", kTRUE);
180 fNOPIDTOF = BookPIDHistogram("TOF signal, all particles", kFALSE);
181 fPIDTOF = BookPIDHistogram("TOF signal, kaons", kFALSE);
182 }
ffef9156 183 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
184 fInvMNP[ptbin] = BookHistogram(Form("NP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
185 fInvMNN[ptbin] = BookHistogram(Form("NN, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
186 fInvMPP[ptbin] = BookHistogram(Form("PP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
5725c281 187 }
188 if(fQA) {
ffef9156 189 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) fPtSpectra[ptbin] = InitPtSpectraHistograms(fPtBins[ptbin], fPtBins[ptbin+1]);
863bc87f 190 fPtP = BookPtHistogram("i^{+}");
191 fPtN = BookPtHistogram("i^{-}");
192 fPtKP = BookPtHistogram("K^{+}");
193 fPtKN = BookPtHistogram("K^{-}");
194 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
195 fOutputList->Add(fPhi);
196 fPt = new TH1F("fPt", "p_{T}", 100, 0, 5.5);
197 fOutputList->Add(fPt);
198 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
199 fOutputList->Add(fEta);
200 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
201 fOutputList->Add(fVZEROA);
202 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
203 fOutputList->Add(fVZEROC);
204 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
205 fOutputList->Add(fTPCM);
206 fDCAXYQA = new TH1F("fDCAXYQA", "fDCAXYQA", 1000, -5, 5);
207 fOutputList->Add(fDCAXYQA);
208 fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
209 fOutputList->Add(fDCAZQA);
6940e1f0 210 if(fCentralityCut2010 || fCentralityCut2011) {
211 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
212 fOutputList->Add(fMultCorAfterCuts);
213 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
214 fOutputList->Add(fMultvsCentr);
215 }
863bc87f 216 }
217 if(fIsMC || fQA) {
863bc87f 218 fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
219 fOutputList->Add(fDCAAll);
863bc87f 220 fDCAPrim = new TH2F("fDCAprim","fDCAprim", 1000, 0, 10, 1000, -5, 5);
221 fOutputList->Add(fDCAPrim);
863bc87f 222 fDCASecondaryWeak = new TH2F("fDCASecondaryWeak","fDCASecondaryWeak", 1000, 0, 10, 1000, -5, 5);
223 fOutputList->Add(fDCASecondaryWeak);
863bc87f 224 fDCAMaterial = new TH2F("fDCAMaterial","fDCAMaterial", 1000, 0, 10, 1000, -5, 5);
225 fOutputList->Add(fDCAMaterial);
226 }
5725c281 227 if(fV0) {
ffef9156 228 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 5, 0, 5);
229 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<#Psi_{a} - #Psi_{b}>");
230 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<#Psi_{a} - #Psi_{c}>");
231 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<#Psi_{b} - #Psi_{c}>");
232 fSubEventDPhiv2->GetXaxis()->SetBinLabel(4, "#sqrt{#frac{<#Psi_{a} - #Psi_{b}><#Psi_{a} - #Psi_{c}>}{<#Psi_{b} - #Psi_{c}>}}");
233 fSubEventDPhiv2->GetXaxis()->SetBinLabel(5, "#sqrt{#frac{<#Psi_{a} - #Psi_{c}><#Psi_{b} - #Psi_{c}>}{<#Psi_{a} - #Psi_{b}>}}");
5725c281 234 fOutputList->Add(fSubEventDPhiv2);
5725c281 235 const char* V0[] = {"V0A", "V0C"};
ffef9156 236 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++)
237 for(Int_t iV0(0); iV0 < 2; iV0++) {
238 fV0Data[ptbin][iV0] = new TProfile(Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), fMassBins, fMinMass, fMaxMass);
239 fOutputList->Add(fV0Data[ptbin][iV0]);
240 }
5725c281 241 }
19d77c98 242}
243//_____________________________________________________________________________
244void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
245{
246 // Create user defined output objects
5725c281 247 if(fDebug > 0) cout << " *** UserCreateOutputObjects() *** " << endl;
b4f29a18 248 fNullCuts = new AliFlowTrackCuts("null_cuts");
249 fBayesianResponse = new AliFlowBayesianPID();
9fb90465 250 // combined pid
251 fPIDCombined = new AliPIDCombined;
252 fPIDCombined->SetDefaultTPCPriors();
253 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
254
9e58d58d 255 // flag to mc
6940e1f0 256 if(fSkipEventSelection || fIsMC) fBayesianResponse->SetMC(kTRUE);
b4f29a18 257 Double_t t(0);
258 for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
ffef9156 259 if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis !!!");
260 if(fNPtBins > 18) AliFatal("Invalid number of pt bins initialied ( > 18 ) -- terminating analysis !!!");
19d77c98 261 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
ffef9156 262 cc->SetNbinsQ(500); cc->SetNbinsPhi(180); cc->SetNbinsMult(10000);
263 cc->SetQMin(0.0); cc->SetPhiMin(0.0); cc->SetMultMin(0);
264 cc->SetQMax(3.0); cc->SetPhiMax(TMath::TwoPi()); cc->SetMultMax(10000);
265 cc->SetNbinsMass(fMassBins); cc->SetNbinsEta(200); (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
266 cc->SetMassMin(fMinMass); cc->SetEtaMin(-5.0); cc->SetPtMin(0);
267 cc->SetMassMax(fMaxMass); cc->SetEtaMax(+5.0); (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
9fb90465 268 fBayesianResponse->SetNewTrackParam();
b4f29a18 269 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
5725c281 270 if (man) {
b4f29a18 271 AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
272 if (inputHandler) fPIDResponse = inputHandler->GetPIDResponse();
273 }
19d77c98 274 // Create all output objects and store them to a list
275 fOutputList = new TList();
b4f29a18 276 fOutputList->SetOwner(kTRUE);
19d77c98 277 // Create and post the Phi identification output objects
278 AddPhiIdentificationOutputObjects();
279 PostData(1, fOutputList);
b4f29a18 280 // create candidate array
281 fCandidates = new TObjArray(1000);
282 fCandidates->SetOwner(kTRUE);
b4f29a18 283 // create and post flowevent
284 fFlowEvent = new AliFlowEvent(10000);
285 PostData(2, fFlowEvent);
5725c281 286 if(fEventMixing) fPoolManager = InitializeEventMixing();
287}
288//_____________________________________________________________________________
289AliEventPoolManager* AliAnalysisTaskPhiFlow::InitializeEventMixing()
290{
291 // initialize event mixing
292 if(fDebug > 0) cout << " *** InitializeEventMixing() *** " << endl;
293 Int_t _c(0), _v(0);
294 for(Int_t i(0); i < 19; i++) {
295 if (fCentralityMixingBins[i+1] < fCentralityMixingBins[i]) { _c = i; break; }
296 else _c = 19;
297 }
298 for(Int_t i(0); i < 19; i++) {
299 if (fVertexMixingBins[i+1] < fVertexMixingBins[i]) { _v = i; break; }
300 else _v = 19;
301 }
302 if(fDebug > 0 ) cout << Form(" --> found %d centrality bins for mixing, %d vertex bins for mixing", _c, _v) << endl;
303 Double_t centralityBins[_c];
304 Double_t vertexBins[_v];
305 for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
306 for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
307 return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
19d77c98 308}
309//_____________________________________________________________________________
48077518 310template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* track1, const T* track2) const
19d77c98 311{
312 // Return the invariant mass of two tracks, assuming both tracks are kaons
5725c281 313 if(fDebug > 1) cout << " *** InvariantMass() *** " << endl;
19d77c98 314 if ((!track2) || (!track1)) return 0.;
19d77c98 315 Double_t masss = TMath::Power(4.93676999999999977e-01, 2);
316 Double_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
317 Double_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
318 Double_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
319 Double_t e1 = TMath::Sqrt(track1->P() * track1->P() + masss);
320 Double_t e2 = TMath::Sqrt(track2->P() * track2->P() + masss);
321 Double_t es = TMath::Power((e1 + e2), 2);
48077518 322 if ((es - (pxs + pys + pzs)) < 0) return 0.;
19d77c98 323 return TMath::Sqrt((es - (pxs + pys + pzs)));
324}
325//_____________________________________________________________________________
9fb90465 326/*
48077518 327template <typename T> Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const T* track1, const T* track2) const
19d77c98 328{
b4f29a18 329 // Calculate the delta dip angle between two particles
5725c281 330 if(fDebug > 1) cout << " *** DeltaDipAngle() *** " << endl;
19d77c98 331 if (track1->P()*track2->P() == 0) return 999;
332 return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
333}
334//_____________________________________________________________________________
48077518 335template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const T* track1, const T* track2) const
19d77c98 336{
337 // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
5725c281 338 if(fDebug > 1) cout << " *** CheckDeltaDipAngle() *** " << endl;
19d77c98 339 if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PhiPt(track1, track2) < fDeltaDipPt)) return kFALSE;
340 return kTRUE;
341}
9fb90465 342*/
19d77c98 343//_____________________________________________________________________________
48077518 344template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
345{
b4f29a18 346 // Check if pair passes eta and pt cut
5725c281 347 if(fDebug > 1) cout << " *** CheckCandidateEtaPtCut() *** " << endl;
b4f29a18 348 if (fCandidateMinPt > PhiPt(track1, track2) || fCandidateMaxPt < PhiPt(track1, track2)) return kFALSE;
349 TVector3 a(track1->Px(), track1->Py(), track1->Pz());
350 TVector3 b(track2->Px(), track2->Py(), track2->Pz());
351 TVector3 c = a + b;
352 if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
353 return kTRUE;
48077518 354}
355//_____________________________________________________________________________
6940e1f0 356template <typename T> Bool_t AliAnalysisTaskPhiFlow::EventCut(T* event)
357{
358 // Impose event cuts
359 if(fDebug > 0) cout << " *** EventCut() *** " << endl;
360 if (!event) return kFALSE;
361 if (fSkipEventSelection) return kTRUE;
362 if (!CheckVertex(event)) return kFALSE;
363 if (!CheckCentrality(event)) return kFALSE;
364 if(fQA) PlotMultiplcities(event);
365 return kTRUE;
366}
367//_____________________________________________________________________________
863bc87f 368template <typename T> void AliAnalysisTaskPhiFlow::PlotMultiplcities(const T* event) const
19d77c98 369{
370 // QA multiplicity plots
5725c281 371 if(fDebug > 1) cout << " *** PlotMultiplcities() *** " << endl;
48077518 372 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
373 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
863bc87f 374 fTPCM->Fill(event->GetNumberOfTracks());
19d77c98 375}
376//_____________________________________________________________________________
6940e1f0 377template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckVertex(const T* event)
378{
379 // Check if event vertex is within given range
380 if(fDebug > 0) cout << " *** CheckVertex() *** " << endl;
381 if (!event->GetPrimaryVertex()) return 0x0;
382 fVertex = event->GetPrimaryVertex()->GetZ();
383 if (TMath::Abs(fVertex) > fVertexRange) return 0x0;
384 return kTRUE;
385}
386//_____________________________________________________________________________
387template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCentrality(T* event)
388{
389 // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
390 if(fDebug > 0) cout << " *** CheckCentrality() *** " << endl;
391 if (!fkCentralityMethodA) AliFatal("No centrality method set! FATAL ERROR!");
392 fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodA);
393 Double_t cenB(-999);
394 // if a second centrality estimator is requited, set it
395 (fkCentralityMethodB) ? cenB = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodB) : cenB = fCentrality;
396 if (TMath::Abs(fCentrality-cenB) > 5 || cenB >= 80 || cenB < 0 || fCentrality <= fCentralityMin || fCentrality > fCentralityMax) {
397 if(fQA) fCentralityNoPass->Fill(fCentrality) ;
398 return kFALSE;
399 }
400 const Int_t nGoodTracks = event->GetNumberOfTracks();
401 if(fCentralityCut2010) { // cut on outliers
402 Float_t multTPC(0.); // tpc mult estimate
403 Float_t multGlob(0.); // global multiplicity
404 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
f15c1f69 405 AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
406 if(!trackAOD) AliFatal("Not a standard AOD");
6940e1f0 407 if (!trackAOD) continue;
408 if (!(trackAOD->TestFilterBit(1))) continue;
409 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
410 multTPC++;
411 }
412 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
f15c1f69 413 AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
414 if(!trackAOD) AliFatal("Not a standard AOD");
6940e1f0 415 if (!trackAOD) continue;
416 if (!(trackAOD->TestFilterBit(16))) continue;
417 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
418 Double_t b[2] = {-99., -99.};
419 Double_t bCov[3] = {-99., -99., -99.};
420 AliAODTrack copy(*trackAOD);
421 if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
422 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
423 multGlob++;
424 } //track loop
425 // printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
426 if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
427 if(fQA) {
428 fMultCorAfterCuts->Fill(multGlob, multTPC);
429 fMultvsCentr->Fill(fCentrality, multTPC);
430 }
431 }
432
433 if(fCentralityCut2011) { // cut on outliers
434 Float_t multTPC(0.); // tpc mult estimate
435 Float_t multGlob(0.); // global multiplicity
436 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
f15c1f69 437 AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
438 if(!trackAOD) AliFatal("Not a standard AOD");
6940e1f0 439 if (!trackAOD) continue;
440 if (!(trackAOD->TestFilterBit(1))) continue;
441 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
442 multTPC++;
443 }
444 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
f15c1f69 445 AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
446 if(!trackAOD) AliFatal("Not a standard AOD");
6940e1f0 447 if (!trackAOD) continue;
448 if (!(trackAOD->TestFilterBit(16))) continue;
449 if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
450 Double_t b[2] = {-99., -99.};
451 Double_t bCov[3] = {-99., -99., -99.};
452 AliAODTrack copy(*trackAOD);
453 if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
454 if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
455 multGlob++;
456 } //track loop
457 //printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
458 if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))) return kFALSE;
459 if(fQA) {
460 fMultCorAfterCuts->Fill(multGlob, multTPC);
461 fMultvsCentr->Fill(fCentrality, multTPC);
462 }
463 }
464
465 fCentralityPass->Fill(fCentrality);
466 return kTRUE;
467}
468//_____________________________________________________________________________
48077518 469void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliAODEvent* event)
470{
471 // Initialize the Bayesian PID object for AOD
5725c281 472 if(fDebug > 0) cout << " *** InitializeBayesianPID() *** " << endl;
48077518 473 fBayesianResponse->SetDetResponse(event, fCentrality);
474}
475//_____________________________________________________________________________
476template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* track) const
477{
478 // Check if the particle passes the TPC TOF bayesian cut.
5725c281 479 if(fDebug > 1) cout << " *** PassesTPCbayesianCut() *** " << endl;
48077518 480 fBayesianResponse->ComputeProb(track);
481 if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
482 Float_t *probabilities = fBayesianResponse->GetProb();
5725c281 483 if (probabilities[3] > fPIDConfig[6]) {
863bc87f 484 if(fQA) {fPhi->Fill(track->Phi()); fPt->Fill(track->Pt()); fEta->Fill(track->Eta());}
48077518 485 return kTRUE;
19d77c98 486 }
48077518 487 return kFALSE;
488}
489//_____________________________________________________________________________
863bc87f 490Bool_t AliAnalysisTaskPhiFlow::PassesDCACut(AliAODTrack* track) const
491{
492 // check if track passes dca cut according to dca routine
493 // setup the routine as follows:
494 // fDCAConfig[0] < -1 no pt dependence
495 // fDCAConfig[0] = 0 do nothing
496 // fDCAConfig[0] > 1 pt dependent dca cut
5725c281 497 if(fDebug > 1) cout << " *** PassesDCACut() *** " << endl;
863bc87f 498 if(fIsMC) return kTRUE;
646a809c 499 AliAODTrack copy(*track);
ffef9156 500 if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
501 if(fQA) {
502 Double_t b[2] = { -99., -99.};
503 Double_t bCov[3] = { -99., -99., -99.};
646a809c 504 if(copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) {
9ea73210 505 fDCAXYQA->Fill(b[0]);
506 fDCAZQA->Fill(b[1]);
507 fDCAPrim->Fill(track->Pt(), b[0]);
508 }
ffef9156 509 }
510 return kTRUE;
511 }
863bc87f 512 Double_t b[2] = { -99., -99.};
513 Double_t bCov[3] = { -99., -99., -99.};
646a809c 514 if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
863bc87f 515 if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
516 if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
517 if(fDCAConfig[0] > .9) {
518 if(fDCAConfig[4] < TMath::Abs(b[1])) return kFALSE;
519 Double_t denom = TMath::Power(track->Pt(), fDCAConfig[3]);
5725c281 520 if( denom < 0.0000001 ) return kFALSE; // avoid division by zero
863bc87f 521 if( (fDCAConfig[1] + fDCAConfig[2] / denom) < TMath::Abs(b[0]) ) return kFALSE;
522 }
523 if(fQA) {
524 fDCAXYQA->Fill(b[0]);
525 fDCAZQA->Fill(b[1]);
526 fDCAPrim->Fill(track->Pt(), b[0]);
527 }
528 return kTRUE;
529}
530//_____________________________________________________________________________
48077518 531Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
532{
b4f29a18 533 // Kaon identification routine, based on multiple detectors and approaches
5725c281 534 if(fDebug > 1) cout << " *** IsKaon() *** " << endl;
9fb90465 535 if(fUsePidResponse) {
536 Double_t prob[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
537 fPIDCombined->ComputeProbabilities(track, fPIDResponse, prob);
538 if(prob[3] > fPIDConfig[6]) return kTRUE;
539 }
863bc87f 540 if(!PassesDCACut(track)) return kFALSE;
541 if(fQA) {fNOPID->Fill(track->P(), track->GetTPCsignal());fNOPIDTOF->Fill(track->P(), track->GetTOFsignal());}
5725c281 542 if(track->Pt() < fPIDConfig[1]) {
543 if(fDebug>1) cout << " ITS received track with p_t " << track->Pt() << endl;
b4f29a18 544 // if tpc control is disabled, pure its pid
5725c281 545 if(fPIDConfig[2] < 0.) {
b4f29a18 546 if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) return kTRUE;
547 return kFALSE;
548 }
549 // else, switch to ITS pid with TPC rejection of protons and pions
550 if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) < fPIDConfig[3]) return kFALSE;
551 else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) < fPIDConfig[3]) return kFALSE;
5725c281 552 else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) {
863bc87f 553 if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
b4f29a18 554 return kTRUE;
555 }
556 return kFALSE;
557 }
5725c281 558 if((track->Pt() > fPIDConfig[1]) && (track->Pt() < fPIDConfig[4])) {
559 if(fDebug>1) cout << " TPC received track with p_t " << track->Pt() << endl;
b4f29a18 560 // if its control is disabled, pure tpc pid
5725c281 561 if(fPIDConfig[5] < 0.) {
b4f29a18 562 if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) return kTRUE;
563 return kFALSE;
564 }
565 // else, switch to TPC pid with ITS rejection of protons and pions
566 if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kProton)) < fPIDConfig[0]) return kFALSE;
567 else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kPion)) < fPIDConfig[0]) return kFALSE;
5725c281 568 else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) {
863bc87f 569 if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
b4f29a18 570 return kTRUE;
571 }
572 return kFALSE;
573 }
5725c281 574 if(fDebug>1) cout << " Bayesian method received track with p_t " << track->Pt() << endl;
b4f29a18 575 // switch to bayesian PID
5725c281 576 if (PassesTPCbayesianCut(track)) {
863bc87f 577 if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal());fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
b4f29a18 578 return kTRUE;
19d77c98 579 }
580 return kFALSE;
581}
582//_____________________________________________________________________________
48077518 583template <typename T> Double_t AliAnalysisTaskPhiFlow::PhiPt(const T* track1, const T* track2) const
19d77c98 584{
5725c281 585 // return p_t of track pair
19d77c98 586 TVector3 a(track1->Px(), track1->Py(), track1->Pz());
587 TVector3 b(track2->Px(), track2->Py(), track2->Pz());
588 TVector3 c = a + b;
589 return c.Pt();
590}
591//_____________________________________________________________________________
48077518 592template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, const T* track1, const T* track2) const
19d77c98 593{
5725c281 594 // plot m_inv spectra of like- and unlike-sign kaon pairs for each pt bin
19d77c98 595 Double_t pt = PhiPt(track1, track2);
5725c281 596 if (tracktype == 0) {
ffef9156 597 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
598 if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
5725c281 599 fInvMNP[ptbin]->Fill(InvariantMass(track1, track2));
600 if(fQA) fPtSpectra[ptbin]->Fill(pt);
601 }
602 }
19d77c98 603 }
5725c281 604 if (tracktype == 1) {
ffef9156 605 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
606 if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
5725c281 607 fInvMPP[ptbin]->Fill(InvariantMass(track1, track2));
608 }
609 }
19d77c98 610 }
5725c281 611 if (tracktype == 2) {
ffef9156 612 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
613 if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
5725c281 614 fInvMNN[ptbin]->Fill(InvariantMass(track1, track2));
615 }
616 }
48077518 617 }
618}
619//_____________________________________________________________________________
6940e1f0 620template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
621{
622 // check if track meets quality cuts
623 if(!track) return kFALSE;
624 return fPOICuts->IsSelected(track);
625}
626//_____________________________________________________________________________
48077518 627template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
19d77c98 628{
629 // Set null cuts
5725c281 630 if (fDebug > 0) cout << " *** SetNullCuts() *** " << fCutsRP << endl;
48077518 631 fCutsRP->SetEvent(event, MCEvent());
19d77c98 632 fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
633 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
634 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
48077518 635 fNullCuts->SetEvent(event, MCEvent());
19d77c98 636}
637//_____________________________________________________________________________
638void AliAnalysisTaskPhiFlow::PrepareFlowEvent(Int_t iMulti)
639{
640 // Prepare flow events
5725c281 641 if (fDebug > 0 ) cout << " *** PrepareFlowEvent() *** " << endl;
b4f29a18 642 fFlowEvent->ClearFast();
643 fFlowEvent->Fill(fCutsRP, fNullCuts);
644 fFlowEvent->SetReferenceMultiplicity(iMulti);
645 fFlowEvent->DefineDeadZone(0, 0, 0, 0);
19d77c98 646}
647//_____________________________________________________________________________
5725c281 648void AliAnalysisTaskPhiFlow::VZEROSubEventAnalysis()
649{
650 // vzero event plane analysis using three subevents
651 if(fDebug > 0) cout << " ** VZEROSubEventAnalysis() *** " << endl;
652 if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
653 if(fDebug > 0 ) cout << "Fatal error: unable to retrieve VZERO task output !!! Exiting ..." << endl;
654 return;
655 }
656 // retrieve data
657 Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
ffef9156 658 for(Int_t i(0); i < 3; i++) if(abcPsi2[i] == 0) {
5725c281 659 if(fDebug > 0) cout << " Warning: I've encountered 0 in one of the symmetry panes (TPC, VZERO) - skipping event !" << endl;
660 return;
661 }
662 // save info for resolution calculations
ffef9156 663 Float_t qaqb = TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1])); // vzeroa - tpc
664 Float_t qaqc = TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2])); // vzeroa - vzeroc
665 Float_t qbqc = TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2])); // tpc - vzeroc
ffef9156 666 fSubEventDPhiv2->Fill(0.5, qaqb);
667 fSubEventDPhiv2->Fill(1.5, qaqc);
668 fSubEventDPhiv2->Fill(2.5, qbqc);
5725c281 669 Float_t minv[31];
ffef9156 670 Float_t _dphi[30][fNPtBins][2]; //minv, pt, v0a-c
671 Int_t occurence[30][fNPtBins]; //minv, pt
5725c281 672 for(Int_t mb(0); mb < 31; mb++) minv[mb] = 0.99 + mb * 0.0034;
ffef9156 673 for(Int_t i(0); i < 30; i++) for (Int_t j(0); j < fNPtBins; j++) for(Int_t k(0); k < 2; k ++) {
674 _dphi[i][j][k] = 0;
675 if(k==0) occurence[i][j] = 0;
5725c281 676 }
677 for(Int_t iCand(0); iCand < fCandidates->GetEntriesFast(); iCand++) { // loop over unlike sign kaon pairs
678 AliFlowTrackSimple *track = dynamic_cast<AliFlowTrackSimple*>(fCandidates->At(iCand));
679 if(!track) {
680 if(fDebug > 1) cout << " --> dynamic_cast returned null-pointer, skipping candidate <-- " << endl;
681 continue;
682 }
683 for(Int_t mb(0); mb < 30; mb++) { // loop over massbands
ffef9156 684 if(track->Mass() <= minv[mb] || track->Mass() > minv[mb+1]) continue;
685 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
686 if(track->Pt() <= fPtBins[ptbin] || track->Pt() > fPtBins[ptbin+1]) continue;
687 _dphi[mb][ptbin][0]+=TMath::Cos(2.*(track->Phi() - abcPsi2[0]));
688 _dphi[mb][ptbin][1]+=TMath::Cos(2.*(track->Phi() - abcPsi2[2]));
5725c281 689 occurence[mb][ptbin]+=1;
690 }
691 }
692 for(Int_t mb(0); mb < 30; mb++) // write vn values to tprofiles
ffef9156 693 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
5725c281 694 if(occurence[mb][ptbin]==0) continue;
ffef9156 695 fV0Data[ptbin][0]->Fill(mb*0.0034+0.99+0.0017, _dphi[mb][ptbin][0]/(Float_t)occurence[mb][ptbin]);
696 fV0Data[ptbin][1]->Fill(mb*0.0034+0.99+0.0017, _dphi[mb][ptbin][1]/(Float_t)occurence[mb][ptbin]);
5725c281 697 }
698 }
699}
700//_____________________________________________________________________________
19d77c98 701void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
702{
b4f29a18 703 // UserExec: called for each event. Commented where necessary
5725c281 704 if(fDebug > 0 ) cout << " *** UserExec() *** " << endl;
705 TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
706 if(fEventMixing) {
707 MixingCandidates = new TObjArray();
708 MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
709 }
710 if (!fPIDResponse) {
711 if(fDebug > 0 ) cout << " Could not get PID response " << endl;
b4f29a18 712 return;
713 }
5725c281 714 fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
715 if (fAOD) {
6940e1f0 716 if (!EventCut(fAOD)) return; // check for event cuts
5725c281 717 InitializeBayesianPID(fAOD); // init event objects
48077518 718 SetNullCuts(fAOD);
719 PrepareFlowEvent(fAOD->GetNumberOfTracks());
b4f29a18 720 fCandidates->SetLast(-1);
5725c281 721 if(fIsMC) IsMC(); // launch mc stuff
863bc87f 722 if(fQA) fEventStats->Fill(0);
48077518 723 Int_t unTracks = fAOD->GetNumberOfTracks();
724 AliAODTrack* un[unTracks];
725 AliAODTrack* up[unTracks];
726 Int_t unp(0);
727 Int_t unn(0);
5725c281 728 for (Int_t iTracks = 0; iTracks < unTracks; iTracks++) { // select analysis candidates
f15c1f69 729 AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
730 if(!track) AliFatal("Not a standard AOD");
6940e1f0 731 if (!PhiTrack(track)) continue;
5725c281 732 if (fQA) {
733 if(track->Charge() > 0) {fEventStats->Fill(1); fPtP->Fill(track->Pt());}
734 if(track->Charge() < 0) {fEventStats->Fill(2); fPtN->Fill(track->Pt());}
48077518 735 }
5725c281 736 if (IsKaon(track)) {
737 if(fEventMixing) MixingCandidates->Add(new AliPhiMesonHelperTrack(track->Eta(), track->Phi(), track->P(),
738 track->Px(), track->Py(), track->Pz(),
739 track->Pt(), track->Charge()));
740 if (track->Charge() > 0) {
48077518 741 up[unp] = track;
742 unp++;
863bc87f 743 if(fQA) {fEventStats->Fill(3);fPtKP->Fill(track->Pt());}
48077518 744 }
5725c281 745 else if (track->Charge() < 0) {
48077518 746 un[unn] = track;
747 unn++;
863bc87f 748 if(fQA) {fEventStats->Fill(4); fPtKN->Fill(track->Pt());}
48077518 749 }
750 }
19d77c98 751 }
5725c281 752 for (Int_t pTracks = 0; pTracks < unp ; pTracks++) { // perform analysis
753 for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
9fb90465 754// if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
48077518 755 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
756 PtSelector(0, up[pTracks], un[nTracks]);
48077518 757 Double_t pt = PhiPt(up[pTracks], un[nTracks]);
758 Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
759 TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
760 TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
761 TVector3 c = a + b;
762 Double_t phi = c.Phi();
763 Double_t eta = c.Eta();
0d91ada7 764 Double_t p = TMath::Sqrt(c.Px()*c.Px()+c.Py()*c.Py()+c.Pz()*c.Pz());
48077518 765 Int_t nIDs[2];
766 nIDs[0] = up[pTracks]->GetID();
767 nIDs[1] = un[nTracks]->GetID();
0d91ada7 768 MakeTrack(mass, pt, phi, eta, 2, nIDs, p, c.Pz());
b4f29a18 769 }
770 }
ffef9156 771 if(fV0) VZEROSubEventAnalysis();
5725c281 772 if (fDebug > 0) printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates into flow events
773 for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
b4f29a18 774 AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
775 if (!cand) continue;
5725c281 776 if (fDebug > 1) printf(" --> Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
777 for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
778 if (fDebug>1) printf(" *** Daughter %d with fID %d ***", iDau, cand->GetIDDaughter(iDau));
779 for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
b4f29a18 780 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
781 if (!iRP) continue;
782 if (!iRP->InRPSelection()) continue;
5725c281 783 if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
784 if (fDebug > 1) printf(" was in RP set");
b4f29a18 785 iRP->SetForRPSelection(kFALSE);
786 fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
48077518 787 }
b4f29a18 788 }
5725c281 789 if (fDebug > 1) printf("\n");
48077518 790 }
b4f29a18 791 cand->SetForPOISelection(kTRUE);
792 fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
19d77c98 793 }
5725c281 794 if (fDebug > 0) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
5725c281 795 if(!fEventMixing) { // combinatorial background
796 for (Int_t pTracks = 0; pTracks < unp ; pTracks++) {
797 for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++) {
9fb90465 798// if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
5725c281 799 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
800 PtSelector(1, up[pTracks], up[nTracks]);
801 }
802 }
803 for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
804 for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++) {
9fb90465 805// if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
5725c281 806 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
807 PtSelector(2, un[nTracks], un[pTracks]);
808 }
809 }
48077518 810 }
5725c281 811 if(fEventMixing) ReconstructionWithEventMixing(MixingCandidates);
48077518 812 PostData(1, fOutputList);
b4f29a18 813 PostData(2, fFlowEvent);
48077518 814 }
19d77c98 815}
816//_____________________________________________________________________________
9e58d58d 817void AliAnalysisTaskPhiFlow::Exec(Option_t* c)
818{
819 // skip the event selection for SE task (e.g. for MC productions)
6940e1f0 820 if(fSkipEventSelection) AliAnalysisTaskPhiFlow::UserExec(c);
9e58d58d 821 // exec of task se will do event selection and call UserExec
822 else AliAnalysisTaskSE::Exec(c);
823}
824//_____________________________________________________________________________
5725c281 825void AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing(TObjArray* MixingCandidates) const
826{
827 // perform phi reconstruction with event mixing
828 if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
829 AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
830 if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
831 if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
832 Int_t nEvents = pool->GetCurrentNEvents();
833 if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
834 for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
835 TObjArray* mixed_candidates = pool->GetEvent(iEvent);
836 if(!mixed_candidates) continue; // this should NEVER happen
837 Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
838 Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
839 if(fDebug > 0) cout << Form(" - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
840 AliPhiMesonHelperTrack* buffer_un[bufferTracks]; // set values for buffered candidates
841 AliPhiMesonHelperTrack* buffer_up[bufferTracks];
842 Int_t buffer_unp(0);
843 Int_t buffer_unn(0);
844 AliPhiMesonHelperTrack* mix_un[candidates];// set values for mixing candidates
845 AliPhiMesonHelperTrack* mix_up[candidates];
846 Int_t mix_unp(0);
847 Int_t mix_unn(0);
848 for (Int_t iTracks = 0; iTracks < candidates; iTracks++) { // distinguish between k+ and k- for mixing candidates
849 AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)MixingCandidates->At(iTracks);
850 if(!track) continue;
851 if (track->Charge() > 0) {
852 mix_up[mix_unp] = track;
853 mix_unp++;
854 }
855 else if (track->Charge() < 0 ) {
856 mix_un[mix_unn] = track;
857 mix_unn++;
858 }
859 }
860 for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) { // distinguish between k+ and k- for buffered candidates
861 AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)mixed_candidates->At(iTracks);
862 if(!track) continue;
863 if (track->Charge() > 0) {
864 buffer_up[buffer_unp] = track;
865 buffer_unp++;
866 }
867 else if (track->Charge() < 0 ) {
868 buffer_un[buffer_unn] = track;
869 buffer_unn++;
870 }
871 }
872 for (Int_t pMix = 0; pMix < mix_unp; pMix++) { // mix k+ (current event) k+ (buffer)
873 if(fDebug > 1 ) cout << Form("mixing current k+ track %d with", pMix);
874 if(!fTypeMixing) { // mix with like-sign kaons
875 for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
876 if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
9fb90465 877// if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_up[pBuf]))) continue;
5725c281 878 if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_up[pBuf]))) continue;
879 PtSelector(1, mix_up[pMix], buffer_up[pBuf]);
880 }
881 }
882 else if(fTypeMixing) { // mix with unlike-sign kaons
883 for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
884 if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
9fb90465 885// if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_un[nBuf]))) continue;
5725c281 886 if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_un[nBuf]))) continue;
887 PtSelector(2, mix_up[pMix], buffer_un[nBuf]);
888 }
889 }
890 }
891 for (Int_t nMix = 0; nMix < mix_unn; nMix++) { // mix k- (current event) k- (buffer)
892 if(fDebug > 1 ) cout << Form("mixing current k- track %d with", nMix);
893 if(!fTypeMixing) { // mix with like-sign kaons
894 for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
895 if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
9fb90465 896// if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_un[nBuf]))) continue;
5725c281 897 if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_un[nBuf]))) continue;
898 PtSelector(2, mix_un[nMix], buffer_un[nBuf]);
899 }
900 }
901 else if(fTypeMixing) { // mix with unlike-sign kaons
902 for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
903 if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
9fb90465 904// if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_up[pBuf]))) continue;
5725c281 905 if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_up[pBuf]))) continue;
906 PtSelector(1, mix_un[nMix], buffer_up[pBuf]);
907 }
908 }
909 }
910 } // end of mixed events loop
911 } // end of checking to see whether pool is filled correctly
912 pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
913 if(fDebug > 0 ) pool->PrintInfo();
914}
915//_____________________________________________________________________________
19d77c98 916void AliAnalysisTaskPhiFlow::Terminate(Option_t *)
917{
5725c281 918 // terminate
919 if(fDebug > 0) cout << " *** Terminate() *** " << endl;
19d77c98 920}
921//______________________________________________________________________________
0d91ada7 922void AliAnalysisTaskPhiFlow::MakeTrack(Double_t mass, Double_t pt, Double_t phi, Double_t eta, Int_t nDau, Int_t iID[], Double_t p, Double_t pz) const
19d77c98 923{
5725c281 924 // Construct Flow Candidate Track from two selected candidates
925 if(fDebug > 1 ) cout << " *** MakeTrack() *** " << endl;
0d91ada7 926 // if requested, check rapidity at this point
927 if(fCandidateYCut) {
928 Double_t y = 0.5*TMath::Log((TMath::Sqrt(mass*mass+p*p)+pz)/(TMath::Sqrt(mass*mass+p*p)-pz));
929 if (y > fCandidateMaxY || y < fCandidateMinY) return;
930 }
b4f29a18 931 Bool_t overwrite = kTRUE;
932 AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
5725c281 933 if (!sTrack) {
b4f29a18 934 sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
935 overwrite = kFALSE;
936 }
937 else sTrack->ClearMe();
19d77c98 938 sTrack->SetMass(mass);
939 sTrack->SetPt(pt);
940 sTrack->SetPhi(phi);
941 sTrack->SetEta(eta);
b4f29a18 942 for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
19d77c98 943 sTrack->SetForPOISelection(kTRUE);
944 sTrack->SetForRPSelection(kFALSE);
b4f29a18 945 if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
946 else fCandidates->AddLast(sTrack);
947 return;
948}
949//_____________________________________________________________________________
863bc87f 950void AliAnalysisTaskPhiFlow::IsMC()
951{
952 // Fill QA histos for MC analysis
953 TClonesArray *arrayMC = 0;
5725c281 954 if(fDebug > 0) cout << " -> Switching to MC mode <- " << endl;
863bc87f 955 // fill array with mc tracks
956 arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
957 if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
5725c281 958 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
f15c1f69 959 AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
960 if(!track) AliFatal("Not a standard AOD");
6940e1f0 961 if(!PhiTrack(track) || !IsKaon(track)) { // check for kaons
5725c281 962 if(fDebug>1) cout << " Rejected track" << endl;
863bc87f 963 continue;
964 }
5725c281 965 if (fDebug>1) cout << " Received MC kaon " << endl;
863bc87f 966 Double_t b[2] = { -99., -99.};
967 Double_t bCov[3] = { -99., -99., -99.};
646a809c 968 AliAODTrack copy(*track);
969 if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return;
863bc87f 970 // find corresponding mc particle
971 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
5725c281 972 if (!partMC) {
973 if(fDebug > 1) cout << "Cannot get MC particle" << endl;
863bc87f 974 continue;
975 }
976 // Check if it is primary, secondary from material or secondary from weak decay
977 Bool_t isPrimary = partMC->IsPhysicalPrimary();
978 Bool_t isSecondaryMaterial = kFALSE;
979 Bool_t isSecondaryWeak = kFALSE;
5725c281 980 if (!isPrimary) {
863bc87f 981 Int_t mfl = -999, codemoth = -999;
5725c281 982 Int_t indexMoth = partMC->GetMother();
983 if (indexMoth >= 0) { // is not fake
863bc87f 984 AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
985 codemoth = TMath::Abs(moth->GetPdgCode());
986 mfl = Int_t(codemoth / TMath::Power(10, Int_t(TMath::Log10(codemoth))));
987 }
5725c281 988 if (mfl == 3) isSecondaryWeak = kTRUE;
863bc87f 989 else isSecondaryMaterial = kTRUE;
990 }
991 if (isPrimary) {
5725c281 992 fDCAPrim->Fill(track->Pt(), b[0]);
863bc87f 993 fDCAXYQA->Fill(b[0]);
994 fDCAZQA->Fill(b[1]);
995 }
996 if (isSecondaryWeak) fDCASecondaryWeak->Fill(track->Pt(), b[0]);
997 if (isSecondaryMaterial) fDCAMaterial->Fill(track->Pt(), b[0]);
998 }
999}
5725c281 1000//_____________________________________________________________________________