AliPhysics  1c9c77b (1c9c77b)
 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 "AliTPCPIDResponse.h"
48 #include "AliAODMCParticle.h"
49 #include "AliAnalysisTaskVnV0.h"
50 #include "AliEventPoolManager.h"
51 #include "AliMultSelection.h"
52 #include "TMatrixDSym.h"
53 #include "AliVVertex.h"
54 #include "AliAODVertex.h"
55 
56 class AliFlowTrackCuts;
57 
58 using std::cout;
59 using std::endl;
60 
63 
65  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), fCentralityEstimator("V0M"), fPOICuts(0), fVertexRange(10.), 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), fPileUp(kTRUE), fLowCut(0), fHighCut(0), fMultTOFLowCut(0), fMultTOFHighCut(0), fHistCentralityWeights(0x0), fCentralityWeight(1.), fVertexZ(0), fHarmonic(2)
66 {
67  // Default constructor
68  for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 1000.;
69  for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
70  for(Int_t i(0); i < 20; i++) {
71  fVertexMixingBins[i] = 0;
72  fCentralityMixingBins[i] = 0;
73  }
74  fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
75  for(Int_t i(0); i < 18; i++) {
76  for(Int_t j(0); j < 2; j++) fV0Data[i][j] = 0;
77  fInvMNP[i] = 0; fInvMNN[i] = 0; fInvMPP[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.;
78  }
79 }
80 //_____________________________________________________________________________
82  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), fCentralityEstimator("V0M"), fPOICuts(0), fVertexRange(10.), 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), fPileUp(kTRUE), fLowCut(0), fHighCut(0), fMultTOFLowCut(0), fMultTOFHighCut(0), fHistCentralityWeights(0x0), fCentralityWeight(1.), fVertexZ(0), fHarmonic(2)
83 {
84  // Constructor
85  for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 1000.;
86  for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
87  for(Int_t i(0); i < 20; i++) {
88  fVertexMixingBins[i] = 0;
89  fCentralityMixingBins[i] = 0;
90  }
91  fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
92  for(Int_t i(0); i < 18; i++) {
93  for(Int_t j(0); j < 2; j++) fV0Data[i][j] = 0;
94  fInvMNP[i] = 0; fInvMNN[i] = 0; fInvMPP[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.;
95  }
96  DefineInput(0, TChain::Class());
97  DefineOutput(1, TList::Class());
98  DefineOutput(2, AliFlowEventSimple::Class());
99  if(fDebug > 0) cout << " === Created instance of AliAnaysisTaskPhiFlow === " << endl;
100 }
101 //_____________________________________________________________________________
103 {
104  // Destructor
105  if (fNullCuts) delete fNullCuts;
106  if (fOutputList) delete fOutputList;
107  if (fCandidates) delete fCandidates;
108  if (fFlowEvent) delete fFlowEvent;
110  if (fEventMixing) delete fPoolManager;
111  if (fPIDCombined) delete fPIDCombined;
112  if (fDebug > 0) cout << " === Deleted instance of AliAnalysisTaskPhiFlow === " << endl;
113 }
114 //_____________________________________________________________________________
116 {
117  // Return a pointer to a TH1 with predefined binning
118  if(fDebug > 0) cout << " *** BookHistogram() *** " << name << endl;
119  TH1F *hist = new TH1F(name, Form("M_{INV} (%s)", name), 60, .99, 1.092);
120  hist->GetXaxis()->SetTitle("M_{INV} (GeV / c^{2})");
121  hist->GetYaxis()->SetTitle("No. of pairs");
122  hist->SetMarkerStyle(kFullCircle);
123  hist->Sumw2();
124  fOutputList->Add(hist);
125  return hist;
126 }
127 //_____________________________________________________________________________
129 {
130  // Return a pointer to a TH2 with predefined binning
131  if(fDebug > 0) cout << " *** BookPIDHisotgram() *** " << endl;
132  TH2F *hist = 0x0;
133  if(TPC) {
134  hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
135  hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
136  }
137  if(!TPC) {
138  hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1.5);
139  hist->GetYaxis()->SetTitle("#beta");
140  }
141  hist->GetXaxis()->SetTitle("P (GeV / c)");
142  fOutputList->Add(hist);
143  return hist;
144 }
145 //_____________________________________________________________________________
147 {
148  // intialize p_t histograms for each p_t bin
149  if(fDebug > 0) cout << " *** InitPtSpectraHistograms() *** " << endl;
150  TH1F* hist = new TH1F(Form("%4.2f p_{t} %4.2f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 60, nmin, nmax);
151  hist->GetXaxis()->SetTitle("p_{T} GeV / c");
152  fOutputList->Add(hist);
153  return hist;
154 }
155 //_____________________________________________________________________________
157 {
158  // Return a pointer to a p_T spectrum histogram
159  if(fDebug > 0) cout << " *** BookPtHistogram() *** " << endl;
160  TH1F* ratio = new TH1F(name, name, 100, 0, 7);
161  ratio->GetXaxis()->SetTitle("p_{T} ( GeV / c^{2} )");
162  ratio->GetYaxis()->SetTitle("No. of events");
163  ratio->Sumw2();
164  fOutputList->Add(ratio);
165  return ratio;
166 }
167 //_____________________________________________________________________________
169 {
170  // Add Phi Identification Output Objects
171  if(fDebug > 0) cout << " ** AddPhiIdentificationOutputObjects() *** " << endl;
172  if(fQA) {
173  fEventStats = new TH1F("fHistStats", "Event Statistics", 18, -.25, 4.25);
174  fEventStats->GetXaxis()->SetTitle("No. of events");
175  fEventStats->GetYaxis()->SetTitle("Statistics");
176  fOutputList->Add(fEventStats);
177  }
178  fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
180  if(fQA) {
181  fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
183  fNOPID = BookPIDHistogram("TPC signal, all particles", kTRUE);
184  fPIDk = BookPIDHistogram("TPC signal, kaons", kTRUE);
185  fNOPIDTOF = BookPIDHistogram("TOF signal, all particles", kFALSE);
186  fPIDTOF = BookPIDHistogram("TOF signal, kaons", kFALSE);
187  }
188  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
189  fInvMNP[ptbin] = BookHistogram(Form("NP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
190  fInvMNN[ptbin] = BookHistogram(Form("NN, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
191  fInvMPP[ptbin] = BookHistogram(Form("PP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
192  }
193  if(fQA) {
194  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) fPtSpectra[ptbin] = InitPtSpectraHistograms(fPtBins[ptbin], fPtBins[ptbin+1]);
195  fPtP = BookPtHistogram("i^{+}");
196  fPtN = BookPtHistogram("i^{-}");
197  fPtKP = BookPtHistogram("K^{+}");
198  fPtKN = BookPtHistogram("K^{-}");
199  fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
200  fOutputList->Add(fPhi);
201  fPt = new TH1F("fPt", "p_{T}", 100, 0, 5.5);
202  fOutputList->Add(fPt);
203  fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
204  fOutputList->Add(fEta);
205  fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
206  fOutputList->Add(fVZEROA);
207  fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
208  fOutputList->Add(fVZEROC);
209  fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
210  fOutputList->Add(fTPCM);
211  fDCAXYQA = new TH1F("fDCAXYQA", "fDCAXYQA", 1000, -5, 5);
212  fOutputList->Add(fDCAXYQA);
213  fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
214  fOutputList->Add(fDCAZQA);
215  fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "centrality correlations;V0 centrality;CL0 centrality", 100, 0, 100, 100, 0, 100);
217  fMultvsCentr = new TH2F("fMultvsCentr", "multiplicty centrality;TPC multiplicity;V0 centrality", 250, 0, 5000, 100, 0, 100);
219  }
220  if(fIsMC || fQA) {
221  fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
222  fOutputList->Add(fDCAAll);
223  fDCAPrim = new TH2F("fDCAprim","fDCAprim", 1000, 0, 10, 1000, -5, 5);
224  fOutputList->Add(fDCAPrim);
225  fDCASecondaryWeak = new TH2F("fDCASecondaryWeak","fDCASecondaryWeak", 1000, 0, 10, 1000, -5, 5);
227  fDCAMaterial = new TH2F("fDCAMaterial","fDCAMaterial", 1000, 0, 10, 1000, -5, 5);
229  }
230  if(fV0) {
231  fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 5, 0, 5);
232  fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<#Psi_{a} - #Psi_{b}>");
233  fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<#Psi_{a} - #Psi_{c}>");
234  fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<#Psi_{b} - #Psi_{c}>");
235  fSubEventDPhiv2->GetXaxis()->SetBinLabel(4, "#sqrt{#frac{<#Psi_{a} - #Psi_{b}><#Psi_{a} - #Psi_{c}>}{<#Psi_{b} - #Psi_{c}>}}");
236  fSubEventDPhiv2->GetXaxis()->SetBinLabel(5, "#sqrt{#frac{<#Psi_{a} - #Psi_{c}><#Psi_{b} - #Psi_{c}>}{<#Psi_{a} - #Psi_{b}>}}");
238  const char* V0[] = {"V0A", "V0C"};
239  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++)
240  for(Int_t iV0(0); iV0 < 2; iV0++) {
241  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);
242  fOutputList->Add(fV0Data[ptbin][iV0]);
243  }
244  }
245 }
246 //_____________________________________________________________________________
248 {
249  // Create user defined output objects
250  if(fDebug > 0) cout << " *** UserCreateOutputObjects() *** " << endl;
251  fNullCuts = new AliFlowTrackCuts("null_cuts");
253  // combined pid
254  fPIDCombined = new AliPIDCombined;
255  fPIDCombined->SetDefaultTPCPriors();
256  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
257 
258  // flag to mc
260  Double_t t(0);
261  for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
262  if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis !!!");
263  if(fNPtBins > 18) AliFatal("Invalid number of pt bins initialied ( > 18 ) -- terminating analysis !!!");
265  cc->SetNbinsQ(500); cc->SetNbinsPhi(180); cc->SetNbinsMult(10000);
266  cc->SetQMin(0.0); cc->SetPhiMin(0.0); cc->SetMultMin(0);
267  cc->SetQMax(3.0); cc->SetPhiMax(TMath::TwoPi()); cc->SetMultMax(10000);
268  cc->SetNbinsMass(fMassBins); cc->SetNbinsEta(200); (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
269  cc->SetMassMin(fMinMass); cc->SetEtaMin(-5.0); cc->SetPtMin(0);
270  cc->SetMassMax(fMaxMass); cc->SetEtaMax(+5.0); (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
272  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
273  if (man) {
274  AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
275  if (inputHandler) fPIDResponse = inputHandler->GetPIDResponse();
276  }
277  // Create all output objects and store them to a list
278  fOutputList = new TList();
279  fOutputList->SetOwner(kTRUE);
280  // Create and post the Phi identification output objects
282  PostData(1, fOutputList);
283  // create candidate array
284  fCandidates = new TObjArray(1000);
285  fCandidates->SetOwner(kTRUE);
286  // create and post flowevent
287  fFlowEvent = new AliFlowEvent(10000);
288  PostData(2, fFlowEvent);
290 
291  fLowCut = new TF1("fLowCut", "[0]+[1]*x - 5.*([2]+[3]*x+[4]*x*x+[5]*x*x*x)", 0, 100);
292  fHighCut = new TF1("fHighCut", "[0]+[1]*x + 5.5*([2]+[3]*x+[4]*x*x+[5]*x*x*x)", 0, 100);
293 
294 
295 
296  fMultTOFLowCut = new TF1("fMultTOFLowCut", "[0]+[1]*x+[2]*x*x+[3]*x*x*x - 4.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x+[9]*x*x*x*x*x)", 0, 10000);
297 
298  fMultTOFHighCut = new TF1("fMultTOFHighCut", "[0]+[1]*x+[2]*x*x+[3]*x*x*x + 4.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x+[9]*x*x*x*x*x)", 0, 10000);
299 
300 
301  fLowCut->SetParameters(0.0157497, 0.973488, 0.673612, 0.0290718, -0.000546728, 5.82749e-06);
302  fHighCut->SetParameters(0.0157497, 0.973488, 0.673612, 0.0290718, -0.000546728, 5.82749e-06);
303 
304  fMultTOFLowCut->SetParameters(-1.0178, 0.333132, 9.10282e-05, -1.61861e-08, 1.47848, 0.0385923, -5.06153e-05, 4.37641e-08, -1.69082e-11, 2.35085e-15);
305  fMultTOFHighCut->SetParameters(-1.0178, 0.333132, 9.10282e-05, -1.61861e-08, 1.47848, 0.0385923, -5.06153e-05, 4.37641e-08, -1.69082e-11, 2.35085e-15);
306 
307 
308 
310  Double_t integral(fHistCentralityWeights->Integral()/fHistCentralityWeights->GetNbinsX());
311  TH1* temp((TH1*)fHistCentralityWeights->Clone("EP_weights_cen"));
312  for(Int_t i(0); i < temp->GetNbinsX(); i++) {
313  temp->SetBinError(1+i, 0.);// unertainty is irrelevant
314  temp->SetBinContent(1+i, integral/fHistCentralityWeights->GetBinContent(1+i));
315  }
316  delete fHistCentralityWeights;
317  fHistCentralityWeights = temp;
318  fOutputList->Add(temp);
319  }
320  fVertexZ = new TH1F("fVertexZ", "vertex Z position:vertex Z position", 60, -15., 15.);
321  fOutputList->Add(fVertexZ);
322 }
323 //_____________________________________________________________________________
325 {
326  // initialize event mixing
327  if(fDebug > 0) cout << " *** InitializeEventMixing() *** " << endl;
328  Int_t _c(0), _v(0);
329  for(Int_t i(0); i < 19; i++) {
330  if (fCentralityMixingBins[i+1] < fCentralityMixingBins[i]) { _c = i; break; }
331  else _c = 19;
332  }
333  for(Int_t i(0); i < 19; i++) {
334  if (fVertexMixingBins[i+1] < fVertexMixingBins[i]) { _v = i; break; }
335  else _v = 19;
336  }
337  if(fDebug > 0 ) cout << Form(" --> found %d centrality bins for mixing, %d vertex bins for mixing", _c, _v) << endl;
338  Double_t centralityBins[_c];
339  Double_t vertexBins[_v];
340  for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
341  for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
342  return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
343 }
344 //_____________________________________________________________________________
345 template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* track1, const T* track2) const
346 {
347  // Return the invariant mass of two tracks, assuming both tracks are kaons
348 // if(fDebug > 1) cout << " *** InvariantMass() *** " << endl;
349  if ((!track2) || (!track1)) return 0.;
350  Double_t masss = TMath::Power(4.93676999999999977e-01, 2);
351  Double_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
352  Double_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
353  Double_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
354  Double_t e1 = TMath::Sqrt(track1->P() * track1->P() + masss);
355  Double_t e2 = TMath::Sqrt(track2->P() * track2->P() + masss);
356  Double_t es = TMath::Power((e1 + e2), 2);
357  if ((es - (pxs + pys + pzs)) < 0) return 0.;
358  return TMath::Sqrt((es - (pxs + pys + pzs)));
359 }
360 //_____________________________________________________________________________
361 /*
362  template <typename T> Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const T* track1, const T* track2) const
363  {
364 // Calculate the delta dip angle between two particles
365 if(fDebug > 1) cout << " *** DeltaDipAngle() *** " << endl;
366 if (track1->P()*track2->P() == 0) return 999;
367 return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
368 }
369 //_____________________________________________________________________________
370 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const T* track1, const T* track2) const
371 {
372 // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
373 if(fDebug > 1) cout << " *** CheckDeltaDipAngle() *** " << endl;
374 if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PhiPt(track1, track2) < fDeltaDipPt)) return kFALSE;
375 return kTRUE;
376 }
377 */
378 //_____________________________________________________________________________
379 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
380 {
381  // Check if pair passes eta and pt cut
382 // if(fDebug > 1) cout << " *** CheckCandidateEtaPtCut() *** " << endl;
383  if (fCandidateMinPt > PhiPt(track1, track2) || fCandidateMaxPt < PhiPt(track1, track2)) return kFALSE;
384  TVector3 a(track1->Px(), track1->Py(), track1->Pz());
385  TVector3 b(track2->Px(), track2->Py(), track2->Pz());
386  TVector3 c = a + b;
387  if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
388  return kTRUE;
389 }
390 //_____________________________________________________________________________
391 template <typename T> Bool_t AliAnalysisTaskPhiFlow::EventCut(T* event)
392 {
393  // Impose event cuts
394 // if(fDebug > 0) cout << " *** EventCut() *** " << endl;
395  if (!event) return kFALSE;
396  // if (fSkipEventSelection) return kTRUE;
397  if (!CheckCentrality(event)) return kFALSE;
398  if(fQA) PlotMultiplcities(event);
399  return kTRUE;
400 }
401 //_____________________________________________________________________________
402 template <typename T> void AliAnalysisTaskPhiFlow::PlotMultiplcities(const T* event) const
403 {
404  // QA multiplicity plots
405  if(fDebug > 1) cout << " *** PlotMultiplcities() *** " << endl;
406  fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A(), fCentralityWeight);
407  fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C(), fCentralityWeight);
408  fTPCM->Fill(event->GetNumberOfTracks(), fCentralityWeight);
409 }
410 //_____________________________________________________________________________
411 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckVertex(const T* event)
412 {
413  // Check if event vertex is within given range
414 // if(fDebug > 0) cout << " *** CheckVertex() *** " << endl;
415  if (!event->GetPrimaryVertex()) return 0x0;
416  fVertex = event->GetPrimaryVertex()->GetZ();
417  if (TMath::Abs(fVertex) > fVertexRange) return 0x0;
418  return kTRUE;
419 }
420 //_____________________________________________________________________________
421 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCentrality(T* event)
422 {
423  // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
424 // if(fDebug > 0) cout << " *** CheckCentrality() *** " << endl;
425  if (!fkCentralityMethodA) AliFatal("No centrality method set! FATAL ERROR!");
426 
427  // check if the AliMultSelection object is present. If so, we should invoke the
428  // new centrality framework
429 
430  AliMultSelection *multSelection = 0x0;
431  multSelection = static_cast<AliMultSelection*>(event->FindListObject("MultSelection"));
432  if(multSelection) {
433  fCentrality = multSelection->GetMultiplicityPercentile(fCentralityEstimator.Data());
434  if(fCentrality > fCentralityMin && fCentrality < fCentralityMax && multSelection->GetMultiplicityPercentile("CL1") < 90) {
437  }
439  return kTRUE;
440  } else {
442  return kFALSE;
443  }
444  }
445 
446  else fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodA);
447  Double_t cenB(-999);
448  // if a second centrality estimator is requited, set it
449  (fkCentralityMethodB) ? cenB = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodB) : cenB = fCentrality;
450  if (TMath::Abs(fCentrality-cenB) > 5 || cenB >= 80 || cenB < 0 || fCentrality <= fCentralityMin || fCentrality > fCentralityMax) {
451  if(fQA) fCentralityNoPass->Fill(fCentrality) ;
452  return kFALSE;
453  }
454  const Int_t nGoodTracks = event->GetNumberOfTracks();
455  if(fCentralityCut2010) { // cut on outliers
456  Float_t multTPC(0.); // tpc mult estimate
457  Float_t multGlob(0.); // global multiplicity
458  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
459  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
460  if(!trackAOD) AliFatal("Not a standard AOD");
461  if (!trackAOD) continue;
462  if (!(trackAOD->TestFilterBit(1))) continue;
463  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;
464  multTPC++;
465  }
466  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
467  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
468  if(!trackAOD) AliFatal("Not a standard AOD");
469  if (!trackAOD) continue;
470  if (!(trackAOD->TestFilterBit(16))) continue;
471  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;
472  Double_t b[2] = {-99., -99.};
473  Double_t bCov[3] = {-99., -99., -99.};
474  AliAODTrack copy(*trackAOD);
475  if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
476  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
477  multGlob++;
478  } //track loop
479  // printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
480  if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
481  if(fQA) {
482  fMultCorAfterCuts->Fill(multGlob, multTPC);
483  fMultvsCentr->Fill(fCentrality, multTPC);
484  }
485  }
486 
487  if(fCentralityCut2011) { // cut on outliers
488  Float_t multTPC(0.); // tpc mult estimate
489  Float_t multGlob(0.); // global multiplicity
490  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
491  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
492  if(!trackAOD) AliFatal("Not a standard AOD");
493  if (!trackAOD) continue;
494  if (!(trackAOD->TestFilterBit(1))) continue;
495  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;
496  multTPC++;
497  }
498  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
499  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
500  if(!trackAOD) AliFatal("Not a standard AOD");
501  if (!trackAOD) continue;
502  if (!(trackAOD->TestFilterBit(16))) continue;
503  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;
504  Double_t b[2] = {-99., -99.};
505  Double_t bCov[3] = {-99., -99., -99.};
506  AliAODTrack copy(*trackAOD);
507  if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
508  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
509  multGlob++;
510  } //track loop
511  //printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
512  if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))) return kFALSE;
513  if(fQA) {
514  fMultCorAfterCuts->Fill(multGlob, multTPC);
515  fMultvsCentr->Fill(fCentrality, multTPC);
516  }
517  }
518 
520  return kTRUE;
521 }
522 //_____________________________________________________________________________
524 {
525  // Initialize the Bayesian PID object for AOD
526  if(fDebug > 0) cout << " *** InitializeBayesianPID() *** " << endl;
528 }
529 //_____________________________________________________________________________
530 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* track) const
531 {
532  // Check if the particle passes the TPC TOF bayesian cut.
533  if(fDebug > 1) cout << " *** PassesTPCbayesianCut() *** " << endl;
535  if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
536  Float_t *probabilities = fBayesianResponse->GetProb();
537  if (probabilities[3] > fPIDConfig[6]) {
538  if(fQA) {fPhi->Fill(track->Phi()); fPt->Fill(track->Pt()); fEta->Fill(track->Eta());}
539  return kTRUE;
540  }
541  return kFALSE;
542 }
543 //_____________________________________________________________________________
545 {
546  // check if track passes dca cut according to dca routine
547  // setup the routine as follows:
548  // fDCAConfig[0] < -1 no pt dependence
549  // fDCAConfig[0] = 0 do nothing
550  // fDCAConfig[0] > 1 pt dependent dca cut
551 // if(fDebug > 1) cout << " *** PassesDCACut() *** " << endl;
552 
553  AliAODTrack copy(*track);
554  if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
555  if(fQA) {
556  Double_t b[2] = { -99., -99.};
557  Double_t bCov[3] = { -99., -99., -99.};
558  if(copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) {
559  fDCAXYQA->Fill(b[0]);
560  fDCAZQA->Fill(b[1]);
561  fDCAPrim->Fill(track->Pt(), b[0]);
562  }
563  }
564  return kTRUE;
565  }
566  Double_t b[2] = { -99., -99.};
567  Double_t bCov[3] = { -99., -99., -99.};
568  if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
569  if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
570  if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
571  if(fDCAConfig[0] > .9) {
572  if(fDCAConfig[4] < TMath::Abs(b[1])) return kFALSE;
573  Double_t denom = TMath::Power(track->Pt(), fDCAConfig[3]);
574  if( denom < 0.0000001 ) return kFALSE; // avoid division by zero
575  if( (fDCAConfig[1] + fDCAConfig[2] / denom) < TMath::Abs(b[0]) ) return kFALSE;
576  }
577  if(fQA) {
578  fDCAXYQA->Fill(b[0]);
579  fDCAZQA->Fill(b[1]);
580  fDCAPrim->Fill(track->Pt(), b[0]);
581  }
582  return kTRUE;
583 }
584 //_____________________________________________________________________________
585 Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
586 {
587  // Kaon identification routine, based on multiple detectors and approaches
588 // if(fDebug > 1) cout << " *** IsKaon() *** " << endl;
589 
590  /* if(fUsePidResponse) {
591  Double_t prob[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
592  fPIDCombined->ComputeProbabilities(track, fPIDResponse, prob);
593  if(prob[3] > fPIDConfig[6]) return kTRUE;
594  }
595  if(!PassesDCACut(track)) return kFALSE;
596  if(fQA) {fNOPID->Fill(track->P(), track->GetTPCsignal());fNOPIDTOF->Fill(track->P(), track->GetTOFsignal());}
597  if(track->Pt() < fPIDConfig[1]) {
598  if(fDebug>1) cout << " ITS received track with p_t " << track->Pt() << endl;
599  // if tpc control is disabled, pure its pid
600  if(fPIDConfig[2] < 0.) {
601  if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) return kTRUE;
602  return kFALSE;
603  }
604  // else, switch to ITS pid with TPC rejection of protons and pions
605  if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) < fPIDConfig[3]) return kFALSE;
606  else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) < fPIDConfig[3]) return kFALSE;
607  else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) {
608  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
609  return kTRUE;
610  }
611  return kFALSE;
612  }
613  if((track->Pt() > fPIDConfig[1]) && (track->Pt() < fPIDConfig[4])) {
614  if(fDebug>1) cout << " TPC received track with p_t " << track->Pt() << endl;
615  // if its control is disabled, pure tpc pid
616  if(fPIDConfig[5] < 0.) {
617  if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) return kTRUE;
618  return kFALSE;
619  }
620  // else, switch to TPC pid with ITS rejection of protons and pions
621  if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kProton)) < fPIDConfig[0]) return kFALSE;
622  else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kPion)) < fPIDConfig[0]) return kFALSE;
623  else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) {
624  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
625  return kTRUE;
626  }
627  return kFALSE;
628  }
629  if(fDebug>1) cout << " Bayesian method received track with p_t " << track->Pt() << endl;
630  // switch to bayesian PID
631  if (PassesTPCbayesianCut(track)) {
632  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal());fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
633  return kTRUE;
634  }*/
635  if(fQA) {
636  Float_t length = track->GetIntegratedLength();
637  Float_t time = track->GetTOFsignal();
638  Double_t beta = -.05;
639 
640  if((length > 0) && (time > 0)) beta = length / 2.99792458e-2 / time;
641 
642  fNOPID->Fill(track->P(), track->GetTPCsignal(), fCentralityWeight);
643  fNOPIDTOF->Fill(track->P(), beta, fCentralityWeight);
644  }
645 
646 
647 
648  //kaon stuff
649  Double_t nSigKTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
650  Double_t nSigKTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
651 
652  Double_t nSigmaK = 999;
653 
654  if(track->Pt() > .4) {
655  nSigmaK = TMath::Sqrt(nSigKTPC*nSigKTPC+nSigKTOF*nSigKTOF);
656  } else {
657  if(track->GetTPCsignal() > 110) nSigmaK = TMath::Abs(nSigKTPC);
658  }
659  // pion
660  Double_t nSigPiTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
661  Double_t nSigPiTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
662 
663  Double_t nSigmaPi = 999;
664 
665  if(track->Pt() > .4) {
666  nSigmaPi = TMath::Sqrt(nSigPiTPC*nSigPiTPC+nSigPiTOF*nSigPiTOF);
667  } else {
668  if(track->GetTPCsignal() < 60) nSigmaPi = TMath::Abs(nSigPiTPC);
669  }
670 
671  //proton
672  Double_t nSigPTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
673  Double_t nSigPTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
674 
675  Double_t nSigmaP = 999;
676 
677  if(track->Pt() > .4) {
678  nSigmaP = TMath::Sqrt(nSigPTPC*nSigPTPC+nSigPTOF*nSigPTOF);
679  } else {
680  if(track->GetTPCsignal() > 110) nSigmaP = TMath::Abs(nSigPTPC);
681  }
682 
683  Short_t minSigma = FindMinNSigma(nSigmaPi, nSigmaK, nSigmaP);
684 
685 
686 
687  if(minSigma == 0) return kFALSE;
688  if((nSigmaK == nSigmaPi) && ( nSigmaK == nSigmaP)) return kFALSE;
689 
690  if(minSigma == 2 &&!GetDoubleCountingK(nSigmaK, minSigma)) {
691 
692  if(fQA) {
693  Float_t length = track->GetIntegratedLength();
694  Float_t time = track->GetTOFsignal();
695  Double_t beta = -.05;
696 
697  if((length > 0) && (time > 0)) beta = length / 2.99792458e-2 / time;
698 
699  // additional selection to bruteforce throw away bad betas
700  if(beta < 0.4) return kFALSE;
701 
702  fPIDk->Fill(track->P(), track->GetTPCsignal(), fCentralityWeight);
703  if(track->Pt() > .4) fPIDTOF->Fill(track->P(), beta, fCentralityWeight);
704  }
705 
706 
707  return kTRUE;
708  }
709 
710  return kFALSE;
711 }
712 //_____________________________________________________________________________
713 template <typename T> Double_t AliAnalysisTaskPhiFlow::PhiPt(const T* track1, const T* track2) const
714 {
715  // return p_t of track pair
716  TVector3 a(track1->Px(), track1->Py(), track1->Pz());
717  TVector3 b(track2->Px(), track2->Py(), track2->Pz());
718  TVector3 c = a + b;
719  return c.Pt();
720 }
721 //_____________________________________________________________________________
722 template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, const T* track1, const T* track2) const
723 {
724  // plot m_inv spectra of like- and unlike-sign kaon pairs for each pt bin
725  Double_t pt = PhiPt(track1, track2);
726  if (tracktype == 0) {
727  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
728  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
729  fInvMNP[ptbin]->Fill(InvariantMass(track1, track2), fCentralityWeight);
730  if(fQA) fPtSpectra[ptbin]->Fill(pt, fCentralityWeight);
731  }
732  }
733  }
734  if (tracktype == 1) {
735  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
736  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
737  fInvMPP[ptbin]->Fill(InvariantMass(track1, track2), fCentralityWeight);
738  }
739  }
740  }
741  if (tracktype == 2) {
742  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
743  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
744  fInvMNN[ptbin]->Fill(InvariantMass(track1, track2), fCentralityWeight);
745  }
746  }
747  }
748 }
749 //_____________________________________________________________________________
750 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
751 {
752  // check if track meets quality cuts
753  if(!track) return kFALSE;
754  return fPOICuts->IsSelected(track);
755 }
756 //_____________________________________________________________________________
757 template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
758 {
759  // Set null cuts
760 // if (fDebug > 0) cout << " *** SetNullCuts() *** " << fCutsRP << endl;
761  fCutsRP->SetEvent(event, MCEvent());
763  fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
764  fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
765  fNullCuts->SetEvent(event, MCEvent());
766 }
767 //_____________________________________________________________________________
769 {
770  // Prepare flow events
771 // if (fDebug > 0 ) cout << " *** PrepareFlowEvent() *** " << endl;
775  fFlowEvent->DefineDeadZone(0, 0, 0, 0);
776 }
777 //_____________________________________________________________________________
779 {
780  // UserExec: called for each event. Commented where necessary
781 // if(fDebug > 0 ) cout << " *** UserExec() *** " << endl;
782  TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
783  if(fEventMixing) {
784  MixingCandidates = new TObjArray();
785  MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
786  }
787  if (!fPIDResponse) {
788  if(fDebug > 0 ) cout << " Could not get PID response " << endl;
789  return;
790  }
791  fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
792  if (fAOD) {
793  if (!EventCut(fAOD)) return; // check for event cuts
794 
795 
796  // add extra pileup cuts for high intensity runs
797  // courtesy of alex dobrin
798  if (fPileUp){
799 
800 
801 
802  // 12 11 centrality cut for pileup
803  //
804  //Centrality
805  Float_t v0Centr = -100.;
806  Float_t cl1Centr = -100.;
807  Float_t cl0Centr = -100.;
808 
809  AliMultSelection* MultSelection = 0x0;
810  MultSelection = (AliMultSelection*) fAOD->FindListObject("MultSelection");
811  if( !MultSelection) {
812  AliWarning("AliMultSelection object not found!");
813  return;
814  } else {
815  v0Centr = MultSelection->GetMultiplicityPercentile("V0M");
816  cl1Centr = MultSelection->GetMultiplicityPercentile("CL1");
817  cl0Centr = MultSelection->GetMultiplicityPercentile("CL0");
818  }
819 
820  if (v0Centr >= 90. || v0Centr < 0)
821  return;
822 
823 
824  const Int_t nTracks = fAOD->GetNumberOfTracks();
825  Int_t multEsd = ((AliAODHeader*)fAOD->GetHeader())->GetNumberOfESDTracks();
826 
827  Int_t multTrk = 0;
828  Int_t multTrkBefC = 0;
829  Int_t multTrkTOFBefC = 0;
830  Int_t multTPC = 0;
831 
832  for (Int_t it = 0; it < nTracks; it++) {
833 
834  AliAODTrack* aodTrk = (AliAODTrack*)fAOD->GetTrack(it);
835 
836  if (!aodTrk){
837  delete aodTrk;
838  continue;
839  }
840 
841  if (aodTrk->TestFilterBit(32)){
842 
843  multTrkBefC++;
844 
845  if ( TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10 && aodTrk->GetTOFsignal() >= 12000 && aodTrk->GetTOFsignal() <= 25000)
846  multTrkTOFBefC++;
847 
848  if ((TMath::Abs(aodTrk->Eta()) < TMath::Abs(fCandidateMinEta)) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= .2) && (aodTrk->Pt() < 20))
849  multTrk++;
850 
851  }
852 
853  if (aodTrk->TestFilterBit(32))
854  multTPC++;
855 
856  }
857 
858 
859  Float_t multTPCn = multTPC;
860  Float_t multEsdn = multEsd;
861  Float_t multESDTPCDif = multEsdn - multTPCn*3.38;
862 
863 
864 
865  if (cl0Centr < fLowCut->Eval(v0Centr))
866  return;
867 
868  if (cl0Centr > fHighCut->Eval(v0Centr))
869  return;
870 
871 
872  if (multESDTPCDif > 15000)
873  return;
874 
875 
876  if (multTrkTOFBefC < fMultTOFLowCut->Eval(Float_t(multTrkBefC)))
877  return;
878 
879  if (multTrkTOFBefC > fMultTOFHighCut->Eval(Float_t(multTrkBefC)))
880  return;
881 
882 
883  if (plpMV(fAOD))
884  return;
885 
886  Short_t isPileup = fAOD->IsPileupFromSPD(3);
887  if (isPileup != 0)
888  return;
889 
890  if (((AliAODHeader*)fAOD->GetHeader())->GetRefMultiplicityComb08() < 0)
891  return;
892 
893 
894  /*
895  Int_t bc2 = ((AliAODHeader*)fAOD->GetHeader())->GetIRInt2ClosestInteractionMap();
896  if (bc2 != 0)
897  return;
898 
899  Int_t bc1 = ((AliAODHeader*)fAOD->GetHeader())->GetIRInt1ClosestInteractionMap();
900  if (bc1 != 0)
901  return;
902  */
903 
904  // add vertexer selection
905  AliAODVertex* vtTrc = fAOD->GetPrimaryVertex();
906  AliAODVertex* vtSPD = fAOD->GetPrimaryVertexSPD();
907  if (vtTrc->GetNContributors()<2 || vtSPD->GetNContributors()<1) return; // one of vertices is missing
908  double covTrc[6],covSPD[6];
909  vtTrc->GetCovarianceMatrix(covTrc);
910  vtSPD->GetCovarianceMatrix(covSPD);
911  double dz = vtTrc->GetZ()-vtSPD->GetZ();
912  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
913  double errTrc = TMath::Sqrt(covTrc[5]);
914  double nsigTot = TMath::Abs(dz)/errTot, nsigTrc = TMath::Abs(dz)/errTrc;
915  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) return; // bad vertexing
916 
917  if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fVertexRange) return;
918  fVertex = fAOD->GetPrimaryVertex()->GetZ();
919  fVertexZ->Fill(fVertex);
920 
921  fMultvsCentr->Fill(multTPCn, v0Centr);
922 
923  fMultCorAfterCuts->Fill(v0Centr,cl0Centr);
924 
925 
926  //new function for 2015 to remove incomplete events
927  if (fAOD->IsIncompleteDAQ()) return;
928  } else {
929  if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fVertexRange) return;
930  fVertex = fAOD->GetPrimaryVertex()->GetZ();
931  fVertexZ->Fill(fVertex);
932  }
933 
934 
935  // InitializeBayesianPID(fAOD); // init event objects
936 
937 
938  double qx(0); // for tpc ep
939  double qy(0);
940  SetNullCuts(fAOD);
941  PrepareFlowEvent(fAOD->GetNumberOfTracks());
942  fCandidates->SetLast(-1);
943  if(fIsMC) IsMC(); // launch mc stuff
944  if(fQA) fEventStats->Fill(0);
945  Int_t unTracks = fAOD->GetNumberOfTracks();
946  AliAODTrack* un[unTracks];
947  AliAODTrack* up[unTracks];
948  Int_t unp(0);
949  Int_t unn(0);
950 // if(fDebug > 1) cout << " started with " << unTracks << " of tracks "<< endl;
951  for (Int_t iTracks = 0; iTracks < unTracks; iTracks++) { // select analysis candidates
952  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
953  if(!track) AliFatal("Not a standard AOD");
954  if (!PhiTrack(track)) continue;
955  qx+=TMath::Cos(fHarmonic*track->Phi());
956  qy+=TMath::Sin(fHarmonic*track->Phi());
957  if (fQA) {
958  if(track->Charge() > 0) {fEventStats->Fill(1, fCentralityWeight); fPtP->Fill(track->Pt(), fCentralityWeight);}
959  if(track->Charge() < 0) {fEventStats->Fill(2, fCentralityWeight); fPtN->Fill(track->Pt(), fCentralityWeight);}
960  }
961  if (IsKaon(track)) {
962  if(fEventMixing) MixingCandidates->Add(new AliPhiMesonHelperTrack(track->Eta(), track->Phi(), track->P(),
963  track->Px(), track->Py(), track->Pz(),
964  track->Pt(), track->Charge()));
965  if (track->Charge() > 0) {
966  up[unp] = track;
967  unp++;
968  if(fQA) {fEventStats->Fill(3, fCentralityWeight);fPtKP->Fill(track->Pt(), fCentralityWeight);}
969  }
970  else if (track->Charge() < 0) {
971  un[unn] = track;
972  unn++;
973  if(fQA) {fEventStats->Fill(4, fCentralityWeight); fPtKN->Fill(track->Pt(), fCentralityWeight);}
974  }
975  }
976  }
977  for (Int_t pTracks = 0; pTracks < unp ; pTracks++) { // perform analysis
978  for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
979  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
980  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
981  PtSelector(0, up[pTracks], un[nTracks]);
982  Double_t pt = PhiPt(up[pTracks], un[nTracks]);
983  Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
984  TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
985  TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
986  TVector3 c = a + b;
987  Double_t phi = c.Phi();
988  Double_t eta = c.Eta();
989  Double_t p = TMath::Sqrt(c.Px()*c.Px()+c.Py()*c.Py()+c.Pz()*c.Pz());
990  Int_t nIDs[2];
991  nIDs[0] = up[pTracks]->GetID();
992  nIDs[1] = un[nTracks]->GetID();
993  MakeTrack(mass, pt, phi, eta, 2, nIDs, p, c.Pz());
994  }
995  }
996 // if (fDebug > 0) printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates into flow events
997  for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
998  AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
999  if (!cand) continue;
1000 // if (fDebug > 1) printf(" --> Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
1001  for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
1002  if (fDebug>1) printf(" *** Daughter %d with fID %d ***", iDau, cand->GetIDDaughter(iDau));
1003  for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
1004  AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
1005  if (!iRP) continue;
1006  if (!iRP->InRPSelection()) continue;
1007  if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
1008 // if (fDebug > 1) printf(" was in RP set");
1009  iRP->SetForRPSelection(kFALSE);
1011  }
1012  }
1013 // if (fDebug > 1) printf("\n");
1014  }
1015  cand->SetForPOISelection(kTRUE);
1016  fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
1017  }
1018 // if (fDebug > 0) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
1019  if(!fEventMixing) { // combinatorial background
1020  for (Int_t pTracks = 0; pTracks < unp ; pTracks++) {
1021  for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++) {
1022  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
1023  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
1024  PtSelector(1, up[pTracks], up[nTracks]);
1025  }
1026  }
1027  for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
1028  for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++) {
1029  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
1030  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
1031  PtSelector(2, un[nTracks], un[pTracks]);
1032  }
1033  }
1034  }
1035  if(fEventMixing) ReconstructionWithEventMixing(MixingCandidates);
1036 
1037 
1038 
1039  // ugly, ugly hack but i want to pass the centrality weight to the sp task
1041  // and second guly hack, save the tpc ep
1042  fFlowEvent->SetPsi2(qx);
1043  fFlowEvent->SetPsi3(qy);
1044  PostData(1, fOutputList);
1045  PostData(2, fFlowEvent);
1046  }
1047 }
1048 //_____________________________________________________________________________
1050 {
1051  // skip the event selection for SE task (e.g. for MC productions)
1053  // exec of task se will do event selection and call UserExec
1054  else AliAnalysisTaskSE::Exec(c);
1055 }
1056 //_____________________________________________________________________________
1058 {
1059  // perform phi reconstruction with event mixing
1060 // if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
1061  AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
1062  if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
1063  if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
1064  Int_t nEvents = pool->GetCurrentNEvents();
1065 // if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
1066  for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
1067  TObjArray* mixed_candidates = pool->GetEvent(iEvent);
1068  if(!mixed_candidates) continue; // this should NEVER happen
1069  Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
1070  Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
1071 // if(fDebug > 0) cout << Form(" - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
1072  AliPhiMesonHelperTrack* buffer_un[bufferTracks]; // set values for buffered candidates
1073  AliPhiMesonHelperTrack* buffer_up[bufferTracks];
1074  Int_t buffer_unp(0);
1075  Int_t buffer_unn(0);
1076  AliPhiMesonHelperTrack* mix_un[candidates];// set values for mixing candidates
1077  AliPhiMesonHelperTrack* mix_up[candidates];
1078  Int_t mix_unp(0);
1079  Int_t mix_unn(0);
1080  for (Int_t iTracks = 0; iTracks < candidates; iTracks++) { // distinguish between k+ and k- for mixing candidates
1081  AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)MixingCandidates->At(iTracks);
1082  if(!track) continue;
1083  if (track->Charge() > 0) {
1084  mix_up[mix_unp] = track;
1085  mix_unp++;
1086  }
1087  else if (track->Charge() < 0 ) {
1088  mix_un[mix_unn] = track;
1089  mix_unn++;
1090  }
1091  }
1092  for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) { // distinguish between k+ and k- for buffered candidates
1093  AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)mixed_candidates->At(iTracks);
1094  if(!track) continue;
1095  if (track->Charge() > 0) {
1096  buffer_up[buffer_unp] = track;
1097  buffer_unp++;
1098  }
1099  else if (track->Charge() < 0 ) {
1100  buffer_un[buffer_unn] = track;
1101  buffer_unn++;
1102  }
1103  }
1104  for (Int_t pMix = 0; pMix < mix_unp; pMix++) { // mix k+ (current event) k+ (buffer)
1105 // if(fDebug > 1 ) cout << Form("mixing current k+ track %d with", pMix);
1106  if(!fTypeMixing) { // mix with like-sign kaons
1107  for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
1108 // if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
1109  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_up[pBuf]))) continue;
1110  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_up[pBuf]))) continue;
1111  PtSelector(1, mix_up[pMix], buffer_up[pBuf]);
1112  }
1113  }
1114  else if(fTypeMixing) { // mix with unlike-sign kaons
1115  for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
1116 // if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
1117  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_un[nBuf]))) continue;
1118  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_un[nBuf]))) continue;
1119  PtSelector(2, mix_up[pMix], buffer_un[nBuf]);
1120  }
1121  }
1122  }
1123  for (Int_t nMix = 0; nMix < mix_unn; nMix++) { // mix k- (current event) k- (buffer)
1124 // if(fDebug > 1 ) cout << Form("mixing current k- track %d with", nMix);
1125  if(!fTypeMixing) { // mix with like-sign kaons
1126  for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
1127 // if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
1128  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_un[nBuf]))) continue;
1129  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_un[nBuf]))) continue;
1130  PtSelector(2, mix_un[nMix], buffer_un[nBuf]);
1131  }
1132  }
1133  else if(fTypeMixing) { // mix with unlike-sign kaons
1134  for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
1135 // if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
1136  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_up[pBuf]))) continue;
1137  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_up[pBuf]))) continue;
1138  PtSelector(1, mix_un[nMix], buffer_up[pBuf]);
1139  }
1140  }
1141  }
1142  } // end of mixed events loop
1143  } // end of checking to see whether pool is filled correctly
1144  pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
1145 // if(fDebug > 0 ) pool->PrintInfo();
1146 }
1147 //_____________________________________________________________________________
1149 {
1150  // terminate
1151  if(fDebug > 0) cout << " *** Terminate() *** " << endl;
1152 }
1153 //______________________________________________________________________________
1155 {
1156  // Construct Flow Candidate Track from two selected candidates
1157 // if(fDebug > 1 ) cout << " *** MakeTrack() *** " << endl;
1158  // if requested, check rapidity at this point
1159  if(fCandidateYCut) {
1160  Double_t y = 0.5*TMath::Log((TMath::Sqrt(mass*mass+p*p)+pz)/(TMath::Sqrt(mass*mass+p*p)-pz));
1161  if (y > fCandidateMaxY || y < fCandidateMinY) return;
1162  }
1163  Bool_t overwrite = kTRUE;
1164  AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
1165  if (!sTrack) {
1166  sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
1167  overwrite = kFALSE;
1168  }
1169  else sTrack->ClearMe();
1170  sTrack->SetMass(mass);
1171  sTrack->SetPt(pt);
1172  sTrack->SetPhi(phi);
1173  sTrack->SetEta(eta);
1174  for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
1175  sTrack->SetForPOISelection(kTRUE);
1176  sTrack->SetForRPSelection(kFALSE);
1177  if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
1178  else fCandidates->AddLast(sTrack);
1179  return;
1180 }
1181 //_____________________________________________________________________________
1183 {
1184  // Fill QA histos for MC analysis
1185  TClonesArray *arrayMC = 0;
1186  if(fDebug > 0) cout << " -> Switching to MC mode <- " << endl;
1187  // fill array with mc tracks
1188  arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1189  if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
1190  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
1191  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
1192  if(!track) AliFatal("Not a standard AOD");
1193  if(!PhiTrack(track) || !IsKaon(track)) { // check for kaons
1194  if(fDebug>1) cout << " Rejected track" << endl;
1195  continue;
1196  }
1197  if (fDebug>1) cout << " Received MC kaon " << endl;
1198  Double_t b[2] = { -99., -99.};
1199  Double_t bCov[3] = { -99., -99., -99.};
1200  AliAODTrack copy(*track);
1201  if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return;
1202  // find corresponding mc particle
1203  AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
1204  if (!partMC) {
1205  if(fDebug > 1) cout << "Cannot get MC particle" << endl;
1206  continue;
1207  }
1208  // Check if it is primary, secondary from material or secondary from weak decay
1209  Bool_t isPrimary = partMC->IsPhysicalPrimary();
1210  Bool_t isSecondaryMaterial = kFALSE;
1211  Bool_t isSecondaryWeak = kFALSE;
1212  if (!isPrimary) {
1213  Int_t mfl = -999, codemoth = -999;
1214  Int_t indexMoth = partMC->GetMother();
1215  if (indexMoth >= 0) { // is not fake
1216  AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
1217  codemoth = TMath::Abs(moth->GetPdgCode());
1218  mfl = Int_t(codemoth / TMath::Power(10, Int_t(TMath::Log10(codemoth))));
1219  }
1220  if (mfl == 3) isSecondaryWeak = kTRUE;
1221  else isSecondaryMaterial = kTRUE;
1222  }
1223  if (isPrimary) {
1224  fDCAPrim->Fill(track->Pt(), b[0]);
1225  fDCAXYQA->Fill(b[0]);
1226  fDCAZQA->Fill(b[1]);
1227  }
1228  if (isSecondaryWeak) fDCASecondaryWeak->Fill(track->Pt(), b[0]);
1229  if (isSecondaryMaterial) fDCAMaterial->Fill(track->Pt(), b[0]);
1230  }
1231 }
1232 
1233 
1234 
1235 //_____________________________________________________________________________
1236 Double_t AliAnalysisTaskPhiFlow::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1237 {
1238 
1239  // calculate sqrt of weighted distance to other vertex
1240  if (!v0 || !v1) {
1241  printf("One of vertices is not valid\n");
1242  return 0;
1243  }
1244  static TMatrixDSym vVb(3);
1245  double dist = -1;
1246  double dx = v0->GetX()-v1->GetX();
1247  double dy = v0->GetY()-v1->GetY();
1248  double dz = v0->GetZ()-v1->GetZ();
1249  double cov0[6],cov1[6];
1250  v0->GetCovarianceMatrix(cov0);
1251  v1->GetCovarianceMatrix(cov1);
1252  //
1253  // fQxavsV0[0] fQxnmV0A
1254  // fQyavsV0[0] fQynmV0A
1255  // fQxavsV0[1] fQxnsV0A
1256  // fQyavsV0[1] fQynsV0A
1257  // fQxcvsV0[0] fQxnmV0C
1258  // fQycvsV0[0] fQynmV0C
1259  // fQxcvsV0[1] fQxnsV0C
1260  // fQycvsV0[1] fQynsV0C vVb(0,0) = cov0[0]+cov1[0];
1261  vVb(1,1) = cov0[2]+cov1[2];
1262  vVb(2,2) = cov0[5]+cov1[5];
1263  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1264  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1265  vVb.InvertFast();
1266  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
1267  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1268  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1269  return dist>0 ? TMath::Sqrt(dist) : -1;
1270 
1271 }
1272 
1273 
1274 
1275 
1276 //_____________________________________________________________________________
1278 {
1279  // check for multi-vertexer pile-up
1280  //
1281  const int kMinPlpContrib = 5;
1282  const double kMaxPlpChi2 = 5.0;
1283  const double kMinWDist = 15;
1284  //
1285  const AliVVertex* vtPrm = 0;
1286  const AliVVertex* vtPlp = 0;
1287  int nPlp = 0;
1288  //
1289  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
1290  vtPrm = aod->GetPrimaryVertex();
1291  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
1292 
1293  //int bcPrim = vtPrm->GetBC();
1294  //
1295  for (int ipl=0;ipl<nPlp;ipl++) {
1296  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1297  //
1298  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1299  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1300  // int bcPlp = vtPlp->GetBC();
1301  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
1302  //
1303  double wDst = GetWDist(vtPrm,vtPlp);
1304  if (wDst<kMinWDist) continue;
1305  //
1306  return kTRUE; // pile-up: well separated vertices
1307  }
1308  //
1309  return kFALSE;
1310  //
1311 }
1312 
1313 //_____________________________________________________________________________
1314 //---------------------------------------------------
1316 {
1317 
1318  Short_t kPID = 0;
1319 
1320  if((nSk == nSpi) && (nSk == nSp))
1321  return kPID;
1322  if((nSk < nSpi) && (nSk < nSp) && (nSk < fPIDConfig[0]))
1323  kPID = 2;
1324 
1325  if((nSpi < nSk) && (nSpi < nSp) && (nSpi < fPIDConfig[0]))
1326  kPID = 1;
1327 
1328  if((nSp < nSk) && (nSp < nSpi) && (nSp < fPIDConfig[0]))
1329  kPID = 3;
1330  return kPID;
1331 
1332 }
1333 
1334 
1335 //_____________________________________________________________________________
1337 {
1338 
1339  Bool_t kDC = kFALSE;
1340 
1341  if (nSk < fPIDConfig[0] && minNSigma != 2)
1342  kDC = kTRUE;
1343 
1344  return kDC;
1345 
1346 }
1347 
1348 
Bool_t plpMV(const AliAODEvent *aod)
const Color_t cc[]
Definition: DrawKs.C:1
Bool_t PassesTPCbayesianCut(T *track) const
TProfile * fV0Data[18][2]
subevent resolution info for v2
Double_t fCentralityMin
QA profile of centralty vs multiplicity.
double Double_t
Definition: External.C:58
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.
Definition: External.C:236
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)
TCanvas * c
Definition: TestFitELoss.C:172
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)
AliESDv0 * fV0
placeholder for the current kink
Double_t GetWDist(const AliVVertex *v0, const AliVVertex *v1)
Bool_t InRPSelection() const
AliFlowBayesianPID * fBayesianResponse
AliFlowTrack * GetTrack(Int_t i)
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
void SetMC(Bool_t flag)
Bool_t fSkipEventSelection
profiles for vzero vn(minv)
void SetPsi5(Double_t gPsi5)
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
virtual void Exec(Option_t *)
void SetNumberOfRPs(Int_t nr)
void Fill(AliFlowTrackCuts *rpCuts, AliFlowTrackCuts *poiCuts)
int Int_t
Definition: External.C:63
Bool_t PassesDCACut(AliAODTrack *track) const
void PrepareFlowEvent(Int_t iMulti)
float Float_t
Definition: External.C:68
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.
TH1F * fEta
QA plot of p_t sectrum of tracks used for event plane estimation.
TH1F * fInvMNN[18]
unlike sign kaon pairs
AliEventPoolManager * InitializeEventMixing()
void SetPsi2(Double_t gPsi2)
Bool_t PhiTrack(T *track) const
AliPIDResponse * fPIDResponse
void SetPhi(Double_t phi)
TH2F * fMultCorAfterCuts
QA histogram of p_t distribution of negative kaons.
short Short_t
Definition: External.C:23
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
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)
Float_t fHarmonic
z vertex position
void InitializeBayesianPID(AliAODEvent *event)
void SetPsi3(Double_t gPsi3)
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)
Bool_t GetDoubleCountingK(Double_t nSk, Short_t minNSigma) 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.
const char Option_t
Definition: External.C:48
Bool_t IsKaon(AliAODTrack *track) const
TH2F * fDCAAll
QA plot TPC multiplicity (tracks used for event plane estimation)
TList * fOutputList
event pool manager
bool Bool_t
Definition: External.C:53
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)
Short_t FindMinNSigma(Double_t nSpi, Double_t nSk, Double_t nSp) const
TH2F * BookPIDHistogram(const char *name, Bool_t TPC)
virtual void UserExec(Option_t *option)
void InsertTrack(AliFlowTrack *)
virtual void ClearFast()
AliEventPoolManager * fPoolManager
AOD oject.
Definition: External.C:196
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