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